diff --git a/avaframe/com6RockAvalanche/scarp.py b/avaframe/com6RockAvalanche/scarp.py index 612d56380..804dc3602 100644 --- a/avaframe/com6RockAvalanche/scarp.py +++ b/avaframe/com6RockAvalanche/scarp.py @@ -112,28 +112,30 @@ def scarpAnalysisMain(cfg, baseDir): ellipsoidsMaxDepth = list(map(float, SHPdata['maxdepth'])) ellipsoidsSemiMajor = list(map(float, SHPdata['semimajor'])) ellipsoidsSemiMinor = list(map(float, SHPdata['semiminor'])) + ellipsoidsTilt = list(map(float, SHPdata['tilt'])) + ellipsoidsDir = list(map(float, SHPdata['direc'])) + ellipsoidsOffset = list(map(float, SHPdata['offset'])) + ellipsoidDip = list(map(float, SHPdata['dip'])) except KeyError as e: - raise ValueError(f"Required attribute '{e.args[0]}' not found in shapefile. Ensure the fields 'maxdepth', 'semimajor', and 'semiminor' exist.") - - if not ( - len(ellipsoidsMaxDepth) - == len(ellipsoidsSemiMajor) - == len(ellipsoidsSemiMinor) - == SHPdata["nFeatures"] - ): - raise ValueError( - "Mismatch between number of shapefile features and ellipsoid parameters in the .ini file." - ) - + raise ValueError(f"Required attribute '{e.args[0]}' not found in shapefile. Ensure the fields 'maxdepth', 'semimajor', 'semiminor', 'tilt', 'dir', 'dip', and 'offset' exist.") + + if not all(len(lst) == SHPdata["nFeatures"] for lst in [ellipsoidsMaxDepth, ellipsoidsSemiMajor, ellipsoidsSemiMinor, ellipsoidsTilt, ellipsoidsDir, ellipsoidsOffset, ellipsoidDip]): + raise ValueError("Mismatch between number of shapefile features and ellipsoid parameters.") + for i in range(SHPdata["nFeatures"]): xCenter = SHPdata["x"][int(SHPdata["Start"][i])] yCenter = SHPdata["y"][int(SHPdata["Start"][i])] maxDepth = ellipsoidsMaxDepth[i] semiMajor = ellipsoidsSemiMajor[i] semiMinor = ellipsoidsSemiMinor[i] - ellipsoidFeatures.extend([xCenter, yCenter, maxDepth, semiMajor, semiMinor]) - + tilt = ellipsoidsTilt[i] + direction = ellipsoidsDir[i] + offset = ellipsoidsOffset[i] + dip = ellipsoidDip[i] + ellipsoidFeatures.extend([xCenter, yCenter, maxDepth, semiMajor, semiMinor, tilt, direction, offset, dip]) + features = ",".join(map(str, ellipsoidFeatures)) + log.debug("Ellipsoid features extracted and combined: %s", features) if method == 'plane': @@ -262,55 +264,104 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes): return scarpData def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids): - """Calculate the scarp using sliding circles (rotational ellipsoids). + """Calculate the scarp using tilted and offset ellipsoids. Parameters ---------- elevData : np.ndarray The elevation data as a 2D array. periData : np.ndarray - The perimeter data as a 2D array, which defines the extent of the scarp. + The perimeter data as a 2D array. elevTransform : Affine - The affine transformation matrix of the raster (used to convert pixel to geographic coordinates). + The affine transformation matrix of the raster. ellipsoids : str - Comma-separated string defining ellipsoids (x_center, y_center, max_depth, semi_major_axis, semi_minor_axis). + Comma-separated string defining ellipsoids with parameters: + (x_center, y_center, max_depth, semi_major, semi_minor, tilt, dir, offset) Returns ------- np.ndarray A 2D array with the calculated scarp values. """ + n, m = elevData.shape scarpData = np.zeros_like(elevData, dtype=np.float32) + # Parse ellipsoid definitions ellipsoids = list(map(float, ellipsoids.split(','))) - nEllipsoids = int(len(ellipsoids) / 5) - - xCenter = [ellipsoids[0]] - yCenter = [ellipsoids[1]] - maxDepth = [ellipsoids[2]] - semiMajor = [ellipsoids[3]] - semiMinor = [ellipsoids[4]] - - for i in range(1, nEllipsoids): - xCenter.append(ellipsoids[5 * i]) - yCenter.append(ellipsoids[5 * i + 1]) - maxDepth.append(ellipsoids[5 * i + 2]) - semiMajor.append(ellipsoids[5 * i + 3]) - semiMinor.append(ellipsoids[5 * i + 4]) + nEllipsoids = int(len(ellipsoids) / 9) + + xCenter, yCenter, maxDepth = [], [], [] + semiMajor, semiMinor = [], [] + tilt, tiltDir, offset, dip = [], [], [], [] + + for i in range(nEllipsoids): + xCenter.append(ellipsoids[9 * i]) + yCenter.append(ellipsoids[9 * i + 1]) + maxDepth.append(ellipsoids[9 * i + 2]) + semiMajor.append(ellipsoids[9 * i + 3]) + semiMinor.append(ellipsoids[9 * i + 4]) + tilt.append(ellipsoids[9 * i + 5]) + tiltDir.append(np.radians(ellipsoids[9 * i + 6])) # tilt direction in radians + offset.append(ellipsoids[9 * i + 7]) + dip.append(np.radians(ellipsoids[9 * i + 8])) # rotation of base ellipse in radians + + # Compute slope direction and magnitude from DEM gradients + grad_y, grad_x = np.gradient(elevData, abs(elevTransform[4]), elevTransform[0]) + slopeDir = np.arctan2(-grad_y, grad_x) # Azimuth of steepest descent + slopeMagnitude = np.sqrt(grad_x**2 + grad_y**2) for row in range(n): for col in range(m): west, north = rasterio.transform.xy(elevTransform, row, col, offset='center') scarpVal = elevData[row, col] + for k in range(nEllipsoids): - x2 = ((west - xCenter[k]) / semiMajor[k]) ** 2 - y2 = ((north - yCenter[k]) / semiMinor[k]) ** 2 - if x2 + y2 < 1: - scarpVal = min(scarpVal, elevData[row, col] - maxDepth[k] * (1 - (x2 + y2))) - if periData[row, col] > 0: - scarpData[row, col] = scarpVal - else: - scarpData[row, col] = elevData[row, col] + center_col, center_row = ~elevTransform * (xCenter[k], yCenter[k]) + center_row = int(round(center_row)) + center_col = int(round(center_col)) + + if 0 <= center_row < n and 0 <= center_col < m: + slopeAzimuth = slopeDir[center_row, center_col] + slopeAzimuth = (slopeAzimuth + 2 * np.pi) % (2 * np.pi) + slopeMag = slopeMagnitude[center_row, center_col] + + if not np.isnan(slopeAzimuth) and slopeMag > 0: + normal_dx = np.cos(slopeAzimuth) + normal_dy = np.sin(slopeAzimuth) + slopeAngle = math.atan(slopeMag) + dxOffset = -offset[k] * normal_dx * math.sin(slopeAngle) + dyOffset = -offset[k] * normal_dy * math.sin(slopeAngle) + dzOffset = -offset[k] * math.cos(slopeAngle) + else: + dxOffset = dyOffset = dzOffset = 0 + else: + dxOffset = dyOffset = dzOffset = 0 + + x0 = xCenter[k] + dxOffset + y0 = yCenter[k] + dyOffset + z0 = dzOffset + + dxPos = west - x0 + dyPos = north - y0 + + # Rotate the position by dip angle + dxRot = dxPos * np.cos(dip[k]) + dyPos * np.sin(dip[k]) + dyRot = -dxPos * np.sin(dip[k]) + dyPos * np.cos(dip[k]) + + # Normalize to ellipsoid axes + xNorm = dxRot / semiMajor[k] + yNorm = dyRot / semiMinor[k] + distance = xNorm**2 + yNorm**2 + + if distance <= 1: + baseDepth = maxDepth[k] * (1 - distance) + tiltEffect = math.tan(math.radians(tilt[k])) * ( + dxRot * np.cos(tiltDir[k]) + dyRot * np.sin(tiltDir[k]) + ) + totalDepth = baseDepth + tiltEffect + z0 + scarpVal = min(scarpVal, elevData[row, col] - totalDepth) + + scarpData[row, col] = scarpVal if periData[row, col] > 0 else elevData[row, col] return scarpData diff --git a/avaframe/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index b18f39d5f..b6493bfc1 100644 --- a/avaframe/in2Trans/shpConversion.py +++ b/avaframe/in2Trans/shpConversion.py @@ -88,6 +88,9 @@ def SHP2Array(infile, defname=None): semiminor_value = None maxdepth_value = None semimajor_value = None + tilt_value = None + direc_value = None + offset_value = None # get coordinate system @@ -111,6 +114,10 @@ def SHP2Array(infile, defname=None): semiminorList = [] semimajorList = [] maxdepthList = [] + tiltList = [] + direcList = [] + offsetList = [] + Length = np.empty((0)) Start = np.empty((0)) Coordx = np.empty((0)) @@ -163,6 +170,12 @@ def SHP2Array(infile, defname=None): maxdepth_value = value if name == "semimajor": semimajor_value = value + if name == "tilt": + tilt_value = value + if name == "direc": + direc_value = value + if name == "offset": + offset_value = value if name == "rho": rho = value if name == "sks": @@ -194,6 +207,9 @@ def SHP2Array(infile, defname=None): semiminorList.append(semiminor_value) semimajorList.append(semimajor_value) maxdepthList.append(maxdepth_value) + tiltList.append(tilt_value) + direcList.append(direc_value) + offsetList.append(offset_value) Start = np.append(Start, start) length = len(pts) @@ -224,6 +240,9 @@ def SHP2Array(infile, defname=None): SHPdata["maxdepth"] = maxdepthList SHPdata["semimajor"] = semimajorList SHPdata["semiminor"] = semiminorList + SHPdata["tilt"] = tiltList + SHPdata["direc"] = direcList + SHPdata["offset"] = offsetList sf.close() diff --git a/docs/moduleCom6RockAvalanche.rst b/docs/moduleCom6RockAvalanche.rst index fe75d019f..087b19159 100644 --- a/docs/moduleCom6RockAvalanche.rst +++ b/docs/moduleCom6RockAvalanche.rst @@ -46,12 +46,27 @@ Input * geometries: a shapefile containing point geometries. These points represent the centers of the ellipsoids or planes. The coordinates (x,y) of these points are used. If the plane method is used, the shape file must contain the Attributes "zseed", "dip" and "slopeangle" as float values. If the ellipsoid method is used, the shape file must - contain the Attributes "maxdepth", "semimajor" and "semiminor". The file must be located in avalancheDir/Inputs/ - POINTS and file name must end with “_coordinates”. If you are using the QGis Connector, the naming and location of the - file is not relevant. + contain the attributes "maxdepth", "semimajor", "semiminor", "tilt", "direc", "dip", "offset" (see below). + The file must be located in avalancheDir/Inputs/POINTS and file name must end with “_coordinates”. + If you are using the QGis Connector, the naming and location of the file is not relevant. + * perimeter: A shapefile that specifies a boundary area. Must be located in avalancheDir/Inputs/POLYGONS and file name must end with “_perimeter”. If you are using the QGis Connector, the naming and location of the file is not relevant. +**Attribute meanings:** + +* zseed: defines z coordinate of plane Center (m) +* dip: direction in which the plane/slope is facing (degree) +* slopeangle: steepness/angle of the slope (degree) + +* maxdepth: maximum depth of the ellipsoid (m) +* semimajor: length of the major axis (m) +* semiminor: length of the minor axis (m) +* tilt: steepness/angle of the slope (degree) +* direc: direction in which the slope is facing (degree) +* dip: direction in which the ellipsoid is facing (degree) +* offset: offset, normal to the DEM slope (m) + Output ~~~~~~ diff --git a/pyproject.toml b/pyproject.toml index 02b82f993..a4210dcb9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,7 +9,7 @@ name = "avaframe" description = "The Open Avalanche Framework" readme = "README.md" dynamic = ["version"] -license = { text = "EUPL" } +license = "EUPL-1.2" authors = [ { name = "AvaFrame Contributors", email = "felix@avaframe.org" } ] @@ -17,7 +17,6 @@ urls = { Homepage = "http://avaframe.org" } classifiers = [ "Development Status :: 5 - Production/Stable", "Intended Audience :: Science/Research", - "License :: OSI Approved :: European Union Public Licence 1.2 (EUPL 1.2)", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", @@ -99,7 +98,8 @@ avaframe = { path = "./", editable = true } pixi-pycharm = "*" pytest = "*" pytest-cov = "*" -pytest-mock ="*" +pytest-mock = "*" +black = "*" #Feature doc [tool.pixi.feature.doc.dependencies]