diff --git a/avaframe/com6RockAvalanche/scarp.py b/avaframe/com6RockAvalanche/scarp.py new file mode 100644 index 000000000..612d56380 --- /dev/null +++ b/avaframe/com6RockAvalanche/scarp.py @@ -0,0 +1,316 @@ +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'])) + except KeyError as e: + raise ValueError(f"Required attribute '{e.args[0]}' not found in shapefile. Ensure the fields 'maxdepth', 'semimajor', and 'semiminor' exist.") + + if not ( + len(ellipsoidsMaxDepth) + == len(ellipsoidsSemiMajor) + == len(ellipsoidsSemiMinor) + == SHPdata["nFeatures"] + ): + raise ValueError( + "Mismatch between number of shapefile features and ellipsoid parameters in the .ini file." + ) + + 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] + ellipsoidFeatures.extend([xCenter, yCenter, maxDepth, semiMajor, semiMinor]) + + 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 sliding circles (rotational ellipsoids). + + 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). + ellipsoids : str + Comma-separated string defining ellipsoids (x_center, y_center, max_depth, semi_major_axis, semi_minor_axis). + + Returns + ------- + np.ndarray + A 2D array with the calculated scarp values. + """ + n, m = elevData.shape + scarpData = np.zeros_like(elevData, dtype=np.float32) + + ellipsoids = list(map(float, ellipsoids.split(','))) + nEllipsoids = int(len(ellipsoids) / 5) + + xCenter = [ellipsoids[0]] + yCenter = [ellipsoids[1]] + maxDepth = [ellipsoids[2]] + semiMajor = [ellipsoids[3]] + semiMinor = [ellipsoids[4]] + + for i in range(1, nEllipsoids): + xCenter.append(ellipsoids[5 * i]) + yCenter.append(ellipsoids[5 * i + 1]) + maxDepth.append(ellipsoids[5 * i + 2]) + semiMajor.append(ellipsoids[5 * i + 3]) + semiMinor.append(ellipsoids[5 * i + 4]) + + 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): + x2 = ((west - xCenter[k]) / semiMajor[k]) ** 2 + y2 = ((north - yCenter[k]) / semiMinor[k]) ** 2 + if x2 + y2 < 1: + scarpVal = min(scarpVal, elevData[row, col] - maxDepth[k] * (1 - (x2 + y2))) + if periData[row, col] > 0: + scarpData[row, col] = scarpVal + else: + scarpData[row, col] = 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 d23ea43d9..b18f39d5f 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,13 @@ 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 + # get coordinate system sks = getSHPProjection(infile) @@ -86,6 +104,13 @@ def SHP2Array(infile, defname=None): idList = [] ci95List = [] layerNameList = [] + zseedList = [] + dipList = [] + slopeList = [] + slopeAngleList = [] + semiminorList = [] + semimajorList = [] + maxdepthList = [] Length = np.empty((0)) Start = np.empty((0)) Coordx = np.empty((0)) @@ -125,6 +150,19 @@ 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 == "rho": rho = value if name == "sks": @@ -149,6 +187,13 @@ 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) Start = np.append(Start, start) length = len(pts) @@ -173,6 +218,12 @@ 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 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