Scarp Update

avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini

@@ -19,8 +19,8 @@ defaultConfig = True
 # desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, particles) - separated by |
 resType = pft|pfv|ppr|FT
 
-#+++++++++SNOW properties
-# density of snow [kg/m³]
+#+++++++++ROCKMASS properties
+# density of rockmass (mobilised rock-avalanche) [kg/m³]
 rho = 2500
 
 #+++++++++++++SPH parameters
@@ -51,13 +51,13 @@ xsivoellmy = 700.
 # expected mesh size [m]
 meshCellSize = 5
 
-# density of entrained snow [kg/m³]
-rhoEnt = 100
+# density of entrained material (soil/soft rock) [kg/m³]
+rhoEnt = 2000
 #+++++Entrainment thickness++++
 # True if entrainment thickness should be read from file (shapefile, raster); if False - entTh read from ini file (only available for shapefile)
 entThFromFile = True
 # if a thickness value is missing for the entrainment feature in the provided shp file this value is used for all features [m]
-entThIfMissingInShp = 0.3
+entThIfMissingInShp = 0
 # VARIATION options only if ENT file is a shapefile --------
 # if a variation on entTh shall be performed add here +- percent and number of steps separated by $
 # for example entThPercentVariation=50$10 [%]
@@ -77,35 +77,8 @@ entThDistVariation =
 entTh =
 
 
-#++++++++++++ Entrainment Erosion Energy
+#++++++++++++ Entrainment Erosion Energy [J/m²]
 # Used to determine speed loss via energy loss due to entrained mass
 entEroEnergy = 5000
 entShearResistance = 0
 entDefResistance = 0
-
-#++++++++++++ Resistance model
-# default setup:
-ResistanceModel = default
-# At each time step, ResistanceModel default applies increased friction and optional detrainment (see below). Only relevant in resistance areas.
-# NOTE: development setup; parameter values need more testing and calibration!!
-
-#  parameter for increased friction in resistance areas
-cResH = 0.01
-
-# Apply detrainment in resistance areas in accordance with the flow thickness and flow velocity thresholds specified below.
-# if False - only increased friction is applied in resistance areas
-detrainment = True
-
-# detrainment parameter defined by Feistl et al. (2014)
-detK = 5
-
-# thresholds if detrainment is set to True
-# FV OR FT below min thresholds: apply only detrainment. no increased friction
-# FV AND FT within min and max thresholds: no detrainment, only apply increased friction
-# FV OR FT above max thresholds: no detrainment and no increased friction applied
-forestVMin = 6.
-forestThMin = 0.6
-forestVMax = 40.
-forestThMax = 10.
-
-

avaframe/com6RockAvalanche/scarp.py

@@ -14,6 +14,70 @@
 logging.basicConfig(level=logging.DEBUG)
 log = logging.getLogger(__name__)
 
