From 7ad240115e1b87cdac054c014707568ac05356ba Mon Sep 17 00:00:00 2001 From: dwolfsch <168710278+dwolfsch@users.noreply.github.com> Date: Wed, 22 Apr 2026 11:04:59 +0200 Subject: [PATCH] Update Scarp.py and config file changes for com6RockAvalanche An Update for the scarp script, so that it does not allow negative values in the generated raster files anymore. Consistency fix fpr scarp.py and config file changes for com6RockAvalanche Small fixes for scarp.py, so that the code is consistent. Options added for the com6RockAvalancheCfg.ini file, so that it includes all parameters, that should be changable --- .../com6RockAvalancheCfg.ini | 61 +++++++++++++++++++ avaframe/com6RockAvalanche/scarp.py | 30 ++++++++- 2 files changed, 88 insertions(+), 3 deletions(-) diff --git a/avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini b/avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini index 852e8fc12..423307056 100644 --- a/avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini +++ b/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. + + diff --git a/avaframe/com6RockAvalanche/scarp.py b/avaframe/com6RockAvalanche/scarp.py index d3ae2f762..dc3b9b938 100644 --- a/avaframe/com6RockAvalanche/scarp.py +++ b/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