[com6] Scarp: FinalUpdate

avaframe/com6RockAvalanche/scarp.py

@@ -83,10 +83,11 @@ 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']))
+            planesDip = list(map(float, SHPdata['dipdir']))
+            planesSlope = 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"]):
             raise ValueError("Mismatch between number of features and extracted plane attributes in the shapefile.")
@@ -112,12 +113,12 @@ 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']))
+            ellipsoidsTilt = list(map(float, SHPdata['dipAngle']))
+            ellipsoidsDir = list(map(float, SHPdata['dipdir']))
             ellipsoidsOffset = list(map(float, SHPdata['offset']))
-            ellipsoidDip = list(map(float, SHPdata['dip']))
+            ellipsoidDip = 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]):
             raise ValueError("Mismatch between number of shapefile features and ellipsoid parameters.")
@@ -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