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.

Author: dwolfsch

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:

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

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.")

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

avaframe/in2Trans/shpConversion.py

@@ -94,6 +94,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
@@ -113,9 +116,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"):
@@ -132,13 +140,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
@@ -173,6 +183,7 @@ def SHP2Array(infile, defname=None):
     SHPdata["layerName"] = layerNameList
     SHPdata["nParts"] = nParts
     SHPdata["nFeatures"] = len(Start)
+    SHPdata["attributes"] = attributes  # Add attributes to SHPdata
 
     sf.close()
 

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))