+def _extract_numeric_attributes(shpData, requiredAttributes):
+    """Validate and extract required numeric shapefile attributes.
+
+    Parameters
+    ----------
+    shpData : dict
+        Shapefile data returned by SHP2Array.
+    requiredAttributes : list of str
+        Attribute names that must be present and contain numeric values.
+
+    Returns
+    -------
+    dict
+        Dictionary mapping each required attribute name to a list of floats.
+
+    Raises
+    ------
+    ValueError
+        If an attribute is missing, set to None, or contains non-numeric values.
+    """
+    missingAttributes = [
+        attr for attr in requiredAttributes
+        if attr not in shpData
+        or shpData.get(attr) is None
+        or any(value is None for value in shpData[attr])
+    ]
+
+    if missingAttributes:
+        raise ValueError(
+            f"Required attribute(s) {missingAttributes} not found in shapefile. "
+            f"Required fields are: {requiredAttributes}."
+        )
+
+    try:
+        return {
+            attr: list(map(float, shpData[attr]))
+            for attr in requiredAttributes
+        }
+    except (TypeError, ValueError) as e:
+        raise ValueError(
+            f"Required shapefile attributes {requiredAttributes} must contain valid numeric values."
+        ) from e
+
+
+
+def _geological_azimuth_components(azimuth_deg):
+    """Return the horizontal unit vector of a geological dip direction.
+
+    Geological azimuths are measured clockwise from north:
+    0° = north, 90° = east, 180° = south, 270° = west.
+
+    Parameters
+    ----------
+    azimuth_deg : float
+        Geological dip-direction azimuth in degrees.
+
+    Returns
+    -------
+    tuple of float
+        Easting and northing components of the downslope unit vector.
+    """
+    azimuth_rad = math.radians(azimuth_deg % 360.0)
+    return math.sin(azimuth_rad), math.cos(azimuth_rad)
+
 
 def scarpAnalysisMain(cfg, baseDir):
     """Run the scarp analysis using parameters from an .ini cfguration file and input from a directory
@@ -80,14 +144,13 @@ def scarpAnalysisMain(cfg, baseDir):
     method = cfg["SETTINGS"]["method"].lower()
 
     if method == "plane":
-        # Read required attributes directly from the shapefile's attribute table
-        try:
-            planesZseed = list(map(float, SHPdata['zseed']))
-            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', 'dipdir', and 'dipangle' fields exist.")
+        # Read and validate required attributes from the shapefile's attribute table
+        requiredAttributes = ["zseed", "dipdir_azi", "dipAngle"]
+        planeAttributes = _extract_numeric_attributes(SHPdata, requiredAttributes)
+
+        planesZseed = planeAttributes["zseed"]
+        planesDipDir = planeAttributes["dipdir_azi"]
+        planesDipAngle = planeAttributes["dipAngle"]
 
         if not (len(planesZseed) == len(planesDipDir) == len(planesDipAngle) == SHPdata["nFeatures"]):
             raise ValueError("Mismatch between number of features and extracted plane attributes in the shapefile.")
@@ -109,16 +172,24 @@ def scarpAnalysisMain(cfg, baseDir):
         log.debug("Plane features extracted and combined: %s", features)
 
     elif method == "ellipsoid":
-        try:
-            ellipsoidsMaxDepth = list(map(float, SHPdata['maxdepth']))
-            ellipsoidsSemiMajor = list(map(float, SHPdata['semimajor']))
-            ellipsoidsSemiMinor = list(map(float, SHPdata['semiminor']))
-            ellipsoidsDipAngle = list(map(float, SHPdata['dipAngle']))
-            ellipsoidsDipDir = list(map(float, SHPdata['dipdir']))
-            ellipsoidsOffset = list(map(float, SHPdata['offset']))
-            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', 'rotangle', 'dipdir', 'dipangle', and 'offset' exist.")
+        requiredAttributes = [
+            "maxdepth",
+            "semimajor",
+            "semiminor",
+            "dipAngle",
+            "dipdir_azi",
+            "offset",
+            "rotAngle",
+        ]
+        ellipsoidAttributes = _extract_numeric_attributes(SHPdata, requiredAttributes)
+
+        ellipsoidsMaxDepth = ellipsoidAttributes["maxdepth"]
+        ellipsoidsSemiMajor = ellipsoidAttributes["semimajor"]
+        ellipsoidsSemiMinor = ellipsoidAttributes["semiminor"]
+        ellipsoidsDipAngle = ellipsoidAttributes["dipAngle"]
+        ellipsoidsDipDir = ellipsoidAttributes["dipdir_azi"]
+        ellipsoidsOffset = ellipsoidAttributes["offset"]
+        ellipsoidsRotAngle = ellipsoidAttributes["rotAngle"]
 
         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.")
@@ -229,7 +300,7 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes):
     elevTransform : Affine
         The affine transformation matrix of the raster (used to convert pixel to geographic coordinates).
     planes : str
-        Comma-separated string defining sliding planes (xseed, yseed, zseed, dip, slope).
+        Comma-separated string defining sliding planes (xseed, yseed, zseed, dip direction azimuth, dip angle).
 
     Returns
     -------
@@ -249,9 +320,11 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes):
     slope = [planes[4]]
 
     slopeRad = math.radians(slope[0])
-    dipRad   = math.radians(dip[0])
-    betaX = [ math.tan(slopeRad) * math.sin(dipRad) ]
-    betaY = [ math.tan(slopeRad) * math.cos(dipRad) ]
+    dipEast, dipNorth = _geological_azimuth_components(dip[0])
+    gradientMagnitude = math.tan(slopeRad)
+    # Elevation must decrease in the geological dip direction.
+    betaX = [-gradientMagnitude * dipEast]
+    betaY = [-gradientMagnitude * dipNorth]
 
     min_clipped = 0.0
 
@@ -263,17 +336,27 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes):
         slope.append(planes[5 * i + 4])
 
         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) )
+        dipEast, dipNorth = _geological_azimuth_components(dip[i])
+        gradientMagnitude = math.tan(slopeRad)
+        betaX.append(-gradientMagnitude * dipEast)
+        betaY.append(-gradientMagnitude * dipNorth)
 
     for row in range(n):
         for col in range(m):
             west, north = rasterio.transform.xy(elevTransform, row, col, offset='center')
 
-            scarpVal = zSeed[0] + (north - ySeed[0]) * betaY[0] - (west - xSeed[0]) * betaX[0]
+            scarpVal = (
+                zSeed[0]
+                + (west - xSeed[0]) * betaX[0]
+                + (north - ySeed[0]) * betaY[0]
+            )
             for k in range(1, nPlanes):
-                scarpVal = max(scarpVal, zSeed[k] + (north - ySeed[k]) * betaY[k] - (west - xSeed[k]) * betaX[k])
+                planeVal = (
+                    zSeed[k]
+                    + (west - xSeed[k]) * betaX[k]
+                    + (north - ySeed[k]) * betaY[k]
+                )
+                scarpVal = max(scarpVal, planeVal)
 
             if periData[row, col] > 0:
                 val = min(elevData[row, col], scarpVal)
@@ -301,7 +384,7 @@ def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids):
         The affine transformation matrix of the raster.
     ellipsoids : str
         Comma-separated string defining ellipsoids with parameters:
-        (x_center, y_center, max_depth, semi_major, semi_minor, tilt, dir, offset)
+        (maxdepth, semimajor, semiminor, dipAngle, dipdir_azi, offset, rotangle)
 
     Returns
     -------
@@ -328,7 +411,8 @@ def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids):
         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
+        # Geological dip-direction azimuth in degrees (0° N, 90° E).
+        tiltDir.append(ellipsoids[9 * i + 6] % 360.0)
         offset.append(ellipsoids[9 * i + 7])
         dip.append(np.radians(ellipsoids[9 * i + 8]))      # rotation of base ellipse in radians
 
@@ -382,8 +466,11 @@ def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids):
 
                 if distance <= 1:
                     baseDepth = maxDepth[k] * (1 - distance)
+                    dipEast, dipNorth = _geological_azimuth_components(tiltDir[k])
+                    # Positive dipAngle increases release depth in the geological
+                    # dip direction. rotAngle only rotates the ellipse footprint.
                     tiltEffect = math.tan(math.radians(tilt[k])) * (
-                        dxRot * np.cos(tiltDir[k]) + dyRot * np.sin(tiltDir[k])
+                        dxPos * dipEast + dyPos * dipNorth
                     )
                     totalDepth = baseDepth + tiltEffect + z0
                     scarpVal = min(scarpVal, elevData[row, col] - totalDepth)

avaframe/in2Trans/shpConversion.py

@@ -59,7 +59,7 @@ 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)
-        dipDir
+        dipDir_azimuth
             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)
@@ -84,7 +84,7 @@ def SHP2Array(infile, defname=None):
     ci95 = None
     layerN = None
     zseed_value = None
-    dipdir_value = None
+    dipdirAzimuth_value = None
     dipAngle_value = None
     semiminor_value = None
     maxdepth_value = None
@@ -108,7 +108,7 @@ def SHP2Array(infile, defname=None):
     ci95List = []
     layerNameList = []
     zseedList = []
-    dipdirList = []
+    dipdirAzimuthList = []
     slopeList = []
     dipAngleList = []
     semiminorList = []
@@ -162,8 +162,8 @@ def SHP2Array(infile, defname=None):
                     dipAngle_value = value
                 if name == "zseed":
                     zseed_value = value
-                if name == "dipdir":
-                    dipdir_value = value
+                if name == "dipdir_azi":
+                    dipdirAzimuth_value = value
                 if name == "semiminor":
                     semiminor_value = value
                 if name == "maxdepth":
@@ -201,7 +201,7 @@ def SHP2Array(infile, defname=None):
         layerNameList.append(layerN)
         idList.append(str(rec.oid))
         zseedList.append(zseed_value)
-        dipdirList.append(dipdir_value)
+        dipdirAzimuthList.append(dipdirAzimuth_value)
         slopeList.append(slope)
         dipAngleList.append(dipAngle_value)
         semiminorList.append(semiminor_value)
@@ -236,7 +236,7 @@ def SHP2Array(infile, defname=None):
     SHPdata["nFeatures"] = len(Start)
     SHPdata["dipAngle"] = dipAngleList
     SHPdata["zseed"] = zseedList
-    SHPdata["dipdir"] = dipdirList
+    SHPdata["dipdir_azi"] = dipdirAzimuthList
     SHPdata["maxdepth"] = maxdepthList
     SHPdata["semimajor"] = semimajorList
     SHPdata["semiminor"] = semiminorList