Update Scarp.py

avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini

@@ -48,3 +48,64 @@ frictModel = Voellmy
 #+++++++++++++Voellmy friction model
 muvoellmy = 0.035
 xsivoellmy = 700.
+# expected mesh size [m]
+meshCellSize = 5
+
+# density of entrained snow [kg/m³]
+rhoEnt = 100
+#+++++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
+# 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 [%]
+entThPercentVariation =
+# if a variation on entTh shall be performed add here +- absolute value and number of steps separated by $
+# for example entThRangeVariation=0.5$10 [m]
+entThRangeVariation =
+# if a variation on entTh shall be performed add here +- ci% value and number of steps separated by $
+# for example entThRangeFromCiVariation= ci95$10
+entThRangeFromCiVariation =
+# if variation on entTh shall be performed using a normal distribution in number of steps,
+# value of buildType (ci95 value), min and max of dist in percent, buildType (ci95 only allowed),
+# support (e.g. 10000) all separated by $: e.g. normaldistribution$numberOfSteps$0.3$95$ci95$10000
+# if entFromShp=True ci95 is read from shp file too
+entThDistVariation =
+# entrainment thickness (only considered if ENT file is shapefile and entThFromFile=False) [m]
+entTh =
+
+
+#++++++++++++ Entrainment Erosion Energy
+# 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

@@ -149,8 +149,14 @@ def scarpAnalysisMain(cfg, baseDir):
         )
     else:
         raise ValueError("Unsupported method. Choose 'plane' or 'ellipsoid'.")
+
 
     hRelData = dem["rasterData"] - scarpData
+    min_hrel = hRelData.min()
+    if min_hrel < 0:
+        log.warning("hRel: most negative clipped value = %.2f m", min_hrel)
+    hRelData = np.maximum(0, hRelData)   #enforce non-negative release
+
 
     #Compute and log excavated volume
     cellArea = abs(dem["header"]["cellsize"] ** 2)
@@ -247,6 +253,8 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes):
     betaX = [ math.tan(slopeRad) * math.sin(dipRad) ]
     betaY = [ math.tan(slopeRad) * math.cos(dipRad) ]
 
+    min_clipped = 0.0
+
     for i in range(1, nPlanes):
         xSeed.append(planes[5 * i])
         ySeed.append(planes[5 * i + 1])
@@ -266,11 +274,17 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes):
             scarpVal = zSeed[0] + (north - ySeed[0]) * betaY[0] - (west - xSeed[0]) * betaX[0]
             for k in range(1, nPlanes):
                 scarpVal = max(scarpVal, zSeed[k] + (north - ySeed[k]) * betaY[k] - (west - xSeed[k]) * betaX[k])
-
+                
             if periData[row, col] > 0:
-                scarpData[row, col] = min(elevData[row, col], scarpVal)
+                val = min(elevData[row, col], scarpVal)
+                if val < 0:
+                    min_clipped = min(min_clipped, val)
+                scarpData[row, col] = max(0, val)
             else:
                 scarpData[row, col] = elevData[row, col]
+
+    if min_clipped < 0:
+        log.warning("Plane scarp: most negative clipped value = %.2f m", min_clipped)
 
     return scarpData
 
@@ -305,6 +319,7 @@ def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids):
     xCenter, yCenter, maxDepth = [], [], []
     semiMajor, semiMinor = [], []
     tilt, tiltDir, offset, dip = [], [], [], []
+    min_clipped = 0.0
 
     for i in range(nEllipsoids):
         xCenter.append(ellipsoids[9 * i])
@@ -373,6 +388,15 @@ def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids):
                     totalDepth = baseDepth + tiltEffect + z0
                     scarpVal = min(scarpVal, elevData[row, col] - totalDepth)
 
-            scarpData[row, col] = scarpVal if periData[row, col] > 0 else elevData[row, col]
+            if periData[row, col] > 0:
+                val = min(elevData[row, col], scarpVal)
+                if val < 0:
+                    min_clipped = min(min_clipped, val)
+                scarpData[row, col] = max(0, val)
+            else:
+                scarpData[row, col] = elevData[row, col]
+
+    if min_clipped < 0:
+        log.warning("Ellipsoid scarp: most negative clipped value = %.2f m", min_clipped)            
 
     return scarpData