diff --git a/avaframe/com6RockAvalanche/scarp.py b/avaframe/com6RockAvalanche/scarp.py index 10eb5d0f0..08f5b7a42 100644 --- a/avaframe/com6RockAvalanche/scarp.py +++ b/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,8 +242,10 @@ 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]) @@ -245,8 +253,11 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes): 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] diff --git a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.dbf b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.dbf index ef0170f04..caca82bce 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.dbf and b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.dbf differ diff --git a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shp b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shp index 0cc7e0016..a1cb3ab00 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shp and b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shp differ diff --git a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shx b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shx index adec32cf5..9f4f0f8bf 100644 Binary files a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shx and b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.shx differ diff --git a/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.dbf b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.dbf index 742c45010..46aebd550 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.shp b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shp index ec5f2194b..a07334ed1 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 c2859ca3a..c11c1b409 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/data/scarpExample/points_coordinates.dbf b/avaframe/data/scarpExample/points_coordinates.dbf new file mode 100644 index 000000000..b8517812c Binary files /dev/null and b/avaframe/data/scarpExample/points_coordinates.dbf differ diff --git a/avaframe/data/scarpExample/points_coordinates.shp b/avaframe/data/scarpExample/points_coordinates.shp new file mode 100644 index 000000000..a1cb3ab00 Binary files /dev/null and b/avaframe/data/scarpExample/points_coordinates.shp differ diff --git a/avaframe/data/scarpExample/points_coordinates.shx b/avaframe/data/scarpExample/points_coordinates.shx new file mode 100644 index 000000000..9f4f0f8bf Binary files /dev/null and b/avaframe/data/scarpExample/points_coordinates.shx differ diff --git a/avaframe/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index 620d81ced..ffebf650d 100644 --- a/avaframe/in2Trans/shpConversion.py +++ b/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