diff --git a/avaframe/com6RockAvalanche/scarp.py b/avaframe/com6RockAvalanche/scarp.py index 10eb5d0f0..d3ae2f762 100644 --- a/avaframe/com6RockAvalanche/scarp.py +++ b/avaframe/com6RockAvalanche/scarp.py @@ -83,15 +83,16 @@ def scarpAnalysisMain(cfg, baseDir): # Read required attributes directly from the shapefile's attribute table try: planesZseed = list(map(float, SHPdata['zseed'])) - planesDip = list(map(float, SHPdata['dip'])) - planesSlope = list(map(float, SHPdata['slopeangle'])) + planesDipDir = list(map(float, SHPdata['dipdir'])) + planesDipAngle = list(map(float, SHPdata['dipAngle'])) + except KeyError as e: - raise ValueError(f"Required attribute '{e.args[0]}' not found in shapefile. Make sure 'zseed', 'dip', and 'slope' fields exist.") + raise ValueError(f"Required attribute '{e.args[0]}' not found in shapefile. Make sure 'zseed', 'dipdir', and 'dipangle' fields exist.") - if not (len(planesZseed) == len(planesDip) == len(planesSlope) == SHPdata["nFeatures"]): + if not (len(planesZseed) == len(planesDipDir) == len(planesDipAngle) == SHPdata["nFeatures"]): raise ValueError("Mismatch between number of features and extracted plane attributes in the shapefile.") - if not (len(planesZseed) == len(planesDip) == len(planesSlope) == SHPdata["nFeatures"]): + if not (len(planesZseed) == len(planesDipDir) == len(planesDipAngle) == SHPdata["nFeatures"]): raise ValueError( "Mismatch between number of shapefile features and plane parameters in the .ini file." ) @@ -100,8 +101,8 @@ def scarpAnalysisMain(cfg, baseDir): xSeed = SHPdata["x"][int(SHPdata["Start"][i])] ySeed = SHPdata["y"][int(SHPdata["Start"][i])] zSeed = planesZseed[i] - dip = planesDip[i] - slopeAngle = planesSlope[i] + dip = planesDipDir[i] + slopeAngle = planesDipAngle[i] planeFeatures.extend([xSeed, ySeed, zSeed, dip, slopeAngle]) features = ",".join(map(str, planeFeatures)) @@ -112,14 +113,14 @@ 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'])) + ellipsoidsDipAngle = list(map(float, SHPdata['dipAngle'])) + ellipsoidsDipDir = list(map(float, SHPdata['dipdir'])) ellipsoidsOffset = list(map(float, SHPdata['offset'])) - ellipsoidDip = list(map(float, SHPdata['dip'])) + ellipsoidsRotAngle = list(map(float, SHPdata['rotAngle'])) except KeyError as e: - raise ValueError(f"Required attribute '{e.args[0]}' not found in shapefile. Ensure the fields 'maxdepth', 'semimajor', 'semiminor', 'tilt', 'dir', 'dip', and 'offset' exist.") + raise ValueError(f"Required attribute '{e.args[0]}' not found in shapefile. Ensure the fields 'maxdepth', 'semimajor', 'semiminor', 'rotangle', 'dipdir', 'dipangle', and 'offset' exist.") - if not all(len(lst) == SHPdata["nFeatures"] for lst in [ellipsoidsMaxDepth, ellipsoidsSemiMajor, ellipsoidsSemiMinor, ellipsoidsTilt, ellipsoidsDir, ellipsoidsOffset, ellipsoidDip]): + if not all(len(lst) == SHPdata["nFeatures"] for lst in [ellipsoidsMaxDepth, ellipsoidsSemiMajor, ellipsoidsSemiMinor, ellipsoidsDipAngle, ellipsoidsDipDir, ellipsoidsOffset, ellipsoidsRotAngle]): raise ValueError("Mismatch between number of shapefile features and ellipsoid parameters.") for i in range(SHPdata["nFeatures"]): @@ -128,10 +129,10 @@ def scarpAnalysisMain(cfg, baseDir): maxDepth = ellipsoidsMaxDepth[i] semiMajor = ellipsoidsSemiMajor[i] semiMinor = ellipsoidsSemiMinor[i] - tilt = ellipsoidsTilt[i] - direction = ellipsoidsDir[i] + tilt = ellipsoidsDipAngle[i] + direction = ellipsoidsDipDir[i] offset = ellipsoidsOffset[i] - dip = ellipsoidDip[i] + dip = ellipsoidsRotAngle[i] ellipsoidFeatures.extend([xCenter, yCenter, maxDepth, semiMajor, semiMinor, tilt, direction, offset, dip]) features = ",".join(map(str, ellipsoidFeatures)) @@ -150,6 +151,11 @@ def scarpAnalysisMain(cfg, baseDir): raise ValueError("Unsupported method. Choose 'plane' or 'ellipsoid'.") hRelData = dem["rasterData"] - scarpData + + #Compute and log excavated volume + cellArea = abs(dem["header"]["cellsize"] ** 2) + volume = np.sum(hRelData[periData > 0]) * cellArea + log.info(f"Excavated volume (within perimeter): {volume:.2f} m³") # create output directory and files outDir = pathlib.Path(baseDir) @@ -236,8 +242,10 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes): dip = [planes[3]] slope = [planes[4]] - betaX = [math.tan(math.radians(slope[0])) * math.cos(math.radians(dip[0]))] - betaY = [math.tan(math.radians(slope[0])) * math.sin(math.radians(dip[0]))] + slopeRad = math.radians(slope[0]) + dipRad = math.radians(dip[0]) + betaX = [ math.tan(slopeRad) * math.sin(dipRad) ] + betaY = [ math.tan(slopeRad) * math.cos(dipRad) ] for i in range(1, nPlanes): xSeed.append(planes[5 * i]) @@ -245,8 +253,11 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes): zSeed.append(planes[5 * i + 2]) dip.append(planes[5 * i + 3]) slope.append(planes[5 * i + 4]) - betaX.append(math.tan(math.radians(slope[i])) * math.cos(math.radians(dip[i]))) - betaY.append(math.tan(math.radians(slope[i])) * math.sin(math.radians(dip[i]))) + + slopeRad = math.radians(slope[i]) + dipRad = math.radians(dip[i]) + betaX.append( math.tan(slopeRad) * math.sin(dipRad) ) + betaY.append( math.tan(slopeRad) * math.cos(dipRad) ) for row in range(n): for col in range(m): @@ -333,6 +344,8 @@ def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids): dxOffset = -offset[k] * normal_dx * math.sin(slopeAngle) dyOffset = -offset[k] * normal_dy * math.sin(slopeAngle) dzOffset = -offset[k] * math.cos(slopeAngle) + #clamp z-offset + dzOffset = max(-maxDepth[k], min(dzOffset, maxDepth[k])) else: dxOffset = dyOffset = dzOffset = 0 else: @@ -344,11 +357,9 @@ def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids): 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]) - + + dxRot = dxPos * np.sin(dip[k]) - dyPos * np.cos(dip[k]) + dyRot = dxPos * np.cos(dip[k]) + dyPos * np.sin(dip[k]) # Normalize to ellipsoid axes xNorm = dxRot / semiMajor[k] yNorm = dyRot / semiMinor[k] diff --git a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.dbf b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.dbf index ef0170f04..b0e1d2757 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.dbf and b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.dbf differ diff --git a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shp b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shp index 0cc7e0016..a1cb3ab00 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shp and b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shp differ diff --git a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shx b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shx index adec32cf5..9f4f0f8bf 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shx and b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shx differ diff --git a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.dbf b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.dbf index 742c45010..46aebd550 100644 Binary files a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.dbf and b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.dbf differ diff --git a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shp b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shp index ec5f2194b..a07334ed1 100644 Binary files a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shp and b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shp differ diff --git a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shx b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shx index c2859ca3a..c11c1b409 100644 Binary files a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shx and b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shx differ diff --git a/avaframe/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index 620d81ced..6eef2d743 100644 --- a/avaframe/in2Trans/shpConversion.py +++ b/avaframe/in2Trans/shpConversion.py @@ -59,9 +59,9 @@ def SHP2Array(infile, defname=None): number of features per line (parts) zseed np array with the height of each scarp plane-feature (as many values as features) - slopeAngle - np array with the slope angle of each scarp plane-feature (as many values as features) - dip + dipDir + np array with the dip direction of each scarp plane-feature (as many values as features) + dipAngle np array with the dip angle of each scarp plane-feature (as many values as features) semiminor np array with the semi-minor axis of each scarp ellipsoid-feature (as many values as features) @@ -84,12 +84,12 @@ def SHP2Array(infile, defname=None): ci95 = None layerN = None zseed_value = None - dip_value = None - slopeAngle_value = None + dipdir_value = None + dipAngle_value = None semiminor_value = None maxdepth_value = None semimajor_value = None - tilt_value = None + rotAngle_value = None direc_value = None offset_value = None @@ -108,13 +108,13 @@ def SHP2Array(infile, defname=None): ci95List = [] layerNameList = [] zseedList = [] - dipList = [] + dipdirList = [] slopeList = [] - slopeAngleList = [] + dipAngleList = [] semiminorList = [] semimajorList = [] maxdepthList = [] - tiltList = [] + rotAngleList = [] direcList = [] offsetList = [] @@ -157,21 +157,21 @@ def SHP2Array(infile, defname=None): if name == "slope": # for dams slope = value - if name == "slopeangle": + if name == "dipangle": # for dams - slopeAngle_value = value + dipAngle_value = value if name == "zseed": zseed_value = value - if name == "dip": - dip_value = value + if name == "dipdir": + dipdir_value = value if name == "semiminor": semiminor_value = value if name == "maxdepth": maxdepth_value = value if name == "semimajor": semimajor_value = value - if name == "tilt": - tilt_value = value + if name == "rotangle": + rotAngle_value = value if name == "direc": direc_value = value if name == "offset": @@ -201,13 +201,13 @@ def SHP2Array(infile, defname=None): layerNameList.append(layerN) idList.append(str(rec.oid)) zseedList.append(zseed_value) - dipList.append(dip_value) + dipdirList.append(dipdir_value) slopeList.append(slope) - slopeAngleList.append(slopeAngle_value) + dipAngleList.append(dipAngle_value) semiminorList.append(semiminor_value) semimajorList.append(semimajor_value) maxdepthList.append(maxdepth_value) - tiltList.append(tilt_value) + rotAngleList.append(rotAngle_value) direcList.append(direc_value) offsetList.append(offset_value) @@ -234,13 +234,13 @@ def SHP2Array(infile, defname=None): SHPdata["layerName"] = layerNameList SHPdata["nParts"] = nParts SHPdata["nFeatures"] = len(Start) - SHPdata["slopeangle"] = slopeAngleList + SHPdata["dipAngle"] = dipAngleList SHPdata["zseed"] = zseedList - SHPdata["dip"] = dipList + SHPdata["dipdir"] = dipdirList SHPdata["maxdepth"] = maxdepthList SHPdata["semimajor"] = semimajorList SHPdata["semiminor"] = semiminorList - SHPdata["tilt"] = tiltList + SHPdata["rotAngle"] = rotAngleList SHPdata["direc"] = direcList SHPdata["offset"] = offsetList diff --git a/avaframe/tests/test_scarp.py b/avaframe/tests/test_scarp.py index cb1b89071..1692b34a3 100644 --- a/avaframe/tests/test_scarp.py +++ b/avaframe/tests/test_scarp.py @@ -121,8 +121,8 @@ def test_plane_parameter_extraction(scarp_test_data): # Extract plane parameters (as done in scarpAnalysisMain) planesZseed = list(map(float, SHPdata["zseed"])) - planesDip = list(map(float, SHPdata["dip"])) - planesSlope = list(map(float, SHPdata["slopeangle"])) + planesDip = list(map(float, SHPdata["dipdir"])) + planesSlope = list(map(float, SHPdata["dipAngle"])) # Assertions assert len(planesZseed) == SHPdata["nFeatures"], ( @@ -132,7 +132,7 @@ def test_plane_parameter_extraction(scarp_test_data): assert len(planesSlope) == SHPdata["nFeatures"], ( "Should have slope for each feature" ) - assert SHPdata["nFeatures"] == 2, "Test data should have 2 features" + assert SHPdata["nFeatures"] == 1, "Test data should have 1 feature" # Build feature string planeFeatures = [] diff --git a/docs/moduleCom6RockAvalanche.rst b/docs/moduleCom6RockAvalanche.rst index 087b19159..f9b76e5ba 100644 --- a/docs/moduleCom6RockAvalanche.rst +++ b/docs/moduleCom6RockAvalanche.rst @@ -45,8 +45,8 @@ Input * elevation: DEM (ASCII), which serves as the basis for calculating the scarps. Must be in avalancheDir/Inputs. * 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", "semiminor", "tilt", "direc", "dip", "offset" (see below). + attributes "zseed", "dipdir" and "dipAngle" as float values. If the ellipsoid method is used, the shape file must + contain the attributes "maxdepth", "semimajor", "semiminor", "dipAngle", "dipdir", "rotAngle", "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. @@ -55,16 +55,16 @@ Input **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) +* zseed: defines z coordinate of plane center (m) +* dipdir: direction in which the plane/slope is facing (degree) +* dipAngle: 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) +* dipAngle: steepness/angle of the ellipsoid tilt (degree) +* dipdir: direction in which the ellipsoid slope is facing (degree) +* rotAngle: rotation angle of the ellipsoid base (degree) * offset: offset, normal to the DEM slope (m) Output