From 7d430d3ede4aa0f4b44a7dbe7b6b069e1aaab028 Mon Sep 17 00:00:00 2001 From: paulasp Date: Wed, 9 Jul 2025 17:57:13 +0200 Subject: [PATCH 1/4] all output rasters have value -9999 in not affected cells little bug --- avaframe/com4FlowPy/com4FlowPy.py | 28 ++++++++++++++++++++++++++++ avaframe/com4FlowPy/flowCore.py | 12 +++++++++++- docs/moduleCom4FlowPy.rst | 1 + 3 files changed, 40 insertions(+), 1 deletion(-) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index 146eaa63d..1fba68d09 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -516,6 +516,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): log.info("-------------------------") # Merge calculated tiles + affectedCells = SPAM.mergeRaster(modelPaths["tempDir"], "res_affected") zDelta = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta") flux = SPAM.mergeRaster(modelPaths["tempDir"], "res_flux") cellCounts = SPAM.mergeRaster(modelPaths["tempDir"], "res_count", method="sum") @@ -554,18 +555,23 @@ def mergeAndWriteResults(modelPaths, modelOptions): output = IOf.writeResultToRaster(outputHeader, flux, modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), flip=True) if 'zDelta' in _outputs: + defineNotAffectedCells(zDelta, affectedCells) output = IOf.writeResultToRaster(outputHeader, zDelta, modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), flip=True) if 'cellCounts' in _outputs: + defineNotAffectedCells(cellCounts, affectedCells) output = IOf.writeResultToRaster(outputHeader, cellCounts, modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), flip=True) if 'zDeltaSum' in _outputs: + defineNotAffectedCells(zDeltaSum, affectedCells) output = IOf.writeResultToRaster(outputHeader, zDeltaSum, modelPaths["resDir"] / "com4_{}_{}_zDeltaSum".format(_uid, _ts), flip=True) if 'routFluxSum' in _outputs: + defineNotAffectedCells(routFluxSum, affectedCells) output = IOf.writeResultToRaster(outputHeader, routFluxSum, modelPaths["resDir"] / "com4_{}_{}_routFluxSum".format(_uid, _ts), flip=True) if 'depFluxSum' in _outputs: + defineNotAffectedCells(depFluxSum, affectedCells) output = IOf.writeResultToRaster(outputHeader, depFluxSum, modelPaths["resDir"] / "com4_{}_{}_depFluxSum".format(_uid, _ts), flip=True) if "fpTravelAngle" in _outputs or "fpTravelAngleMax" in _outputs: @@ -701,3 +707,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 diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index 731195178..888284119 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -273,6 +273,7 @@ def run(optTuple): slTravelAngleArray = np.ones_like(dem, dtype=np.float32) * -9999 travelLengthMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 + affectedCellsArray = np.zeros_like(dem, dtype=np.int32) if forestInteraction: forestIntArray = np.ones_like(dem, dtype=np.float32) * -9999 @@ -289,6 +290,7 @@ def run(optTuple): slTravelAngleList = [] travelLengthMaxList = [] travelLengthMinList = [] + affectedCellsList = [] if forestInteraction: forestIntList = [] @@ -308,11 +310,13 @@ def run(optTuple): fpTravelAngleMinList.append(res[9]) routFluxSumList.append(res[10]) depFluxSumList.append(res[11]) + affectedCellsList.append(res[12]) if forestInteraction: - forestIntList.append(res[12]) + forestIntList.append(res[13]) logging.info("Calculation finished, getting results.") for i in range(len(zDeltaList)): + affectedCellsArray = np.maximum(affectedCellsArray, affectedCellsList[i]) zDeltaArray = np.maximum(zDeltaArray, zDeltaList[i]) fluxArray = np.maximum(fluxArray, fluxList[i]) countArray += ccList[i] @@ -344,6 +348,7 @@ def run(optTuple): np.maximum(forestIntArray, forestIntList[i])) # Save Calculated tiles + np.save(tempDir / ("res_affected_%s_%s" % (optTuple[0], optTuple[1])), affectedCellsArray) np.save(tempDir / ("res_z_delta_%s_%s" % (optTuple[0], optTuple[1])), zDeltaArray) np.save(tempDir / ("res_z_delta_sum_%s_%s" % (optTuple[0], optTuple[1])), zDeltaSumArray) np.save(tempDir / ("res_rout_flux_sum_%s_%s" % (optTuple[0], optTuple[1])), routFluxSumArray) @@ -480,6 +485,8 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 travelLengthMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 + affectedCellsArray = np.zeros_like(dem, dtype=np.int32) + if infraBool: backcalc = np.ones_like(dem, dtype=np.int32) * -9999 else: @@ -610,6 +617,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): forestParams=forestParams, )) + affectedCellsArray[cell.rowindex, cell.colindex] = 1 zDeltaArray[cell.rowindex, cell.colindex] = max(zDeltaArray[cell.rowindex, cell.colindex], cell.z_delta) fluxArray[cell.rowindex, cell.colindex] = max(fluxArray[cell.rowindex, cell.colindex], cell.flux) routFluxSumArray[cell.rowindex, cell.colindex] += cell.flux @@ -698,6 +706,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): fpTravelAngleMinArray, routFluxSumArray, depFluxSumArray, + affectedCellsArray, forestIntArray, ) else: @@ -714,6 +723,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): fpTravelAngleMinArray, routFluxSumArray, depFluxSumArray, + affectedCellsArray, ) diff --git a/docs/moduleCom4FlowPy.rst b/docs/moduleCom4FlowPy.rst index 28514dcff..335186fae 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -224,6 +224,7 @@ Output All outputs are written in *'.tif'* or in *'.asc'* raster format (controlable via the ``outputFileFormat`` option in ``(local_)com4FlowPyCfg.ini``, default is *'.tif'*) in the same resolution and extent as the input raster layers. You can customize which output rasters are written at the end of the model run by selecting the desired output files through the ``outputFiles`` option in ``(local_)com4FlowPyCfg.ini``. +In the output rasters, the cells, that are not affected by the process, are specified by the value -9999. By default the following four output layers are written to disk at the end of the model run: From d2cb0c0dc13da7f1deef7380fd6228feb91df014 Mon Sep 17 00:00:00 2001 From: paulasp Date: Wed, 6 Aug 2025 15:04:45 +0200 Subject: [PATCH 2/4] new parameter to define value for not affected cells --- avaframe/com4FlowPy/com4FlowPy.py | 28 +++++++++++++++++++++------ avaframe/com4FlowPy/com4FlowPyCfg.ini | 5 +++++ avaframe/runCom4FlowPy.py | 4 +++- 3 files changed, 30 insertions(+), 7 deletions(-) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index 1fba68d09..626bf33a3 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -92,6 +92,7 @@ def com4FlowPyMain(cfgPath, cfgSetup): modelPaths["uid"] = cfgPath["uid"] modelPaths["timeString"] = cfgPath["timeString"] modelPaths["outputFileList"] = cfgPath["outputFiles"].split('|') + modelPaths["outputNoDataValue"] = cfgPath["outputNoDataValue"] modelPaths["outputFileFormat"] = cfgPath["outputFileFormat"] if modelPaths["outputFileFormat"] in [".asc", ".ASC"]: @@ -175,7 +176,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 +512,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): """ _uid = modelPaths["uid"] _outputs = set(modelPaths['outputFileList']) + _outputNoDataValue = modelPaths["outputNoDataValue"] log.info(" merging results ...") log.info("-------------------------") @@ -552,29 +554,31 @@ def mergeAndWriteResults(modelPaths, modelOptions): outputHeader["driver"] = "GTiff" if 'flux' in _outputs: + flux = defineNotAffectedCells(flux, affectedCells, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, flux, modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), flip=True) if 'zDelta' in _outputs: - defineNotAffectedCells(zDelta, affectedCells) + zDelta = defineNotAffectedCells(zDelta, affectedCells, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, zDelta, modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), flip=True) if 'cellCounts' in _outputs: - defineNotAffectedCells(cellCounts, affectedCells) + cellCounts = defineNotAffectedCells(cellCounts, affectedCells, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, cellCounts, modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), flip=True) if 'zDeltaSum' in _outputs: - defineNotAffectedCells(zDeltaSum, affectedCells) + zDeltaSum = defineNotAffectedCells(zDeltaSum, affectedCells, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, zDeltaSum, modelPaths["resDir"] / "com4_{}_{}_zDeltaSum".format(_uid, _ts), flip=True) if 'routFluxSum' in _outputs: - defineNotAffectedCells(routFluxSum, affectedCells) + routFluxSum = defineNotAffectedCells(routFluxSum, affectedCells, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, routFluxSum, modelPaths["resDir"] / "com4_{}_{}_routFluxSum".format(_uid, _ts), flip=True) if 'depFluxSum' in _outputs: - defineNotAffectedCells(depFluxSum, affectedCells) + depFluxSum = defineNotAffectedCells(depFluxSum, affectedCells, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, depFluxSum, modelPaths["resDir"] / "com4_{}_{}_depFluxSum".format(_uid, _ts), flip=True) if "fpTravelAngle" in _outputs or "fpTravelAngleMax" in _outputs: + fpTaMax = defineNotAffectedCells(fpTaMax, affectedCells, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, fpTaMax, @@ -582,6 +586,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, ) if "fpTravelAngleMin" in _outputs: + fpTaMin = defineNotAffectedCells(fpTaMin, affectedCells, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, fpTaMin, @@ -589,9 +594,13 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, ) if 'slTravelAngle' in _outputs: + slTa = defineNotAffectedCells(slTa, affectedCells, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, slTa, modelPaths["resDir"] / "com4_{}_{}_slTravelAngle".format(_uid, _ts), flip=True) if "travelLength" in _outputs or "travelLengthMax" in _outputs: + travelLengthMax = defineNotAffectedCells( + travelLengthMax, affectedCells, noDataValue=_outputNoDataValue + ) output = IOf.writeResultToRaster( outputHeader, travelLengthMax, @@ -599,6 +608,9 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, ) if "travelLengthMin" in _outputs: + travelLengthMin = defineNotAffectedCells( + travelLengthMin, affectedCells, noDataValue=_outputNoDataValue + ) output = IOf.writeResultToRaster( outputHeader, travelLengthMin, @@ -611,9 +623,13 @@ def mergeAndWriteResults(modelPaths, modelOptions): # 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 + backcalc = defineNotAffectedCells(backcalc, affectedCells, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, backcalc, modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts), flip=True) if modelOptions["forestInteraction"]: + forestInteraction = defineNotAffectedCells( + forestInteraction, affectedCells, noDataValue=_outputNoDataValue + ) output = IOf.writeResultToRaster(outputHeader, forestInteraction, modelPaths["resDir"] / "com4_{}_{}_forestInteraction".format(_uid, _ts), flip=True) del output diff --git a/avaframe/com4FlowPy/com4FlowPyCfg.ini b/avaframe/com4FlowPy/com4FlowPyCfg.ini index 19e17d887..54baab0a9 100644 --- a/avaframe/com4FlowPy/com4FlowPyCfg.ini +++ b/avaframe/com4FlowPy/com4FlowPyCfg.ini @@ -204,6 +204,11 @@ maxChunks = 500 # then you should choose '.asc' here! outputFileFormat = .tif +# 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: diff --git a/avaframe/runCom4FlowPy.py b/avaframe/runCom4FlowPy.py index 3fbde5a2d..b709874a4 100644 --- a/avaframe/runCom4FlowPy.py +++ b/avaframe/runCom4FlowPy.py @@ -111,6 +111,7 @@ def main(avalancheDir=''): cfgPath["uid"] = uid cfgPath["timeString"] = timeString cfgPath["outputFiles"] = cfgCustomPaths["outputFiles"] + cfgPath["outputNoDataValue"] = cfgCustomPaths["outputNoDataValue"] com4FlowPy.com4FlowPyMain(cfgPath, cfgSetup) @@ -169,6 +170,7 @@ def main(avalancheDir=''): cfgPath["deleteTemp"] = cfgCustomPaths["deleteTempFolder"] cfgPath["outputFileFormat"] = cfgCustomPaths["outputFileFormat"] cfgPath["outputFiles"] = cfgCustomPaths["outputFiles"] + cfgPath["outputNoDataValue"] = cfgCustomPaths.getfloat("outputNoDataValue") cfgSetup["cpuCount"] = str(cfgUtils.getNumberOfProcesses(cfgMain, 9999)) cfgPath["customDirs"] = cfgCustomPaths["useCustomPaths"] @@ -369,4 +371,4 @@ def writeCfgJSON(cfg, uid, workDir): if __name__ == "__main__": - main() \ No newline at end of file + main() From 90faf626bfa577d42091e27ad2a6947962bfe8f1 Mon Sep 17 00:00:00 2001 From: paulasp Date: Wed, 6 Aug 2025 15:21:18 +0200 Subject: [PATCH 3/4] adapt test Update doc compute affected cells using cellCounts fixed typo note for tilesize --- avaframe/com4FlowPy/com4FlowPy.py | 33 +++++++++++++------------------ avaframe/com4FlowPy/flowCore.py | 14 ++----------- avaframe/tests/test_com4FlowPy.py | 12 +++++------ docs/moduleCom4FlowPy.rst | 5 +++-- 4 files changed, 25 insertions(+), 39 deletions(-) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index 626bf33a3..3b360fcde 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -518,7 +518,6 @@ def mergeAndWriteResults(modelPaths, modelOptions): log.info("-------------------------") # Merge calculated tiles - affectedCells = SPAM.mergeRaster(modelPaths["tempDir"], "res_affected") zDelta = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta") flux = SPAM.mergeRaster(modelPaths["tempDir"], "res_flux") cellCounts = SPAM.mergeRaster(modelPaths["tempDir"], "res_count", method="sum") @@ -554,31 +553,31 @@ def mergeAndWriteResults(modelPaths, modelOptions): outputHeader["driver"] = "GTiff" if 'flux' in _outputs: - flux = defineNotAffectedCells(flux, affectedCells, noDataValue=_outputNoDataValue) + flux = defineNotAffectedCells(flux, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, flux, modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), flip=True) if 'zDelta' in _outputs: - zDelta = defineNotAffectedCells(zDelta, affectedCells, noDataValue=_outputNoDataValue) + zDelta = defineNotAffectedCells(zDelta, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, zDelta, modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), flip=True) if 'cellCounts' in _outputs: - cellCounts = defineNotAffectedCells(cellCounts, affectedCells, noDataValue=_outputNoDataValue) + cellCounts = defineNotAffectedCells(cellCounts, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, cellCounts, modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), flip=True) if 'zDeltaSum' in _outputs: - zDeltaSum = defineNotAffectedCells(zDeltaSum, affectedCells, noDataValue=_outputNoDataValue) + zDeltaSum = defineNotAffectedCells(zDeltaSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, zDeltaSum, modelPaths["resDir"] / "com4_{}_{}_zDeltaSum".format(_uid, _ts), flip=True) if 'routFluxSum' in _outputs: - routFluxSum = defineNotAffectedCells(routFluxSum, affectedCells, noDataValue=_outputNoDataValue) + routFluxSum = defineNotAffectedCells(routFluxSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, routFluxSum, modelPaths["resDir"] / "com4_{}_{}_routFluxSum".format(_uid, _ts), flip=True) if 'depFluxSum' in _outputs: - depFluxSum = defineNotAffectedCells(depFluxSum, affectedCells, noDataValue=_outputNoDataValue) + depFluxSum = defineNotAffectedCells(depFluxSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, depFluxSum, modelPaths["resDir"] / "com4_{}_{}_depFluxSum".format(_uid, _ts), flip=True) if "fpTravelAngle" in _outputs or "fpTravelAngleMax" in _outputs: - fpTaMax = defineNotAffectedCells(fpTaMax, affectedCells, noDataValue=_outputNoDataValue) + fpTaMax = defineNotAffectedCells(fpTaMax, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, fpTaMax, @@ -586,7 +585,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, ) if "fpTravelAngleMin" in _outputs: - fpTaMin = defineNotAffectedCells(fpTaMin, affectedCells, noDataValue=_outputNoDataValue) + fpTaMin = defineNotAffectedCells(fpTaMin, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, fpTaMin, @@ -594,13 +593,11 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, ) if 'slTravelAngle' in _outputs: - slTa = defineNotAffectedCells(slTa, affectedCells, noDataValue=_outputNoDataValue) + slTa = defineNotAffectedCells(slTa, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, slTa, modelPaths["resDir"] / "com4_{}_{}_slTravelAngle".format(_uid, _ts), flip=True) if "travelLength" in _outputs or "travelLengthMax" in _outputs: - travelLengthMax = defineNotAffectedCells( - travelLengthMax, affectedCells, noDataValue=_outputNoDataValue - ) + travelLengthMax = defineNotAffectedCells(travelLengthMax, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, travelLengthMax, @@ -608,9 +605,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, ) if "travelLengthMin" in _outputs: - travelLengthMin = defineNotAffectedCells( - travelLengthMin, affectedCells, noDataValue=_outputNoDataValue - ) + travelLengthMin = defineNotAffectedCells(travelLengthMin, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, travelLengthMin, @@ -623,12 +618,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): # 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 - backcalc = defineNotAffectedCells(backcalc, affectedCells, noDataValue=_outputNoDataValue) + backcalc = defineNotAffectedCells(backcalc, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster(outputHeader, backcalc, modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts), flip=True) if modelOptions["forestInteraction"]: forestInteraction = defineNotAffectedCells( - forestInteraction, affectedCells, noDataValue=_outputNoDataValue + forestInteraction, cellCounts, noDataValue=_outputNoDataValue ) output = IOf.writeResultToRaster(outputHeader, forestInteraction, modelPaths["resDir"] / "com4_{}_{}_forestInteraction".format(_uid, _ts), flip=True) @@ -743,5 +738,5 @@ def defineNotAffectedCells(raster, affectedCells, noDataValue=-9999): raster: np. array raster with not affected cells have the value noDataValue """ - raster[affectedCells == 0] = noDataValue + raster[affectedCells <= 0] = noDataValue return raster diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index 888284119..a9a2315f7 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -273,7 +273,6 @@ def run(optTuple): slTravelAngleArray = np.ones_like(dem, dtype=np.float32) * -9999 travelLengthMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 - affectedCellsArray = np.zeros_like(dem, dtype=np.int32) if forestInteraction: forestIntArray = np.ones_like(dem, dtype=np.float32) * -9999 @@ -290,7 +289,6 @@ def run(optTuple): slTravelAngleList = [] travelLengthMaxList = [] travelLengthMinList = [] - affectedCellsList = [] if forestInteraction: forestIntList = [] @@ -310,13 +308,11 @@ def run(optTuple): fpTravelAngleMinList.append(res[9]) routFluxSumList.append(res[10]) depFluxSumList.append(res[11]) - affectedCellsList.append(res[12]) if forestInteraction: - forestIntList.append(res[13]) + forestIntList.append(res[12]) logging.info("Calculation finished, getting results.") for i in range(len(zDeltaList)): - affectedCellsArray = np.maximum(affectedCellsArray, affectedCellsList[i]) zDeltaArray = np.maximum(zDeltaArray, zDeltaList[i]) fluxArray = np.maximum(fluxArray, fluxList[i]) countArray += ccList[i] @@ -329,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]), ) @@ -348,7 +344,6 @@ def run(optTuple): np.maximum(forestIntArray, forestIntList[i])) # Save Calculated tiles - np.save(tempDir / ("res_affected_%s_%s" % (optTuple[0], optTuple[1])), affectedCellsArray) np.save(tempDir / ("res_z_delta_%s_%s" % (optTuple[0], optTuple[1])), zDeltaArray) np.save(tempDir / ("res_z_delta_sum_%s_%s" % (optTuple[0], optTuple[1])), zDeltaSumArray) np.save(tempDir / ("res_rout_flux_sum_%s_%s" % (optTuple[0], optTuple[1])), routFluxSumArray) @@ -485,8 +480,6 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 travelLengthMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 - affectedCellsArray = np.zeros_like(dem, dtype=np.int32) - if infraBool: backcalc = np.ones_like(dem, dtype=np.int32) * -9999 else: @@ -617,7 +610,6 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): forestParams=forestParams, )) - affectedCellsArray[cell.rowindex, cell.colindex] = 1 zDeltaArray[cell.rowindex, cell.colindex] = max(zDeltaArray[cell.rowindex, cell.colindex], cell.z_delta) fluxArray[cell.rowindex, cell.colindex] = max(fluxArray[cell.rowindex, cell.colindex], cell.flux) routFluxSumArray[cell.rowindex, cell.colindex] += cell.flux @@ -706,7 +698,6 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): fpTravelAngleMinArray, routFluxSumArray, depFluxSumArray, - affectedCellsArray, forestIntArray, ) else: @@ -723,7 +714,6 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): fpTravelAngleMinArray, routFluxSumArray, depFluxSumArray, - affectedCellsArray, ) diff --git a/avaframe/tests/test_com4FlowPy.py b/avaframe/tests/test_com4FlowPy.py index aabd8c4d1..8c4321149 100644 --- a/avaframe/tests/test_com4FlowPy.py +++ b/avaframe/tests/test_com4FlowPy.py @@ -32,7 +32,7 @@ def test_reverseTopology(): 4: [], 5: [6], 6: []} - + testGraphReverse = {0: [], 1: [0], 2: [0,1,3], @@ -40,7 +40,7 @@ def test_reverseTopology(): 4: [2,1], 5: [2], 6: [5]} - + reverseGraphCalc = flowCore.reverseTopology(testGraph) for key, item in reverseGraphCalc.items(): @@ -48,7 +48,7 @@ def test_reverseTopology(): 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,7 +76,7 @@ def test_reverseTopology(): 10:[5,6], 11:[7], 12:[10]} - + reverseGraphCalc = flowCore.reverseTopology(testGraph) for key, item in reverseGraphCalc.items(): @@ -84,7 +84,7 @@ def test_reverseTopology(): setTestChildren = set(item) setCalcChildren = set(testGraphReverse[key]) assert setTestChildren == setCalcChildren - + def test_backTracking(): ''' diff --git a/docs/moduleCom4FlowPy.rst b/docs/moduleCom4FlowPy.rst index 335186fae..c9f45570e 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -162,7 +162,8 @@ If the model extent (i.e. number of cells and/or rows in the input layers) is la will automatically split the input layers into different tiles (these are pickled to ``.npy`` files inside ``\temp`` folder). Each (quadratic) tile will then be consecutively calculated using all CPUs as defined by ``nCPU`` in ``avaframeCfg.ini``. The ``tileOverlap`` option defines by which margins the tiles overlap; in overlapping parts of the model domain the outputs -of the single tiles are combined (maximum, sum - depending on output variable). +of the single tiles are combined (maximum, sum - depending on output variable). Note that if a single tile does not cover +the entire process path, the path will be cut off, leading to incorrect results. The default settings provide reasonable performance on standard machines/model domains - however for special applications (e.g. modeling over large areas or on HPC hardware, **different raster resolution**) tweaking parameters might improve model performance. @@ -224,7 +225,7 @@ Output All outputs are written in *'.tif'* or in *'.asc'* raster format (controlable via the ``outputFileFormat`` option in ``(local_)com4FlowPyCfg.ini``, default is *'.tif'*) in the same resolution and extent as the input raster layers. You can customize which output rasters are written at the end of the model run by selecting the desired output files through the ``outputFiles`` option in ``(local_)com4FlowPyCfg.ini``. -In the output rasters, the cells, that are not affected by the process, are specified by the value -9999. +In the output rasters, the cells, that are not affected by the process, are specified by ``outputNoDataValue`` (default is -9999). By default the following four output layers are written to disk at the end of the model run: From 23b3e9c936e0ab7482142f9f004384ef25e803d4 Mon Sep 17 00:00:00 2001 From: paulasp Date: Tue, 18 Nov 2025 18:24:22 +0100 Subject: [PATCH 4/4] fix bug; written tif rasters are compressed per default fix bugs for tif compression --- avaframe/com4FlowPy/com4FlowPy.py | 97 +++++++++++++++++++++------ avaframe/com4FlowPy/com4FlowPyCfg.ini | 3 + avaframe/runCom4FlowPy.py | 20 +++--- 3 files changed, 91 insertions(+), 29 deletions(-) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index 3b360fcde..b448049ed 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -93,6 +93,7 @@ def com4FlowPyMain(cfgPath, cfgSetup): 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"]: @@ -545,37 +546,69 @@ 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: flux = defineNotAffectedCells(flux, cellCounts, noDataValue=_outputNoDataValue) - output = IOf.writeResultToRaster(outputHeader, flux, - modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), flip=True) + output = IOf.writeResultToRaster( + outputHeader, + flux, + modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) if 'zDelta' in _outputs: zDelta = defineNotAffectedCells(zDelta, cellCounts, noDataValue=_outputNoDataValue) - output = IOf.writeResultToRaster(outputHeader, zDelta, - modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), flip=True) + output = IOf.writeResultToRaster( + outputHeader, + zDelta, + modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) if 'cellCounts' in _outputs: cellCounts = defineNotAffectedCells(cellCounts, cellCounts, noDataValue=_outputNoDataValue) - output = IOf.writeResultToRaster(outputHeader, cellCounts, - modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), flip=True) + output = IOf.writeResultToRaster( + outputHeader, + cellCounts, + modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) if 'zDeltaSum' in _outputs: zDeltaSum = defineNotAffectedCells(zDeltaSum, cellCounts, noDataValue=_outputNoDataValue) - output = IOf.writeResultToRaster(outputHeader, zDeltaSum, - modelPaths["resDir"] / "com4_{}_{}_zDeltaSum".format(_uid, _ts), flip=True) + output = IOf.writeResultToRaster( + outputHeader, + zDeltaSum, + modelPaths["resDir"] / "com4_{}_{}_zDeltaSum".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) if 'routFluxSum' in _outputs: routFluxSum = defineNotAffectedCells(routFluxSum, cellCounts, noDataValue=_outputNoDataValue) - output = IOf.writeResultToRaster(outputHeader, routFluxSum, - modelPaths["resDir"] / "com4_{}_{}_routFluxSum".format(_uid, _ts), flip=True) + output = IOf.writeResultToRaster( + outputHeader, + routFluxSum, + modelPaths["resDir"] / "com4_{}_{}_routFluxSum".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) if 'depFluxSum' in _outputs: depFluxSum = defineNotAffectedCells(depFluxSum, cellCounts, noDataValue=_outputNoDataValue) - output = IOf.writeResultToRaster(outputHeader, depFluxSum, - modelPaths["resDir"] / "com4_{}_{}_depFluxSum".format(_uid, _ts), flip=True) + 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( @@ -583,6 +616,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): fpTaMax, modelPaths["resDir"] / "com4_{}_{}_fpTravelAngleMax".format(_uid, _ts), flip=True, + useCompression=useCompression, ) if "fpTravelAngleMin" in _outputs: fpTaMin = defineNotAffectedCells(fpTaMin, cellCounts, noDataValue=_outputNoDataValue) @@ -591,11 +625,17 @@ def mergeAndWriteResults(modelPaths, modelOptions): fpTaMin, modelPaths["resDir"] / "com4_{}_{}_fpTravelAngleMin".format(_uid, _ts), flip=True, + useCompression=useCompression, ) if 'slTravelAngle' in _outputs: slTa = defineNotAffectedCells(slTa, cellCounts, noDataValue=_outputNoDataValue) - output = IOf.writeResultToRaster(outputHeader, slTa, - modelPaths["resDir"] / "com4_{}_{}_slTravelAngle".format(_uid, _ts), flip=True) + 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( @@ -603,6 +643,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): travelLengthMax, modelPaths["resDir"] / "com4_{}_{}_travelLengthMax".format(_uid, _ts), flip=True, + useCompression=useCompression, ) if "travelLengthMin" in _outputs: travelLengthMin = defineNotAffectedCells(travelLengthMin, cellCounts, noDataValue=_outputNoDataValue) @@ -611,6 +652,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): travelLengthMin, modelPaths["resDir"] / "com4_{}_{}_travelLengthMin".format(_uid, _ts), flip=True, + useCompression=useCompression, ) # NOTE: @@ -619,14 +661,24 @@ def mergeAndWriteResults(modelPaths, modelOptions): # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("z_delta_sum%s" %(output_format)),z_delta_sum) if modelOptions["infraBool"]: # if infra backcalc = defineNotAffectedCells(backcalc, cellCounts, noDataValue=_outputNoDataValue) - output = IOf.writeResultToRaster(outputHeader, backcalc, - modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts), flip=True) + output = IOf.writeResultToRaster( + outputHeader, + backcalc, + modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) if modelOptions["forestInteraction"]: forestInteraction = defineNotAffectedCells( forestInteraction, cellCounts, noDataValue=_outputNoDataValue ) - output = IOf.writeResultToRaster(outputHeader, forestInteraction, - modelPaths["resDir"] / "com4_{}_{}_forestInteraction".format(_uid, _ts), flip=True) + output = IOf.writeResultToRaster( + outputHeader, + forestInteraction, + modelPaths["resDir"] / "com4_{}_{}_forestInteraction".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) del output @@ -664,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"] diff --git a/avaframe/com4FlowPy/com4FlowPyCfg.ini b/avaframe/com4FlowPy/com4FlowPyCfg.ini index 54baab0a9..f5be66d51 100644 --- a/avaframe/com4FlowPy/com4FlowPyCfg.ini +++ b/avaframe/com4FlowPy/com4FlowPyCfg.ini @@ -204,6 +204,9 @@ 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) diff --git a/avaframe/runCom4FlowPy.py b/avaframe/runCom4FlowPy.py index b709874a4..c8af4737a 100644 --- a/avaframe/runCom4FlowPy.py +++ b/avaframe/runCom4FlowPy.py @@ -111,7 +111,8 @@ def main(avalancheDir=''): cfgPath["uid"] = uid cfgPath["timeString"] = timeString cfgPath["outputFiles"] = cfgCustomPaths["outputFiles"] - cfgPath["outputNoDataValue"] = cfgCustomPaths["outputNoDataValue"] + cfgPath["outputNoDataValue"] = cfgCustomPaths.getfloat("outputNoDataValue") + cfgPath["useCompression"] = cfgCustomPaths.getboolean("useCompression") com4FlowPy.com4FlowPyMain(cfgPath, cfgSetup) @@ -171,6 +172,7 @@ def main(avalancheDir=''): 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"] @@ -213,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: @@ -239,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) @@ -252,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) @@ -265,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) @@ -279,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)