From ed6bc89f93da12235ae302e4a324e9a1701b7606 Mon Sep 17 00:00:00 2001 From: dwolfsch <168710278+dwolfsch@users.noreply.github.com> Date: Thu, 23 Jan 2025 11:50:13 +0100 Subject: [PATCH 1/5] spatial Voellmy friction model: add VariableVoellmyShapeToRaster to generate rasters from shapes Adds functionality for generating raster files for mu and xsi values from a DEM and shapefiles including: A main script to rasterize mu and xsi based on polygon shapefiles and attribute fields. A corresponding run script and configuration file for user customization. --- avaframe/avaframeCfg.ini | 2 +- .../readMeVariableVoellmyShapeToRaster.txt | 23 +++++ .../variableVoellmyShapeToRaster.py | 86 +++++++++++++++++++ .../variableVoellmyShapeToRasterCfg.ini | 6 ++ avaframe/in2Trans/shpConversion.py | 13 ++- avaframe/runVariableVoellmyShapeToRaster.py | 60 +++++++++++++ 6 files changed, 188 insertions(+), 2 deletions(-) create mode 100644 avaframe/com6RockAvalanche/readMeVariableVoellmyShapeToRaster.txt create mode 100644 avaframe/com6RockAvalanche/variableVoellmyShapeToRaster.py create mode 100644 avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini create mode 100644 avaframe/runVariableVoellmyShapeToRaster.py diff --git a/avaframe/avaframeCfg.ini b/avaframe/avaframeCfg.ini index 4572b2567..ea96547ff 100644 --- a/avaframe/avaframeCfg.ini +++ b/avaframe/avaframeCfg.ini @@ -4,7 +4,7 @@ [MAIN] # Path to avalanche directory -avalancheDir = data/avaParabola +avalancheDir = # number of CPU cores to use for the computation of com1DFA # possible values are: diff --git a/avaframe/com6RockAvalanche/readMeVariableVoellmyShapeToRaster.txt b/avaframe/com6RockAvalanche/readMeVariableVoellmyShapeToRaster.txt new file mode 100644 index 000000000..761fe6719 --- /dev/null +++ b/avaframe/com6RockAvalanche/readMeVariableVoellmyShapeToRaster.txt @@ -0,0 +1,23 @@ +The VariableVoellmyShapeToRaster.py script allows the user to define spatially different values for the voellmy Parameters mu and xsi, with the use of polygon shapefiles. For the extent of a DEM raster, all the areas that are not covered by a polygon get assigned a default mu or xsi value. The script then converts this Information into a raster mu and a raster xsi file, which can then be used in Avaframe Simulation runs, using the "spatialVoellmy" friction model. + +First, set up the Config File and provide inputs: +•Config File: + oIn the first step, the Config File needs to be configured and all input files have to be provided + Main Config (avaframeCfg.ini): + •Set the path to the avalanche directory +oInputs: + All the Input Files are automatically fetched through the set avalanche directory. It is not necessary to provide a file path. + dem: DEM Raster that is later on used for the avaframe simulation. This is needed, because the mu and xsi output rasters need to be the exact same size. Has to lie in avadir/Inputs. + mu_shapefile: Mu shapefile, that is then converted to a raster file. Be aware, that the attribute has to be named “mu” and the file name has to end with “_mu”. Has to lie in avadir/Inputs/POLYGONS. + xsi_shapefile: Xsi shapefile, that is then converted to a raster file. Be aware, that the attribute has to be named “xsi” and the file name has to end with “_xsi”. Has to lie in avadir/Inputs/POLYGONS. +oDefaults: + default_mu: this is the default mu value, that gets assigned to all areas in the raster, that are not covered by shapefile-polygons + default_xsi: this is the default xsi value, that gets assigned to all areas in the raster, that are not covered by shapefile-polygons +oOutputs: + For the variable Voellmy calculations in the com1DFA algorithm to work, it is mandatory, that the files are stored in: avaframe\data\*yourAvalancheDir*\Inputs\RASTERS\ + mu_raster: Output for the generated mu raster file stored as *_mu.asc + xsi_raster: Output for the generated xsi raster file stored as *_xi.asc + +•RunScript: + oOnce everything is set up, run the script “runVariableVoellmyShapeToRaster.py” + oIf libraries are missing use: pip install *name of missing library diff --git a/avaframe/com6RockAvalanche/variableVoellmyShapeToRaster.py b/avaframe/com6RockAvalanche/variableVoellmyShapeToRaster.py new file mode 100644 index 000000000..79337f589 --- /dev/null +++ b/avaframe/com6RockAvalanche/variableVoellmyShapeToRaster.py @@ -0,0 +1,86 @@ + +import rasterio +import numpy as np +import pathlib +from rasterio.features import rasterize +from shapely.geometry import shape, mapping +from in2Trans.shpConversion import SHP2Array +from in1Data.getInput import getAndCheckInputFiles +import logging + +# Configure logging +logging.basicConfig(level=logging.DEBUG) +log = logging.getLogger(__name__) + +def generateMuXsiRasters(avadir, variableVoellmyCfg): + """ + Generate raster files for \u03bc and \u03be based on input DEM and shapefiles. + + Parameters + ---------- + avadir : str + Path to the avalanche directory. + variableVoellmyCfg : Config Parser Object + variableVoellmyCfg Configuration File + + Returns + ------- + None + """ + avadir = pathlib.Path(avadir) + + config = variableVoellmyCfg # Directly use the ConfigParser object + + inputDir = avadir / "Inputs" + outputDir = avadir / "Inputs" # Output directory is Inputs, because Outputs of this Script will be used as Inputs for AvaFrame + + demPath, _ = getAndCheckInputFiles(inputDir, '', 'DEM', fileExt='asc') + muShapefile, _ = getAndCheckInputFiles(inputDir, 'POLYGONS', '\u03bc Shapefile', fileExt='shp', fileSuffix='_mu') + xsiShapefile, _ = getAndCheckInputFiles(inputDir, 'POLYGONS', '\u03be Shapefile', fileExt='shp', fileSuffix='_xsi') + + muOutputPath = outputDir / "RASTERS" / "raster_mu.asc" + xsiOutputPath = outputDir / "RASTERS" /"raster_xi.asc" + + defaultMu = float(config['DEFAULTS']['default_mu']) + defaultXsi = float(config['DEFAULTS']['default_xsi']) + + # Read DEM + with rasterio.open(demPath) as demSrc: + demData = demSrc.read(1) + demTransform = demSrc.transform + demCrs = demSrc.crs + demShape = demData.shape + + def rasterizeShapefile(shapefilePath, defaultValue, attributeName): + if not shapefilePath: + return np.full(demShape, defaultValue, dtype=np.float32) + + shpData = SHP2Array(shapefilePath) + 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]}) + value = shpData['attributes'][i][attributeName] + shapes.append((mapping(poly), value)) + + return rasterize(shapes, out_shape=demShape, transform=demTransform, fill=defaultValue, all_touched=True, dtype=np.float32) + + log.info("Rasterizing \u03bc shapefile.") + muRaster = rasterizeShapefile(muShapefile, defaultMu, "mu") + + log.info("Rasterizing \u03be shapefile.") + xsiRaster = rasterizeShapefile(xsiShapefile, defaultXsi, "xsi") + + def saveRaster(outputPath, data): + with rasterio.open(outputPath, 'w', driver='GTiff', height=data.shape[0], width=data.shape[1], count=1, dtype=data.dtype, crs=demCrs, transform=demTransform) as dst: + dst.write(data, 1) + + log.info("Saving \u03bc raster to %s", muOutputPath) + saveRaster(muOutputPath, muRaster) + + log.info("Saving \u03be raster to %s", xsiOutputPath) + saveRaster(xsiOutputPath, xsiRaster) + + log.info("Raster generation completed.") diff --git a/avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini b/avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini new file mode 100644 index 000000000..8ef1f0139 --- /dev/null +++ b/avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini @@ -0,0 +1,6 @@ +[DEFAULTS] +# Default \u03bc value for areas not covered by shapefiles +default_mu = 0.1 + +# Default \u03be value for areas not covered by shapefiles +default_xsi = 300.0 diff --git a/avaframe/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index 6eef2d743..1d95418e3 100644 --- a/avaframe/in2Trans/shpConversion.py +++ b/avaframe/in2Trans/shpConversion.py @@ -126,6 +126,9 @@ def SHP2Array(infile, defname=None): start = 0 nParts = [] + # New: Create an empty list to store attributes + attributes = [] + for n, (item, rec) in enumerate(zip(shps, sf.records())): pts = item.points # if feature has no points - ignore @@ -145,9 +148,14 @@ def SHP2Array(infile, defname=None): # check if records are available and extract if records: # loop through fields + # Extract attributes for the feature + attr_dict = {} for (name, typ, size, deci), value in zip(sf.fields[1:], records[n].record): # get entity name name = name.lower() + attr_dict[name] = value # Store attributes in dictionary + + # Specific field handling (existing code) if name == "name": layername = str(value) if (name == "thickness") or (name == "d0"): @@ -183,13 +191,15 @@ def SHP2Array(infile, defname=None): if name == "iso": iso = value if name == "layer": - layerN = value + layerN = value # if name is still empty go through file again and take Layer instead if (type(layername) is bytes) or (layername is None): for (name, typ, size, deci), value in zip(sf.fields[1:], records[n].record): if name == "Layer": layername = value + attributes.append(attr_dict) # Add the attribute dictionary to the list + # if layer still not defined, use generic if layername is None: layername = defname @@ -243,6 +253,7 @@ def SHP2Array(infile, defname=None): SHPdata["rotAngle"] = rotAngleList SHPdata["direc"] = direcList SHPdata["offset"] = offsetList + SHPdata["attributes"] = attributes # Add attributes to SHPdata sf.close() diff --git a/avaframe/runVariableVoellmyShapeToRaster.py b/avaframe/runVariableVoellmyShapeToRaster.py new file mode 100644 index 000000000..6f97a6eb2 --- /dev/null +++ b/avaframe/runVariableVoellmyShapeToRaster.py @@ -0,0 +1,60 @@ + +import argparse +import pathlib +import time +from avaframe.in3Utils import cfgUtils +from avaframe.in3Utils import logUtils +import avaframe.in3Utils.initializeProject as initProj +from com6RockAvalanche import variableVoellmyShapeToRaster +from com6RockAvalanche.variableVoellmyShapeToRaster import generateMuXsiRasters + +def runMuXsiWorkflow(avadir=''): + """ + Run the workflow to generate \u03bc and \u03be rasters. + + Parameters + ---------- + avadir : str + Path to the avalanche directory containing input and output folders. + + Returns + ------- + None + """ + startTime = time.time() + logName = 'runMuXsi' + + # Load general configuration file + cfgMain = cfgUtils.getGeneralConfig() + if avadir: + cfgMain['MAIN']['avalancheDir'] = avadir + else: + avadir = cfgMain['MAIN']['avalancheDir'] + + avadir = pathlib.Path(avadir) + + # Start logging + log = logUtils.initiateLogger(avadir, logName) + log.info('MAIN SCRIPT') + log.info('Using avalanche directory: %s', avadir) + + # Clean input directory(ies) of old work files + initProj.cleanSingleAvaDir(avadir, deleteOutput=False) + + # Load module-specific configuration for Variable Voellmy + variableVoellmyCfg = cfgUtils.getModuleConfig(variableVoellmyShapeToRaster) + + # Run the raster generation process + generateMuXsiRasters(avadir, variableVoellmyCfg) + + endTime = time.time() + log.info("Took %6.1f seconds to calculate.", (endTime - startTime)) + log.info('Workflow completed successfully.') + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Run \u03bc and \u03be raster generation workflow') + parser.add_argument('avadir', metavar='a', type=str, nargs='?', default='', + help='Path to the avalanche directory') + + args = parser.parse_args() + runMuXsiWorkflow(str(args.avadir)) From ac487461dcce959183209d9c9446af526bb806fb Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Mon, 17 Nov 2025 15:06:37 +0100 Subject: [PATCH 2/5] refactor(config): update default configuration values and comments - Updated `variableVoellmyShapeToRasterCfg.ini` to use `mu` and `xsi` instead of Unicode characters for clarity. - Modified `avaframeCfg.ini` to set a default path for `avalancheDir`. --- avaframe/avaframeCfg.ini | 2 +- .../com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/avaframe/avaframeCfg.ini b/avaframe/avaframeCfg.ini index ea96547ff..4572b2567 100644 --- a/avaframe/avaframeCfg.ini +++ b/avaframe/avaframeCfg.ini @@ -4,7 +4,7 @@ [MAIN] # Path to avalanche directory -avalancheDir = +avalancheDir = data/avaParabola # number of CPU cores to use for the computation of com1DFA # possible values are: diff --git a/avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini b/avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini index 8ef1f0139..ff2fea660 100644 --- a/avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini +++ b/avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini @@ -1,6 +1,6 @@ [DEFAULTS] -# Default \u03bc value for areas not covered by shapefiles +# Default mu value for areas not covered by shapefiles default_mu = 0.1 -# Default \u03be value for areas not covered by shapefiles +# Default xsi value for areas not covered by shapefiles default_xsi = 300.0 From bcf635804696d326948e2e830ed28f32a24240c3 Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Tue, 16 Jun 2026 09:55:43 +0200 Subject: [PATCH 3/5] refactor(spatialVoellmyInputs): replace and relocate `VariableVoellmyShapeToRaster` workflow - Removed `variableVoellmyShapeToRaster` scripts and related configurations. - Introduced `spatialVoellmyInputs` for handling spatial Voellmy friction parameter generation. - Migrated default input values to `spatialVoellmyInputsCfg.ini` with updated defaults. - Replaced inline raster generation logic with standardized utility functions. --- .../readMeVariableVoellmyShapeToRaster.txt | 23 -- .../variableVoellmyShapeToRaster.py | 86 -------- avaframe/in2Trans/shpConversion.py | 13 +- avaframe/in3Utils/spatialVoellmyInputs.py | 136 ++++++++++++ .../spatialVoellmyInputsCfg.ini} | 4 +- .../runScripts/runSpatialVoellmyInputs.py | 55 +++++ avaframe/runVariableVoellmyShapeToRaster.py | 60 ------ avaframe/tests/test_spatialVoellmyInputs.py | 204 ++++++++++++++++++ 8 files changed, 398 insertions(+), 183 deletions(-) delete mode 100644 avaframe/com6RockAvalanche/readMeVariableVoellmyShapeToRaster.txt delete mode 100644 avaframe/com6RockAvalanche/variableVoellmyShapeToRaster.py create mode 100644 avaframe/in3Utils/spatialVoellmyInputs.py rename avaframe/{com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini => in3Utils/spatialVoellmyInputsCfg.ini} (75%) create mode 100644 avaframe/runScripts/runSpatialVoellmyInputs.py delete mode 100644 avaframe/runVariableVoellmyShapeToRaster.py create mode 100644 avaframe/tests/test_spatialVoellmyInputs.py diff --git a/avaframe/com6RockAvalanche/readMeVariableVoellmyShapeToRaster.txt b/avaframe/com6RockAvalanche/readMeVariableVoellmyShapeToRaster.txt deleted file mode 100644 index 761fe6719..000000000 --- a/avaframe/com6RockAvalanche/readMeVariableVoellmyShapeToRaster.txt +++ /dev/null @@ -1,23 +0,0 @@ -The VariableVoellmyShapeToRaster.py script allows the user to define spatially different values for the voellmy Parameters mu and xsi, with the use of polygon shapefiles. For the extent of a DEM raster, all the areas that are not covered by a polygon get assigned a default mu or xsi value. The script then converts this Information into a raster mu and a raster xsi file, which can then be used in Avaframe Simulation runs, using the "spatialVoellmy" friction model. - -First, set up the Config File and provide inputs: -•Config File: - oIn the first step, the Config File needs to be configured and all input files have to be provided - Main Config (avaframeCfg.ini): - •Set the path to the avalanche directory -oInputs: - All the Input Files are automatically fetched through the set avalanche directory. It is not necessary to provide a file path. - dem: DEM Raster that is later on used for the avaframe simulation. This is needed, because the mu and xsi output rasters need to be the exact same size. Has to lie in avadir/Inputs. - mu_shapefile: Mu shapefile, that is then converted to a raster file. Be aware, that the attribute has to be named “mu” and the file name has to end with “_mu”. Has to lie in avadir/Inputs/POLYGONS. - xsi_shapefile: Xsi shapefile, that is then converted to a raster file. Be aware, that the attribute has to be named “xsi” and the file name has to end with “_xsi”. Has to lie in avadir/Inputs/POLYGONS. -oDefaults: - default_mu: this is the default mu value, that gets assigned to all areas in the raster, that are not covered by shapefile-polygons - default_xsi: this is the default xsi value, that gets assigned to all areas in the raster, that are not covered by shapefile-polygons -oOutputs: - For the variable Voellmy calculations in the com1DFA algorithm to work, it is mandatory, that the files are stored in: avaframe\data\*yourAvalancheDir*\Inputs\RASTERS\ - mu_raster: Output for the generated mu raster file stored as *_mu.asc - xsi_raster: Output for the generated xsi raster file stored as *_xi.asc - -•RunScript: - oOnce everything is set up, run the script “runVariableVoellmyShapeToRaster.py” - oIf libraries are missing use: pip install *name of missing library diff --git a/avaframe/com6RockAvalanche/variableVoellmyShapeToRaster.py b/avaframe/com6RockAvalanche/variableVoellmyShapeToRaster.py deleted file mode 100644 index 79337f589..000000000 --- a/avaframe/com6RockAvalanche/variableVoellmyShapeToRaster.py +++ /dev/null @@ -1,86 +0,0 @@ - -import rasterio -import numpy as np -import pathlib -from rasterio.features import rasterize -from shapely.geometry import shape, mapping -from in2Trans.shpConversion import SHP2Array -from in1Data.getInput import getAndCheckInputFiles -import logging - -# Configure logging -logging.basicConfig(level=logging.DEBUG) -log = logging.getLogger(__name__) - -def generateMuXsiRasters(avadir, variableVoellmyCfg): - """ - Generate raster files for \u03bc and \u03be based on input DEM and shapefiles. - - Parameters - ---------- - avadir : str - Path to the avalanche directory. - variableVoellmyCfg : Config Parser Object - variableVoellmyCfg Configuration File - - Returns - ------- - None - """ - avadir = pathlib.Path(avadir) - - config = variableVoellmyCfg # Directly use the ConfigParser object - - inputDir = avadir / "Inputs" - outputDir = avadir / "Inputs" # Output directory is Inputs, because Outputs of this Script will be used as Inputs for AvaFrame - - demPath, _ = getAndCheckInputFiles(inputDir, '', 'DEM', fileExt='asc') - muShapefile, _ = getAndCheckInputFiles(inputDir, 'POLYGONS', '\u03bc Shapefile', fileExt='shp', fileSuffix='_mu') - xsiShapefile, _ = getAndCheckInputFiles(inputDir, 'POLYGONS', '\u03be Shapefile', fileExt='shp', fileSuffix='_xsi') - - muOutputPath = outputDir / "RASTERS" / "raster_mu.asc" - xsiOutputPath = outputDir / "RASTERS" /"raster_xi.asc" - - defaultMu = float(config['DEFAULTS']['default_mu']) - defaultXsi = float(config['DEFAULTS']['default_xsi']) - - # Read DEM - with rasterio.open(demPath) as demSrc: - demData = demSrc.read(1) - demTransform = demSrc.transform - demCrs = demSrc.crs - demShape = demData.shape - - def rasterizeShapefile(shapefilePath, defaultValue, attributeName): - if not shapefilePath: - return np.full(demShape, defaultValue, dtype=np.float32) - - shpData = SHP2Array(shapefilePath) - 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]}) - value = shpData['attributes'][i][attributeName] - shapes.append((mapping(poly), value)) - - return rasterize(shapes, out_shape=demShape, transform=demTransform, fill=defaultValue, all_touched=True, dtype=np.float32) - - log.info("Rasterizing \u03bc shapefile.") - muRaster = rasterizeShapefile(muShapefile, defaultMu, "mu") - - log.info("Rasterizing \u03be shapefile.") - xsiRaster = rasterizeShapefile(xsiShapefile, defaultXsi, "xsi") - - def saveRaster(outputPath, data): - with rasterio.open(outputPath, 'w', driver='GTiff', height=data.shape[0], width=data.shape[1], count=1, dtype=data.dtype, crs=demCrs, transform=demTransform) as dst: - dst.write(data, 1) - - log.info("Saving \u03bc raster to %s", muOutputPath) - saveRaster(muOutputPath, muRaster) - - log.info("Saving \u03be raster to %s", xsiOutputPath) - saveRaster(xsiOutputPath, xsiRaster) - - log.info("Raster generation completed.") diff --git a/avaframe/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index 1d95418e3..6eef2d743 100644 --- a/avaframe/in2Trans/shpConversion.py +++ b/avaframe/in2Trans/shpConversion.py @@ -126,9 +126,6 @@ def SHP2Array(infile, defname=None): start = 0 nParts = [] - # New: Create an empty list to store attributes - attributes = [] - for n, (item, rec) in enumerate(zip(shps, sf.records())): pts = item.points # if feature has no points - ignore @@ -148,14 +145,9 @@ def SHP2Array(infile, defname=None): # check if records are available and extract if records: # loop through fields - # Extract attributes for the feature - attr_dict = {} for (name, typ, size, deci), value in zip(sf.fields[1:], records[n].record): # get entity name name = name.lower() - attr_dict[name] = value # Store attributes in dictionary - - # Specific field handling (existing code) if name == "name": layername = str(value) if (name == "thickness") or (name == "d0"): @@ -191,15 +183,13 @@ def SHP2Array(infile, defname=None): if name == "iso": iso = value if name == "layer": - layerN = value + layerN = value # if name is still empty go through file again and take Layer instead if (type(layername) is bytes) or (layername is None): for (name, typ, size, deci), value in zip(sf.fields[1:], records[n].record): if name == "Layer": layername = value - attributes.append(attr_dict) # Add the attribute dictionary to the list - # if layer still not defined, use generic if layername is None: layername = defname @@ -253,7 +243,6 @@ def SHP2Array(infile, defname=None): SHPdata["rotAngle"] = rotAngleList SHPdata["direc"] = direcList SHPdata["offset"] = offsetList - SHPdata["attributes"] = attributes # Add attributes to SHPdata sf.close() diff --git a/avaframe/in3Utils/spatialVoellmyInputs.py b/avaframe/in3Utils/spatialVoellmyInputs.py new file mode 100644 index 000000000..1b1d4b83c --- /dev/null +++ b/avaframe/in3Utils/spatialVoellmyInputs.py @@ -0,0 +1,136 @@ +""" +Functions for generating spatial Voellmy friction raster inputs. +""" + +import pathlib +import logging +import numpy as np +import shapefile +from rasterio.features import rasterize +from shapely.geometry import shape, mapping + +from avaframe.in1Data.getInput import getAndCheckInputFiles, getDEMPath +from avaframe.in2Trans.rasterUtils import readRasterHeader, writeResultToRaster + +log = logging.getLogger(__name__) + + +def generateMuXsiRasters(avaDir, cfg): + """Generate mu and xi raster files from polygon shapefiles. + + Reads polygon shapefiles with "mu" and "xsi" attribute fields, + rasterizes the attribute values onto a grid matching the DEM extent + and resolution, and writes the rasters to Inputs/RASTERS/. + + Parameters + ---------- + avaDir : pathlib.Path + Path to avalanche directory containing Inputs/DEM and + Inputs/POLYGONS/ with *_mu.shp and *_xsi.shp shapefiles. + cfg : configparser.ConfigParser + Configuration with [DEFAULTS] section containing + default_mu and default_xsi values for uncovered areas. + """ + avaDir = pathlib.Path(avaDir) + inputDir = avaDir / "Inputs" + outDir = inputDir / "RASTERS" + outDir.mkdir(parents=True, exist_ok=True) + + # Find DEM + demPath = getDEMPath(avaDir) + demSuffix = demPath.suffix + + # Find shapefiles + muShpPath, muAvailable, _ = getAndCheckInputFiles( + inputDir, "POLYGONS", "mu shapefile", fileExt="shp", fileSuffix="_mu" + ) + if muAvailable == "No": + raise FileNotFoundError("No *_mu.shp found in %s/POLYGONS/" % inputDir) + xsiShpPath, xsiAvailable, _ = getAndCheckInputFiles( + inputDir, "POLYGONS", "xsi shapefile", fileExt="shp", fileSuffix="_xsi" + ) + if xsiAvailable == "No": + raise FileNotFoundError("No *_xsi.shp found in %s/POLYGONS/" % inputDir) + + # Read DEM header + demHeader = readRasterHeader(demPath) + demTransform = demHeader["transform"] + demCrs = demHeader["crs"] + demShape = (demHeader["nrows"], demHeader["ncols"]) + + defaultMu = float(cfg["DEFAULTS"]["default_mu"]) + defaultXsi = float(cfg["DEFAULTS"]["default_xsi"]) + + # Rasterize mu + log.info("Rasterizing mu shapefile: %s", muShpPath) + muRaster = _rasterizeShapefile(muShpPath, defaultMu, "mu", demShape, demTransform) + + # Rasterize xsi + log.info("Rasterizing xsi shapefile: %s", xsiShpPath) + xsiRaster = _rasterizeShapefile(xsiShpPath, defaultXsi, "xsi", demShape, demTransform) + + # Determine output driver + if demSuffix == ".asc": + driver = "AAIGrid" + else: + driver = "GTiff" + + # Write output + outHeader = { + "driver": driver, + "crs": demCrs, + "transform": demTransform, + "nodata_value": None, + } + log.info("Writing mu raster") + writeResultToRaster(outHeader, muRaster, outDir / "raster_mu") + log.info("Writing xsi raster") + writeResultToRaster(outHeader, xsiRaster, outDir / "raster_xi") + log.info("Raster generation completed.") + + +def _rasterizeShapefile(shpPath, defaultValue, fieldName, demShape, demTransform): + """Rasterize a polygon shapefile attribute field onto a DEM-matching grid. + + Parameters + ---------- + shpPath : pathlib.Path + Path to shapefile. + defaultValue : float + Fill value for cells not covered by any polygon. + fieldName : str + Attribute field name to extract from each feature. + demShape : tuple + (height, width) of the output raster. + demTransform : affine.Affine + Geotransform of the DEM. + + Returns + ------- + raster : numpy.ndarray + Rasterized array with shape demShape. + """ + with shapefile.Reader(str(shpPath)) as sf: + fieldNames = [f[0].lower() for f in sf.fields[1:]] + if fieldName not in fieldNames: + raise KeyError( + "Field '%s' not found in %s. Available fields: %s" % (fieldName, shpPath.name, fieldNames) + ) + fieldIdx = fieldNames.index(fieldName) + + shapes = [] + for rec in sf.shapeRecords(): + geom = rec.shape.__geo_interface__ + poly = shape(geom) + value = rec.record[fieldIdx] + shapes.append((mapping(poly), value)) + + raster = rasterize( + shapes, + out_shape=demShape, + transform=demTransform, + fill=defaultValue, + all_touched=True, + dtype=np.float32, + ) + return raster diff --git a/avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini b/avaframe/in3Utils/spatialVoellmyInputsCfg.ini similarity index 75% rename from avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini rename to avaframe/in3Utils/spatialVoellmyInputsCfg.ini index ff2fea660..1fb927e33 100644 --- a/avaframe/com6RockAvalanche/variableVoellmyShapeToRasterCfg.ini +++ b/avaframe/in3Utils/spatialVoellmyInputsCfg.ini @@ -1,6 +1,6 @@ [DEFAULTS] # Default mu value for areas not covered by shapefiles -default_mu = 0.1 +default_mu = 0.155 # Default xsi value for areas not covered by shapefiles -default_xsi = 300.0 +default_xsi = 4000. diff --git a/avaframe/runScripts/runSpatialVoellmyInputs.py b/avaframe/runScripts/runSpatialVoellmyInputs.py new file mode 100644 index 000000000..f02c89ad8 --- /dev/null +++ b/avaframe/runScripts/runSpatialVoellmyInputs.py @@ -0,0 +1,55 @@ +""" +Run script for generating spatial Voellmy friction raster inputs. +""" + +import argparse +import pathlib +import time + +from avaframe.in3Utils import cfgUtils +from avaframe.in3Utils import logUtils +from avaframe.in3Utils import initializeProject as initProj +from avaframe.in3Utils import spatialVoellmyInputs + + +def runSpatialVoellmyInputs(avaDir=""): + """Run generation of mu and xi rasters from shapefiles. + + Parameters + ---------- + avaDir : str + Path to the avalanche directory. If empty, read from + avaframeCfg.ini. + """ + startTime = time.time() + cfgMain = cfgUtils.getGeneralConfig() + if avaDir: + avaDir = pathlib.Path(avaDir) + else: + avaDir = pathlib.Path(cfgMain["MAIN"]["avalancheDir"]) + + logName = "runSpatialVoellmyInputs" + log = logUtils.initiateLogger(avaDir, logName) + log.info("MAIN SCRIPT") + log.info("Avalanche directory: %s", avaDir) + + initProj.cleanSingleAvaDir(avaDir, deleteOutput=False) + + cfg = cfgUtils.getModuleConfig(spatialVoellmyInputs) + spatialVoellmyInputs.generateMuXsiRasters(avaDir, cfg) + + endTime = time.time() + log.info("Took %6.1f seconds to calculate.", endTime - startTime) + log.info("Workflow completed successfully.") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="Generate spatial Voellmy mu/xi rasters from shapefiles") + parser.add_argument( + "avaDir", + nargs="?", + default="", + help="Path to avalanche directory", + ) + args = parser.parse_args() + runSpatialVoellmyInputs(str(args.avaDir)) diff --git a/avaframe/runVariableVoellmyShapeToRaster.py b/avaframe/runVariableVoellmyShapeToRaster.py deleted file mode 100644 index 6f97a6eb2..000000000 --- a/avaframe/runVariableVoellmyShapeToRaster.py +++ /dev/null @@ -1,60 +0,0 @@ - -import argparse -import pathlib -import time -from avaframe.in3Utils import cfgUtils -from avaframe.in3Utils import logUtils -import avaframe.in3Utils.initializeProject as initProj -from com6RockAvalanche import variableVoellmyShapeToRaster -from com6RockAvalanche.variableVoellmyShapeToRaster import generateMuXsiRasters - -def runMuXsiWorkflow(avadir=''): - """ - Run the workflow to generate \u03bc and \u03be rasters. - - Parameters - ---------- - avadir : str - Path to the avalanche directory containing input and output folders. - - Returns - ------- - None - """ - startTime = time.time() - logName = 'runMuXsi' - - # Load general configuration file - cfgMain = cfgUtils.getGeneralConfig() - if avadir: - cfgMain['MAIN']['avalancheDir'] = avadir - else: - avadir = cfgMain['MAIN']['avalancheDir'] - - avadir = pathlib.Path(avadir) - - # Start logging - log = logUtils.initiateLogger(avadir, logName) - log.info('MAIN SCRIPT') - log.info('Using avalanche directory: %s', avadir) - - # Clean input directory(ies) of old work files - initProj.cleanSingleAvaDir(avadir, deleteOutput=False) - - # Load module-specific configuration for Variable Voellmy - variableVoellmyCfg = cfgUtils.getModuleConfig(variableVoellmyShapeToRaster) - - # Run the raster generation process - generateMuXsiRasters(avadir, variableVoellmyCfg) - - endTime = time.time() - log.info("Took %6.1f seconds to calculate.", (endTime - startTime)) - log.info('Workflow completed successfully.') - -if __name__ == '__main__': - parser = argparse.ArgumentParser(description='Run \u03bc and \u03be raster generation workflow') - parser.add_argument('avadir', metavar='a', type=str, nargs='?', default='', - help='Path to the avalanche directory') - - args = parser.parse_args() - runMuXsiWorkflow(str(args.avadir)) diff --git a/avaframe/tests/test_spatialVoellmyInputs.py b/avaframe/tests/test_spatialVoellmyInputs.py new file mode 100644 index 000000000..6b6f0b9ae --- /dev/null +++ b/avaframe/tests/test_spatialVoellmyInputs.py @@ -0,0 +1,204 @@ +"""Tests for module spatialVoellmyInputs""" + +import pathlib +import tempfile +import shutil +import numpy as np +import rasterio +import shapefile +import pytest +from avaframe.in3Utils import spatialVoellmyInputs + + +def _makeSyntheticDEM(tmpDir, suffix=".asc"): + """Create a small DEM raster for testing. + + Returns path to DEM, its transform, crs, and shape. + """ + demPath = tmpDir / f"DEM{suffix}" + data = np.arange(100, dtype=np.float32).reshape(10, 10) + # North-up: 10x10 grid from (0,0) to (10,10) + transform = rasterio.transform.from_bounds(0, 0, 10, 10, 10, 10) + crs = "EPSG:32633" + driver = "AAIGrid" if suffix == ".asc" else "GTiff" + with rasterio.open( + demPath, + "w", + driver=driver, + height=10, + width=10, + count=1, + dtype=data.dtype, + crs=crs, + transform=transform, + ) as dst: + dst.write(data, 1) + return demPath, transform, crs, data.shape + + +def _makeSyntheticShapefile(shpPath, fieldName, featureCoordsValues): + """Create a shapefile with a single field and polygon features. + + featureCoordsValues: list of (coords_list, field_value) tuples. + coords_list is list of (x, y) tuples forming a clockwise ring. + """ + with shapefile.Writer(shpPath, shapeType=shapefile.POLYGON) as w: + w.field(fieldName, "F", decimal=6) + for coords, value in featureCoordsValues: + w.poly([coords]) + w.record(value) + + +def test_generateMuXsiRasters_asc(): + """Test raster generation with .asc DEM and shapefiles.""" + tmpDir = pathlib.Path(tempfile.mkdtemp()) + try: + # Setup: DEM + demPath, transform, crs, demShape = _makeSyntheticDEM(tmpDir, ".asc") + inputsDir = tmpDir / "Inputs" + inputsDir.mkdir() + shutil.move(str(demPath), str(inputsDir / "DEM.asc")) + demPath = inputsDir / "DEM.asc" + + # Setup: mu shapefile with two polygons + # Polygon 1: geographic (2..5, 7..9) -> rows 1-3, cols 2-5 + # Polygon 2: geographic (6..9, 2..5) -> rows 5-8, cols 6-9 + polyDir = inputsDir / "POLYGONS" + polyDir.mkdir() + muShp = polyDir / "zones_mu.shp" + _makeSyntheticShapefile( + muShp, + "mu", + [ + ([(2, 7), (5, 7), (5, 9), (2, 9), (2, 7)], 0.300), + ([(6, 2), (9, 2), (9, 5), (6, 5), (6, 2)], 0.500), + ], + ) + + # Setup: xsi shapefile (same geometry, different values) + xsiShp = polyDir / "zones_xsi.shp" + _makeSyntheticShapefile( + xsiShp, + "xsi", + [ + ([(2, 7), (5, 7), (5, 9), (2, 9), (2, 7)], 3000.0), + ([(6, 2), (9, 2), (9, 5), (6, 5), (6, 2)], 5000.0), + ], + ) + + # Setup: config + import configparser + + cfg = configparser.ConfigParser() + cfg["DEFAULTS"] = {"default_mu": "0.155", "default_xsi": "4000."} + + # Run + spatialVoellmyInputs.generateMuXsiRasters(tmpDir, cfg) + + # Assert output files exist + rastersDir = inputsDir / "RASTERS" + muRasterPath = rastersDir / "raster_mu.asc" + xsiRasterPath = rastersDir / "raster_xi.asc" + assert muRasterPath.exists() + assert xsiRasterPath.exists() + + # Assert mu raster values + with rasterio.open(muRasterPath) as src: + muData = src.read(1) + assert src.transform == transform + assert src.shape == demShape + # Covered cells (interior of polygons) + assert muData[2, 3] == pytest.approx(0.300) # poly 1: row 2, col 3 + assert muData[6, 7] == pytest.approx(0.500) # poly 2: row 6, col 7 + # Uncovered cell should have default + assert muData[9, 0] == pytest.approx(0.155) # row 9, col 0 outside + + # Assert xsi raster values + with rasterio.open(xsiRasterPath) as src: + xsiData = src.read(1) + assert xsiData[2, 3] == pytest.approx(3000.0) + assert xsiData[6, 7] == pytest.approx(5000.0) + assert xsiData[9, 0] == pytest.approx(4000.0) + + finally: + shutil.rmtree(tmpDir) + + +def test_generateMuXsiRasters_tif(): + """Test raster generation with .tif DEM -- output should be .tif.""" + tmpDir = pathlib.Path(tempfile.mkdtemp()) + try: + demPath, transform, crs, demShape = _makeSyntheticDEM(tmpDir, ".tif") + inputsDir = tmpDir / "Inputs" + inputsDir.mkdir() + shutil.move(str(demPath), str(inputsDir / "DEM.tif")) + + polyDir = inputsDir / "POLYGONS" + polyDir.mkdir() + muShp = polyDir / "zones_mu.shp" + _makeSyntheticShapefile(muShp, "mu", [([(2, 7), (5, 7), (5, 9), (2, 9), (2, 7)], 0.300)]) + xsiShp = polyDir / "zones_xsi.shp" + _makeSyntheticShapefile(xsiShp, "xsi", [([(2, 7), (5, 7), (5, 9), (2, 9), (2, 7)], 3000.0)]) + + import configparser + + cfg = configparser.ConfigParser() + cfg["DEFAULTS"] = {"default_mu": "0.155", "default_xsi": "4000."} + + spatialVoellmyInputs.generateMuXsiRasters(tmpDir, cfg) + + rastersDir = inputsDir / "RASTERS" + muPath = rastersDir / "raster_mu.tif" + xsiPath = rastersDir / "raster_xi.tif" + assert muPath.exists() + assert xsiPath.exists() + assert muPath.suffix == ".tif" + assert xsiPath.suffix == ".tif" + + finally: + shutil.rmtree(tmpDir) + + +def test_missingMuFieldRaises(): + """Test that missing 'mu' field in shapefile raises clear error.""" + tmpDir = pathlib.Path(tempfile.mkdtemp()) + try: + demPath, _, _, _ = _makeSyntheticDEM(tmpDir, ".asc") + inputsDir = tmpDir / "Inputs" + inputsDir.mkdir() + shutil.move(str(demPath), str(inputsDir / "DEM.asc")) + polyDir = inputsDir / "POLYGONS" + polyDir.mkdir() + # Shapefile with wrong field name + _makeSyntheticShapefile( + polyDir / "zones_mu.shp", "friction_mu", [([(2, 7), (5, 7), (5, 9), (2, 9), (2, 7)], 0.3)] + ) + _makeSyntheticShapefile( + polyDir / "zones_xsi.shp", "xsi", [([(2, 7), (5, 7), (5, 9), (2, 9), (2, 7)], 3000.0)] + ) + + import configparser + + cfg = configparser.ConfigParser() + cfg["DEFAULTS"] = {"default_mu": "0.1", "default_xsi": "300."} + + with pytest.raises(KeyError, match="mu"): + spatialVoellmyInputs.generateMuXsiRasters(tmpDir, cfg) + finally: + shutil.rmtree(tmpDir) + + +def test_missingDEMRaises(): + """Test that missing DEM raises an error.""" + tmpDir = pathlib.Path(tempfile.mkdtemp()) + try: + inputsDir = tmpDir / "Inputs" + inputsDir.mkdir() + import configparser + + cfg = configparser.ConfigParser() + cfg["DEFAULTS"] = {"default_mu": "0.1", "default_xsi": "300."} + with pytest.raises(FileNotFoundError): + spatialVoellmyInputs.generateMuXsiRasters(tmpDir, cfg) + finally: + shutil.rmtree(tmpDir) From 171fc17c9511a6c5a5df34824eb04c843800b2c5 Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Tue, 16 Jun 2026 10:06:54 +0200 Subject: [PATCH 4/5] docs(spatialVoellmyInputs): add documentation for generating rasters from shapefiles --- docs/moduleCom1DFA.rst | 1 + docs/moduleIn3Utils.rst | 31 +++++++++++++++++++++++++++++++ 2 files changed, 32 insertions(+) diff --git a/docs/moduleCom1DFA.rst b/docs/moduleCom1DFA.rst index 31ad44898..b38798909 100644 --- a/docs/moduleCom1DFA.rst +++ b/docs/moduleCom1DFA.rst @@ -107,6 +107,7 @@ at least two results are generated: the *null* variant and the variant with entr - only one file per parameter allowed - if ``meshCellSize`` is different from simulation ``meshCellSize`` fields will be remeshed - only used if ``frictionModel`` is set to ``spatialVoellmy`` + - to generate these rasters from polygon shapefiles, see :ref:`moduleIn3Utils:Spatial Voellmy inputs` * **one ``_cropshape.shp`` shape file (in Inputs/POLYGONS)** diff --git a/docs/moduleIn3Utils.rst b/docs/moduleIn3Utils.rst index e16e51620..4c5d38286 100644 --- a/docs/moduleIn3Utils.rst +++ b/docs/moduleIn3Utils.rst @@ -10,6 +10,37 @@ comparison or interpolation on rasters, lines, points... Further information about the available functions can be found in :py:mod:`in3Utils.geoTrans` +Spatial Voellmy inputs +====================== + +The :py:mod:`in3Utils.spatialVoellmyInputs` module generates raster files +for the Voellmy friction parameters :math:`\mu` and :math:`\xi` from polygon +shapefiles. This is required when using the ``spatialVoellmy`` friction model +in com1DFA (see :ref:`moduleCom1DFA:Input`). + +Provide polygon shapefiles with ``mu`` and ``xsi`` attribute fields in +``Inputs/POLYGONS/``, with file names ending in ``_mu.shp`` and ``_xsi.shp``. +The DEM must be placed in ``Inputs/``. The generated rasters are written to +``Inputs/RASTERS/`` with the same file format as the DEM. + +Default values for areas not covered by the polygon shapefiles are set in +``avaframe/in3Utils/spatialVoellmyInputsCfg.ini``. + + +To run +------ + +* first go to ``AvaFrame/avaframe`` +* copy ``in3Utils/spatialVoellmyInputsCfg.ini`` to + ``in3Utils/local_spatialVoellmyInputsCfg.ini`` and set desired + ``default_mu`` and ``default_xsi`` values (if not, the default values + are used) +* ensure the DEM and shapefiles are in the avalanche directory as + described above +* run:: + + python3 runScripts/runSpatialVoellmyInputs.py + Generate Topography =================== From e55a2b0888aa94c536dae559c845fbf0f429344f Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Tue, 16 Jun 2026 13:16:05 +0200 Subject: [PATCH 5/5] refactor(spatialVoellmyInputs): use `getfloat` for configuration parsing and improve docs - Replaced `float` casting with `getfloat` for cleaner configuration parsing. - Updated `_rasterizeShapefile` docstring to clarify default cell value behavior. --- avaframe/in3Utils/spatialVoellmyInputs.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/avaframe/in3Utils/spatialVoellmyInputs.py b/avaframe/in3Utils/spatialVoellmyInputs.py index 1b1d4b83c..4ffad9347 100644 --- a/avaframe/in3Utils/spatialVoellmyInputs.py +++ b/avaframe/in3Utils/spatialVoellmyInputs.py @@ -58,8 +58,8 @@ def generateMuXsiRasters(avaDir, cfg): demCrs = demHeader["crs"] demShape = (demHeader["nrows"], demHeader["ncols"]) - defaultMu = float(cfg["DEFAULTS"]["default_mu"]) - defaultXsi = float(cfg["DEFAULTS"]["default_xsi"]) + defaultMu = cfg["DEFAULTS"].getfloat("default_mu") + defaultXsi = cfg["DEFAULTS"].getfloat("default_xsi") # Rasterize mu log.info("Rasterizing mu shapefile: %s", muShpPath) @@ -92,6 +92,8 @@ def generateMuXsiRasters(avaDir, cfg): def _rasterizeShapefile(shpPath, defaultValue, fieldName, demShape, demTransform): """Rasterize a polygon shapefile attribute field onto a DEM-matching grid. + All cells not covered by any polygon are filled with defaultValue. + Parameters ---------- shpPath : pathlib.Path