[com1]: add saving option for resistance fields

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,38 +1325,40 @@ 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,
                 )
                 log.info(
                     "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,
         )
         log.info(
             "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
 
 

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)

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

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,
     }