diff --git a/avaframe/ana3AIMEC/aimecTools.py b/avaframe/ana3AIMEC/aimecTools.py index 501b43d0d..9aac43451 100644 --- a/avaframe/ana3AIMEC/aimecTools.py +++ b/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 + _ = 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 diff --git a/avaframe/ana3AIMEC/ana3AIMEC.py b/avaframe/ana3AIMEC/ana3AIMEC.py index 383dd9bc3..8af04b01c 100644 --- a/avaframe/ana3AIMEC/ana3AIMEC.py +++ b/avaframe/ana3AIMEC/ana3AIMEC.py @@ -128,6 +128,7 @@ 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 @@ -135,7 +136,7 @@ def mainAIMEC(pathDict, inputsDF, cfg): 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 diff --git a/avaframe/data/avaGar/Inputs/LINES/path_aimec.dbf b/avaframe/data/avaGar/Inputs/LINES/path_aimec.dbf index e064fbbb0..875ae8e8d 100644 Binary files a/avaframe/data/avaGar/Inputs/LINES/path_aimec.dbf and b/avaframe/data/avaGar/Inputs/LINES/path_aimec.dbf differ diff --git a/avaframe/data/avaGar/Inputs/LINES/path_aimec.shp b/avaframe/data/avaGar/Inputs/LINES/path_aimec.shp index 00c38c040..9c46cd2f6 100644 Binary files a/avaframe/data/avaGar/Inputs/LINES/path_aimec.shp and b/avaframe/data/avaGar/Inputs/LINES/path_aimec.shp differ diff --git a/avaframe/data/avaMal/Inputs/LINES/path_aimec.dbf b/avaframe/data/avaMal/Inputs/LINES/path_aimec.dbf index 03360e0e4..46e1726f9 100644 Binary files a/avaframe/data/avaMal/Inputs/LINES/path_aimec.dbf and b/avaframe/data/avaMal/Inputs/LINES/path_aimec.dbf differ diff --git a/avaframe/data/avaMal/Inputs/LINES/path_aimec.shp b/avaframe/data/avaMal/Inputs/LINES/path_aimec.shp index 357a0844e..19ff4fa39 100644 Binary files a/avaframe/data/avaMal/Inputs/LINES/path_aimec.shp and b/avaframe/data/avaMal/Inputs/LINES/path_aimec.shp differ diff --git a/avaframe/data/avaMal/Inputs/LINES/path_aimec.shx b/avaframe/data/avaMal/Inputs/LINES/path_aimec.shx index 378f16046..efb5fc666 100644 Binary files a/avaframe/data/avaMal/Inputs/LINES/path_aimec.shx and b/avaframe/data/avaMal/Inputs/LINES/path_aimec.shx differ diff --git a/avaframe/in3Utils/geoTrans.py b/avaframe/in3Utils/geoTrans.py index 81f221dfe..eb415640e 100644 --- a/avaframe/in3Utils/geoTrans.py +++ b/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) + diff --git a/avaframe/out3Plot/outAIMEC.py b/avaframe/out3Plot/outAIMEC.py index c08fb7d26..d26b33f47 100644 --- a/avaframe/out3Plot/outAIMEC.py +++ b/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) diff --git a/avaframe/tests/test_aimecTools.py b/avaframe/tests/test_aimecTools.py index a6f6819d0..58be31f68 100644 --- a/avaframe/tests/test_aimecTools.py +++ b/avaframe/tests/test_aimecTools.py @@ -314,3 +314,103 @@ 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 + + +def test_checkOverlapDBXYWithData(): + """test checkOverlapDBXYWithData - detect intersection points in coordinate grid""" + from shapely import geometry as shp + + # Test case 1: No intersection - parallel lines + x1 = np.arange(0, 5, 1) + y1 = np.arange(0, 4, 1) + X, Y = np.meshgrid(x1, y1) + rasterTransfo = {"gridx": X, "gridy": Y} + pointTolerance = 0.01 + + intPointsArray = aT.checkOverlapDBXYWithData(rasterTransfo, pointTolerance) + + # Verify output is boolean array with correct shape + assert intPointsArray.dtype == bool + assert intPointsArray.shape == X.shape + # No intersections should be found + assert np.sum(intPointsArray) == 0 + + # Test case 2: Lines with intersection + # Create a grid where two columns cross each other + # Column 0: vertical line at x=0 + # Column 1: diagonal line from (1,0) through (0,2) to (-1,4) + # These should intersect at approximately (0.5, 1) + X = np.array([[0, 1, 0.5, 2], [0, 0.5, 1, 2], [0, 0, 1.5, 2], [0, -0.5, 2, 2], [0, -1, 2.5, 2]]) + Y = np.array([[0, 0, 0, 0], [1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]]) + rasterTransfo = {"gridx": X, "gridy": Y} + + intPointsArray = aT.checkOverlapDBXYWithData(rasterTransfo, pointTolerance) + + # Verify output is boolean array with correct shape + assert intPointsArray.dtype == bool + assert intPointsArray.shape == X.shape + # With crossing lines, intersections may or may not be found depending on geometry + # The function should at least run without errors + assert isinstance(intPointsArray, np.ndarray) + + +def test_findIntSectCoors(): + """test findIntSectCoors - find indices of intersection points in coordinate arrays""" + from shapely import geometry as shp + + # Setup test data + x = np.array([[0, 1, 2], [0, 1, 2], [0, 1, 2]]) + y = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]]) + intPointsArray = np.zeros((3, 3)) + pointTolerance = 0.01 + + # Test case 1: Single point intersection at (1, 1) + intersectionPoint = shp.Point(1.0, 1.0) + + intPointsArray = aT.findIntSectCoors(intersectionPoint, x, y, intPointsArray, pointTolerance) + + # Verify that the point at (1,1) is marked + assert intPointsArray[1, 1] == 1 + # Verify only one point is marked + assert np.sum(intPointsArray) == 1 + + # Test case 2: Point with tolerance + intPointsArray = np.zeros((3, 3)) + # Point slightly off from (2, 2) but within tolerance + intersectionPoint = shp.Point(2.005, 2.005) + + intPointsArray = aT.findIntSectCoors(intersectionPoint, x, y, intPointsArray, pointTolerance) + + # Verify that the point at (2,2) is marked despite slight offset + assert intPointsArray[2, 2] == 1 + assert np.sum(intPointsArray) == 1 + + # Test case 3: Point outside grid + intPointsArray = np.zeros((3, 3)) + intersectionPoint = shp.Point(5.0, 5.0) + + intPointsArray = aT.findIntSectCoors(intersectionPoint, x, y, intPointsArray, pointTolerance) + + # Verify no points are marked + assert np.sum(intPointsArray) == 0 diff --git a/avaframe/tests/test_ana3AIMEC.py b/avaframe/tests/test_ana3AIMEC.py index a727ae0f2..2ac45ef43 100644 --- a/avaframe/tests/test_ana3AIMEC.py +++ b/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 diff --git a/docs/moduleAna3AIMEC.rst b/docs/moduleAna3AIMEC.rst index 182582c6a..f90fce799 100644 --- a/docs/moduleAna3AIMEC.rst +++ b/docs/moduleAna3AIMEC.rst @@ -30,14 +30,23 @@ Inputs * DEM (digital elevation model) as raster file with either `ESRI grid format `_ 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)