feat(rasterUtils): add optional compression parameter for raster writing

- Introduced `useCompression` parameter to `writeResultToRaster` for configurable compression of `GTiff` files.
- Defaulted the parameter to `True` for LZW compression

refactor(rasterUtils): simplify raster file writing logic and enforce driver validation; fixes #1163

- Added validation to ensure only supported drivers (`AAIGrid`, `GTiff`) are used.
- Utilized `with` statement to handle file writing safely and cleanly.

- Implemented LZW compression for `GTiff` raster outputs to optimize file size.
- Refactored file writing logic to support driver-specific configurations.

Author: fso42

avaframe/com1DFA/com1DFA.py

@@ -2017,8 +2017,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 +2058,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 +2073,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 +2090,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 +2188,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 +2317,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 +2374,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":
@@ -2968,12 +2971,18 @@ 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)
+        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 "
@@ -2985,11 +2994,6 @@ def exportFields(
             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)
-            )
 
 
 def prepareVarSimDict(standardCfg, inputSimFiles, variationDict, simNameExisting="", module=com1DFA):

avaframe/com1DFA/com1DFACfg.ini

@@ -13,6 +13,7 @@ modelType = dfa
 #+++++++++++++ Output++++++++++++
 # desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, dmDet, sfcChange, demAdapted, particles) - separated by |
 resType = ppr|pft|pfv
+
 # saving time step, i.e.  time in seconds (first and last time step are always saved)
 # option 1: give an interval with start:interval in seconds (tStep = 0:5 - this will save desired results every 5 seconds for the full simulation)
 # option 2: explicitly list all desired time steps (closest to actual computational time step) separated by | (example tSteps = 1|50.2|100)

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,41 @@ 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,
+    }
+
+    if useCompression:
+        extraKwargs = {"compress": "lzw"} if driver == "GTiff" else {}
+
+    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: