From 33482a16087dbb70907aef86986b4c27992b8aff Mon Sep 17 00:00:00 2001 From: PaulaSp3 Date: Fri, 29 Aug 2025 13:18:09 +0200 Subject: [PATCH 1/3] 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 --- avaframe/com4FlowPy/com4FlowPy.py | 207 ++++++++++++++++++++------ avaframe/com4FlowPy/com4FlowPyCfg.ini | 7 + avaframe/com4FlowPy/flowClass.py | 34 ++++- avaframe/com4FlowPy/flowCore.py | 196 ++++++++++++++++++------ avaframe/com4FlowPy/splitAndMerge.py | 118 +++++++++++++++ avaframe/runCom4FlowPy.py | 45 +++++- docs/moduleCom4FlowPy.rst | 6 + 7 files changed, 505 insertions(+), 108 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 a9a2315f7..822d9bc23 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -13,6 +13,7 @@ import gc import psutil import time +import pickle if os.name == "nt": from multiprocessing.pool import Pool as Pool @@ -146,6 +147,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) @@ -199,6 +202,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, @@ -207,6 +220,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)): @@ -342,7 +371,26 @@ def run(optTuple): forestIntArray = np.where((forestIntArray >= 0) & (forestIntList[i] >= 0), 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) @@ -355,6 +403,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: @@ -388,6 +438,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 ----------- @@ -455,6 +506,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] @@ -480,6 +535,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: @@ -492,12 +550,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: @@ -525,18 +586,37 @@ 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, - dem_ng, cellsize, - 1, 0, None, - alpha, exp, flux_threshold, max_z_delta, - startcell=True, fluxDistOldVersionBool=fluxDistOldVersionBool, + row_idx, + col_idx, + dem_ng, + cellsize, + 1, + 0, + None, + alpha, + exp, + flux_threshold, + max_z_delta, + startcell=True, + 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) @@ -545,6 +625,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 @@ -566,6 +654,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) @@ -599,16 +689,26 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): else: processedCells[(row[k], col[k])] = 1 - cell_list.append(Cell( - row[k], col[k], - dem_ng, cellsize, - flux[k], z_delta[k], - cell, - alpha, exp, flux_threshold, max_z_delta, - startcell, fluxDistOldVersionBool=fluxDistOldVersionBool, - FSI=forestArray[row[k], col[k]] if isinstance(forestArray, np.ndarray) else None, - forestParams=forestParams, - )) + cell_list.append( + Cell( + row[k], + col[k], + dem_ng, + cellsize, + flux[k], + z_delta[k], + cell, + alpha, + exp, + flux_threshold, + max_z_delta, + startcell, + fluxDistOldVersionBool=fluxDistOldVersionBool, + FSI=forestArray[row[k], col[k]] if isinstance(forestArray, np.ndarray) else None, + forestParams=forestParams, + startcellVol=startcellVol, + ) + ) 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) @@ -643,6 +743,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) @@ -653,6 +754,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: # if 'infraBool' is True - i.e. calculation is performed with infrastructure information @@ -683,38 +797,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..245f7fefb 100644 --- a/avaframe/com4FlowPy/splitAndMerge.py +++ b/avaframe/com4FlowPy/splitAndMerge.py @@ -11,6 +11,10 @@ import gc import numpy as np import avaframe.in2Trans.rasterUtils as IOf +import shapely +import shapely.ops +import geopandas as gpd +from shapely.geometry import Polygon # create local logger log = logging.getLogger(__name__) @@ -233,3 +237,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/docs/moduleCom4FlowPy.rst b/docs/moduleCom4FlowPy.rst index c9f45570e..52a35bc49 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -191,6 +191,7 @@ the following folder structure inside the ``avalancheDir`` directory inside whic Inputs/ 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) @@ -206,6 +207,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``) @@ -243,6 +245,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``): From 120cc0e78b81281550e1bdd6dcda3ed8e7fab630 Mon Sep 17 00:00:00 2001 From: paulasp Date: Tue, 25 Nov 2025 11:25:10 +0100 Subject: [PATCH 2/3] adapt pytest --- avaframe/tests/test_com4FlowPy.py | 30 +++++++++++++++++++++++++++--- 1 file changed, 27 insertions(+), 3 deletions(-) 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) From e75f10d3eb93a98f58cec833aa088976b9cdca16 Mon Sep 17 00:00:00 2001 From: paulasp Date: Tue, 25 Nov 2025 11:35:59 +0100 Subject: [PATCH 3/3] delete unused package --- avaframe/com4FlowPy/splitAndMerge.py | 1 - 1 file changed, 1 deletion(-) diff --git a/avaframe/com4FlowPy/splitAndMerge.py b/avaframe/com4FlowPy/splitAndMerge.py index 245f7fefb..7ba57bc6f 100644 --- a/avaframe/com4FlowPy/splitAndMerge.py +++ b/avaframe/com4FlowPy/splitAndMerge.py @@ -14,7 +14,6 @@ import shapely import shapely.ops import geopandas as gpd -from shapely.geometry import Polygon # create local logger log = logging.getLogger(__name__)