[ana3;in3] add check for aimec new coordinate system

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,11 @@ 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
+    intPoints = checkOverlapDBXY(
+        rasterTransfo,
+    )
+
     # add surface parallel coordinate (sParallel)
     rasterTransfo = addSurfaceParalleCoord(rasterTransfo)
 
@@ -760,8 +765,8 @@ def getSArea(rasterTransfo, dem):
     return rasterTransfo
 
 
-def transform(data, name, rasterTransfo, interpMethod):
-    """ Transfer data from old raster to new raster
+def transform(data, name, rasterTransfo, interpMethod, dem=False):
+    """Transfer data from old raster to new raster
 
     Assign value to the points of the new raster (after domain transormation)
 
@@ -775,6 +780,8 @@ def transform(data, name, rasterTransfo, interpMethod):
         transformation information
     interpMethod: str
         interpolation method to chose between 'nearest' and 'bilinear'
+    dem: bool
+        indicator if field to be transformed is a DEM
 
     Returns
     -------
@@ -798,6 +805,21 @@ def transform(data, name, rasterTransfo, interpMethod):
     log.debug('Data-file: %s - %d raster values transferred - %d out of original raster'
               'bounds!' % (name, iib-ioob, ioob))
 
+    # # TODO: only required if check of folding on coordinate system where data is to be analysed
+    # # check if data overlap in s, l coordinate system
+    # # first set no data (nan) values to zero for next check
+    # newDataInt = np.where(np.isnan(newData), 0, newData)
+    # newDataMask = np.where(rasterTransfo["intersectionPoints"], newDataInt, 0)
+    # if np.any(newDataMask) and not dem:
+    #     message = (
+    #         "New coordinate system has overlaps - "
+    #         "The provided path has too much curvature. "
+    #         "Try: (1) smoothing the path, (2) reducing the domain width, or (3) using fewer points."
+    #     )
+    #     outAimec.plotOverlap(rasterTransfo, newData, name, rasterTransfo["avaDir"])
+    #     log.error(message)
+    #     raise AssertionError(message)
+
     return newData
 
 
