diff --git a/.github/workflows/runTestSinglePython.yml b/.github/workflows/runTestSinglePython.yml index 9f7eeab54..104bc4fa3 100644 --- a/.github/workflows/runTestSinglePython.yml +++ b/.github/workflows/runTestSinglePython.yml @@ -54,16 +54,6 @@ jobs: - name: Test with pytest run: | pytest -ra --cov --cov-report=xml --cov-report lcov:cov.info --cov-config=.coveragerc - - name: Test & publish code coverage - uses: paambaati/codeclimate-action@v9.0.0 - env: - CC_TEST_REPORTER_ID: ${{ secrets.CODECLIMATE_ID }} - # - name: Upload coverage to Codecov - # uses: codecov/codecov-action@v3 - # with: - # fail_ci_if_error: false - # name: codecov-umbrella - # verbose: true - uses: qltysh/qlty-action/coverage@main with: coverage-token: ${{secrets.QLTY_COVERAGE_TOKEN}} diff --git a/avaframe/avaframeCfg.ini b/avaframe/avaframeCfg.ini index 4572b2567..2c07a426e 100644 --- a/avaframe/avaframeCfg.ini +++ b/avaframe/avaframeCfg.ini @@ -4,7 +4,8 @@ [MAIN] # Path to avalanche directory -avalancheDir = data/avaParabola +avalancheDir = data/scarpExample + # number of CPU cores to use for the computation of com1DFA # possible values are: diff --git a/avaframe/com6RockAvalanche/ScarpReadMe.txt b/avaframe/com6RockAvalanche/ScarpReadMe.txt new file mode 100644 index 000000000..4f6084487 --- /dev/null +++ b/avaframe/com6RockAvalanche/ScarpReadMe.txt @@ -0,0 +1,32 @@ +Prepare files and config file (avaframeCfg.ini): +- Set path to avalancheDir in avaframeCfg.ini and create it +- Inputs: + o All input files are automatically read from the set avalancheDir. No file paths need to be specified + o elevation: DEM (ASCII), which serves as the basis for calculating the scarps. Must be in avalancheDir/Inputs. + o shapefile: 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", "dip" and "slopeangle" as float values. If the ellipsoid method is used, the shape file must contain the Attributes "maxdepth", "semimajor", "semiminor", "tilt", "direc", "dip", "offset". The file must be located in avalancheDir/Inputs/POINTS and file name must end with “_coordinates”. + o perimeter_shapefile: A shapefile that specifies a boundary area. Must be located in avalancheDir/Inputs/POLYGONS and file name must end with “_perimeter”. + +-Output: + o elevscarp: Output DGM (ASCII or GeoTIFF), which maps the input DGM minus the calculated scarp. Is saved under “raster_scarp.asc” in avalancheDir/Inputs/RASTERS. + o hrelease: File path to the output DGM (ASCII or GeoTIFF), which represents the calculated scarp volumes.Is saved under “raster_rel.asc” in avalancheDir/Inputs/REL. + +Prepare the config file (scarpCfg.ini): +- Input: + o set useShapefiles = True +- Settings: + o method: Here you specify whether the plane or the ellipsoid method should be used + +If all the data is provided successfully, start the script by running runCom6Scarp.py + +Attribute meanings: +zseed: defines z coordinate of plane Center (m) +dip: direction in which the plane/slope is facing (degree) +slopeangle: steepness/angle of the slope (degree) + +maxdepth: maximum depth of the ellipsoid (m) +semimajor: length of the major axis (m) +semiminor: length of the minor axis (m) +tilt: steepness/angle of the slope (degree) +direc: direction in which the slope is facing (degree) +dip: direction in which the ellipsoid is facing (degree) +Offset: offset, normal to the DEM slope (m) \ No newline at end of file diff --git a/avaframe/com6RockAvalanche/scarp.py b/avaframe/com6RockAvalanche/scarp.py new file mode 100644 index 000000000..804dc3602 --- /dev/null +++ b/avaframe/com6RockAvalanche/scarp.py @@ -0,0 +1,367 @@ +import rasterio +import numpy as np +import math +import logging +import pathlib +from rasterio.features import rasterize +from shapely.geometry import shape, mapping +from avaframe.in2Trans.shpConversion import SHP2Array +from avaframe.in1Data import getInput as gI +from avaframe.in3Utils import fileHandlerUtils as fU +import avaframe.in2Trans.rasterUtils as IOf + +# Configure logging +logging.basicConfig(level=logging.DEBUG) +log = logging.getLogger(__name__) + + +def scarpAnalysisMain(cfg, baseDir): + """Run the scarp analysis using parameters from an .ini cfguration file and input from a directory + + Parameters + ---------- + cfg : cfgparser + with all required fields in scarpCfg.ini + baseDir: str + path to directory of data to analyze + + Returns + ------- + None + """ + + # Get input from baseDir and check if all files exist + if cfg["INPUT"].getboolean("useShapefiles"): + + # Set directories for inputs, outputs and current work + inputDir = pathlib.Path(baseDir, "Inputs") + + # get the path to the perimeter shapefile + perimeterShapefilePath, periFile = gI.getAndCheckInputFiles( + inputDir, "POLYGONS", "scarp perimeter", fileExt="shp", fileSuffix="_perimeter" + ) + if periFile: + log.info("Perimeterfile is: %s" % perimeterShapefilePath) + else: + log.error("Perimeter shapefile not found in %s", inputDir) + + # get the path to the coordinates shapefile + coordinatesShapefilePath, coordFile = gI.getAndCheckInputFiles( + inputDir, "POINTS", "scarp coordinates", fileExt="shp", fileSuffix="_coordinates" + ) + if coordFile: + log.info("Coordinate shapefile is: %s" % coordinatesShapefilePath) + else: + log.error("Coordinate shapefile not found in %s", inputDir) + + else: + log.error( + "Shapefile option not selected. Please set 'useShapefiles' to True in the configuration file. A " + "raster version is currently not implemented." + ) + + # Initialize feature parameters + planeFeatures = [] + ellipsoidFeatures = [] + + # Read the DEM + dem = gI.readDEM(baseDir) + dem["rasterData"] = np.flipud(dem["rasterData"]) + + n = dem["header"]["nrows"] + m = dem["header"]["ncols"] + + periData = readPerimeterSHP(perimeterShapefilePath, dem["header"]["transform"], (n, m)) + + SHPdata = SHP2Array(coordinatesShapefilePath) + + log.debug("Feature shapefile data loaded: %s", SHPdata) + + 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'])) + planesDip = list(map(float, SHPdata['dip'])) + planesSlope = list(map(float, SHPdata['slopeangle'])) + except KeyError as e: + raise ValueError(f"Required attribute '{e.args[0]}' not found in shapefile. Make sure 'zseed', 'dip', and 'slope' 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.") + + if not (len(planesZseed) == len(planesDip) == len(planesSlope) == SHPdata["nFeatures"]): + raise ValueError( + "Mismatch between number of shapefile features and plane parameters in the .ini file." + ) + + for i in range(SHPdata["nFeatures"]): + xSeed = SHPdata["x"][int(SHPdata["Start"][i])] + ySeed = SHPdata["y"][int(SHPdata["Start"][i])] + zSeed = planesZseed[i] + dip = planesDip[i] + slopeAngle = planesSlope[i] + planeFeatures.extend([xSeed, ySeed, zSeed, dip, slopeAngle]) + + features = ",".join(map(str, planeFeatures)) + 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'])) + ellipsoidsTilt = list(map(float, SHPdata['tilt'])) + ellipsoidsDir = list(map(float, SHPdata['direc'])) + ellipsoidsOffset = list(map(float, SHPdata['offset'])) + ellipsoidDip = list(map(float, SHPdata['dip'])) + 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.") + + 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.") + + for i in range(SHPdata["nFeatures"]): + xCenter = SHPdata["x"][int(SHPdata["Start"][i])] + yCenter = SHPdata["y"][int(SHPdata["Start"][i])] + maxDepth = ellipsoidsMaxDepth[i] + semiMajor = ellipsoidsSemiMajor[i] + semiMinor = ellipsoidsSemiMinor[i] + tilt = ellipsoidsTilt[i] + direction = ellipsoidsDir[i] + offset = ellipsoidsOffset[i] + dip = ellipsoidDip[i] + ellipsoidFeatures.extend([xCenter, yCenter, maxDepth, semiMajor, semiMinor, tilt, direction, offset, dip]) + + features = ",".join(map(str, ellipsoidFeatures)) + + log.debug("Ellipsoid features extracted and combined: %s", features) + + if method == 'plane': + scarpData = calculateScarpWithPlanes( + dem["rasterData"], periData, dem["header"]["transform"], features + ) + elif method == 'ellipsoid': + scarpData = calculateScarpWithEllipsoids( + dem["rasterData"], periData, dem["header"]["transform"], features + ) + else: + raise ValueError("Unsupported method. Choose 'plane' or 'ellipsoid'.") + + hRelData = dem["rasterData"] - scarpData + + # create output directory and files + outDir = pathlib.Path(baseDir) + outDir = outDir / "Outputs" / "com6RockAvalanche" / "scarp" + fU.makeADir(outDir) + + # Set file names and write + elevationOut = outDir / "scarpElevation" + hRelOut = outDir / "scarpHRel" + IOf.writeResultToRaster(dem["header"], scarpData, elevationOut) + IOf.writeResultToRaster(dem["header"], hRelData, hRelOut) + +def readPerimeterSHP(perimeterShapefilePath, elevTransform, elevShape): + """Read perimeter from shapefile and return as numpy array. + + Parameters + ---------- + perimeterShapefilePath : str + Path to the perimeter shapefile. + elevTransform : Affine + transformation of the elevation raster. + elevShape : tuple + Shape (height, width) of the elevation raster. + + Returns + ------- + periData : np.ndarray + 2D array representing the perimeter mask. + """ + + log.debug("Processing perimeter shapefile: %s", perimeterShapefilePath) + SHPdata = SHP2Array(perimeterShapefilePath) + log.debug("Perimeter shapefile data loaded: %s", SHPdata) + + shapes = [] + for i in range(SHPdata["nFeatures"]): + start = int(SHPdata["Start"][i]) + length = int(SHPdata["Length"][i]) + coords = [(SHPdata["x"][j], SHPdata["y"][j]) for j in range(start, start + length)] + poly = shape({"type": "Polygon", "coordinates": [coords]}) + shapes.append(poly) + + periData = rasterize( + [(mapping(poly), 1) for poly in shapes], + out_shape=elevShape, + transform=elevTransform, + fill=0, + all_touched=True, + dtype=np.uint8, + ) + log.debug("Perimeter shapefile rasterized.") + + return periData + + +def calculateScarpWithPlanes(elevData, periData, elevTransform, planes): + """Calculate the scarp using sliding planes. + + Parameters + ---------- + elevData : np.ndarray + The elevation data as a 2D array. + periData : np.ndarray + The perimeter data as a 2D array, which defines the extent of the scarp. + 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). + + Returns + ------- + np.ndarray + A 2D array with the calculated scarp values. + """ + n, m = elevData.shape + scarpData = np.zeros_like(elevData, dtype=np.float32) + + planes = list(map(float, planes.split(','))) + nPlanes = int(len(planes) / 5) + + xSeed = [planes[0]] + ySeed = [planes[1]] + zSeed = [planes[2]] + 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]))] + + for i in range(1, nPlanes): + xSeed.append(planes[5 * i]) + ySeed.append(planes[5 * i + 1]) + 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]))) + + 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] + 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) + else: + scarpData[row, col] = elevData[row, col] + + return scarpData + +def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids): + """Calculate the scarp using tilted and offset ellipsoids. + + Parameters + ---------- + elevData : np.ndarray + The elevation data as a 2D array. + periData : np.ndarray + The perimeter data as a 2D array. + elevTransform : Affine + 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) + + Returns + ------- + np.ndarray + A 2D array with the calculated scarp values. + """ + + n, m = elevData.shape + scarpData = np.zeros_like(elevData, dtype=np.float32) + + # Parse ellipsoid definitions + ellipsoids = list(map(float, ellipsoids.split(','))) + nEllipsoids = int(len(ellipsoids) / 9) + + xCenter, yCenter, maxDepth = [], [], [] + semiMajor, semiMinor = [], [] + tilt, tiltDir, offset, dip = [], [], [], [] + + for i in range(nEllipsoids): + xCenter.append(ellipsoids[9 * i]) + yCenter.append(ellipsoids[9 * i + 1]) + maxDepth.append(ellipsoids[9 * i + 2]) + 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 + offset.append(ellipsoids[9 * i + 7]) + dip.append(np.radians(ellipsoids[9 * i + 8])) # rotation of base ellipse in radians + + # Compute slope direction and magnitude from DEM gradients + grad_y, grad_x = np.gradient(elevData, abs(elevTransform[4]), elevTransform[0]) + slopeDir = np.arctan2(-grad_y, grad_x) # Azimuth of steepest descent + slopeMagnitude = np.sqrt(grad_x**2 + grad_y**2) + + for row in range(n): + for col in range(m): + west, north = rasterio.transform.xy(elevTransform, row, col, offset='center') + scarpVal = elevData[row, col] + + for k in range(nEllipsoids): + center_col, center_row = ~elevTransform * (xCenter[k], yCenter[k]) + center_row = int(round(center_row)) + center_col = int(round(center_col)) + + if 0 <= center_row < n and 0 <= center_col < m: + slopeAzimuth = slopeDir[center_row, center_col] + slopeAzimuth = (slopeAzimuth + 2 * np.pi) % (2 * np.pi) + slopeMag = slopeMagnitude[center_row, center_col] + + if not np.isnan(slopeAzimuth) and slopeMag > 0: + normal_dx = np.cos(slopeAzimuth) + normal_dy = np.sin(slopeAzimuth) + slopeAngle = math.atan(slopeMag) + dxOffset = -offset[k] * normal_dx * math.sin(slopeAngle) + dyOffset = -offset[k] * normal_dy * math.sin(slopeAngle) + dzOffset = -offset[k] * math.cos(slopeAngle) + else: + dxOffset = dyOffset = dzOffset = 0 + else: + dxOffset = dyOffset = dzOffset = 0 + + x0 = xCenter[k] + dxOffset + y0 = yCenter[k] + dyOffset + z0 = dzOffset + + 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]) + + # Normalize to ellipsoid axes + xNorm = dxRot / semiMajor[k] + yNorm = dyRot / semiMinor[k] + distance = xNorm**2 + yNorm**2 + + if distance <= 1: + baseDepth = maxDepth[k] * (1 - distance) + tiltEffect = math.tan(math.radians(tilt[k])) * ( + dxRot * np.cos(tiltDir[k]) + dyRot * np.sin(tiltDir[k]) + ) + totalDepth = baseDepth + tiltEffect + z0 + scarpVal = min(scarpVal, elevData[row, col] - totalDepth) + + scarpData[row, col] = scarpVal if periData[row, col] > 0 else elevData[row, col] + + return scarpData diff --git a/avaframe/com6RockAvalanche/scarpCfg.ini b/avaframe/com6RockAvalanche/scarpCfg.ini new file mode 100644 index 000000000..d9efb72f0 --- /dev/null +++ b/avaframe/com6RockAvalanche/scarpCfg.ini @@ -0,0 +1,11 @@ +[INPUT] + +# Use shapefiles for feature coordinates and perimeter features +# They have to be saved as .shp file in Inputs/POINTS/*_coordinates.* +# and INPUT/POLYGONS/*_perimeter.* +useShapefiles = True + + +[SETTINGS] +# Choose the method: 'plane' or 'ellipsoid' +method = plane diff --git a/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.dbf b/avaframe/data/scarpExample/Inputs/POINTS/points_coordinates.dbf new file mode 100644 index 000000000..ef0170f04 Binary files /dev/null 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 new file mode 100644 index 000000000..0cc7e0016 Binary files /dev/null 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 new file mode 100644 index 000000000..adec32cf5 Binary files /dev/null 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 new file mode 100644 index 000000000..742c45010 Binary files /dev/null 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 new file mode 100644 index 000000000..ec5f2194b Binary files /dev/null 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 new file mode 100644 index 000000000..c2859ca3a Binary files /dev/null and b/avaframe/data/scarpExample/Inputs/POLYGONS/scarpFluchthorn_perimeter.shx differ diff --git a/avaframe/data/scarpExample/Inputs/fluchthorn.tif b/avaframe/data/scarpExample/Inputs/fluchthorn.tif new file mode 100644 index 000000000..7cc119c9d Binary files /dev/null and b/avaframe/data/scarpExample/Inputs/fluchthorn.tif differ diff --git a/avaframe/in1Data/getInput.py b/avaframe/in1Data/getInput.py index 710a1aa10..cd7b3611d 100644 --- a/avaframe/in1Data/getInput.py +++ b/avaframe/in1Data/getInput.py @@ -29,7 +29,7 @@ def readDEM(avaDir): - """read the ascii DEM file from a provided avalanche directory + """read the DEM (ascii or tif) file from a provided avalanche directory Parameters ---------- @@ -292,7 +292,7 @@ def getAndCheckInputFiles(inputDir, folder, inputType, fileExt="shp", fileSuffix folder : str subfolder name where the shape file should be located (SECREL, ENT or RES) inputType : str - type of input (used for the logging messages). Secondary release or Entrainment or Resistance + type of input (used for the logging messages). fileExt: str file extension e.g. shp, asc, tif - optional; default is shp fileSuffix: str @@ -340,7 +340,6 @@ def getAndCheckInputFiles(inputDir, folder, inputType, fileExt="shp", fileSuffix log.error(message) raise AssertionError(message) - return OutputFile, available diff --git a/avaframe/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index f86d2ad33..958453acc 100644 --- a/avaframe/in2Trans/shpConversion.py +++ b/avaframe/in2Trans/shpConversion.py @@ -56,7 +56,18 @@ def SHP2Array(infile, defname=None): list of parts of polygon (added the total number of points as list item, so if multiple parts len>2) nFeatures: int 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 + 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) + maxdepth + np array with the masimum depth of each scarp ellipsoid-feature (as many values as features) + semimajor + np array with the semi-major axis of each scarp ellipsoid-feature (as many values as features) """ # Input shapefile sf = shapefile.Reader(str(infile)) @@ -71,6 +82,16 @@ def SHP2Array(infile, defname=None): id = None ci95 = None layerN = None + zseed_value = None + dip_value = None + slopeAngle_value = None + semiminor_value = None + maxdepth_value = None + semimajor_value = None + tilt_value = None + direc_value = None + offset_value = None + # get coordinate system sks = getSHPProjection(infile) @@ -86,6 +107,17 @@ def SHP2Array(infile, defname=None): idList = [] ci95List = [] layerNameList = [] + zseedList = [] + dipList = [] + slopeList = [] + slopeAngleList = [] + semiminorList = [] + semimajorList = [] + maxdepthList = [] + tiltList = [] + direcList = [] + offsetList = [] + Length = np.empty((0)) Start = np.empty((0)) Coordx = np.empty((0)) @@ -125,6 +157,25 @@ def SHP2Array(infile, defname=None): if name == "slope": # for dams slope = value + if name == "slopeangle": + # for dams + slopeAngle_value = value + if name == "zseed": + zseed_value = value + if name == "dip": + dip_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 == "direc": + direc_value = value + if name == "offset": + offset_value = value if name == "rho": rho = value if name == "sks": @@ -149,6 +200,16 @@ def SHP2Array(infile, defname=None): ci95List.append(str(ci95)) layerNameList.append(layerN) idList.append(str(rec.oid)) + zseedList.append(zseed_value) + dipList.append(dip_value) + slopeList.append(slope) + slopeAngleList.append(slopeAngle_value) + semiminorList.append(semiminor_value) + semimajorList.append(semimajor_value) + maxdepthList.append(maxdepth_value) + tiltList.append(tilt_value) + direcList.append(direc_value) + offsetList.append(offset_value) Start = np.append(Start, start) length = len(pts) @@ -173,6 +234,15 @@ def SHP2Array(infile, defname=None): SHPdata["layerName"] = layerNameList SHPdata["nParts"] = nParts SHPdata["nFeatures"] = len(Start) + SHPdata["slopeangle"] = slopeAngleList + SHPdata["zseed"] = zseedList + SHPdata["dip"] = dipList + SHPdata["maxdepth"] = maxdepthList + SHPdata["semimajor"] = semimajorList + SHPdata["semiminor"] = semiminorList + SHPdata["tilt"] = tiltList + SHPdata["direc"] = direcList + SHPdata["offset"] = offsetList sf.close() diff --git a/avaframe/runCom6Scarp.py b/avaframe/runCom6Scarp.py new file mode 100644 index 000000000..92680325f --- /dev/null +++ b/avaframe/runCom6Scarp.py @@ -0,0 +1,67 @@ +""" + Run script for scarp analysis + In this runscript, the path to the scarp configuartion File has to be defined in the main config file +""" + +import pathlib +import argparse + +# Local imports +from avaframe.com6RockAvalanche import scarp +from avaframe.in3Utils import cfgUtils +from avaframe.in3Utils import logUtils +from avaframe.in3Utils import initializeProject as initProj + + +def runScarpAnalysisWorkflow(inputDir=""): + """Run the scarp analysis workflow. + + Parameters + ---------- + inputDir: str + Path to the input directory containing DEM, coordinates and perimeter files. + + Returns + ------- + None + """ + + logName = 'runScarpAnalysis' + + # Load general configuration file + cfgMain = cfgUtils.getGeneralConfig() + + # Load configuration file path from general config if not provided + if inputDir != "": + cfgMain["MAIN"]["avalancheDir"] = inputDir + else: + inputDir = cfgMain["MAIN"]["avalancheDir"] + + # Start logging + log = logUtils.initiateLogger(inputDir, logName) + log.info("MAIN SCRIPT") + log.info("Using inputDir: %s", inputDir) + + # ---------------- + # Clean input directory(ies) of old work files + initProj.cleanSingleAvaDir(inputDir, deleteOutput=False) + + # load scarp config + scarpCfg = cfgUtils.getModuleConfig(scarp) + + # Run the scarp analysis + scarp.scarpAnalysisMain(scarpCfg, str(inputDir)) + + log.info('Scarp analysis completed successfully.') + + return + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Run scarp analysis workflow') + parser.add_argument( + "inputDir", metavar="inputDir", type=str, nargs="?", default="", help="the input directory" + ) + + args = parser.parse_args() + runScarpAnalysisWorkflow(str(args.inputDir)) diff --git a/docs/moduleCom6RockAvalanche.rst b/docs/moduleCom6RockAvalanche.rst index f8f121233..97dd06dff 100644 --- a/docs/moduleCom6RockAvalanche.rst +++ b/docs/moduleCom6RockAvalanche.rst @@ -27,5 +27,50 @@ To run * run: :: - python3 runCom6RockAvalanche.py + python runCom6RockAvalanche.py + +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 +~~~~~ + +* 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. +* shapefile: 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", "dip" and "slopeangle" as float values. If the ellipsoid method is used, the shape file must + contain the Attributes "maxdepth", "semimajor" and "semiminor". The file must be located in avalancheDir/Inputs/ + POINTS and file name must end with “_coordinates”. +* perimeter_shapefile: A shapefile that specifies a boundary area. Must be located in avalancheDir/Inputs/POLYGONS and file name must end with “_perimeter”. + +Output +~~~~~~ + +* elevscarp: Output DGM (ASCII or GeoTIFF), which maps the input DGM minus the calculated scarp. Is saved under + ``raster_scarp.asc`` in ``avalancheDir/Outputs/com6RockAvalanche/scarp``. +* hrelease: File path to the output DGM (ASCII or GeoTIFF), which represents the calculated scarp volumes.Is saved + under ``raster_rel.asc`` in ``avalancheDir/Outputs/com6RockAvalanche/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 + +If all the data is provided successfully, start the script by running:: + + python runCom6Scarp.py \ No newline at end of file