[com1;in2] Add tiff compression; adjust particles export logic

avaframe/com1DFA/com1DFA.py

@@ -1184,7 +1184,14 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     if cfg["EXPORTS"].getboolean("exportRasters"):
         outDir = pathlib.Path(cfgGen["avalancheDir"], "Outputs", "internalRasters")
         fU.makeADir(outDir)
-        IOf.writeResultToRaster(dem["originalHeader"], relRaster, (outDir / "releaseRaster"), flip=True)
+        useCompression = cfg["EXPORTS"].getboolean("useCompression")
+        IOf.writeResultToRaster(
+            dem["originalHeader"],
+            relRaster,
+            (outDir / "releaseRaster"),
+            flip=True,
+            useCompression=useCompression,
+        )
         log.info(
             "Release area raster derived from %s saved to %s"
             % (releaseLine["initializedFrom"], str(outDir / "releaseRaster"))
@@ -1268,11 +1275,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")
+                useCompression = cfg["EXPORTS"].getboolean("useCompression")
                 IOf.writeResultToRaster(
                     dem["originalHeader"],
                     secRelRaster,
                     (outDir / ("secondaryReleaseRaster_%d" % secIndex)),
                     flip=True,
+                    useCompression=useCompression,
                 )
                 log.info(
                     "SecondaryRelease area raster derived from %s saved to %s"
@@ -1284,11 +1293,13 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, 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")
+        useCompression = cfg["EXPORTS"].getboolean("useCompression")
         IOf.writeResultToRaster(
             dem["originalHeader"],
             entrMassRaster / cfg["GENERAL"].getfloat("rhoEnt"),
             (outDir / "entrainmentRaster"),
             flip=True,
+            useCompression=useCompression,
         )
         log.info(
             "Entrainment area raster derived from %s saved to %s"
@@ -2017,8 +2028,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     # make sure to save all desiered resuts for first and last time step for
     # the report
     resTypesReport = fU.splitIniValueToArraySteps(cfg["REPORT"]["plotFields"])
-    # always add particles to first and last time step
-    resTypesLast = list(set(resTypes + resTypesReport + ["particles"]))
+    resTypesLast = list(set(resTypes + resTypesReport))
     # derive friction type
     # turn friction model into integer
     frictModelsList = [
@@ -2059,7 +2069,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     # setup a result fields info data frame to save max values of fields and avalanche front
     resultsDF = setupresultsDF(resTypesLast, cfg["VISUALISATION"].getboolean("createRangeTimeDiagram"))
 
-    # TODO: add here different time stepping options
+    # Add different time stepping options here
     log.debug("Use standard time stepping")
     # Initialize time and counters
     nSave = 1
@@ -2074,6 +2084,12 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     # export initial time step
     if cfg["EXPORTS"].getboolean("exportData"):
         exportFields(cfg, t, fields, dem, outDir, cuSimName, TSave="initial")
+
+        if "particles" in resTypes:
+            outDirData = outDir / "particles"
+            fU.makeADir(outDirData)
+            savePartToPickle(particles, outDirData, cuSimName)
+
     # export particles properties for visulation
     if cfg["VISUALISATION"].getboolean("writePartToCSV"):
         particleTools.savePartToCsv(
@@ -2085,10 +2101,6 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
         countParticleCsv = countParticleCsv + 1
 
     # export particles dictionaries of saving time steps
-    # (if particles is not in resType, only first and last time step are saved)
-    outDirData = outDir / "particles"
-    fU.makeADir(outDirData)
-    savePartToPickle(particles, outDirData, cuSimName)
 
     zPartArray0 = copy.deepcopy(particles["z"])
 
@@ -2187,7 +2199,8 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
                 exportFields(cfg, t, fields, dem, outDir, cuSimName, TSave="intermediate")
 
                 # export particles dictionaries of saving time steps
-                savePartToPickle(particles, outDirData, cuSimName)
+                if "particles" in resTypes:
+                    savePartToPickle(particles, outDirData, cuSimName)
 
             # export particles properties for visulation
             if cfg["VISUALISATION"].getboolean("writePartToCSV"):
@@ -2315,7 +2328,8 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
         exportFields(cfg, t, fields, dem, outDir, cuSimName, TSave="final")
 
         # export particles dictionaries of saving time steps
-        savePartToPickle(particles, outDirData, cuSimName)
+        if "particles" in resTypes:
+            savePartToPickle(particles, outDirData, cuSimName)
     else:
         # fetch contourline info
         contourDictXY = outCom1DFA.fetchContCoors(
@@ -2371,7 +2385,7 @@ def setupresultsDF(resTypes, cfgRangeTime):
     resultsDF: dataframe
         data frame with on line for iniital time step and max and mean values of fields
     """
-
+    # TODO catch empty resTypes
     resDict = {"timeStep": [0.0]}
     for resT in resTypes:
         if resT != "particles" and resT != "FTDet":
@@ -2566,7 +2580,7 @@ def computeEulerTimeStep(
         cfg["ResistanceModel"].lower() == "default" and cfg.getboolean("detrainment")
     ):
         # update resistance area fields using thresholds
-        fields = com1DFATools.updateResCoeffFields(fields, cfg, float(particles["t"]), dem)
+        fields = com1DFATools.updateResCoeffFields(fields, cfg)
         if debugPlot:
             outCom1DFA.plotResFields(fields, cfg, particles["tPlot"], dem, particles["mTot"])
 
@@ -2603,7 +2617,6 @@ def computeEulerTimeStep(
     particles = DFAfunC.updatePositionC(cfg, particles, dem, force, fields, typeStop=0)
     tCPUPos = time.time() - startTime
     tCPU["timePos"] = tCPU["timePos"] + tCPUPos
-
     # Split particles
     if cfg.getint("splitOption") == 0:
         # split particles with too much mass
@@ -2968,12 +2981,21 @@ def exportFields(
             # convert from J/cell to kJ/m²
             # (by dividing the peak kinetic energy per cell by the real area of the cell)
             resField = resField * 0.001 / dem["areaRaster"]
+
         dataName = cuSimName + "_" + resType + "_" + "t%.2f" % (timeStep)
         # create directory
         outDirPeak = outDir / "peakFiles" / "timeSteps"
         fU.makeADir(outDirPeak)
         outFile = outDirPeak / dataName
-        IOf.writeResultToRaster(dem["originalHeader"], resField, outFile, flip=True)
+        useCompression = cfg["EXPORTS"].getboolean("useCompression")
+        IOf.writeResultToRaster(
+            dem["originalHeader"], resField, outFile, flip=True, useCompression=useCompression
+        )
+        log.debug(
+            "Results parameter: %s has been exported to Outputs/peakFiles for time step: %.2f "
+            % (resType, timeStep)
+        )
+
         if TSave == "final":
             log.debug(
                 "Results parameter: %s exported to Outputs/peakFiles for time step: %.2f - FINAL time step "
@@ -2984,11 +3006,9 @@ def exportFields(
             outDirPeakAll = outDir / "peakFiles"
             fU.makeADir(outDirPeakAll)
             outFile = outDirPeakAll / dataName
-            IOf.writeResultToRaster(dem["originalHeader"], resField, outFile, flip=True)
-        else:
-            log.debug(
-                "Results parameter: %s has been exported to Outputs/peakFiles for time step: %.2f "
-                % (resType, timeStep)
+            useCompression = cfg["EXPORTS"].getboolean("useCompression")
+            IOf.writeResultToRaster(
+                dem["originalHeader"], resField, outFile, flip=True, useCompression=useCompression
             )
 
 

avaframe/com1DFA/com1DFACfg.ini

@@ -545,4 +545,6 @@ unitpfv = ms-1
 exportData = True
 # export release and optional entrainment raster files derived from shp files saved to Outputs/com1DFA/internalRasters
 exportRasters = False
+# use LZW compression when writing TIFF raster files
+useCompression = True
 

avaframe/com1DFA/com1DFATools.py

@@ -399,7 +399,7 @@ def checkCfgInfoType(cfgInfo):
     return typeCfgInfo
 
 
-def updateResCoeffFields(fields, cfg, t, dem):
+def updateResCoeffFields(fields, cfg):
     """update fields of cRes and detK, coefficients of resistance parameter and detrainment parameter
     according to the thresholds of FV and FT
     if FV OR FT below min thresholds -> apply detrainment in that area
@@ -412,10 +412,6 @@ def updateResCoeffFields(fields, cfg, t, dem):
         dictionary with cResRasterOrig, cResRaster and detRasterOrig, detRaster fields
     cfg: configparser object
         configuration of com1DFA, thresholds
-    t: float
-        time step
-    dem: dict
-        dictionary with info on DEM
 
     Returns
     ------------
@@ -453,10 +449,6 @@ def updateResCoeffFields(fields, cfg, t, dem):
             "Resistance area removed %d cells because FV or FT exceeded %.2f ms-1, %.2f m"
             % (lTh, vMax, thMax)
         )
-        outFileRes = pathlib.Path(
-            cfg["avalancheDir"], "Outputs", "com1DFA", "reports", ("resArea_t%.2f" % t)
-        )
-        IOf.writeResultToRaster(dem["originalHeader"], fields["cResRasterOrig"], outFileRes, flip=True)
 
     # update fields dictionary
     fields["cResRaster"] = cResRaster

avaframe/com1DFA/deriveParameterSet.py

@@ -1017,7 +1017,10 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType):
             raise FileExistsError(message)
 
         # write raster to file
-        outFile = IOf.writeResultToRaster(dem["header"], inputField["rasterData"], outFile, flip=True)
+        useCompression = cfg["EXPORTS"].getboolean("useCompression")
+        outFile = IOf.writeResultToRaster(
+            dem["header"], inputField["rasterData"], outFile, flip=True, useCompression=useCompression
+        )
         log.info("Saved remeshed raster to %s" % outFile)
         returnStr = str(pathlib.Path("remeshedRasters", outFile.name))
         remeshedFlag = "Yes"

avaframe/in2Trans/rasterUtils.py

@@ -1,5 +1,5 @@
 """
-    Raster (ascii and tif) file reader and handler
+Raster (ascii and tif) file reader and handler
 
 """
 
@@ -150,7 +150,7 @@ def isEqualASCheader(headerA, headerB):
     )
 
 
-def writeResultToRaster(header, resultArray, outFileName, flip=False):
+def writeResultToRaster(header, resultArray, outFileName, flip=False, useCompression=True):
     """Write 2D array to a raster file with header and save to location of outFileName
 
     Parameters
@@ -165,39 +165,42 @@ class with methods that give cellsize, nrows, ncols, xllcenter
     flip: boolean
         if True, flip the rows of the resultArray when writing. AF considers the first line in a data array to be the
         southernmost one. Some formats (e.g. tif) have the northernmost line first
+    useCompression: boolean
+        True if compression should be used on writing tiff files (lzw)
+
 
     Returns
     -------
     outFile: path
         to file being written
     """
 
-    if header["driver"] == "AAIGrid":
-        outFile = outFileName.parent / (outFileName.name + ".asc")
-    elif header["driver"] == "GTiff":
-        outFile = outFileName.parent / (outFileName.name + ".tif")
-
-    # try:
-    rasterOut = rasterio.open(
-        outFile,
-        "w",
-        driver=header["driver"],
-        crs=header["crs"],
-        nodata=header["nodata_value"],
-        transform=header["transform"],
-        height=resultArray.shape[0],
-        width=resultArray.shape[1],
-        count=1,
-        dtype=resultArray.dtype,
-        # decimal_precision=3,
-    )
-    if flip:
-        rasterOut.write(np.flipud(resultArray), 1)
-    else:
-        rasterOut.write(resultArray, 1)
-    rasterOut.close()
-    # except:
-    #     log.error("could not write {} to {}".format(resultArray, outFileName))
+    driver = header["driver"]
+    if driver not in ("AAIGrid", "GTiff"):
+        raise ValueError(f"Unsupported driver: {driver}")
+
+    extMap = {"AAIGrid": ".asc", "GTiff": ".tif"}
+    outFile = outFileName.parent / (outFileName.name + extMap[driver])
+
+    commonKwargs = {
+        "driver": driver,
+        "crs": header["crs"],
+        "nodata": header["nodata_value"],
+        "transform": header["transform"],
+        "height": resultArray.shape[0],
+        "width": resultArray.shape[1],
+        "count": 1,
+        "dtype": resultArray.dtype,
+    }
+
+    extraKwargs = {}
+    if useCompression and driver == "GTiff":
+        extraKwargs = {"compress": "lzw"}
+
+    with rasterio.open(outFile, "w", **commonKwargs, **extraKwargs) as rasterOut:
+        data = np.flipud(resultArray) if flip else resultArray
+        rasterOut.write(data, 1)
+
     return outFile
 
 

avaframe/in3Utils/fileHandlerUtils.py

@@ -368,7 +368,7 @@ def splitTimeValueToArrayInterval(cfgValues, endTime):
     Returns
     --------
     items : 1D numpy array
-        time step values as 1D numpy array
+        sorted time step values as 1D numpy array
     """
 
     if ":" in cfgValues:

avaframe/tests/test_com1DFA.py

@@ -1746,6 +1746,7 @@ def test_exportFields(tmp_path):
     cfg = configparser.ConfigParser()
     cfg["GENERAL"] = {"resType": "ppr|pft|FT"}
     cfg["REPORT"] = {"plotFields": "ppr|pft|pfv|pke"}
+    cfg["EXPORTS"] = {"useCompression": "True"}
     Tsave = [0, 10, 15, 25, 40]
     demHeader = {}
     demHeader["cellsize"] = 1

avaframe/tests/test_com1DFATools.py

@@ -119,7 +119,7 @@ def test_updateResCoeffFields(tmp_path):
     dem["originalHeader"]["transform"] = transform
     dem["originalHeader"]["crs"] = rasterio.crs.CRS()
 
-    fields = com1DFATools.updateResCoeffFields(fields, cfg["GENERAL"], 0.0, dem)
+    fields = com1DFATools.updateResCoeffFields(fields, cfg["GENERAL"])
 
     assert fields["cResRaster"].any() == False
 
@@ -135,7 +135,7 @@ def test_updateResCoeffFields(tmp_path):
         "FT": FT,
     }
 
-    fields = com1DFATools.updateResCoeffFields(fields, cfg["GENERAL"], 0.0, dem)
+    fields = com1DFATools.updateResCoeffFields(fields, cfg["GENERAL"])
 
     assert fields["cResRaster"].any() == True
     assert fields["cResRaster"][0, 4] == 1

avaframe/tests/test_deriveParameterSet.py

@@ -608,6 +608,7 @@ def test_checkExtentAndCellSize(tmp_path):
     cfg = configparser.ConfigParser()
     cfg["GENERAL"] = {"resizeThreshold": 3.0, "meshCellSize": 1.0, "meshCellSizeThreshold": 0.0001}
     cfg["GENERAL"]["avalancheDir"] = str(testDir)
+    cfg["EXPORTS"] = {"useCompression": "True"}
     inDir = testDir / "Inputs"
     inDirR = inDir / "RASTERS"
     fU.makeADir(inDirR)