diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index 146eaa63d..b448049ed 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/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,61 +546,113 @@ 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: @@ -605,11 +660,25 @@ 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 - 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 diff --git a/avaframe/com4FlowPy/com4FlowPyCfg.ini b/avaframe/com4FlowPy/com4FlowPyCfg.ini index 19e17d887..f5be66d51 100644 --- a/avaframe/com4FlowPy/com4FlowPyCfg.ini +++ b/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: diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index 731195178..a9a2315f7 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/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]), ) diff --git a/avaframe/runCom4FlowPy.py b/avaframe/runCom4FlowPy.py index 3fbde5a2d..c8af4737a 100644 --- a/avaframe/runCom4FlowPy.py +++ b/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() \ No newline at end of file + main() 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 28514dcff..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,6 +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 ``outputNoDataValue`` (default is -9999). By default the following four output layers are written to disk at the end of the model run: