From 14fbf31aa54cbd9628974e8b8b7080092302b874 Mon Sep 17 00:00:00 2001 From: Anna Wirbel Date: Mon, 8 Jun 2026 11:20:03 +0200 Subject: [PATCH] add saving option for resistance fields add simname to name move outside if, add log.info change name from orig to track adjust tests test(com1DFA): add resistance initialization tests with raster export - Added comprehensive tests for initializing resistance with raster inputs, including raster export validation. - Verified creation of output raster files with appropriate naming conventions. --- avaframe/com1DFA/com1DFA.py | 61 ++++++++++++--- avaframe/com1DFA/com1DFATools.py | 11 +-- avaframe/tests/test_com1DFA.py | 115 ++++++++++++++++++++++++++++ avaframe/tests/test_com1DFATools.py | 8 +- 4 files changed, 174 insertions(+), 21 deletions(-) diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 50d34963f..b9cb10aa5 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -1225,19 +1225,19 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): inputSimLines["releaseLine"]["header"] = dem["originalHeader"] # export release area raster to file if cfg["EXPORTS"].getboolean("exportRasters"): - outDir = pathlib.Path(cfgGen["avalancheDir"], "Outputs", "internalRasters") - fU.makeADir(outDir) + outDirRasters = outDir / "internalRasters" + fU.makeADir(outDirRasters) useCompression = cfg["EXPORTS"].getboolean("useCompression") IOf.writeResultToRaster( dem["originalHeader"], relRaster, - (outDir / "releaseRaster"), + (outDirRasters / ("releaseRaster_%s" % logName)), flip=True, useCompression=useCompression, ) log.info( "Release area raster derived from %s saved to %s" - % (releaseLine["initializedFrom"], str(outDir / "releaseRaster")) + % (releaseLine["initializedFrom"], str(outDirRasters / ("releaseRaster_%s" % logName))) ) particles = initializeParticles( cfgGen, @@ -1325,12 +1325,13 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): ) # export secondary release raster used for computations (after cutting potential overlap with release) if cfg["EXPORTS"].getboolean("exportRasters"): - outDir = pathlib.Path(cfg["GENERAL"]["avalancheDir"], "Outputs", "internalRasters") + outDirRasters = outDir / "internalRasters" + fU.makeADir(outDirRasters) useCompression = cfg["EXPORTS"].getboolean("useCompression") IOf.writeResultToRaster( dem["originalHeader"], secRelRaster, - (outDir / ("secondaryReleaseRaster_%d" % secIndex)), + (outDirRasters / ("secondaryReleaseRaster_%d_%s" % (secIndex, logName))), flip=True, useCompression=useCompression, ) @@ -1338,17 +1339,18 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): "SecondaryRelease area raster derived from %s saved to %s" % ( inputSimLines["entResInfo"]["secondaryRelThFileType"], - str(outDir / ("secondaryReleaseRaster_%d" % secIndex)), + str(outDirRasters / ("secondaryReleaseRaster_%d_%s" % (secIndex, logName))), ) ) # export entrainment raster used for computations (after cutting potential overlap with release or secondary release) if cfg["EXPORTS"].getboolean("exportRasters"): - outDir = pathlib.Path(cfg["GENERAL"]["avalancheDir"], "Outputs", "internalRasters") + outDirRasters = outDir / "internalRasters" + fU.makeADir(outDirRasters) useCompression = cfg["EXPORTS"].getboolean("useCompression") IOf.writeResultToRaster( dem["originalHeader"], entrMassRaster / cfg["GENERAL"].getfloat("rhoEnt"), - (outDir / "entrainmentRaster"), + (outDirRasters / ("entrainmentRaster_%s" % logName)), flip=True, useCompression=useCompression, ) @@ -1356,7 +1358,7 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): "Entrainment area raster derived from %s saved to %s" % ( inputSimLines["entResInfo"]["entThFileType"], - str(outDir / "entrainmentRaster"), + str(outDirRasters / ("entrainmentRaster_%s" % logName)), ) ) @@ -1377,8 +1379,27 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): ) fields["cResRaster"] = cResRaster fields["detRaster"] = detRaster - fields["cResRasterOrig"] = cResRaster - fields["detRasterOrig"] = detRaster + fields["cResRasterTrack"] = cResRaster + fields["detRasterTrack"] = detRaster + # export resistance raster used for computations + if (cfgGen["simTypeActual"] in ["entres", "res"]) and cfg["EXPORTS"].getboolean("exportRasters"): + outDirRasters = outDir / "internalRasters" + fU.makeADir(outDirRasters) + useCompression = cfg["EXPORTS"].getboolean("useCompression") + IOf.writeResultToRaster( + dem["originalHeader"], + fields["cResRaster"], + (outDirRasters / ("resistanceRaster_%s" % logName)), + flip=True, + useCompression=useCompression, + ) + log.info( + "Resistance area raster derived from %s saved to %s" + % ( + inputSimLines["resLine"]["fileName"], + str(outDirRasters / ("resistanceRaster_%s" % logName)), + ) + ) for fric in ["mu", "xi"]: if (inputSimLines[fric + "File"] == None) or ( @@ -2448,6 +2469,22 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si % particles["nExitedParticles"] ) + # save final cRes raster + if (cfgGen["simTypeActual"] in ["entres", "res"]) and cfg["EXPORTS"].getboolean("exportRasters"): + outDirRes = outDir / "internalRasters" + fU.makeADir(outDirRes) + IOf.writeResultToRaster( + dem["originalHeader"], + fields["cResRasterTrack"], + (outDirRes / ("cResRaster_Final_%s" % cuSimName)), + flip=True, + useCompression=cfg["EXPORTS"].getboolean("useCompression"), + ) + log.info( + "Resistance area raster (final state) saved to %s" + % (str(outDirRes / ("cResFinal_%s" % cuSimName))) + ) + return Tsave, infoDict, contourDictXY diff --git a/avaframe/com1DFA/com1DFATools.py b/avaframe/com1DFA/com1DFATools.py index 5abd0d857..9c7f73f8c 100644 --- a/avaframe/com1DFA/com1DFATools.py +++ b/avaframe/com1DFA/com1DFATools.py @@ -9,6 +9,7 @@ import math import pathlib import numpy as np +import copy as cp from deepdiff import DeepDiff @@ -366,7 +367,7 @@ def updateResCoeffFields(fields, cfg): Parameters ------------ fields: dict - dictionary with cResRasterOrig, cResRaster and detRasterOrig, detRaster fields + dictionary with cResRasterTrack, cResRaster and detRasterTrack, detRaster fields cfg: configparser object configuration of com1DFA, thresholds @@ -377,8 +378,8 @@ def updateResCoeffFields(fields, cfg): """ # fetch cRes and detK raster and thresholds for FV and FT - cResOrig = fields["cResRasterOrig"].copy() - detOrig = fields["detRasterOrig"].copy() + cResOrig = cp.deepcopy(fields["cResRasterTrack"]) + detOrig = cp.deepcopy(fields["detRasterTrack"]) vMin = cfg.getfloat("forestVMin") thMin = cfg.getfloat("forestThMin") vMax = cfg.getfloat("forestVMax") @@ -400,8 +401,8 @@ def updateResCoeffFields(fields, cfg): if lTh > 0: cResOrig = np.where(((fields["FV"] > vMax) | (fields["FT"] > thMax)), 0, cResOrig) detOrig = np.where(((fields["FV"] > vMax) | (fields["FT"] > thMax)), 0, detOrig) - fields["cResRasterOrig"] = cResOrig - fields["detRasterOrig"] = detOrig + fields["cResRasterTrack"] = cResOrig + fields["detRasterTrack"] = detOrig log.debug( "Resistance area removed %d cells because FV or FT exceeded %.2f ms-1, %.2f m" % (lTh, vMax, thMax) diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py index 70592941f..17f1cd8cd 100644 --- a/avaframe/tests/test_com1DFA.py +++ b/avaframe/tests/test_com1DFA.py @@ -2533,6 +2533,7 @@ def test_initializeSimulation(tmp_path): demHeader["nodata_value"] = -9999 demHeader["nrows"] = 12 demHeader["ncols"] = 12 + demHeader["driver"] = "AAIGrid" demData = np.ones((12, 12)) demOri = {"header": demHeader, "rasterData": demData} @@ -2780,6 +2781,120 @@ def test_initializeSimulation(tmp_path): assert np.isin(np.round(fields4["pfv"]), [0.0, 10.0]).all() assert np.all(particles4["ux"] == 0.0) + # test resistance initialization with rasters export + cfg = configparser.ConfigParser() + cfg["REPORT"] = {} + cfg["GENERAL"] = { + "methodMeshNormal": "1", + "thresholdPointInPoly": "0.001", + "useRelThFromIni": "False", + "resType": "ppr|pft|pfv", + "relTh": "1.0", + "useEntThFromIni": "False", + "meshCellSizeThreshold": "0.0001", + "meshCellSize": "1.", + "simTypeActual": "entres", + "rhoEnt": "100.", + "entTh": "0.3", + "rho": "200.", + "gravAcc": "9.81", + "massPerParticleDeterminationMethod": "MPPDH", + "interpOption": "2", + "sphKernelRadius": "1", + "deltaTh": "0.25", + "seed": "12345", + "initPartDistType": "uniform", + "thresholdPointInPoly": "0.001", + "avalancheDir": "data/avaTest", + "cRes": "0.003", + "initialiseParticlesFromFile": "False", + "entTempRef": "-10.", + "cpIce": "2050.", + "TIni": "-10.", + "ResistanceModel": "default", + "cResH": "0.003", + "detK": "0.05", + "detrainment": "False", + "restitutionCoefficient": 1, + "nIterDam": 1, + } + cfg["EXPORTS"] = {"exportRasters": "True"} + + # create dem with full header fields needed for raster export + demHeaderRes = demHeader.copy() + demHeaderRes["transform"] = transformFromASCHeader(demHeaderRes) + demHeaderRes["crs"] = rasterio.crs.CRS.from_epsg(31287) + demOriRes = {"header": demHeaderRes, "rasterData": demData} + + resRasterData = np.zeros((12, 12)) + resRasterData[5:8, 2:4] = 1 + + # write a real raster file so plotReleaseScenarioView can read it + resRasterPath = tmp_path / "resRasterTest" + IOf.writeResultToRaster(demHeaderRes, resRasterData, resRasterPath, flip=False) + + resLine = { + "fileName": str(resRasterPath) + ".asc", + "Name": ["testRes"], + "initializedFrom": "raster", + "rasterData": resRasterData, + } + + releaseLine = { + "x": np.asarray([6.9, 8.5, 8.5, 6.9, 6.9]), + "y": np.asarray([7.9, 7.9, 9.5, 9.5, 7.9]), + "Start": np.asarray([0]), + "Length": np.asarray([5]), + "Name": [""], + "thickness": [1.0], + "thicknessSource": ["ini File"], + "type": "release", + "file": relFileTest, + "initializedFrom": "shapefile", + } + + inputSimLines = { + "releaseLine": releaseLine, + "entResInfo": {"flagSecondaryRelease": "No", "entThFileType": "shp file"}, + "entLine": { + "fileName": (avaDir / "ENT" / "entAlr.shp"), + "Name": ["testEnt"], + "Start": np.asarray([0.0]), + "thickness": [0.3, 0.3], + "thicknessSource": ["shp file", "shp file"], + "Length": np.asarray([5]), + "x": np.asarray([4, 5.0, 5.0, 4.0, 4.0]), + "type": "entrainment", + "y": np.asarray([4.0, 4.0, 5.0, 5.0, 4.0]), + "initializedFrom": "shapefile", + }, + "secondaryReleaseLine": None, + "resLine": resLine, + "relThFile": "", + "entThFile": "", + "relThField": "", + "damLine": None, + "muFile": None, + "xiFile": None, + } + logName = "simLog" + + particles5, fields5, dem5, reportAreaInfo5 = com1DFA.initializeSimulation( + cfg, outDir, demOriRes, inputSimLines, logName + ) + + assert reportAreaInfo5["resistance"] == "Yes" + assert "cResRasterTrack" in fields5 + assert "detRasterTrack" in fields5 + assert np.array_equal(fields5["cResRaster"], fields5["cResRasterTrack"]) + assert np.array_equal(fields5["detRaster"], fields5["detRasterTrack"]) + + # check that raster files are created with logName in filename + rastersDir = outDir / "internalRasters" + assert (rastersDir / ("releaseRaster_%s.asc" % logName)).is_file() + assert (rastersDir / ("resistanceRaster_%s.asc" % logName)).is_file() + assert (rastersDir / ("entrainmentRaster_%s.asc" % logName)).is_file() + def test_runCom1DFA(tmp_path, caplog): """Check that runCom1DFA produces the good outputs""" diff --git a/avaframe/tests/test_com1DFATools.py b/avaframe/tests/test_com1DFATools.py index 98cc5074e..13c5d1b70 100644 --- a/avaframe/tests/test_com1DFATools.py +++ b/avaframe/tests/test_com1DFATools.py @@ -99,8 +99,8 @@ def test_updateResCoeffFields(tmp_path): FT = np.zeros((4, 10)) fields = { - "cResRasterOrig": cResRasterOrig, - "detRasterOrig": detRasterOrig, + "cResRasterTrack": cResRasterOrig, + "detRasterTrack": detRasterOrig, "FV": FV, "FT": FT, } @@ -129,8 +129,8 @@ def test_updateResCoeffFields(tmp_path): FV[1, 3] = 41 FT[1, 1] = 11 fields = { - "cResRasterOrig": cResRasterOrig, - "detRasterOrig": detRasterOrig, + "cResRasterTrack": cResRasterOrig, + "detRasterTrack": detRasterOrig, "FV": FV, "FT": FT, }