From ae290650063c085e6bf2e95ed9ec7b8ac0ee0a63 Mon Sep 17 00:00:00 2001 From: PaulaSp3 Date: Fri, 29 Aug 2025 13:18:09 +0200 Subject: [PATCH] compute startCellID per cell still in progress relVolMin and relVolMax as output correct position for relId delete print small change fix bug pathPolygons as output rename function tidy up try to keep RAM as small as possible print a warning if the output in cfg is not completly valid small corrections correct bugs and doc update docu adapt pytest delete unused package --- avaframe/com4FlowPy/com4FlowPy.py | 207 ++++++++++++++++++++------ avaframe/com4FlowPy/com4FlowPyCfg.ini | 7 + avaframe/com4FlowPy/flowClass.py | 34 ++++- avaframe/com4FlowPy/flowCore.py | 149 ++++++++++++++---- avaframe/com4FlowPy/splitAndMerge.py | 117 +++++++++++++++ avaframe/runCom4FlowPy.py | 45 +++++- avaframe/tests/test_com4FlowPy.py | 30 +++- docs/moduleCom4FlowPy.rst | 6 + 8 files changed, 499 insertions(+), 96 deletions(-) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index b448049ed..6170631ba 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -158,6 +158,18 @@ def com4FlowPyMain(cfgPath, cfgSetup): else: modelPaths["varExponentPath"] = "" + if "relIdPolygon" in modelPaths["outputFileList"] or "relIdCount" in modelPaths["outputFileList"]: + modelPaths["relIdPath"] = cfgPath["relIdPath"] + modelParameters["outputRelIdBool"] = True + else: + modelPaths["relIdPath"] = "" + modelParameters["outputRelIdBool"] = False + + if "relVolMin" in modelPaths["outputFileList"] or "relVolMax" in modelPaths["outputFileList"]: + modelParameters["outputRelVolBool"] = True + else: + modelParameters["outputRelVolBool"] = False + # TODO: provide some kind of check for the model Parameters # i.e. * sensible value ranges # * contradicting options ... @@ -452,6 +464,8 @@ def tileInputLayers(modelParameters, modelPaths, rasterAttributes, tilingParamet SPAM.tileRaster(modelPaths["varExponentPath"], "varExponent", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) if modelParameters["forestBool"]: SPAM.tileRaster(modelPaths["forestPath"], "forest", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) + if modelParameters["outputRelIdBool"]: + SPAM.tileRaster(modelPaths["relIdPath"], "relId", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) log.info("Finished Tiling All Input Rasters.") log.info("==================================") @@ -515,56 +529,40 @@ def mergeAndWriteResults(modelPaths, modelOptions): _outputs = set(modelPaths['outputFileList']) _outputNoDataValue = modelPaths["outputNoDataValue"] - log.info(" merging results ...") - log.info("-------------------------") - - # Merge calculated tiles - zDelta = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta") - flux = SPAM.mergeRaster(modelPaths["tempDir"], "res_flux") - cellCounts = SPAM.mergeRaster(modelPaths["tempDir"], "res_count", method="sum") - zDeltaSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta_sum", method="sum") - routFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_rout_flux_sum", method="sum") - depFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_dep_flux_sum", method="sum") - fpTaMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp_max") - fpTaMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp_min", method="min") - slTa = SPAM.mergeRaster(modelPaths["tempDir"], "res_sl") - travelLengthMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length_max") - travelLengthMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length_min", method="min") - - if modelOptions["infraBool"]: - backcalc = SPAM.mergeRaster(modelPaths["tempDir"], "res_backcalc") - - if modelOptions["forestInteraction"]: - forestInteraction = SPAM.mergeRaster(modelPaths["tempDir"], "res_forestInt", method='min') - - # Write Output Files to Disk - log.info("-------------------------") - log.info(" writing output files ...") - log.info("-------------------------") - _oF = modelPaths["outputFileFormat"] - _ts = modelPaths["timeString"] - demHeader = IOf.readRasterHeader(modelPaths["demPath"]) outputHeader = demHeader.copy() outputHeader["nodata_value"] = _outputNoDataValue + _oF = modelPaths["outputFileFormat"] + _ts = modelPaths["timeString"] + + useCompression = modelPaths["useCompression"] if _oF == ".asc": outputHeader["driver"] = "AAIGrid" elif _oF == ".tif": outputHeader["driver"] = "GTiff" - useCompression = modelPaths["useCompression"] + log.info(" merging and writing results ...") + log.info("-------------------------") - if 'flux' in _outputs: - flux = defineNotAffectedCells(flux, cellCounts, noDataValue=_outputNoDataValue) + # Merge calculated tiles + # compute cellCounts and don't delete because it is used for defining not affected cells + # other rasters (and polygons) are deleted after writing (to reduce computation time and used RAM) + cellCounts = SPAM.mergeRaster(modelPaths["tempDir"], "res_count", method="sum") + if 'cellCounts' in _outputs: + cellCounts = defineNotAffectedCells(cellCounts, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, - flux, - modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), + cellCounts, + modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), flip=True, useCompression=useCompression, - ) + ) + del output + log.info("com4_{}_{}_cellCounts is written".format(_uid, _ts)) + if 'zDelta' in _outputs: + zDelta = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta") zDelta = defineNotAffectedCells(zDelta, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -572,17 +570,27 @@ def mergeAndWriteResults(modelPaths, modelOptions): modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), flip=True, useCompression=useCompression, - ) - if 'cellCounts' in _outputs: - cellCounts = defineNotAffectedCells(cellCounts, cellCounts, noDataValue=_outputNoDataValue) + ) + del zDelta + del output + log.info("com4_{}_{}_zdelta is written".format(_uid, _ts)) + + if 'flux' in _outputs: + flux = SPAM.mergeRaster(modelPaths["tempDir"], "res_flux") + flux = defineNotAffectedCells(flux, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, - cellCounts, - modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), + flux, + modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), flip=True, useCompression=useCompression, - ) + ) + del flux + del output + log.info("com4_{}_{}_flux is written".format(_uid, _ts)) + if 'zDeltaSum' in _outputs: + zDeltaSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta_sum", method="sum") zDeltaSum = defineNotAffectedCells(zDeltaSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -591,7 +599,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del zDeltaSum + del output + log.info("com4_{}_{}_zDeltaSum is written".format(_uid, _ts)) + if 'routFluxSum' in _outputs: + routFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_rout_flux_sum", method="sum") routFluxSum = defineNotAffectedCells(routFluxSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -600,7 +613,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del routFluxSum + del output + log.info("com4_{}_{}_routFluxSum is written".format(_uid, _ts)) + if 'depFluxSum' in _outputs: + depFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_dep_flux_sum", method="sum") depFluxSum = defineNotAffectedCells(depFluxSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -609,7 +627,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del depFluxSum + del output + log.info("com4_{}_{}_depFluxSum is written".format(_uid, _ts)) + if "fpTravelAngle" in _outputs or "fpTravelAngleMax" in _outputs: + fpTaMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp_max") fpTaMax = defineNotAffectedCells(fpTaMax, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -618,7 +641,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del fpTaMax + del output + log.info("com4_{}_{}_fpTravelAngleMax is written".format(_uid, _ts)) + if "fpTravelAngleMin" in _outputs: + fpTaMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp_min", method="min") fpTaMin = defineNotAffectedCells(fpTaMin, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -627,7 +655,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del fpTaMin + del output + log.info("com4_{}_{}_fpTravelAngleMin is written".format(_uid, _ts)) + if 'slTravelAngle' in _outputs: + slTa = SPAM.mergeRaster(modelPaths["tempDir"], "res_sl") slTa = defineNotAffectedCells(slTa, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -636,7 +669,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del slTa + del output + log.info("com4_{}_{}_slTravelAngle is written".format(_uid, _ts)) + if "travelLength" in _outputs or "travelLengthMax" in _outputs: + travelLengthMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length_max") travelLengthMax = defineNotAffectedCells(travelLengthMax, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -645,7 +683,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del travelLengthMax + del output + log.info("com4_{}_{}_travelLengthMax is written".format(_uid, _ts)) + if "travelLengthMin" in _outputs: + travelLengthMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length_min", method="min") travelLengthMin = defineNotAffectedCells(travelLengthMin, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -654,12 +697,13 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del travelLengthMin + del output + log.info("com4_{}_{}_travelLengthMin is written".format(_uid, _ts)) - # 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 + + if modelOptions["infraBool"]: + backcalc = SPAM.mergeRaster(modelPaths["tempDir"], "res_backcalc") backcalc = defineNotAffectedCells(backcalc, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -667,8 +711,14 @@ def mergeAndWriteResults(modelPaths, modelOptions): modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts), flip=True, useCompression=useCompression, - ) + ) + del backcalc + del output + log.info("com4_{}_{}_backcalculation is written".format(_uid, _ts)) + + if modelOptions["forestInteraction"]: + forestInteraction = SPAM.mergeRaster(modelPaths["tempDir"], "res_forestInt", method='min') forestInteraction = defineNotAffectedCells( forestInteraction, cellCounts, noDataValue=_outputNoDataValue ) @@ -678,8 +728,69 @@ def mergeAndWriteResults(modelPaths, modelOptions): modelPaths["resDir"] / "com4_{}_{}_forestInteraction".format(_uid, _ts), flip=True, useCompression=useCompression, + ) + del forestInteraction + del output + log.info("com4_{}_{}_forestInteraction is written".format(_uid, _ts)) + + + if "relIdPolygon" in _outputs: + pathPolygons = SPAM.mergeDictToPolygon(modelPaths["tempDir"], "res_startCellIdDict", outputHeader) + pathPolygons.to_file( + modelPaths["resDir"] / "com4_{}_{}_pathPolygons.geojson".format(_uid, _ts), driver="GeoJSON" + ) + del pathPolygons + log.info("com4_{}_{}_pathPolygons is written".format(_uid, _ts)) + + + if "relIdCount" in _outputs: + countRelId = SPAM.mergeDictToRaster(modelPaths["tempDir"], "res_startCellIdDict") + countRelId = defineNotAffectedCells(countRelId, cellCounts, noDataValue=_outputNoDataValue) + output = IOf.writeResultToRaster( + outputHeader, + countRelId, + modelPaths["resDir"] / "com4_{}_{}_countRelId".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) + del countRelId + del output + log.info("com4_{}_{}_countRelId is written".format(_uid, _ts)) + + + if "relVolMin" in _outputs: + relVolMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_relVol_min", method="min") + relVolMin = defineNotAffectedCells(relVolMin, cellCounts, noDataValue=_outputNoDataValue) + output = IOf.writeResultToRaster( + outputHeader, + relVolMin, + modelPaths["resDir"] / "com4_{}_{}_relVolMin".format(_uid, _ts), + flip=True, + useCompression=useCompression, ) - del output + del relVolMin + del output + log.info("com4_{}_{}_relVolMin is written".format(_uid, _ts)) + + + if "relVolMax" in _outputs: + relVolMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_relVol_max") + relVolMax = defineNotAffectedCells(relVolMax, cellCounts, noDataValue=_outputNoDataValue) + output = IOf.writeResultToRaster( + outputHeader, + relVolMax, + modelPaths["resDir"] / "com4_{}_{}_relVolMax".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) + del relVolMax + del output + log.info("com4_{}_{}_relVolMax is written".format(_uid, _ts)) + + # 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) def checkConvertReleaseShp2Tif(modelPaths): diff --git a/avaframe/com4FlowPy/com4FlowPyCfg.ini b/avaframe/com4FlowPy/com4FlowPyCfg.ini index c4920e7d8..a56d56b39 100644 --- a/avaframe/com4FlowPy/com4FlowPyCfg.ini +++ b/avaframe/com4FlowPy/com4FlowPyCfg.ini @@ -227,8 +227,14 @@ outputNoDataValue = -9999 # depFluxSum # travelLengthMin # fpTravelAngleMin +# relVolMin +# relVolMax +# relIdPolygon +# relIdCount # if forestInteraction: forestInteraction is automatically added to outputs # if infra: backCalculation is automatically added to output +# if relVolMin or relVolMax is in outputFiles, the Volume of the PRA should be provided in the raster file in the REL folder +# if relIdCount or relIdPolygon is in outputFiles, the ids of the PRAs should be provided in the raster file in the RELID folder outputFiles = zDelta|cellCounts|travelLengthMax|fpTravelAngleMax #++++++++++++ Custom paths True/False @@ -256,6 +262,7 @@ useCustomPathDEM = False workDir = demPath = releasePath = +releaseIdPath = infraPath = forestPath = varUmaxPath = diff --git a/avaframe/com4FlowPy/flowClass.py b/avaframe/com4FlowPy/flowClass.py index 94cfa2dd6..c4e21e3e7 100644 --- a/avaframe/com4FlowPy/flowClass.py +++ b/avaframe/com4FlowPy/flowClass.py @@ -10,14 +10,25 @@ class Cell: + def __init__( self, - rowindex, colindex, - dem_ng, cellsize, - flux, z_delta, parent, - alpha, exp, flux_threshold, max_z_delta, - startcell, fluxDistOldVersionBool=False, - FSI=None, forestParams=None, + rowindex, + colindex, + dem_ng, + cellsize, + flux, + z_delta, + parent, + alpha, + exp, + flux_threshold, + max_z_delta, + startcell, + fluxDistOldVersionBool=False, + FSI=None, + forestParams=None, + startcellVol=None, ): """constructor for the Cell class that describes a raster cell that is hit by the GMF. the constructor function is called every time a new instance of type 'Cell' is @@ -62,6 +73,9 @@ def __init__( self._SQRT2 = np.sqrt(2.0) self._RAD90 = np.deg2rad(90.0) + self.startcellVolMin = startcellVol + self.startcellVolMax = startcellVol + # NOTE: Forest Interaction included here # if FSI != None AND forestParams != None - then self.ForestBool = True and forestParams and # FSI are accordingly initialized @@ -177,6 +191,10 @@ def add_parent(self, parent): if parent.forestIntCount < (self.forestIntCount - self.isForest): self.forestIntCount = parent.forestIntCount + self.isForest + def calc_startCellVol(self, startcellVolNew): + self.startcellVolMin = min(self.startcellVolMin, startcellVolNew) + self.startcellVolMax = max(self.startcellVolMax, startcellVolNew) + def calcDistMin(self, calc3D=False): """ function calculates the projected horizontal (self.min_distance) and 3D (self.minDistXYZ) length @@ -416,7 +434,7 @@ def calc_distribution(self): self.dist[self.dist < threshold] = 0 if np.sum(self.dist) != self.flux and count > 0: # correction/flux conservation for potential rounding losses or gains - # (self.flux - np.sum(self.dist)) will either be negative or positive + # (self.flux - np.sum(self.dist)) will either be negative or positive # depending on the direction of the rounding error self.dist[self.dist >= threshold] += (self.flux - np.sum(self.dist)) / count if count == 0: @@ -446,4 +464,4 @@ def forest_detrainment(self): # rise over run (should be negative slope) slope = (_rest - self.minAddedDetrainmentForest) / (0 - _noDetrainmentEffectZdelta) # y=mx+b, where zDelta is x - self.detrainment = max(self.minAddedDetrainmentForest, slope * self.z_delta + _rest) \ No newline at end of file + self.detrainment = max(self.minAddedDetrainmentForest, slope * self.z_delta + _rest) diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index bd2bc9c89..b004f9954 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -13,6 +13,7 @@ import gc import psutil import time +import pickle from multiprocessing import Pool @@ -143,6 +144,8 @@ def run(optTuple): varAlphaBool = optTuple[2]["varAlphaBool"] varExponentBool = optTuple[2]["varExponentBool"] fluxDistOldVersionBool = optTuple[2]["fluxDistOldVersionBool"] + relIdBool = optTuple[2]["outputRelIdBool"] + relVolBool = optTuple[2]["outputRelVolBool"] previewMode = optTuple[2]["previewMode"] # Temp-Dir (all input files are located here and results are written back in here) @@ -196,6 +199,16 @@ def run(optTuple): else: varExponentArray = None + if relIdBool: + relIdArray = np.load(tempDir / ("relId_%s_%s.npy" % (optTuple[0], optTuple[1]))) + else: + relIdArray = None + + if relVolBool: + relVolArray = release.copy() + else: + relVolArray = None + varParams = { "varUmaxBool": varUmaxBool, "varUmaxArray": varUmaxArray, @@ -204,6 +217,12 @@ def run(optTuple): "varExponentBool": varExponentBool, "varExponentArray": varExponentArray, } + relOutputParams = { + "relIdBool": relIdBool, + "relIdArray": relIdArray, + "relVolBool": relVolBool, + "relVolArray": relVolArray, + } # convert release areas to binary (0: no release areas, 1: release areas) # every positive value >0 is interpreted as release area @@ -250,6 +269,7 @@ def run(optTuple): forestArray, forestParams, outputs, + relOutputParams ] for release_sub in release_list ], @@ -275,6 +295,9 @@ def run(optTuple): travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 if forestInteraction: forestIntArray = np.ones_like(dem, dtype=np.float32) * -9999 + relVolMinArray = np.ones_like(dem, dtype=np.float32) * -9999 + relVolMaxArray = np.zeros_like(dem, dtype=np.float32) + processedStartCellIdDict = {} zDeltaList = [] fluxList = [] @@ -289,6 +312,9 @@ def run(optTuple): slTravelAngleList = [] travelLengthMaxList = [] travelLengthMinList = [] + processedStartCellIdList = [] + relVolMinList = [] + relVolMaxList = [] if forestInteraction: forestIntList = [] @@ -308,8 +334,11 @@ def run(optTuple): fpTravelAngleMinList.append(res[9]) routFluxSumList.append(res[10]) depFluxSumList.append(res[11]) + processedStartCellIdList.append(res[12]) if forestInteraction: - forestIntList.append(res[12]) + forestIntList.append(res[15]) + relVolMinList.append(res[13]) + relVolMaxList.append(res[14]) logging.info("Calculation finished, getting results.") for i in range(len(zDeltaList)): @@ -344,7 +373,26 @@ def run(optTuple): np.minimum(forestIntArray, forestIntList[i]), np.maximum(forestIntArray, forestIntList[i]), ) + if "relVolMin" in outputs: + relVolMinArray = np.where( + (relVolMinArray >= 0) & (relVolMinList[i] >= 0), + np.minimum(relVolMinArray, relVolMinList[i]), + np.maximum(relVolMinArray, relVolMinList[i]), + ) + if "relVolMax" in outputs: + relVolMaxArray = np.maximum(relVolMaxArray, relVolMaxList[i]) + if "relIdPolygon" in outputs or "relIdCount" in outputs: + for key in processedStartCellIdList[i]: + if key in processedStartCellIdDict: + ids = np.append(processedStartCellIdList[i][key], processedStartCellIdDict[key]) + processedStartCellIdDict[key] = np.unique(ids) + else: + processedStartCellIdDict[key] = processedStartCellIdList[i][key] + saveDict = open(tempDir / ("res_startCellIdDict_%s_%s.pickle" % (optTuple[0], optTuple[1])), "wb") + pickle.dump(processedStartCellIdDict, saveDict) + saveDict.close() + del processedStartCellIdDict # Save Calculated tiles 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) @@ -357,6 +405,8 @@ def run(optTuple): np.save(tempDir / ("res_sl_%s_%s" % (optTuple[0], optTuple[1])), slTravelAngleArray) np.save(tempDir / ("res_travel_length_max_%s_%s" % (optTuple[0], optTuple[1])), travelLengthMaxArray) np.save(tempDir / ("res_travel_length_min_%s_%s" % (optTuple[0], optTuple[1])), travelLengthMinArray) + np.save(tempDir / ("res_relVol_max_%s_%s" % (optTuple[0], optTuple[1])), relVolMaxArray) + np.save(tempDir / ("res_relVol_min_%s_%s" % (optTuple[0], optTuple[1])), relVolMinArray) if infraBool: np.save(tempDir / ("res_backcalc_%s_%s" % (optTuple[0], optTuple[1])), backcalc) if forestInteraction: @@ -390,6 +440,7 @@ def calculation(args): - args[14] (numpy array) - contains forest information (None if forestBool=False) - args[15] (dict) - contains parameters for forest interaction models (None if forestBool=False) - args[16] (list) - output names + - args[17] (dict) - contains flags and rasters for release - information outputs Returns ----------- @@ -457,6 +508,10 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): fluxDistOldVersionBool = args[12] previewMode = args[13] outputs = args[16] + relIdArray = args[17]["relIdArray"] + relVolArray = args[17]["relVolArray"] + relIdBool = args[17]["relIdBool"] + relVolBool = args[17]["relVolBool"] if forestBool: forestArray = args[14] @@ -482,6 +537,9 @@ 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 + relVolMinArray = np.ones_like(dem, dtype=np.float32) * -9999 + relVolMaxArray = np.zeros_like(dem, dtype=np.float32) + if infraBool: backcalc = np.ones_like(dem, dtype=np.int32) * -9999 else: @@ -494,12 +552,15 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): if forestInteraction: forestIntArray = np.ones_like(dem, dtype=np.float32) * -9999 + else: + forestIntArray = None # Core # NOTE-TODO: row_list, col_list are tuples - rethink variable naming row_list, col_list = get_start_idx(dem, release) startcell_idx = 0 + startCellIdDict = {} while startcell_idx < len(row_list): if infraBool: @@ -527,6 +588,15 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): startcell_idx += 1 continue + if relIdBool: + startcellId = relIdArray[row_idx, col_idx] + else: + startcellId = None + if relVolBool: + startcellVol = relVolArray[row_idx, col_idx] + else: + startcellVol = None + startcell = Cell( row_idx, col_idx, @@ -543,10 +613,12 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): fluxDistOldVersionBool=fluxDistOldVersionBool, FSI=forestArray[row_idx, col_idx] if isinstance(forestArray, np.ndarray) else None, forestParams=forestParams, + startcellVol=startcellVol, ) # dictionary of all the cells that have been processed and the number of times the cell has been visited processedCells[(startcell.rowindex, startcell.colindex)] = 1 + # list of flowClass.Cell() Objects that is contains the "path" for each release-cell cell_list.append(startcell) @@ -555,6 +627,14 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): updateInfraDirGraph(startcell.rowindex, startcell.colindex) for idx, cell in enumerate(cell_list): + if relIdBool: + if (cell.rowindex, cell.colindex) in startCellIdDict: + startcellIdList = np.append( + startCellIdDict[(cell.rowindex, cell.colindex)], startcellId) + startCellIdDict[(cell.rowindex, cell.colindex)] = np.unique( + startcellIdList) + else: + startCellIdDict[(cell.rowindex, cell.colindex)] = np.array([startcellId]) # calculate flux, z_delta from current cell (cell) to child-cells # lenght of row, col, flux, and z_delta vectors correspond to @@ -576,6 +656,8 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): if row[k] == cell_list[i].rowindex and col[k] == cell_list[i].colindex: cell_list[i].add_os(flux[k]) cell_list[i].add_parent(cell) + if relVolBool: + cell_list[i].calc_startCellVol(startcellVol) if infraBool: updateInfraDirGraph(row[k], col[k], cell.rowindex, cell.colindex) @@ -626,6 +708,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): fluxDistOldVersionBool=fluxDistOldVersionBool, FSI=forestArray[row[k], col[k]] if isinstance(forestArray, np.ndarray) else None, forestParams=forestParams, + startcellVol=startcellVol, ) ) @@ -667,6 +750,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): travelLengthMinArray[cell.rowindex, cell.colindex] = max( travelLengthMinArray[cell.rowindex, cell.colindex], cell.min_distance ) + if processedCells[(cell.rowindex, cell.colindex)] == 1: countArray[cell.rowindex, cell.colindex] += int(1) @@ -678,6 +762,19 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): else: forestIntArray[cell.rowindex, cell.colindex] = max( forestIntArray[cell.rowindex, cell.colindex], cell.forestIntCount + ) + if "relVolMax" in outputs: + relVolMaxArray[cell.rowindex, cell.colindex] = max( + relVolMaxArray[cell.rowindex, cell.colindex], cell.startcellVolMax + ) + if "relVolMin" in outputs: + if relVolMinArray[cell.rowindex, cell.colindex] >= 0 and cell.startcellVolMin >= 0: + relVolMinArray[cell.rowindex, cell.colindex] = min( + relVolMinArray[cell.rowindex, cell.colindex], cell.startcellVolMin + ) + else: + relVolMinArray[cell.rowindex, cell.colindex] = max( + relVolMinArray[cell.rowindex, cell.colindex], cell.startcellVolMin ) if infraBool: @@ -713,38 +810,24 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): zDeltaSumArray += zDeltaPathArray gc.collect() - - if forestInteraction: - return ( - zDeltaArray, - fluxArray, - countArray, - zDeltaSumArray, - backcalc, - fpTravelAngleMaxArray, - slTravelAngleArray, - travelLengthMaxArray, - travelLengthMinArray, - fpTravelAngleMinArray, - routFluxSumArray, - depFluxSumArray, - forestIntArray, - ) - else: - return ( - zDeltaArray, - fluxArray, - countArray, - zDeltaSumArray, - backcalc, - fpTravelAngleMaxArray, - slTravelAngleArray, - travelLengthMaxArray, - travelLengthMinArray, - fpTravelAngleMinArray, - routFluxSumArray, - depFluxSumArray, - ) + return ( + zDeltaArray, + fluxArray, + countArray, + zDeltaSumArray, + backcalc, + fpTravelAngleMaxArray, + slTravelAngleArray, + travelLengthMaxArray, + travelLengthMinArray, + fpTravelAngleMinArray, + routFluxSumArray, + depFluxSumArray, + startCellIdDict, + relVolMinArray, + relVolMaxArray, + forestIntArray, + ) def enoughMemoryAvailable(limit=0.05): diff --git a/avaframe/com4FlowPy/splitAndMerge.py b/avaframe/com4FlowPy/splitAndMerge.py index a4a82c246..7ba57bc6f 100644 --- a/avaframe/com4FlowPy/splitAndMerge.py +++ b/avaframe/com4FlowPy/splitAndMerge.py @@ -11,6 +11,9 @@ import gc import numpy as np import avaframe.in2Trans.rasterUtils as IOf +import shapely +import shapely.ops +import geopandas as gpd # create local logger log = logging.getLogger(__name__) @@ -233,3 +236,117 @@ def mergeRaster(inDirPath, fName, method='max'): del smallRas log.info("appended result %s_%i_%i", fName, i, j) return mergedRas + + +def mergeDict(inDirPath, fName): + """ + Merges the dictionary-results for each tile to one dictionary + + Parameters + ---------- + inDirPath: str + Path to the temporary files, that are results for each tile + fName : str + file name of the parameter which should be merged from tile-results + + Returns + ------- + mergedDict: dict + contains all + """ + nTiles = pickle.load(open(inDirPath / "nTiles", "rb")) + mergedDict = {} + + for i in range(nTiles[0] + 1): + for j in range(nTiles[1] + 1): + pos = pickle.load(open(inDirPath / ("ext_%i_%i" % (i, j)), "rb")) + with open(inDirPath / ("%s_%i_%i.pickle" % (fName, i, j)), "rb") as file: + smallDict = pickle.load(file) + if bool(smallDict): + for cellindSmall in smallDict: + cellind = (cellindSmall[0] + pos[0][0], cellindSmall[1] + pos[1][0]) + if cellind in mergedDict: + mergedDict[cellind] = np.append(smallDict[cellindSmall], mergedDict[cellind]) + else: + mergedDict[cellind] = smallDict[cellindSmall] + mergedDict[cellind] = np.unique(mergedDict[cellind]) + log.info("appended result %s_%i_%i", fName, i, j) + return mergedDict + + +def mergeDictToRaster(inDirPath, fName): + """ + Merges the dictionary-results for each tile to one array using + the length of the array assigned to each cell + + Parameters + ---------- + inDirPath: str + Path to the temporary files, that are results for each tile + fName : str + file name of the parameter which should be merged from tile-results + + Returns + ------- + mergedRas : numpy array + merged raster + """ + extL = pickle.load(open(inDirPath / "extentLarge", "rb")) + mergedRas = np.zeros((extL[0], extL[1])) + + mergedDict = mergeDict(inDirPath, fName) + for cellind in mergedDict: + mergedRas[cellind] = len(np.unique(mergedDict[cellind])) + del mergedDict + return mergedRas + + +def mergeDictToPolygon(inDirPath, fName, demHeader): + """ + Merges the dictionary-results for each tile to polygons for every path per PRA ID + + Parameters + ---------- + inDirPath: str + Path to the temporary files, that are results for each tile + fName : str + file name of the parameter which should be merged from tile-results + demHeader: dict + header of DEM raster + + Returns + ------- + gdfPathPolygons: GeoDataFrame + polygons per path for every PRA ID + """ + # get path polygons for every PRA ID + mergedDict = mergeDict(inDirPath, fName) + cellsize = demHeader["cellsize"] + pathPolygons = {} + + for (row, col), praIds in mergedDict.items(): + # get a polygon around every cell contained in mergedDict + xmin = col * cellsize + demHeader["xllcenter"] - cellsize / 2 + ymin = row * cellsize + demHeader["yllcenter"] - cellsize / 2 + xmax = xmin + cellsize + ymax = ymin + cellsize + cellsPoly = shapely.geometry.box(xmin, ymin, xmax, ymax) + + # reorder the dictionary: keys: PRA ID, values: list of polygons around every cell + for pid in praIds: + if pid not in pathPolygons: + pathPolygons[pid] = [] + pathPolygons[pid].append(cellsPoly) + del mergedDict + + for pid, polys in pathPolygons.items(): + # merge all polygons belonging to a PRA ID + pathPolygons[pid] = shapely.ops.unary_union(polys) + + gdfPathPolygons = gpd.GeoDataFrame( + {"PRA_id": list(pathPolygons.keys())}, + geometry=list(pathPolygons.values()), + crs=demHeader["crs"], + ) + del pathPolygons + return gdfPathPolygons diff --git a/avaframe/runCom4FlowPy.py b/avaframe/runCom4FlowPy.py index c8af4737a..aab8ca789 100644 --- a/avaframe/runCom4FlowPy.py +++ b/avaframe/runCom4FlowPy.py @@ -163,6 +163,7 @@ def main(avalancheDir=''): cfgPath["tempDir"] = pathlib.Path(temp_dir) cfgPath["demPath"] = pathlib.Path(cfgCustomPaths["demPath"]) cfgPath["releasePath"] = pathlib.Path(cfgCustomPaths["releasePath"]) + cfgPath["relIdPath"] = pathlib.Path(cfgCustomPaths["relIdPath"]) cfgPath["infraPath"] = pathlib.Path(cfgCustomPaths["infraPath"]) cfgPath["forestPath"] = pathlib.Path(cfgCustomPaths["forestPath"]) cfgPath["varUmaxPath"] = pathlib.Path(cfgCustomPaths["varUmaxPath"]) @@ -229,7 +230,7 @@ def readFlowPyinputs(avalancheDir, cfgFlowPy, log): # TODO: also use the getAndCheckInputFiles to get the paths for the following files? # read infra area if cfgFlowPy.getboolean("GENERAL", "infra") is True: - infraPath, available = gI.getAndCheckInputFiles(inputDir, "INFRA", "Infra", fileExt="raster") + infraPath, available, _ = gI.getAndCheckInputFiles(inputDir, "INFRA", "Infra", fileExt="raster") if available == "No": message = f"There is no infra file in supported format provided in {avalancheDir}/INFRA" log.error(message) @@ -291,6 +292,20 @@ def readFlowPyinputs(avalancheDir, cfgFlowPy, log): forestPath = "" cfgPath["forestPath"] = forestPath + # read release ID raster + if "relIdPolygon" in cfgFlowPy["PATHS"]["outputFiles"].split("|") or "relIdCount" in cfgFlowPy["PATHS"][ + "outputFiles" + ].split("|"): + relIdPath, available, _ = gI.getAndCheckInputFiles(inputDir, "RELID", "release ID", fileExt="raster") + if available == "No": + message = f"There is no release id file in supported format provided in {avalancheDir}/RELID" + log.error(message) + raise AssertionError(message) + log.info("Release ID file is: %s" % relIdPath) + else: + relIdPath = "" + cfgPath["relIdPath"] = relIdPath + # read DEM if cfgFlowPy.getboolean("PATHS", "useCustomPathDEM") is False: demPath = gI.getDEMPath(avalancheDir) @@ -325,11 +340,33 @@ def checkOutputFilesFormat(strOutputFiles): try: setA = set(strOutputFiles.split('|')) - setB = set(['zDelta', 'cellCounts', 'fpTravelAngle', 'travelLength', 'fpTravelAngleMax', - 'fpTravelAngleMin', 'travelLengthMax', 'travelLengthMin', - 'slTravelAngle', 'flux', 'zDeltaSum', 'routFluxSum', 'depFluxSum']) + setB = set( + [ + "zDelta", + "cellCounts", + "fpTravelAngle", + "travelLength", + "fpTravelAngleMax", + "fpTravelAngleMin", + "travelLengthMax", + "travelLengthMin", + "slTravelAngle", + "flux", + "zDeltaSum", + "routFluxSum", + "depFluxSum", + "relVolMin", + "relVolMax", + "relIdPolygon", + "relIdCount", + ] + ) # if there is at least 1 correct outputfile defined, we use the string provided in the .ini file if (setA & setB): + outNotValid = setA - setB + if outNotValid: + for outFile in outNotValid: + print("WARNING! - {} is not a valid output file and is not computed".format(outFile)) return strOutputFiles else: raise ValueError('outputFiles defined in .ini have wrong format - using default settings') diff --git a/avaframe/tests/test_com4FlowPy.py b/avaframe/tests/test_com4FlowPy.py index 8c4321149..9029fee0b 100644 --- a/avaframe/tests/test_com4FlowPy.py +++ b/avaframe/tests/test_com4FlowPy.py @@ -197,8 +197,32 @@ def test_calculation(): outputs = ['travelLengthMin', 'flux'] forestArray = None forestParams = None - args = [dem, infra, pra, alpha, exp, fluxTh, zDeltaMax, nodata, cellsize, infraBool, forestBool, variableParameters, - fluxDistOldVersionBool, previewMode, forestArray, forestParams, outputs] + relOutputParams = { + "relIdBool": False, + "relIdArray": None, + "relVolBool": False, + "relVolArray": None, + } + args = [ + dem, + infra, + pra, + alpha, + exp, + fluxTh, + zDeltaMax, + nodata, + cellsize, + infraBool, + forestBool, + variableParameters, + fluxDistOldVersionBool, + previewMode, + forestArray, + forestParams, + outputs, + relOutputParams, + ] flux = ones_like(dem) * -9999. flux[[1, 2, 3], [2, 2, 2]] = 1 @@ -210,7 +234,7 @@ def test_calculation(): travelLengthMin[3, 2] = 2 * np.sqrt(cellsize ** 2) results = flowCore.calculation(args) - assert len(results) == 12 + assert len(results) == 16 assert np.all(results[1] == flux) assert np.all(results[10] == routFluxSum) assert np.all(results[11] == depFluxSum) diff --git a/docs/moduleCom4FlowPy.rst b/docs/moduleCom4FlowPy.rst index f5a97ae04..e64b4c105 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -192,6 +192,7 @@ the following folder structure inside the ``avalancheDir`` directory inside whic CFGs/ - expert configuration files (optional) ElevationModel - digital elevation model (.asc) REL/ - release area file (can be either .asc, .tif, or .shp) + RELID/ - release ID file (can be either .asc, .tif) RES/ - forest structure information (FSI) (.asc or .tif) INFRA/ - infrastructure layer (.asc or .tif) ALPHA/ - variable alpha angle layer (.tif) @@ -207,6 +208,7 @@ inputs and working directories/model outputs in different places, which might be - ``workDir`` :math:`\ldots` working directory (a ``temp/`` folder, model log and model results will be written here) - ``demPath`` :math:`\ldots` path to input DEM (must be ``.asc`` currently) - ``releasePath`` :math:`\ldots` path to release area raster (``.asc, .tif``) +- ``releaseIdPath`` :math:`\ldots` path to release ID raster (``.asc, .tif``) - ``infraPath`` :math:`\ldots` path to infrastructure raster (``.asc, .tif``) (required if ``infra = True``) - ``forestPath`` :math:`\ldots` path to forest (FSI) raster (``.asc, .tif``) (required if ``forest = True``) - ``varAlphaPath`` :math:`\ldots` path to variable alpha angle raster (``.tif``) (required if ``variableAlpha = True``) @@ -244,6 +246,10 @@ In addition these output layers are also available: - ``depFluxSum``: deposited flux summed up over all paths - ``travelLengthMin``: the travel length along the flow path (the minimum value of all paths for every raster cell) - ``fpTravelAngleMin``: the gamma angle along the flow path (the minimum value of all paths for every raster cell) +- ``relVolMin``: the minimum of the tracked release volume that route through the raster cell (the volume is provided in the release file in the REL folder) +- ``relVolMax``: the maximum of the tracked release volume that route through the raster cell (the volume is provided in the release file in the REL folder) +- ``relIdPolygon``: polygons (*.geojson) that cover the affected process belonging to one release ID (can include multiple release cells; release IDs are provided in the RELID folder) +- ``relIdCount``: number of paths belonging to different release IDs that route flux through a raster cell (release IDs are provided in the RELID folder) If ``forestInteraction = True`` this layer will be written automatically (no need to separately define in ``outputFiles``):