diff --git a/avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini b/avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini index 423307056..4b5b7abe8 100644 --- a/avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini +++ b/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. - - diff --git a/avaframe/com6RockAvalanche/scarp.py b/avaframe/com6RockAvalanche/scarp.py index dc3b9b938..f790963ea 100644 --- a/avaframe/com6RockAvalanche/scarp.py +++ b/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) diff --git a/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.dbf b/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.dbf index 937a33eac..8a9686bd6 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.dbf and b/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.dbf differ diff --git a/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.prj b/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.prj deleted file mode 100644 index f45cbadf0..000000000 --- a/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.prj +++ /dev/null @@ -1 +0,0 @@ -GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] \ No newline at end of file diff --git a/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.shp b/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.shp index e1b59674d..08d71efbf 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.shp and b/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.shp differ diff --git a/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.shx b/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.shx index e13af6299..fe60c7979 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.shx and b/avaframe/data/scarpExample/Inputs/POINTS/EllipsoidMethodPoints_coordinates.shx differ diff --git a/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.dbf b/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.dbf index 0fcf199fa..743ab5d99 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.dbf and b/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.dbf differ diff --git a/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.prj b/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.prj deleted file mode 100644 index f45cbadf0..000000000 --- a/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.prj +++ /dev/null @@ -1 +0,0 @@ -GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] \ No newline at end of file diff --git a/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.shp b/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.shp index a5689bf17..fd8869249 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.shp and b/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.shp differ diff --git a/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.shx b/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.shx index f84bbc10f..198e1d2a9 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.shx and b/avaframe/data/scarpExample/Inputs/POINTS_plane/PlaneMethodPoints_coordinates.shx differ diff --git a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.dbf b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.dbf index 1d9fbd3ca..ec9b562cd 100644 Binary files a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.dbf and b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.dbf differ diff --git a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.prj b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.prj deleted file mode 100644 index f45cbadf0..000000000 --- a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.prj +++ /dev/null @@ -1 +0,0 @@ -GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] \ No newline at end of file diff --git a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shp b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shp index e633df8e2..91eab622d 100644 Binary files a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shp and b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shp differ diff --git a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shx b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shx index ef40ee0cb..69d2dbc14 100644 Binary files a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shx and b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shx differ diff --git a/avaframe/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index 6eef2d743..57f26ef1c 100644 --- a/avaframe/in2Trans/shpConversion.py +++ b/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 diff --git a/avaframe/tests/test_scarp.py b/avaframe/tests/test_scarp.py index 79031d4f5..19fcde7fb 100644 --- a/avaframe/tests/test_scarp.py +++ b/avaframe/tests/test_scarp.py @@ -116,7 +116,7 @@ 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["dipdir"])) + planesDip = list(map(float, SHPdata["dipdir_azi"])) planesSlope = list(map(float, SHPdata["dipAngle"])) # Assertions @@ -150,22 +150,22 @@ def test_plane_geometry_calculations(): dip = 45.0 # degrees # Expected calculations - expected_betaX = math.tan(math.radians(slope)) * math.cos(math.radians(dip)) - expected_betaY = math.tan(math.radians(slope)) * math.sin(math.radians(dip)) + expected_betaX = -math.tan(math.radians(slope)) * math.sin(math.radians(dip)) + expected_betaY = -math.tan(math.radians(slope)) * math.cos(math.radians(dip)) # Assertions - these are the formulas used in calculateScarpWithPlanes - assert abs(expected_betaX - 0.408248) < 0.001, "betaX calculation should be correct" - assert abs(expected_betaY - 0.408248) < 0.001, "betaY calculation should be correct" + assert abs(expected_betaX + 0.408248) < 0.001, "betaX calculation should be correct" + assert abs(expected_betaY + 0.408248) < 0.001, "betaY calculation should be correct" # Test plane equation xSeed, ySeed, zSeed = 100.0, 200.0, 1000.0 west, north = 150.0, 250.0 # Point coordinates - # Plane equation: z = zSeed + (north - ySeed) * betaY - (west - xSeed) * betaX - scarpVal = zSeed + (north - ySeed) * expected_betaY - (west - xSeed) * expected_betaX + # Plane equation: z = zSeed + (west - xSeed) * betaX + (north - ySeed) * betaY + scarpVal = zSeed + (west - xSeed) * expected_betaX + (north - ySeed) * expected_betaY # Manual calculation - expected_scarpVal = 1000.0 + (50.0 * expected_betaY) - (50.0 * expected_betaX) + expected_scarpVal = 1000.0 + (50.0 * expected_betaX) + (50.0 * expected_betaY) assert abs(scarpVal - expected_scarpVal) < 0.001, "Plane equation should be correct" @@ -390,58 +390,20 @@ def test_scarpAnalysisMain_invalid_method(scarp_test_data, scarp_config, tmp_pat scarp.scarpAnalysisMain(scarp_config, str(test_dir)) -def test_scarpAnalysisMain_missing_required_attributes(scarp_test_data, tmp_path): +def test_scarpAnalysisMain_missing_required_attributes(): """Test that missing required plane attributes raises ValueError""" - # This test checks the error path when shapefile is missing required attributes - # We test this by examining that the code properly validates attribute existence - - # For this test, we'd need to create a shapefile with missing attributes, - # which is complex. The code path is covered by the KeyError handling at lines 88-89 - # in scarp.py. We verify the error message is descriptive. - - # Create config - cfg = configparser.ConfigParser() - cfg["INPUT"] = {"useShapefiles": "True"} - cfg["SETTINGS"] = {"method": "plane"} - - # The test data has correct attributes, so we can't test the error path easily - # without creating invalid shapefiles. We document this limitation. - assert True, "Error path for missing attributes tested through code inspection" + shpData = {"zseed": ["1000"], "dipdir_azi": [None], "dipAngle": ["30"]} + with pytest.raises(ValueError, match="dipdir_azi"): + scarp._extract_numeric_attributes(shpData, ["zseed", "dipdir_azi", "dipAngle"]) def test_error_message_attribute_names(): - """Test that error messages reference correct attribute names""" - # This test verifies that error messages in scarp.py reference - # the same attribute names that are actually used in the code - + """Test that scarp.py references correct attribute names""" import avaframe.com6RockAvalanche.scarp as scarp_module - # Read the scarp.py file to check error messages scarp_path = pathlib.Path(scarp_module.__file__) scarp_content = scarp_path.read_text() - # Check line 90 error message - line_90_match = None - for i, line in enumerate(scarp_content.split("\\n"), 1): - if i == 90: - line_90_match = line - break - - # The error message should reference 'dipAngle' not 'dipangle' - if line_90_match: - assert ( - "'dipAngle'" in line_90_match or "'dipangle'" in line_90_match - ), f"Line 90 should reference dipAngle attribute: {line_90_match}" - - # Check line 121 error message - line_121_match = None - for i, line in enumerate(scarp_content.split("\\n"), 1): - if i == 121: - line_121_match = line - break - - # The error message should reference 'rotAngle' not 'rotangle' - if line_121_match: - assert ( - "'rotAngle'" in line_121_match or "'rotangle'" in line_121_match - ), f"Line 121 should reference rotAngle attribute: {line_121_match}" + assert "dipdir_azi" in scarp_content, "scarp.py should reference dipdir_azi" + assert "dipAngle" in scarp_content, "scarp.py should reference dipAngle" + assert "rotAngle" in scarp_content, "scarp.py should reference rotAngle" diff --git a/docs/_static/com6_ellipsoid_crosssection.png b/docs/_static/com6_ellipsoid_crosssection.png new file mode 100644 index 000000000..b0ec7b844 Binary files /dev/null and b/docs/_static/com6_ellipsoid_crosssection.png differ diff --git a/docs/_static/com6_ellipsoid_offset.png b/docs/_static/com6_ellipsoid_offset.png new file mode 100644 index 000000000..1749e3fb2 Binary files /dev/null and b/docs/_static/com6_ellipsoid_offset.png differ diff --git a/docs/_static/com6_ellipsoid_topview.png b/docs/_static/com6_ellipsoid_topview.png new file mode 100644 index 000000000..53291b133 Binary files /dev/null and b/docs/_static/com6_ellipsoid_topview.png differ diff --git a/docs/_static/com6_plane_topview.png b/docs/_static/com6_plane_topview.png new file mode 100644 index 000000000..882fe4e7e Binary files /dev/null and b/docs/_static/com6_plane_topview.png differ diff --git a/docs/_static/com6_planes_crosssection.png b/docs/_static/com6_planes_crosssection.png new file mode 100644 index 000000000..95595d1da Binary files /dev/null and b/docs/_static/com6_planes_crosssection.png differ diff --git a/docs/moduleCom6RockAvalanche.rst b/docs/moduleCom6RockAvalanche.rst index 6f763e23b..6c3ad0ef3 100644 --- a/docs/moduleCom6RockAvalanche.rst +++ b/docs/moduleCom6RockAvalanche.rst @@ -6,23 +6,48 @@ com6RockAvalanche: Rock Avalanche The com6RockAvalanche computational module provides an override setting for com1DFA targeting the simulation of rock avalanches. +Tips +---- + +* Download the default configuration for each module via + ``OpenNHM > AvaFrame_Experimental > Get default module ini`` (see :ref:`connector:Experimental`). +* After a run, restart the ``Rock Avalanche (com6)`` tool before rerunning with changed parameters; a fresh instance + is often more reliable than reusing the old one. +* Parameter meanings can be looked up in the documentation, e.g. ``massPerPart``: + https://docs.avaframe.org/en/latest/com1DFAAlgorithm.html#initialize-particles + Input ------- -The standard inputs required to perform a simulation run using :py:mod:`com1DFA` +The standard inputs required to perform a simulation run using :py:mod:`com1DFA` can be found here: :ref:`moduleCom1DFA:Input`. However there is one main difference: com6RockAvalanche NEEDS a release thickness raster file. This file has to have -the exact same dimensions as the topography file. +the exact same dimensions and resolution as the topography file. There is a run script to perform a rock avalanche com1DFA run: :py:mod:`runCom6RockAvalanche.py`, and the configuration settings can be found in ``com6RockAvalanche/com6RockAvalancheCfg.ini``. +The following files are required: + +* a DEM from which the release volume has already been removed + (see :ref:`moduleCom6RockAvalanche:Scarp Calculation`) +* a release thickness raster with the same extent and resolution as the DEM + +The following are optional: + +* an entrainment polygon shape file (raster support may follow) +* a configuration file that overrides the default settings (see :ref:`connector:Experimental` for how to obtain it) + To run ------ * first go to ``AvaFrame/avaframe`` * copy ``avaframeCfg.ini`` to ``local_avaframeCfg.ini`` and set your desired avalanche directory name -* create an avalanche directory with required input files - for this task you can use :ref:`moduleIn3Utils:Initialize Project` -* copy ``com6RockAvalanche/com6RockAvalancheCfg.ini`` to ``com6RockAvalanche/local_com6RockAvalancheCfg.ini`` and if desired change configuration settings +* create an avalanche directory with required input files - for this task you can use + :ref:`moduleIn3Utils:Initialize Project` +* copy ``com6RockAvalanche/com6RockAvalancheCfg.ini`` to + ``com6RockAvalanche/local_com6RockAvalancheCfg.ini`` and if desired change configuration settings +* if you are on a develop installation, make sure you have an updated compilation, see + :ref:`complexUsage:Update AvaFrame` * optionally, the ``spatialVoellmy`` friction model can be selected with ``--friction_calibration spatialVoellmy``. When using the QGis Connector, a shapefile with ``mu`` and ``xi`` attributes @@ -36,64 +61,240 @@ To run pixi run python runCom6RockAvalanche.py +Run via QGis Connector +---------------------- + +Alternatively, run from QGIS via the OpenNHM connector: open the OpenNHM tool and select ``Rock Avalanche (com6)``. +Choose a suitable DEM and release layer, and optionally an entrainment layer. A path to a configuration file with +additional settings can also be provided. + +It is strongly recommended to set ``meshCellSize`` to the cell size of the input rasters +(see :ref:`moduleCom6RockAvalanche:Remeshing`); otherwise remeshing artifacts can occur. For example, the Fluchthorn +dataset uses 10 m. + +.. Warning:: Depending on the DEM, release and settings, runs can take very long (up to several hours). + +Remeshing +--------- + +The default configuration remeshes the DEM and release raster to 5 m unless their cell size already equals 5 m. This +occasionally introduces no-data values, which then cause errors. + +Either set the mesh size in the configuration file, or reproject the inputs (DEM and release raster) to the target +cell size beforehand. In QGIS use ``Raster -> Projections -> Warp (Reproject)`` and set the output resolution, then +repeat for the release raster. To adjust the configuration instead, add or change:: + + # expected mesh size [m]; use the cell size of your raster + meshCellSize = 10 + +Entrainment +----------- + +A full description of the entrainment model, including formulas, is available here: +https://docs.avaframe.org/en/latest/theoryCom1DFA.html#entrainment +and an introduction on how to use entrainment is available here: +https://docs.avaframe.org/en/latest/moduleCom1DFA.html#input + +Entrainment is the uptake of material by the rock avalanche during its flow. Two main processes are distinguished: + +* plowing: uptake of material at the rock avalanche front +* erosion: uptake of material at the rock avalanche base + +com1DFA uses basal erosion by default; plowing is disabled. The entrainment rate strongly influences rock avalanche +mass and dynamics. Entrainment can be provided as a shape or raster file, but the QGIS connector currently only supports +shapes. Use a (multi-)polygon shape file that: + +* defines the areas where entrainment occurs (areas should not overlap) +* carries a ``thickness`` attribute (entrainment material thickness, measured normal to the slope) and contains no + holes or rings + +If no ``thickness`` attribute is present, the default ``entThIfMissingInShp`` from the configuration file is used. +The flag ``THICKNESSFromFile`` (i.e. ``relThFromFile``, ``entThFromFile``, ``secondaryRelThFromFile``) must be set to +True in the configuration file (default is True). + +Thickness variations for parameter studies can be defined as follows: + +* ``entThPercentVariation``: ``+-percentage$numberOfSteps``; ``+`` gives a positive variation, ``-`` a negative one, + no sign gives both directions +* ``entThRangeVariation``: ``+-range$numberOfSteps``, same sign convention +* ``entThRangeFromCiVariation``: ``ci95$numberOfSteps``, varies the thickness within +- the 95% confidence interval + read from a ``ci95`` attribute in the shape file + +Entrainment parameters +---------------------- + +.. list-table:: + :header-rows: 1 + :widths: 25 35 40 + + * - Parameter + - Meaning + - Notes / effect + * - ``entThFromFile`` + - Read entrainment areas from a shape file (True) or use the global ``entTh`` value (False). + - If True, expects a file such as ``Inputs/ENT/entrainment.shp``. + * - ``entThIfMissingInShp`` + - Fallback thickness [m] if the shape file has no thickness attribute. + - Default 0 m; used for features without a valid value. + * - ``entThPercentVariation``, ``entThRangeVariation``, ``entThRangeFromCiVariation``, ``entThDistVariation`` + - Options for parameter studies (entrainment thickness variations). + - E.g. ``entThRangeVariation = 0.5$10`` gives a +-0.5 m variation in 10 steps. + * - ``entTh`` + - Fixed entrainment thickness [m], only relevant if ``entThFromFile = False``. + - E.g. ``entTh = 0.2``. + * - ``entEroEnergy`` + - Erosion energy [J/m2]; controls the energy loss from mass uptake. + - Higher values lose more energy (and thus speed) during entrainment. + * - ``entShearResistance`` + - Shear resistance of the entrained material. + - 0 = no additional resistance. + * - ``entDefResistance`` + - Deformation resistance of the material. + - 0 = no deformation resistance. + +Common issues +------------- + +* Remeshing or other processing errors produce no-data values. Set ``meshCellSize`` to the raster cell size + (see :ref:`moduleCom6RockAvalanche:Remeshing`). +* Negative values in the release thickness raster cause errors. The + :ref:`moduleCom6RockAvalanche:Scarp Calculation` step clamps negative thicknesses to 0. + Scarp Calculation ----------------- - * first go to ``AvaFrame/avaframe`` * copy ``avaframeCfg.ini`` to ``local_avaframeCfg.ini`` and set your desired avalanche directory name * create an avalanche directory - for this task you can use :ref:`moduleIn3Utils:Initialize Project` -Input -~~~~~ +Scarp Input +~~~~~~~~~~~ * all input files are automatically read from the set avalancheDir. No file paths need to be specified * 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", "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. + attributes ``zseed``, ``dipdir_azi`` and ``dipAngle`` as float values. If the ellipsoid method is used, the shape + file must contain the attributes ``maxdepth``, ``semimajor``, ``semiminor``, ``dipAngle``, ``dipdir_azi``, + ``rotAngle`` and ``offset`` (see below). The file must be located in avalancheDir/Inputs/POINTS and the file name + must end with ``_coordinates``. If you are using the QGis Connector, the naming and location of the file is not + relevant. +* perimeter: a shapefile that defines the spatial extent within which the scarp geometry is applied. It is rasterized + to a binary mask and used to clip the scarp surface to a predefined area: inside the perimeter the calculated scarp + elevation replaces the DEM where it is lower, outside the perimeter the original DEM is kept. This allows to limit + the scarp to a geologically meaningful release area, prevent artificial terrain modification outside the intended + scarp, and combine multiple scarp elements without affecting the surrounding topography. The file must be located + in avalancheDir/Inputs/POLYGONS and the file name must end with ``_perimeter``. If you are using the QGis + Connector, the naming and location of the file is not relevant. + +Attribute meanings +~~~~~~~~~~~~~~~~~~ + +Plane: + +* ``zseed``: z coordinate of the plane center (m). Usually the pre-event terrain elevation at the scarp initiation + point, but any value can be set. +* ``dipdir_azi``: azimuth, the direction the plane faces (degree). +* ``dipAngle``: tilt angle of the plane (degree). + +Ellipsoid: + +* ``maxdepth``: maximum depth of the untilted ellipsoid (m). E.g. 50 m means the geometric center of the untilted + ellipsoid lies 50 m below the surface in the vertical direction. +* ``semimajor``: half length of the major axis (m). +* ``semiminor``: half length of the minor axis (m). +* ``dipAngle``: tilt angle of the ellipsoid, i.e. inclination of its x axis (degree). +* ``dipdir_azi``: azimuth, the direction the ellipsoid faces (degree). +* ``rotAngle``: rotation angle of the ellipsoid base (degree). +* ``offset``: offset normal to the DEM slope (m). + +Running Scarp +~~~~~~~~~~~~~ + +From QGIS, open the OpenNHM tool and select ``Scarp (com6)``. Select the input DEM, perimeter shape and coordinate +file, and choose the method. The chosen method must match the attributes created in the coordinate file. + +.. Note:: Very large DEMs can lead to long runtimes depending on the machine; keep the input files as small as + possible, and avoid leaving or interrupting QGIS during the computation. + +Alternatively, run from the command line (see :ref:`moduleCom6RockAvalanche:Scarp Config`): + +:: -* perimeter: A shapefile that specifies a boundary area. Must be located in avalancheDir/Inputs/POLYGONS and file name - must end with “_perimeter”. If you are using the QGis Connector, the naming and location of the file is not relevant. + pixi run python runCom6Scarp.py -**Attribute meanings:** +Scarp parameter sketches +~~~~~~~~~~~~~~~~~~~~~~~~ -* 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) +The following sketches illustrate the correct use of the Scarp parameters. -* maxdepth: maximum depth of the ellipsoid (m) -* semimajor: length of the major axis (m) -* semiminor: length of the minor axis (m) -* 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) +.. figure:: /_static/com6_ellipsoid_topview.png + :width: 70% + :alt: Top view of a rotated but untilted ellipsoid -Output -~~~~~~ + Top view of a rotated but untilted ellipsoid (drawn unrotated to avoid distortion). The perimeter boundary does + not cut anything from the ellipsoid. -* elevscarp: Output DGM (ASCII or GeoTIFF), which maps the input DGM minus the calculated scarp. Is saved under +.. figure:: /_static/com6_ellipsoid_crosssection.png + :width: 70% + :alt: Cross-section of a tilted ellipsoid + + Cross-section of an ellipsoid tilted by 30 degrees but not rotated. The maximum depth refers to the untilted + ellipsoid, so ``maxdepth = 50 m`` does not necessarily mean a 50 m cut depth; it depends on the tilt angles. + +.. figure:: /_static/com6_ellipsoid_offset.png + :width: 70% + :alt: Cross-section of an ellipsoid with slope-normal offset + + Cross-section of an ellipsoid with a positive slope-normal offset of the center. + +.. figure:: /_static/com6_plane_topview.png + :width: 70% + :alt: Top view of the plane method + + Top view of the plane method. P1 and P2 define two points/planes dipping in the same direction. Different dip + angles combined with the perimeter clip produce a failure body (see next figure). The grey area is the top-view + failure surface clipped by the perimeter boundary. + +.. figure:: /_static/com6_planes_crosssection.png + :width: 70% + :alt: Cross-section of the planes + + Cross-section of the planes. The planes defined by P1 and P2 dip in the same direction with different dip + angles, forming a failure body (orange area). + +Scarp Output +~~~~~~~~~~~~ + +The Scarp step produces a DEM from which the release area has been cut out, plus a raster file with the release +thickness. Both rasters have the same resolution and extent as the original DEM. Negative thicknesses are set to 0 to +avoid errors in the rock avalanche module; the largest (absolute) negative value is logged to indicate how much the +script had to correct. + +* elevscarp: Output DGM (ASCII or GeoTIFF), which maps the input DGM minus the calculated scarp. It is saved under ``scarpElevation.(asc/tif)`` in ``avalancheDir/Outputs/com6RockAvalanche/scarp``. -* hrelease: File path to the output DGM (ASCII or GeoTIFF), which represents the calculated scarp volumes.Is saved - under ``scarpHRel.(asc/tif)`` in ``avalancheDir/Outputs/com6RockAvalanche/scarp``. +* hrelease: File path to the output DGM (ASCII or GeoTIFF), which represents the calculated scarp volumes. It is + saved under ``scarpHRel.(asc/tif)`` in ``avalancheDir/Outputs/com6RockAvalanche/scarp``. -Config -~~~~~~ +Scarp Config +~~~~~~~~~~~~ Prepare the config file (scarpCfg.ini): * copy ``com6RockAvalanche/scarpCfg.ini`` to ``com6RockAvalanche/local_scarpCfg.ini`` and if desired change -configuration settings - -* Input: - o set useShapefiles = True -* Settings: - o method: Here you specify whether the plane or the ellipsoid method should be used + configuration settings +* Input: set ``useShapefiles = True`` +* Settings: ``method`` specifies whether the plane or the ellipsoid method is used If all the data is provided successfully, start the script by running:: - pixi run python runCom6Scarp.py \ No newline at end of file + pixi run python runCom6Scarp.py + +Scarp common issues +~~~~~~~~~~~~~~~~~~~ + +* Projections: mismatched projections usually result in an input DEM plus a release file with 0 m thickness + everywhere. Use one consistent projection for all inputs. +* Attribute names: incorrect attribute names (e.g. wrong case) fail silently. Copy the names from the attribute + list above. +* Empty attribute fields: every field must contain a value; enter 0 where a field is not applicable.