diff --git a/avaframe/in3Utils/spatialVoellmyInputs.py b/avaframe/in3Utils/spatialVoellmyInputs.py index 4ffad9347..491d0224b 100644 --- a/avaframe/in3Utils/spatialVoellmyInputs.py +++ b/avaframe/in3Utils/spatialVoellmyInputs.py @@ -15,10 +15,10 @@ 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/. @@ -26,10 +26,10 @@ def generateMuXsiRasters(avaDir, cfg): ---------- 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,15 +57,24 @@ 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": @@ -75,6 +82,13 @@ def generateMuXsiRasters(avaDir, cfg): 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.") diff --git a/avaframe/in3Utils/spatialVoellmyInputsCfg.ini b/avaframe/in3Utils/spatialVoellmyInputsCfg.ini index 1fb927e33..5873282dc 100644 --- a/avaframe/in3Utils/spatialVoellmyInputsCfg.ini +++ b/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. diff --git a/avaframe/runCom6RockAvalanche.py b/avaframe/runCom6RockAvalanche.py index 0d670beea..34f02e243 100644 --- a/avaframe/runCom6RockAvalanche.py +++ b/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)) diff --git a/avaframe/runScripts/runSpatialVoellmyInputs.py b/avaframe/runScripts/runSpatialVoellmyInputs.py index f02c89ad8..55fef7e9a 100644 --- a/avaframe/runScripts/runSpatialVoellmyInputs.py +++ b/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) diff --git a/avaframe/tests/test_spatialVoellmyInputs.py b/avaframe/tests/test_spatialVoellmyInputs.py index 6b6f0b9ae..c700ac08f 100644 --- a/avaframe/tests/test_spatialVoellmyInputs.py +++ b/avaframe/tests/test_spatialVoellmyInputs.py @@ -49,8 +49,22 @@ def _makeSyntheticShapefile(shpPath, fieldName, featureCoordsValues): w.record(value) -def test_generateMuXsiRasters_asc(): - """Test raster generation with .asc DEM and shapefiles.""" +def _makeSpatialVoellmyShapefile(shpPath, features): + """Create a *_spatialVoellmy.shp with 'mu' and 'xi' fields. + + features: list of (coords_list, mu_value, xi_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("mu", "F", decimal=6) + w.field("xi", "F", decimal=6) + for coords, muVal, xiVal in features: + w.poly([coords]) + w.record(muVal, xiVal) + + +def test_generateMuXiRasters_asc(): + """Test raster generation with .asc DEM and shapefile.""" tmpDir = pathlib.Path(tempfile.mkdtemp()) try: # Setup: DEM @@ -60,29 +74,17 @@ def test_generateMuXsiRasters_asc(): shutil.move(str(demPath), str(inputsDir / "DEM.asc")) demPath = inputsDir / "DEM.asc" - # Setup: mu shapefile with two polygons + # Setup: spatialVoellmy shapefile with two polygons, each with mu and xi # 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", + shpPath = polyDir / "zones_spatialVoellmy.shp" + _makeSpatialVoellmyShapefile( + shpPath, [ - ([(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), + ([(2, 7), (5, 7), (5, 9), (2, 9), (2, 7)], 0.300, 3000.0), + ([(6, 2), (9, 2), (9, 5), (6, 5), (6, 2)], 0.500, 5000.0), ], ) @@ -90,17 +92,17 @@ def test_generateMuXsiRasters_asc(): import configparser cfg = configparser.ConfigParser() - cfg["DEFAULTS"] = {"default_mu": "0.155", "default_xsi": "4000."} + cfg["DEFAULTS"] = {"default_mu": "0.155", "default_xi": "4000."} # Run - spatialVoellmyInputs.generateMuXsiRasters(tmpDir, cfg) + spatialVoellmyInputs.generateMuXiRasters(tmpDir, cfg) # Assert output files exist rastersDir = inputsDir / "RASTERS" muRasterPath = rastersDir / "raster_mu.asc" - xsiRasterPath = rastersDir / "raster_xi.asc" + xiRasterPath = rastersDir / "raster_xi.asc" assert muRasterPath.exists() - assert xsiRasterPath.exists() + assert xiRasterPath.exists() # Assert mu raster values with rasterio.open(muRasterPath) as src: @@ -113,18 +115,18 @@ def test_generateMuXsiRasters_asc(): # 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) + # Assert xi raster values + with rasterio.open(xiRasterPath) as src: + xiData = src.read(1) + assert xiData[2, 3] == pytest.approx(3000.0) + assert xiData[6, 7] == pytest.approx(5000.0) + assert xiData[9, 0] == pytest.approx(4000.0) finally: shutil.rmtree(tmpDir) -def test_generateMuXsiRasters_tif(): +def test_generateMuXiRasters_tif(): """Test raster generation with .tif DEM -- output should be .tif.""" tmpDir = pathlib.Path(tempfile.mkdtemp()) try: @@ -135,25 +137,26 @@ def test_generateMuXsiRasters_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)]) + shpPath = polyDir / "zones_spatialVoellmy.shp" + _makeSpatialVoellmyShapefile( + shpPath, + [([(2, 7), (5, 7), (5, 9), (2, 9), (2, 7)], 0.300, 3000.0)], + ) import configparser cfg = configparser.ConfigParser() - cfg["DEFAULTS"] = {"default_mu": "0.155", "default_xsi": "4000."} + cfg["DEFAULTS"] = {"default_mu": "0.155", "default_xi": "4000."} - spatialVoellmyInputs.generateMuXsiRasters(tmpDir, cfg) + spatialVoellmyInputs.generateMuXiRasters(tmpDir, cfg) rastersDir = inputsDir / "RASTERS" muPath = rastersDir / "raster_mu.tif" - xsiPath = rastersDir / "raster_xi.tif" + xiPath = rastersDir / "raster_xi.tif" assert muPath.exists() - assert xsiPath.exists() + assert xiPath.exists() assert muPath.suffix == ".tif" - assert xsiPath.suffix == ".tif" + assert xiPath.suffix == ".tif" finally: shutil.rmtree(tmpDir) @@ -169,21 +172,21 @@ def test_missingMuFieldRaises(): 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)] - ) + # Shapefile with wrong field name instead of 'mu' + shpPath = polyDir / "zones_spatialVoellmy.shp" + with shapefile.Writer(shpPath, shapeType=shapefile.POLYGON) as w: + w.field("friction_mu", "F", decimal=6) + w.field("xi", "F", decimal=6) + w.poly([[(2, 7), (5, 7), (5, 9), (2, 9), (2, 7)]]) + w.record(0.3, 3000.0) import configparser cfg = configparser.ConfigParser() - cfg["DEFAULTS"] = {"default_mu": "0.1", "default_xsi": "300."} + cfg["DEFAULTS"] = {"default_mu": "0.1", "default_xi": "300."} with pytest.raises(KeyError, match="mu"): - spatialVoellmyInputs.generateMuXsiRasters(tmpDir, cfg) + spatialVoellmyInputs.generateMuXiRasters(tmpDir, cfg) finally: shutil.rmtree(tmpDir) @@ -197,8 +200,8 @@ def test_missingDEMRaises(): import configparser cfg = configparser.ConfigParser() - cfg["DEFAULTS"] = {"default_mu": "0.1", "default_xsi": "300."} + cfg["DEFAULTS"] = {"default_mu": "0.1", "default_xi": "300."} with pytest.raises(FileNotFoundError): - spatialVoellmyInputs.generateMuXsiRasters(tmpDir, cfg) + spatialVoellmyInputs.generateMuXiRasters(tmpDir, cfg) finally: shutil.rmtree(tmpDir) diff --git a/docs/moduleCom6RockAvalanche.rst b/docs/moduleCom6RockAvalanche.rst index 421314d4f..6f763e23b 100644 --- a/docs/moduleCom6RockAvalanche.rst +++ b/docs/moduleCom6RockAvalanche.rst @@ -23,6 +23,12 @@ To run * copy ``avaframeCfg.ini`` to ``local_avaframeCfg.ini`` and set your desired avalanche directory name * create an avalanche directory with required input files - for this task you can use :ref:`moduleIn3Utils:Initialize Project` * copy ``com6RockAvalanche/com6RockAvalancheCfg.ini`` to ``com6RockAvalanche/local_com6RockAvalancheCfg.ini`` and if desired change configuration settings +* optionally, the ``spatialVoellmy`` friction model can be selected with + ``--friction_calibration spatialVoellmy``. When using the QGis Connector, + a shapefile with ``mu`` and ``xi`` attributes + can be provided and the required rasters will + be generated automatically. See + :ref:`moduleIn3Utils:Spatial Voellmy inputs` for details. * if you are on a develop installation, make sure you have an updated compilation, see :ref:`complexUsage:Update AvaFrame` * run: :: diff --git a/docs/moduleIn3Utils.rst b/docs/moduleIn3Utils.rst index 4c5d38286..4915cff76 100644 --- a/docs/moduleIn3Utils.rst +++ b/docs/moduleIn3Utils.rst @@ -14,16 +14,17 @@ 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`). +for the Voellmy friction parameters :math:`\mu` and :math:`\xi` from a +polygon shapefile. 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. +Provide a polygon shapefile with ``mu`` and ``xi`` attribute fields in +``Inputs/POLYGONS/``, with a file name ending in ``_spatialVoellmy.shp``. +Each polygon feature must include both a ``mu`` and an ``xi`` attribute +value. 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 +Default values for areas not covered by the polygons are set in ``avaframe/in3Utils/spatialVoellmyInputsCfg.ini``. @@ -33,9 +34,9 @@ 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 + ``default_mu`` and ``default_xi`` values (if not, the default values are used) -* ensure the DEM and shapefiles are in the avalanche directory as +* ensure the DEM and shapefile are in the avalanche directory as described above * run::