[com4]: all output rasters have value -9999 in not affected cells

avaframe/com4FlowPy/com4FlowPy.py

@@ -92,6 +92,8 @@ def com4FlowPyMain(cfgPath, cfgSetup):
     modelPaths["uid"] = cfgPath["uid"]
     modelPaths["timeString"] = cfgPath["timeString"]
     modelPaths["outputFileList"] = cfgPath["outputFiles"].split('|')
+    modelPaths["outputNoDataValue"] = cfgPath["outputNoDataValue"]
+    modelPaths["useCompression"] = cfgPath["useCompression"]
 
     modelPaths["outputFileFormat"] = cfgPath["outputFileFormat"]
     if modelPaths["outputFileFormat"] in [".asc", ".ASC"]:
@@ -175,7 +177,7 @@ def com4FlowPyMain(cfgPath, cfgSetup):
 
     # get information on cellsize and nodata value from demHeader
     rasterAttributes = {}
-    
+
     demHeader = IOf.readRasterHeader(modelPaths["demPath"])
     rasterAttributes["nodata"] = demHeader["nodata_value"]
     rasterAttributes["cellsize"] = demHeader["cellsize"]
@@ -511,6 +513,7 @@ def mergeAndWriteResults(modelPaths, modelOptions):
     """
     _uid = modelPaths["uid"]
     _outputs = set(modelPaths['outputFileList'])
+    _outputNoDataValue = modelPaths["outputNoDataValue"]
 
     log.info(" merging results ...")
     log.info("-------------------------")
@@ -543,73 +546,139 @@ def mergeAndWriteResults(modelPaths, modelOptions):
 
     demHeader = IOf.readRasterHeader(modelPaths["demPath"])
     outputHeader = demHeader.copy()
-    outputHeader["nodata_value"]=-9999
+    outputHeader["nodata_value"] = _outputNoDataValue
 
     if _oF == ".asc":
         outputHeader["driver"] = "AAIGrid"
     elif _oF == ".tif":
         outputHeader["driver"] = "GTiff"
 
+    useCompression = modelPaths["useCompression"]
+
     if 'flux' in _outputs:
-        output = IOf.writeResultToRaster(outputHeader, flux,
-                                         modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), flip=True)
+        flux = defineNotAffectedCells(flux, cellCounts, noDataValue=_outputNoDataValue)
+        output = IOf.writeResultToRaster(
+            outputHeader,
+            flux,
+            modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts),
+            flip=True,
+            useCompression=useCompression,
+        )
     if 'zDelta' in _outputs:
-        output = IOf.writeResultToRaster(outputHeader, zDelta,
-                                         modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), flip=True)
+        zDelta = defineNotAffectedCells(zDelta, cellCounts, noDataValue=_outputNoDataValue)
+        output = IOf.writeResultToRaster(
+            outputHeader,
+            zDelta,
+            modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts),
+            flip=True,
+            useCompression=useCompression,
+        )
     if 'cellCounts' in _outputs:
-        output = IOf.writeResultToRaster(outputHeader, cellCounts,
-                                         modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), flip=True)
+        cellCounts = defineNotAffectedCells(cellCounts, cellCounts, noDataValue=_outputNoDataValue)
+        output = IOf.writeResultToRaster(
+            outputHeader,
+            cellCounts,
+            modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts),
+            flip=True,
+            useCompression=useCompression,
+        )
     if 'zDeltaSum' in _outputs:
-        output = IOf.writeResultToRaster(outputHeader, zDeltaSum,
-                                         modelPaths["resDir"] / "com4_{}_{}_zDeltaSum".format(_uid, _ts), flip=True)
+        zDeltaSum = defineNotAffectedCells(zDeltaSum, cellCounts, noDataValue=_outputNoDataValue)
+        output = IOf.writeResultToRaster(
+            outputHeader,
+            zDeltaSum,
+            modelPaths["resDir"] / "com4_{}_{}_zDeltaSum".format(_uid, _ts),
+            flip=True,
+            useCompression=useCompression,
+        )
     if 'routFluxSum' in _outputs:
-        output = IOf.writeResultToRaster(outputHeader, routFluxSum,
-                                         modelPaths["resDir"] / "com4_{}_{}_routFluxSum".format(_uid, _ts), flip=True)
+        routFluxSum = defineNotAffectedCells(routFluxSum, cellCounts, noDataValue=_outputNoDataValue)
+        output = IOf.writeResultToRaster(
+            outputHeader,
+            routFluxSum,
+            modelPaths["resDir"] / "com4_{}_{}_routFluxSum".format(_uid, _ts),
+            flip=True,
+            useCompression=useCompression,
+        )
     if 'depFluxSum' in _outputs:
-        output = IOf.writeResultToRaster(outputHeader, depFluxSum,
-                                         modelPaths["resDir"] / "com4_{}_{}_depFluxSum".format(_uid, _ts), flip=True)
+        depFluxSum = defineNotAffectedCells(depFluxSum, cellCounts, noDataValue=_outputNoDataValue)
+        output = IOf.writeResultToRaster(
+            outputHeader,
+            depFluxSum,
+            modelPaths["resDir"] / "com4_{}_{}_depFluxSum".format(_uid, _ts),
+            flip=True,
+            useCompression=useCompression,
+        )
     if "fpTravelAngle" in _outputs or "fpTravelAngleMax" in _outputs:
+        fpTaMax = defineNotAffectedCells(fpTaMax, cellCounts, noDataValue=_outputNoDataValue)
         output = IOf.writeResultToRaster(
             outputHeader,
             fpTaMax,
             modelPaths["resDir"] / "com4_{}_{}_fpTravelAngleMax".format(_uid, _ts),
             flip=True,
+            useCompression=useCompression,
         )
     if "fpTravelAngleMin" in _outputs:
+        fpTaMin = defineNotAffectedCells(fpTaMin, cellCounts, noDataValue=_outputNoDataValue)
         output = IOf.writeResultToRaster(
             outputHeader,
             fpTaMin,
             modelPaths["resDir"] / "com4_{}_{}_fpTravelAngleMin".format(_uid, _ts),
             flip=True,
+            useCompression=useCompression,
         )
     if 'slTravelAngle' in _outputs:
-        output = IOf.writeResultToRaster(outputHeader, slTa,
-                                         modelPaths["resDir"] / "com4_{}_{}_slTravelAngle".format(_uid, _ts), flip=True)
+        slTa = defineNotAffectedCells(slTa, cellCounts, noDataValue=_outputNoDataValue)
+        output = IOf.writeResultToRaster(
+            outputHeader,
+            slTa,
+            modelPaths["resDir"] / "com4_{}_{}_slTravelAngle".format(_uid, _ts),
+            flip=True,
+            useCompression=useCompression,
+        )
     if "travelLength" in _outputs or "travelLengthMax" in _outputs:
+        travelLengthMax = defineNotAffectedCells(travelLengthMax, cellCounts, noDataValue=_outputNoDataValue)
         output = IOf.writeResultToRaster(
             outputHeader,
             travelLengthMax,
             modelPaths["resDir"] / "com4_{}_{}_travelLengthMax".format(_uid, _ts),
             flip=True,
+            useCompression=useCompression,
         )
     if "travelLengthMin" in _outputs:
+        travelLengthMin = defineNotAffectedCells(travelLengthMin, cellCounts, noDataValue=_outputNoDataValue)
         output = IOf.writeResultToRaster(
             outputHeader,
             travelLengthMin,
             modelPaths["resDir"] / "com4_{}_{}_travelLengthMin".format(_uid, _ts),
             flip=True,
+            useCompression=useCompression,
         )
 
     # NOTE:
     # if not modelOptions["infraBool"]:  # if no infra
     # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("cell_counts%s" %(output_format)),cell_counts)
     # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("z_delta_sum%s" %(output_format)),z_delta_sum)
     if modelOptions["infraBool"]:  # if infra
-        output = IOf.writeResultToRaster(outputHeader, backcalc,
-                                         modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts), flip=True)
+        backcalc = defineNotAffectedCells(backcalc, cellCounts, noDataValue=_outputNoDataValue)
+        output = IOf.writeResultToRaster(
+            outputHeader,
+            backcalc,
+            modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts),
+            flip=True,
+            useCompression=useCompression,
+        )
     if modelOptions["forestInteraction"]:
-        output = IOf.writeResultToRaster(outputHeader, forestInteraction,
-                                         modelPaths["resDir"] / "com4_{}_{}_forestInteraction".format(_uid, _ts), flip=True)
+        forestInteraction = defineNotAffectedCells(
+            forestInteraction, cellCounts, noDataValue=_outputNoDataValue
+        )
+        output = IOf.writeResultToRaster(
+            outputHeader,
+            forestInteraction,
+            modelPaths["resDir"] / "com4_{}_{}_forestInteraction".format(_uid, _ts),
+            flip=True,
+            useCompression=useCompression,
+        )
     del output
 
 
@@ -647,7 +716,12 @@ def checkConvertReleaseShp2Tif(modelPaths):
         # NOTE: nodata_value does not matter anyways - since gT.prepareArea() returns a
         # raster with values of '0' indicating 'no release area' - so there actually should
         # not be any 'no_data' values in the created release raster file
-        IOf.writeResultToRaster(demHeader, releaseArea, modelPaths["workDir"] / "release")
+        IOf.writeResultToRaster(
+            demHeader,
+            releaseArea,
+            modelPaths["workDir"] / "release",
+            useCompression=modelPaths["useCompression"],
+        )
         del releaseLine
     else:
         modelPaths["releasePathWork"] = modelPaths["releasePath"]
@@ -701,3 +775,25 @@ def deleteTempFolder(tempFolderPath):
         log.info(" isDir:{} isTemp:{}}".format(isDir, validTemp))
 
     log.info("+++++++++++++++++++++++")
+
+
+def defineNotAffectedCells(raster, affectedCells, noDataValue=-9999):
+    """
+    define not affected cells as -9999
+
+    Parameters
+    -----------
+    raster: np.array
+        raster whose not affected cells are specified
+    affectedCells: np.array
+        mask for affected cells
+    noDataValue: float
+        value for not affected cells (default: -9999)
+
+    Returns
+    -----------
+    raster: np. array
+        raster with not affected cells have the value noDataValue
+    """
+    raster[affectedCells <= 0] = noDataValue
+    return raster

avaframe/com4FlowPy/com4FlowPyCfg.ini

@@ -204,6 +204,14 @@ maxChunks = 500
 # then you should choose '.asc' here!
 outputFileFormat = .tif
 
+# use LZW compression when writing tif raster files
+useCompression = True
+
+# define noData value that is assigned in output rasters in cells that are not affected by the process (default: -9999)
+# when changing this value BE AWARE that noData values can have same values as a affected cells
+# (e.g., when using outputNoDataValue = 0)
+outputNoDataValue = -9999
+
 # define the different output files that are written to disk
 # default = 'zDelta|cellCounts|travelLengthMax|fpTravelAngleMax'
 # additional options: 

avaframe/com4FlowPy/flowCore.py

@@ -325,7 +325,7 @@ def run(optTuple):
             fpTravelAngleMaxArray = np.maximum(fpTravelAngleMaxArray, fpTravelAngleMaxList[i])
         if "fpTravelAngleMin" in outputs:
             fpTravelAngleMinArray = np.where(
-                (travelLengthMinArray >= 0) & (fpTravelAngleMinList[i] >= 0),
+                (fpTravelAngleMinArray >= 0) & (fpTravelAngleMinList[i] >= 0),
                 np.minimum(fpTravelAngleMinArray, fpTravelAngleMinList[i]),
                 np.maximum(fpTravelAngleMinArray, fpTravelAngleMinList[i]),
             )

avaframe/runCom4FlowPy.py

@@ -111,6 +111,8 @@ def main(avalancheDir=''):
         cfgPath["uid"] = uid
         cfgPath["timeString"] = timeString
         cfgPath["outputFiles"] = cfgCustomPaths["outputFiles"]
+        cfgPath["outputNoDataValue"] = cfgCustomPaths.getfloat("outputNoDataValue")
+        cfgPath["useCompression"] = cfgCustomPaths.getboolean("useCompression")
 
         com4FlowPy.com4FlowPyMain(cfgPath, cfgSetup)
 
@@ -169,6 +171,8 @@ def main(avalancheDir=''):
         cfgPath["deleteTemp"] = cfgCustomPaths["deleteTempFolder"]
         cfgPath["outputFileFormat"] = cfgCustomPaths["outputFileFormat"]
         cfgPath["outputFiles"] = cfgCustomPaths["outputFiles"]
+        cfgPath["outputNoDataValue"] = cfgCustomPaths.getfloat("outputNoDataValue")
+        cfgPath["useCompression"] = cfgCustomPaths.getboolean("useCompression")
 
         cfgSetup["cpuCount"] = str(cfgUtils.getNumberOfProcesses(cfgMain, 9999))
         cfgPath["customDirs"] = cfgCustomPaths["useCustomPaths"]
@@ -211,17 +215,17 @@ def readFlowPyinputs(avalancheDir, cfgFlowPy, log):
     avalancheDir = pathlib.Path(avalancheDir)
 
     inputDir = avalancheDir / "Inputs"
-    relFile, available = gI.getAndCheckInputFiles(inputDir, "REL", "Release Area", fileExt="shp")
-    
+    relFile, available, _ = gI.getAndCheckInputFiles(inputDir, "REL", "Release Area", fileExt="shp")
+
     if available == "No":
-        relFile, available = gI.getAndCheckInputFiles(inputDir, "REL", "Release Area", fileExt="raster")
+        relFile, available, _ = gI.getAndCheckInputFiles(inputDir, "REL", "Release Area", fileExt="raster")
     if available == "No":
         message = f"There is no release area file in supported format provided in {avalancheDir}/REL"
         log.error(message)
         raise AssertionError(message)
     log.info("Release area file is: %s" % relFile)
     cfgPath["releasePath"] = relFile
-    
+
     # TODO: also use the getAndCheckInputFiles to get the paths for the following files?
     # read infra area
     if cfgFlowPy.getboolean("GENERAL", "infra") is True:
@@ -237,7 +241,7 @@ def readFlowPyinputs(avalancheDir, cfgFlowPy, log):
 
     # read uMax Limit Raster
     if cfgFlowPy.getboolean("GENERAL", "variableUmaxLim") is True:
-        varUmaxPath, available = gI.getAndCheckInputFiles(inputDir, "UMAX", "Umax", fileExt="raster")
+        varUmaxPath, available, _ = gI.getAndCheckInputFiles(inputDir, "UMAX", "Umax", fileExt="raster")
         if available == "No":
             message = f"There is no variable UMAX file in supported format provided in {avalancheDir}/UMAX"
             log.error(message)
@@ -250,7 +254,7 @@ def readFlowPyinputs(avalancheDir, cfgFlowPy, log):
     # read variable Alpha Angle Raster
 
     if cfgFlowPy.getboolean("GENERAL", "variableAlpha") is True:
-        varAlphaPath, available = gI.getAndCheckInputFiles(inputDir, "ALPHA", "alpha", fileExt="raster")
+        varAlphaPath, available, _ = gI.getAndCheckInputFiles(inputDir, "ALPHA", "alpha", fileExt="raster")
         if available == "No":
             message = f"There is no variable ALPHA file in supported format provided in {avalancheDir}/ALPHA"
             log.error(message)
@@ -263,7 +267,7 @@ def readFlowPyinputs(avalancheDir, cfgFlowPy, log):
     # read variable Exponent Raster
 
     if cfgFlowPy.getboolean("GENERAL", "variableExponent") is True:
-        varExponentPath, available = gI.getAndCheckInputFiles(inputDir, "EXP", "exp", fileExt="raster")
+        varExponentPath, available, _ = gI.getAndCheckInputFiles(inputDir, "EXP", "exp", fileExt="raster")
         if available == "No":
             message = f"There is no variable EXPONENT file in supported format provided in {avalancheDir}/EXP"
             log.error(message)
@@ -277,7 +281,7 @@ def readFlowPyinputs(avalancheDir, cfgFlowPy, log):
     # TODO: should we also allow forest as shp - file?
 
     if cfgFlowPy.getboolean("GENERAL", "forest") is True:
-        forestPath, available = gI.getAndCheckInputFiles(inputDir, "RES", "forest", fileExt="raster")
+        forestPath, available, _ = gI.getAndCheckInputFiles(inputDir, "RES", "forest", fileExt="raster")
         if available == "No":
             message = f"There is no forest file in supported format provided in {avalancheDir}/RES"
             log.error(message)
@@ -369,4 +373,4 @@ def writeCfgJSON(cfg, uid, workDir):
 
 
 if __name__ == "__main__":
-    main()
+    main()

avaframe/tests/test_com4FlowPy.py

@@ -32,23 +32,23 @@ def test_reverseTopology():
                       4: [],
                       5: [6],
                       6: []}
-    
+
     testGraphReverse = {0: [],
                         1: [0],
                         2: [0,1,3],
                         3: [0],
                         4: [2,1],
                         5: [2],
                         6: [5]}
-    
+
     reverseGraphCalc = flowCore.reverseTopology(testGraph)
 
     for key, item in reverseGraphCalc.items():
         assert key in testGraphReverse.keys()
         setTestChildren = set(item)
         setCalcChildren = set(testGraphReverse[key])
         assert setTestChildren == setCalcChildren
-    
+
     testGraph = {0:[1,2,3],
                  1:[4],
                  2:[5,6],
@@ -62,7 +62,7 @@ def test_reverseTopology():
                  10:[12],
                  11:[],
                  12:[]}
-    
+
     testGraphReverse = {0:[],
                         1:[0],
                         2:[0],
@@ -76,15 +76,15 @@ def test_reverseTopology():
                         10:[5,6],
                         11:[7],
                         12:[10]}
-    
+
     reverseGraphCalc = flowCore.reverseTopology(testGraph)
 
     for key, item in reverseGraphCalc.items():
         assert key in testGraphReverse.keys()
         setTestChildren = set(item)
         setCalcChildren = set(testGraphReverse[key])
         assert setTestChildren == setCalcChildren
-    
+
 
 def test_backTracking():
     '''