[com6]: add support for spatialVoellmy friction model

avaframe/in3Utils/spatialVoellmyInputs.py

@@ -15,21 +15,21 @@
 log = logging.getLogger(__name__)
 
 
-def generateMuXsiRasters(avaDir, cfg):
-    """Generate mu and xi raster files from polygon shapefiles.
+def generateMuXiRasters(avaDir, cfg):
+    """Generate mu and xi raster files from a polygon shapefile.
 
-    Reads polygon shapefiles with "mu" and "xsi" attribute fields,
+    Reads a polygon shapefile with "mu" and "xi" 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.
+        Inputs/POLYGONS/ with *_spatialVoellmy.shp shapefile.
     cfg : configparser.ConfigParser
         Configuration with [DEFAULTS] section containing
-        default_mu and default_xsi values for uncovered areas.
+        default_mu and default_xi values for uncovered areas.
     """
     avaDir = pathlib.Path(avaDir)
     inputDir = avaDir / "Inputs"
@@ -40,17 +40,15 @@ def generateMuXsiRasters(avaDir, cfg):
     demPath = getDEMPath(avaDir)
     demSuffix = demPath.suffix
 
-    # Find shapefiles
-    muShpPath, muAvailable, _ = getAndCheckInputFiles(
-        inputDir, "POLYGONS", "mu shapefile", fileExt="shp", fileSuffix="_mu"
+    # Find shapefile
+    shpPath, shpAvailable, _ = getAndCheckInputFiles(
+        inputDir, "POLYGONS", "spatialVoellmy shapefile", fileExt="shp",
+        fileSuffix="_spatialVoellmy"
     )
-    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)
+    if shpAvailable == "No":
+        raise FileNotFoundError(
+            "No *_spatialVoellmy.shp found in %s/POLYGONS/" % inputDir
+        )
 
     # Read DEM header
     demHeader = readRasterHeader(demPath)
@@ -59,22 +57,38 @@ def generateMuXsiRasters(avaDir, cfg):
     demShape = (demHeader["nrows"], demHeader["ncols"])
 
     defaultMu = cfg["DEFAULTS"].getfloat("default_mu")
-    defaultXsi = cfg["DEFAULTS"].getfloat("default_xsi")
+    defaultXi = cfg["DEFAULTS"].getfloat("default_xi")
 
-    # Rasterize mu
-    log.info("Rasterizing mu shapefile: %s", muShpPath)
-    muRaster = _rasterizeShapefile(muShpPath, defaultMu, "mu", demShape, demTransform)
+    # Validate required fields
+    with shapefile.Reader(str(shpPath)) as sf:
+        fieldNames = [f[0].lower() for f in sf.fields[1:]]
+    for field in ["mu", "xi"]:
+        if field not in fieldNames:
+            raise KeyError(
+                "Field '%s' not found in %s. Available fields: %s"
+                % (field, shpPath.name, fieldNames)
+            )
+
+    # Rasterize mu and xi from the same shapefile
+    log.info("Rasterizing mu from: %s", shpPath)
+    muRaster = _rasterizeShapefile(shpPath, defaultMu, "mu", demShape, demTransform)
 
-    # Rasterize xsi
-    log.info("Rasterizing xsi shapefile: %s", xsiShpPath)
-    xsiRaster = _rasterizeShapefile(xsiShpPath, defaultXsi, "xsi", demShape, demTransform)
+    log.info("Rasterizing xi from: %s", shpPath)
+    xiRaster = _rasterizeShapefile(shpPath, defaultXi, "xi", demShape, demTransform)
 
     # Determine output driver
     if demSuffix == ".asc":
         driver = "AAIGrid"
     else:
         driver = "GTiff"
 
+    # Check if any mu or xi raster files already exist
+    existing = sorted(p.name for p in outDir.glob("*_mu.*")) + sorted(p.name for p in outDir.glob("*_xi.*"))
+    if existing:
+        raise FileExistsError(
+            "Output file(s) already exist in %s: %s" % (outDir, ", ".join(existing))
+        )
+
     # Write output
     outHeader = {
         "driver": driver,
@@ -84,8 +98,8 @@ def generateMuXsiRasters(avaDir, cfg):
     }
     log.info("Writing mu raster")
     writeResultToRaster(outHeader, muRaster, outDir / "raster_mu")
-    log.info("Writing xsi raster")
-    writeResultToRaster(outHeader, xsiRaster, outDir / "raster_xi")
+    log.info("Writing xi raster")
+    writeResultToRaster(outHeader, xiRaster, outDir / "raster_xi")
     log.info("Raster generation completed.")
 
 

avaframe/in3Utils/spatialVoellmyInputsCfg.ini

@@ -2,5 +2,5 @@
 # Default mu value for areas not covered by shapefiles
 default_mu = 0.155
 
-# Default xsi value for areas not covered by shapefiles
-default_xsi = 4000.
+# Default xi value for areas not covered by shapefiles
+default_xi = 4000.

avaframe/runCom6RockAvalanche.py

@@ -15,15 +15,18 @@
 
 # import computation modules
 from avaframe.com6RockAvalanche import com6RockAvalanche
+from avaframe.in3Utils import spatialVoellmyInputs
 
 
-def runCom6RockAvalanche(avalancheDir=""):
-    """Run com1DFA with rock avalanche parameters with only an avalanche directory as input
+def runCom6RockAvalanche(avalancheDir="", calibration="voellmy"):
+    """Run com1DFA with rock avalanche parameters
 
     Parameters
     ----------
     avalancheDir: str
         path to avalanche directory (setup e.g. with init scripts)
+    calibration: str
+        friction model: voellmy (default) or spatialVoellmy
 
     Returns
     -------
@@ -54,14 +57,51 @@ def runCom6RockAvalanche(avalancheDir=""):
     # Clean input directory(ies) of old work files
     initProj.cleanSingleAvaDir(avalancheDir, deleteOutput=False)
 
+    # pathlib version of avalanche dir
+    avaDir = pathlib.Path(avalancheDir)
+
     # load rock avalanche config
     rockAvalancheCfg = cfgUtils.getModuleConfig(com6RockAvalanche, avalancheDir)
 
+    # override friction model if spatialVoellmy calibration is requested
+    if calibration == "spatialVoellmy":
+        rockAvalancheCfg["com1DFA_com1DFA_override"]["frictModel"] = "spatialVoellmy"
+
+        muRasters = list((avaDir / "Inputs" / "RASTERS").glob("*_mu.*"))
+        xiRasters = list((avaDir / "Inputs" / "RASTERS").glob("*_xi.*"))
+        spatialShps = list(
+            (avaDir / "Inputs" / "POLYGONS").glob("*_spatialVoellmy.shp")
+        )
+
+        rastersExist = bool(muRasters and xiRasters)
+        shpExists = bool(spatialShps)
+
+        if rastersExist and shpExists:
+            raise RuntimeError(
+                "spatialVoellmy friction model: both rasters in Inputs/RASTERS/"
+                " and *_spatialVoellmy.shp in Inputs/POLYGONS/ found"
+                " - ambiguous input"
+            )
+        elif shpExists:
+            spatialVoellmyCfg = cfgUtils.getModuleConfig(spatialVoellmyInputs)
+            # set default fill values from rock avalanche Voellmy defaults
+            overrideParams = rockAvalancheCfg["com1DFA_com1DFA_override"]
+            spatialVoellmyCfg["DEFAULTS"]["default_mu"] = overrideParams["muvoellmy"]
+            spatialVoellmyCfg["DEFAULTS"]["default_xi"] = overrideParams["xsivoellmy"]
+            spatialVoellmyInputs.generateMuXiRasters(avaDir, spatialVoellmyCfg)
+        elif rastersExist:
+            log.info("spatialVoellmy: using existing mu/xi rasters from Inputs/RASTERS/")
+        else:
+            raise FileNotFoundError(
+                "spatialVoellmy friction model: no *_mu and *_xi rasters found in"
+                " Inputs/RASTERS/ and no *_spatialVoellmy.shp found in"
+                " Inputs/POLYGONS/"
+            )
+
     # perform com1DFA simulation with rock avalanche settings
     _, plotDict, reportDictList, _ = com6RockAvalanche.com6RockAvalancheMain(cfgMain, rockAvalancheCfg)
 
     # Get peakfiles to return to QGIS
-    avaDir = pathlib.Path(avalancheDir)
     inputDir = avaDir / "Outputs" / "com1DFA" / "peakFiles"
     peakFilesDF = fU.makeSimDF(inputDir, avaDir=avaDir)
 
@@ -77,6 +117,14 @@ def runCom6RockAvalanche(avalancheDir=""):
     parser.add_argument(
         "avadir", metavar="avadir", type=str, nargs="?", default="", help="the avalanche directory"
     )
+    parser.add_argument(
+        "-fc",
+        "--friction_calibration",
+        choices=["voellmy", "spatialVoellmy"],
+        type=str,
+        default="voellmy",
+        help="friction model: voellmy (default) or spatialVoellmy",
+    )
 
     args = parser.parse_args()
-    runCom6RockAvalanche(str(args.avadir))
+    runCom6RockAvalanche(str(args.avadir), str(args.friction_calibration))

avaframe/runScripts/runSpatialVoellmyInputs.py

@@ -13,7 +13,7 @@
 
 
 def runSpatialVoellmyInputs(avaDir=""):
-    """Run generation of mu and xi rasters from shapefiles.
+    """Run generation of mu and xi rasters from shapefile.
 
     Parameters
     ----------
@@ -36,7 +36,7 @@ def runSpatialVoellmyInputs(avaDir=""):
     initProj.cleanSingleAvaDir(avaDir, deleteOutput=False)
 
     cfg = cfgUtils.getModuleConfig(spatialVoellmyInputs)
-    spatialVoellmyInputs.generateMuXsiRasters(avaDir, cfg)
+    spatialVoellmyInputs.generateMuXiRasters(avaDir, cfg)
 
     endTime = time.time()
     log.info("Took %6.1f seconds to calculate.", endTime - startTime)