[com6] Scarp: FinalUpdate

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,17 +242,22 @@ 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])
         ySeed.append(planes[5 * i + 1])
         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]

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
 

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 = []

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