add check for aimec new coordinate system

remove unused imports

add pytest

remove unused return variable

Author: awirb

avaframe/ana3AIMEC/aimecTools.py

@@ -7,7 +7,8 @@
 import numpy as np
 import pandas as pd
 import copy
-
+import shapely as shp
+import itertools as iterT
 
 # Local imports
 import avaframe.in2Trans.shpConversion as shpConv
@@ -21,7 +22,6 @@
 import avaframe.out3Plot.outAIMEC as outAimec
 import avaframe.out3Plot.plotUtils as pU
 
-
 # create local logger
 log = logging.getLogger(__name__)
 
@@ -522,6 +522,9 @@ def makeDomainTransfo(pathDict, dem, refCellSize, cfgSetup):
     # get z coordinate of the center line
     rasterTransfo, _ = geoTrans.projectOnRaster(dem, rasterTransfo)
 
+    # check if any coordinate lines of new coordinate system are overlapping
+    _ = checkOverlapDBXY(rasterTransfo)
+
     # add surface parallel coordinate (sParallel)
     rasterTransfo = addSurfaceParalleCoord(rasterTransfo)
 
@@ -1837,3 +1840,41 @@ def analyzeDiffsRunoutLines(cfgSetup, runoutLine, refDataTransformed, resAnalysi
             log.info('For reference data type %s, runout line comparison is not available' % refLine['type'])
 
     return resAnalysisDF
+
+
+def checkOverlapDBXY(rasterTransfo):
+    """check if x, y coordinates of new coordinate system overlap
+
+    Parameters
+    -----------
+    rasterTransfo: dict
+        dict with gridx, gridy locations of new coordinates in old coordinate system
+
+    """
+
+    x = rasterTransfo["gridx"]
+    y = rasterTransfo["gridy"]
+
+    # create lineStrings from coordinate points
+    multiLine = []
+    for i in range(x.shape[1]):
+        lineArray = np.zeros((x.shape[0], 2))
+        lineArray[:, 0] = x[:, i]
+        lineArray[:, 1] = y[:, i]
+        multiLine.append(shp.LineString(lineArray))
+
+    # check for intersections of the individual lines in the multiline
+    for line1, line2 in iterT.combinations(multiLine, 2):
+        if line1.crosses(line2):
+            intersectionPoints = line1.intersection(line2)
+            if isinstance(intersectionPoints, shp.Point) or isinstance(intersectionPoints, shp.MultiPoint):
+                message = (
+                    "New coordinate system has overlaps - first intersection at %s. "
+                    "The provided path has too much curvature. "
+                    "Try: (1) smoothing the path, (2) reducing the domain width, or (3) using fewer points."
+                    % intersectionPoints
+                )
+                log.error(message)
+                raise AssertionError(message)
+
+    return True

avaframe/in3Utils/geoTrans.py

@@ -1128,6 +1128,9 @@ def path2domain(xyPath, rasterTransfo):
     DBYl = np.array((y + w * np.sin(d)))
     DBYr = np.array((y + w * np.sin(d + math.pi)))
 
+    # check if left or right domain boundary line is selfintersecting
+    checkDBOverlap(DBXl, DBXr, DBYl, DBYr)
+
     rasterTransfo["DBXl"] = DBXl
     rasterTransfo["DBXr"] = DBXr
     rasterTransfo["DBYl"] = DBYl
@@ -2094,3 +2097,31 @@ def cellsAffectedLine(header, pointsXY, typePoints):
                 mask[ind, ind2] = 1.0
 
     return mask, xx, yy
+
+
+def checkDBOverlap(DBXl, DBXr, DBYl, DBYr):
+    """check if lines spanned by DBX and DBY do intersect themselves; if selfintersecting line error
+
+    Parameters
+    -----------
+    DBXl, DBXr, DBYl, DBYr: numpy nd array
+        coordinates of lines
+
+    """
+
+    # create left and right domain boundar lineString
+    DBr = np.zeros((len(DBXr), 2))
+    DBr[:, 0] = DBXr
+    DBr[:, 1] = DBYr
+    DBrLine = shp.LineString(DBr)
+
+    DBl = np.zeros((len(DBXl), 2))
+    DBl[:, 0] = DBXl
+    DBl[:, 1] = DBYl
+    DBlLine = shp.LineString(DBl)
+
+    # check if either of the left or right domain boundary lineString is selfintersecting
+    if not DBrLine.is_simple or not DBlLine.is_simple:
+        message = "Domain transformation not applicable for given line - curvature of provided line would lead to folding"
+        log.error(message)
+        raise AssertionError(message)

avaframe/tests/test_aimecTools.py

@@ -314,3 +314,24 @@ def test_computeRunoutLine(tmp_path):
     assert ("sRunout" in runoutLine.keys()) is False
     assert ("lRunout" in runoutLine.keys()) is False
     assert ("xRunout" in runoutLine.keys()) is False
+
+
+def test_checkOverlapDBXY():
+    """check if lines along coordinate grid do intersect"""
+
+    # setup required input
+    x1 = np.arange(0, 10, 1)
+    y1 = np.arange(2, 10, 1)
+    X, Y = np.meshgrid(x1, y1)
+    rasterTransfo = {"gridx": X, "gridy": Y}
+
+    flagOverlap = aT.checkOverlapDBXY(rasterTransfo)
+
+    assert flagOverlap is True
+
+    X[:, 0] = np.arange(0, 8, 1)
+    rasterTransfo = {"gridx": X, "gridy": Y}
+
+    with pytest.raises(AssertionError) as e:
+        assert aT.checkOverlapDBXY(rasterTransfo)
+    assert "New coordinate system has overlaps - first intersection at POINT (1 3)." in str(e.value)