@@ -1837,3 +1859,130 @@ 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
+
+    Returns
+    ---------
+    flagOverlap: bool
+        False if no folding of coordinate system
+    """
+
+    x = rasterTransfo["gridx"]
+    y = rasterTransfo["gridy"]
+
+    flagOverlap = False
+
+    # 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.warning(message)
+                flagOverlap = True
+                break
+
+    return flagOverlap
+
+
+## CODE CURRENTLY NOT USED: MORE ADVANCED CHECKS
+def checkOverlapDBXYWithData(rasterTransfo, pointTolerance):
+    """check if x, y coordinates of new coordinate system overlap
+
+    TODO: this check is required if folding of cooridnate system only where there is data to be analysed shall we checked
+
+    Parameters
+    -----------
+    rasterTransfo: dict
+        dict with gridx, gridy locations of new coordinates in old coordinate system
+    pointTolerance: float
+        defines the absolute tolerance used to check whether coordinate locations are affected by the intersection
+        of coordinate grid lines (to check if overlap/folding in transformed coordinate system)
+
+    Returns
+    ---------
+    intPointsArray: numpy nd array
+        boolean array used as index array for coordinate arrays (rasterTransfo: gridx, gridy) where intersections occur
+        new coordinate system (s, l) has overlaps
+    """
+
+    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))
+
+    intPointsArray = np.zeros((x.shape[0], x.shape[1]))
+    # 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):
+                intPointsArray = findIntSectCoors(intersectionPoints, x, y, intPointsArray, pointTolerance)
+
+            elif isinstance(intersectionPoints, shp.MultiPoint):
+                for intPoint in intersectionPoints.geoms:
+                    intPointsArray = findIntSectCoors(intPoint, x, y, intPointsArray, pointTolerance)
+
+    # convert to boolean array for use as index array
+    intPointsArray = np.array(intPointsArray, dtype=bool)
+
+    return intPointsArray
+
+
+def findIntSectCoors(intersectionPoint, x, y, intPointsArray, pointTolerance):
+    """find in coordinate arrays (x, y) indices of points in intersectionPoints
+
+    Parameters
+    -----------
+    intersectionPoint: shapely Point
+        point where intersection is found
+    x, y: numpy nd array
+        x, y coordinates of s,l coordinate system
+    intPointsArray: numpy nd array
+        array with 1, 0 if point is intersecting
+
+    Returns
+    --------
+    intPointsArray: numpy nd array
+        updated array with 1, 0 if point is intersecting
+    """
+
+    # get x and y coordinate of intersection
+    y1 = intersectionPoint.y
+    x1 = intersectionPoint.x
+
+    # find corresponding coordinate of x, y arrays of coordinates in s, l system
+    yLoc = np.where(np.isclose(y, y1, atol=pointTolerance, rtol=0.0), y, np.nan)
+    xLoc = np.where(np.isclose(x, x1, atol=pointTolerance, rtol=0.0), x, np.nan)
+    indexArray = np.where(((np.isnan(yLoc) == False) & (np.isnan(xLoc) == False)))
+
+    intPointsArray[indexArray] = 1
+
+    return intPointsArray

avaframe/ana3AIMEC/ana3AIMEC.py

@@ -128,14 +128,15 @@ def mainAIMEC(pathDict, inputsDF, cfg):
     log.debug("Creating new deskewed raster and preparing new raster assignment function")
     log.debug('Analyze and make domain transformation on Data-file %s' % demSource)
     rasterTransfo = aimecTools.makeDomainTransfo(pathDict, dem, refHeader['cellsize'], cfgSetup)
+    rasterTransfo["avaDir"] = pathDict["avalancheDir"]
 
     # ####################################################
     # visualisation
     # TODO: needs to be moved somewhere else
     slRaster = aimecTools.transform(refRaster, refResultSource, rasterTransfo, interpMethod)
     newRasters = {}
     log.debug("Assigning dem data to deskewed raster")
-    newRasters['newRasterDEM'] = aimecTools.transform(dem, demSource, rasterTransfo, interpMethod)
+    newRasters["newRasterDEM"] = aimecTools.transform(dem, demSource, rasterTransfo, interpMethod, dem=True)
 
     inputData = {}
     inputData['slRaster'] = slRaster

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 for given path_aimec - curvature of provided line leads to folding"
+        log.warning(message)
+

avaframe/out3Plot/outAIMEC.py

@@ -20,6 +20,7 @@
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 from matplotlib.cm import ScalarMappable
 from matplotlib.colors import Normalize
+import pathlib
 
 # Local imports
 import avaframe.out3Plot.plotUtils as pU
@@ -2020,3 +2021,63 @@ def defineColorcodingValues(dataDF, paraVar, nSamples):
         values = None
 
     return colorFlag, colorSuffix, minVal, maxVal, values, cmapSCVals
+
+
+def plotOverlap(rasterTransfo, newData, name, avaDir):
+    """plot where s,l coordinate system has overlaps and where this coincides with data to be analysed
+
+    Parameters
+    -----------
+    rasterTransfo: dict
+        dictionary with all the information about the coordinate transformation
+    newData: numpy nd array
+        current data field to be analysed and already transformed in s,l coordinate system
+    """
+
+    l = rasterTransfo["l"]
+    s = rasterTransfo["s"]
+
+    cmap1, _, ticks1, norm1 = pU.makeColorMap(
+        pU.colorMaps["prob"],
+        np.nanmin(newData),
+        np.nanmax(newData),
+        continuous=pU.contCmap,
+    )
+
+    fig, ax = plt.subplots(nrows=3, ncols=1, figsize=(pU.figW * 2, pU.figH * 3))
+
+    ax[0].set_title("Data field")
+    ref0, im = pU.NonUnifIm(
+        ax[0],
+        s,
+        l,
+        np.transpose(np.where((newData > 0), newData, np.nan)),
+        "$S_{XY}$ (thalweg) [m]",
+        "$L_{XY}$ (thalweg) [m]",
+        extent=[s.min(), s.max(), l.min(), l.max()],
+    )
+    ax[1].set_title("Overlap in coordinate system")
+    ref1, im1 = pU.NonUnifIm(
+        ax[1],
+        s,
+        l,
+        np.transpose(np.where(rasterTransfo["intersectionPoints"], 1, np.nan)),
+        "$S_{XY}$ (thalweg) [m]",
+        "$L_{XY}$ (thalweg) [m]",
+        extent=[s.min(), s.max(), l.min(), l.max()],
+    )
+
+    ax[2].set_title("Overlap coinciding with analysed data")
+    fieldOverlap = np.where(rasterTransfo["intersectionPoints"], newData, np.nan)
+    ref2, im2 = pU.NonUnifIm(
+        ax[2],
+        s,
+        l,
+        np.transpose(np.where(fieldOverlap > 0, newData, np.nan)),
+        "$S_{XY}$ (thalweg) [m]",
+        "$L_{XY}$ (thalweg) [m]",
+        extent=[s.min(), s.max(), l.min(), l.max()],
+    )
+
+    outFileName = "ErrorPlot_domainOverlap%s" % (name.stem)
+    pU.saveAndOrPlot({"pathResult": pathlib.Path(avaDir, "Outputs", "ana3AIMEC")}, outFileName, fig)

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 not flagOverlap
+
+    X[:, 0] = np.arange(0, 8, 1)
+    rasterTransfo = {"gridx": X, "gridy": Y}
+
+    flagOverlap = aT.checkOverlapDBXY(rasterTransfo)
+
+    assert flagOverlap

avaframe/tests/test_ana3AIMEC.py

@@ -92,6 +92,7 @@ def test_analyzeArea(tmp_path):
     rasterTransfo["rasterArea"] = np.ones((500, 100))
     rasterTransfo["indStartOfRunout"] = 400
     rasterTransfo["startOfRunoutAreaAngle"] = 10
+    rasterTransfo["intersectionPoints"] = np.array(np.zeros((gridx.shape[0], gridy.shape[1])), dtype=bool)
     contourDict = {}
 
     # Analyze data

docs/moduleAna3AIMEC.rst

@@ -30,14 +30,23 @@ Inputs
 *  DEM  (digital elevation model) as raster file with either `ESRI grid format <https://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/esri-ascii-raster-format.htm>`_
    or GeoTIFF format. The format of the DEM determines
 which format is used for the output.
+
 .. Note:: The spatial resolution of the DEM and its extent can differ from the result raster data.
           Spatial resolution can also differ between simulations. If this is the case, the spatial
           resolution of the reference simulation results raster is used (default) or, if provided,
           the resolution specified in the configuration file (``cellSizeSL``) is used.
           This is done to ensure that all simulations will be transformed and analyzed using the
           same spatial resolution.
+
 *  avalanche thalweg in LINES (as a shapefile named ``NameOfAvalanche/Inputs/LINES/path_aimec.shp``), the line needs to cover the entire affected area but is not allowed
    to exceed the DEM extent
+
+.. Note:: If the thalweg is strongly curved this can lead to overlaps in the transformed coordinate system. If this affects
+          areas where there is data to be analysed, this will lead to wrong/distorted results in the analysis.
+          This is the case if any of the lines normal to the thalweg in the domain transformation figure,
+          e.g. :numref:`fig-aimec-domain-transfo`, cross. If so, computations will still be performed but a
+          a warning will be prompted in the log.
+
 *  results from avalanche simulation (when using results from com1DFA,
    the helper function :py:func:`ana3AIMEC.dfa2Aimec.mainDfa2Aimec` in
    :py:mod:`ana3AIMEC.dfa2Aimec` fetches and prepares the input for Aimec)