From a6e6c390d9466f2f113aa32f92d055696dcfc64a Mon Sep 17 00:00:00 2001 From: PaulaSp3 Date: Fri, 29 Aug 2025 13:18:09 +0200 Subject: [PATCH 1/6] 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 dfc8e18a0..ed1b0e745 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -191,6 +191,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) @@ -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 04bbefc61219aa684c89039bb3955969fbc3b956 Mon Sep 17 00:00:00 2001 From: ahuber-bfw <112931921+ahuber-bfw@users.noreply.github.com> Date: Fri, 8 May 2026 13:10:36 +0200 Subject: [PATCH 2/6] Update com4FlowPy.py deleteTempFolder - small BugFix and adaptation for releaseID functionality --- avaframe/com4FlowPy/com4FlowPy.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index 6170631ba..935bbc7d1 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -841,11 +841,11 @@ def checkConvertReleaseShp2Tif(modelPaths): def deleteTempFolder(tempFolderPath): - """delete tempFolder containing the pickled np.arrays of the input data and output data tiles. + """delete tempFolder containing the pickled np.arrays of the input and output data tiles or other .pickle files. should be called after all merged model results are written to disk. performs a few checks to make sure the folder is indeed a com4FlowPy tempFolder, i.e. - does not contain subfolders - - no other file-extensions than '.npy' and '' + - no other file-extensions than '.npy', '.pickle' and '' Parameters ----------- @@ -866,9 +866,9 @@ def deleteTempFolder(tempFolderPath): if f.is_dir(): validTemp = False break - # check if all files are either .npy or start with "ext, "nTi" + # check if all files are either .npy, .pickle or start with "ext, "nTi" elif f.is_file(): - if not f.path.endswith(".npy"): + if not (f.path.endswith(".npy") or f.path.endswith(".pickle")): if not f.name[:3] in ["ext", "nTi"]: validTemp = False break @@ -883,7 +883,7 @@ def deleteTempFolder(tempFolderPath): print("Error: %s : %s" % (tempFolderPath, e.strerror)) else: log.info("deletion of temp folder {} failed".format(tempFolderPath)) - log.info(" isDir:{} isTemp:{}}".format(isDir, validTemp)) + log.info(" isDir:{} isTemp:{}".format(isDir, validTemp)) log.info("+++++++++++++++++++++++") From 42fbd86d1274bfb99a70b942401c57306bdc954d Mon Sep 17 00:00:00 2001 From: Paula Date: Wed, 13 May 2026 17:46:46 +0200 Subject: [PATCH 3/6] add tests for tiling, merging of tiles fixed bug extend split and merge test add test for dictToPolygon fix bug format and test fix --- avaframe/com4FlowPy/flowCore.py | 54 +- avaframe/com4FlowPy/splitAndMerge.py | 39 +- avaframe/in2Trans/rasterUtils.py | 2 +- .../tests/data/testCom4/refPolygon.geojson | 10 + .../testCom4/refPolygon_manipulated.geojson | 10 + avaframe/tests/data/testCom4/testRaster.tif | Bin 0 -> 1184 bytes avaframe/tests/test_com4FlowPy.py | 541 +++++++++++++----- 7 files changed, 484 insertions(+), 172 deletions(-) create mode 100644 avaframe/tests/data/testCom4/refPolygon.geojson create mode 100644 avaframe/tests/data/testCom4/refPolygon_manipulated.geojson create mode 100644 avaframe/tests/data/testCom4/testRaster.tif diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index b004f9954..318fd00f3 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -269,7 +269,7 @@ def run(optTuple): forestArray, forestParams, outputs, - relOutputParams + relOutputParams, ] for release_sub in release_list ], @@ -521,6 +521,10 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): forestInteraction = False forestArray = None forestParams = None + if infraBool: + # initialize infrastructure array + # TODO: check if this can be simplified + infraArr = infra # infrastructure array (input file) zDeltaArray = np.zeros_like(dem, dtype=np.float32) zDeltaSumArray = np.zeros_like(dem, dtype=np.float32) @@ -529,27 +533,43 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): depFluxSumArray = np.zeros_like(dem, dtype=np.float32) fluxArray = np.ones_like(dem, dtype=np.float32) * -9999 countArray = np.zeros_like(dem, dtype=np.int32) - - fpTravelAngleMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 # fp = Flow Path - fpTravelAngleMinArray = np.ones_like(dem, dtype=np.float32) * -9999 # fp = Flow Path slTravelAngleArray = np.ones_like(dem, dtype=np.float32) * -9999 # sl = Straight Line - travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 - travelLengthMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 + if "fpTravelAngleMax" in outputs or "fpTravelAngle" in outputs: + fpTravelAngleMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 # fp = Flow Path + else: + fpTravelAngleMaxArray = None - relVolMinArray = np.ones_like(dem, dtype=np.float32) * -9999 - relVolMaxArray = np.zeros_like(dem, dtype=np.float32) + if "fpTravelAngleMin" in outputs: + fpTravelAngleMinArray = np.ones_like(dem, dtype=np.float32) * -9999 # fp = Flow Path + else: + fpTravelAngleMinArray = None + + if "travelLengthMin" in outputs: + travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 + else: + travelLengthMinArray = None + + if "travelLengthMax" in outputs or "travelLength" in outputs: + travelLengthMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 + else: + travelLengthMaxArray = None + + if "relVolMin" in outputs: + relVolMinArray = np.ones_like(dem, dtype=np.float32) * -9999 + else: + relVolMinArray = None + + if "relVolMax" in outputs: + relVolMaxArray = np.zeros_like(dem, dtype=np.float32) + else: + relVolMaxArray = None if infraBool: backcalc = np.ones_like(dem, dtype=np.int32) * -9999 else: backcalc = None - if infraBool: - # initialize infrastructure array - # TODO: check if this can be simplified - infraArr = infra # infrastructure array (input file) - if forestInteraction: forestIntArray = np.ones_like(dem, dtype=np.float32) * -9999 else: @@ -629,10 +649,8 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): 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) + 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]) @@ -762,7 +780,7 @@ 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 diff --git a/avaframe/com4FlowPy/splitAndMerge.py b/avaframe/com4FlowPy/splitAndMerge.py index 7ba57bc6f..5a65100a5 100644 --- a/avaframe/com4FlowPy/splitAndMerge.py +++ b/avaframe/com4FlowPy/splitAndMerge.py @@ -2,8 +2,8 @@ # -*- coding: utf-8 -*- """ - Functions to handle the raster tiles. - Tiling is intended to manage processing of large(r) computational domains. +Functions to handle the raster tiles. + Tiling is intended to manage processing of large(r) computational domains. """ import logging @@ -25,11 +25,11 @@ def tileRaster(fNameIn, fNameOut, dirName, xDim, yDim, U, isInit=False): Parameters ----------- - fNameIn : str + fNameIn : pathlib path path to raster that is tiled fNameOut: str name of saved raster file - dirName: str + dirName: pathlib path path to folder, where tiled raster is saved (temp - folder) xDim: int size of one tile in x dimension (number of raster columns) @@ -179,14 +179,14 @@ def tileRaster(fNameIn, fNameOut, dirName, xDim, yDim, U, isInit=False): # return largeRaster -def mergeRaster(inDirPath, fName, method='max'): +def mergeRaster(inDirPath, fName, method="max"): """ Merges the results for each tile to one array using the method provided through the function parameters Parameters ---------- - inDirPath: str + inDirPath: pathlib path 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 @@ -210,7 +210,7 @@ def mergeRaster(inDirPath, fName, method='max'): mergedRas = np.zeros((extL[0], extL[1])) # create Raster with original size - if method != 'sum': + if method != "sum": mergedRas[:, :] = np.nan for i in range(nTiles[0] + 1): @@ -220,18 +220,19 @@ def mergeRaster(inDirPath, fName, method='max'): pos = pickle.load(open(inDirPath / ("ext_%i_%i" % (i, j)), "rb")) # print pos - if method == 'max': - mergedRas[pos[0][0]: pos[0][1], pos[1][0]: pos[1][1]] = np.fmax( - mergedRas[pos[0][0]: pos[0][1], pos[1][0]: pos[1][1]], smallRas + if method == "max": + mergedRas[pos[0][0] : pos[0][1], pos[1][0] : pos[1][1]] = np.fmax( + mergedRas[pos[0][0] : pos[0][1], pos[1][0] : pos[1][1]], smallRas + ) + elif method == "min": + mergedRas[pos[0][0] : pos[0][1], pos[1][0] : pos[1][1]] = np.where( + (mergedRas[pos[0][0] : pos[0][1], pos[1][0] : pos[1][1]] >= 0) & (smallRas >= 0), + np.fmin(mergedRas[pos[0][0] : pos[0][1], pos[1][0] : pos[1][1]], smallRas), + np.fmax(mergedRas[pos[0][0] : pos[0][1], pos[1][0] : pos[1][1]], smallRas), ) - elif method == 'min': - mergedRas[pos[0][0]:pos[0][1], pos[1][0]:pos[1][1]] =\ - np.where((mergedRas[pos[0][0]:pos[0][1], pos[1][0]:pos[1][1]] >= 0) & (smallRas >= 0), - np.fmin(mergedRas[pos[0][0]:pos[0][1], pos[1][0]:pos[1][1]], smallRas), - np.fmax(mergedRas[pos[0][0]:pos[0][1], pos[1][0]:pos[1][1]], smallRas)) - if method == 'sum': - mergedRas[pos[0][0]: pos[0][1], pos[1][0]: pos[1][1]] = np.add( - mergedRas[pos[0][0]: pos[0][1], pos[1][0]: pos[1][1]], smallRas + if method == "sum": + mergedRas[pos[0][0] : pos[0][1], pos[1][0] : pos[1][1]] = np.add( + mergedRas[pos[0][0] : pos[0][1], pos[1][0] : pos[1][1]], smallRas ) del smallRas log.info("appended result %s_%i_%i", fName, i, j) @@ -244,7 +245,7 @@ def mergeDict(inDirPath, fName): Parameters ---------- - inDirPath: str + inDirPath: pathlib path 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 diff --git a/avaframe/in2Trans/rasterUtils.py b/avaframe/in2Trans/rasterUtils.py index a66d1b780..5e1f188d7 100644 --- a/avaframe/in2Trans/rasterUtils.py +++ b/avaframe/in2Trans/rasterUtils.py @@ -160,7 +160,7 @@ class with methods that give cellsize, nrows, ncols, xllcenter yllcenter, nodata_value, driver, transfrom, crs resultArray : 2D numpy array 2D numpy array of values that shall be written to file - outFileName : str + outFileName : pathlib path path incl. name of file to be written flip: boolean if True, flip the rows of the resultArray when writing. AF considers the first line in a data array to be the diff --git a/avaframe/tests/data/testCom4/refPolygon.geojson b/avaframe/tests/data/testCom4/refPolygon.geojson new file mode 100644 index 000000000..bbfab2e48 --- /dev/null +++ b/avaframe/tests/data/testCom4/refPolygon.geojson @@ -0,0 +1,10 @@ +{ +"type": "FeatureCollection", +"name": "test", +"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }, +"features": [ +{ "type": "Feature", "properties": { "PRA_id": 3.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 5.0, -5.0 ], [ 5.0, 5.0 ], [ -5.0, 5.0 ], [ -5.0, -5.0 ], [ 5.0, -5.0 ] ] ] } }, +{ "type": "Feature", "properties": { "PRA_id": 1.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 15.0, 35.0 ], [ 15.0, 45.0 ], [ 25.0, 45.0 ], [ 35.0, 45.0 ], [ 45.0, 45.0 ], [ 45.0, 35.0 ], [ 45.0, 25.0 ], [ 45.0, 15.0 ], [ 35.0, 15.0 ], [ 25.0, 15.0 ], [ 15.0, 15.0 ], [ 15.0, 25.0 ], [ 15.0, 35.0 ] ] ] } }, +{ "type": "Feature", "properties": { "PRA_id": 2.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 45.0, 65.0 ], [ 45.0, 75.0 ], [ 45.0, 85.0 ], [ 55.0, 85.0 ], [ 65.0, 85.0 ], [ 75.0, 85.0 ], [ 85.0, 85.0 ], [ 85.0, 75.0 ], [ 85.0, 65.0 ], [ 85.0, 55.0 ], [ 75.0, 55.0 ], [ 75.0, 45.0 ], [ 65.0, 45.0 ], [ 55.0, 45.0 ], [ 45.0, 45.0 ], [ 45.0, 55.0 ], [ 45.0, 65.0 ] ] ] } } +] +} diff --git a/avaframe/tests/data/testCom4/refPolygon_manipulated.geojson b/avaframe/tests/data/testCom4/refPolygon_manipulated.geojson new file mode 100644 index 000000000..caca8e148 --- /dev/null +++ b/avaframe/tests/data/testCom4/refPolygon_manipulated.geojson @@ -0,0 +1,10 @@ +{ +"type": "FeatureCollection", +"name": "test", +"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }, +"features": [ +{ "type": "Feature", "properties": { "PRA_id": 3.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 5.0, -5.0 ], [ 5.0, 5.0 ], [ -5.0, 5.0 ], [ -5.0, -5.0 ], [ 5.0, -5.0 ] ] ] } }, +{ "type": "Feature", "properties": { "PRA_id": 1.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 15.0, 35.0 ], [ 15.0, 45.0 ], [ 25.0, 45.0 ], [ 25.0, 55.0 ], [ 35.0, 55.0 ], [ 45.0, 55.0 ], [ 45.0, 45.0 ], [ 45.0, 35.0 ], [ 45.0, 25.0 ], [ 45.0, 15.0 ], [ 35.0, 15.0 ], [ 25.0, 15.0 ], [ 15.0, 15.0 ], [ 15.0, 25.0 ], [ 15.0, 35.0 ] ] ] } }, +{ "type": "Feature", "properties": { "PRA_id": 2.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 45.0, 65.0 ], [ 45.0, 75.0 ], [ 45.0, 85.0 ], [ 55.0, 85.0 ], [ 65.0, 85.0 ], [ 75.0, 85.0 ], [ 85.0, 85.0 ], [ 85.0, 75.0 ], [ 85.0, 65.0 ], [ 85.0, 55.0 ], [ 75.0, 55.0 ], [ 75.0, 45.0 ], [ 65.0, 45.0 ], [ 55.0, 45.0 ], [ 45.0, 45.0 ], [ 35.0, 45.0 ], [ 25.0, 45.0 ], [ 25.0, 55.0 ], [ 35.0, 55.0 ], [ 45.0, 55.0 ], [ 45.0, 65.0 ] ] ] } } +] +} diff --git a/avaframe/tests/data/testCom4/testRaster.tif b/avaframe/tests/data/testCom4/testRaster.tif new file mode 100644 index 0000000000000000000000000000000000000000..a47562f837d4cf1594f821813e8dde9087a0a73d GIT binary patch literal 1184 zcmebD)MDUZU|-pT0L7W1Y>+xOB(@+U3s`RhP(l<*Tnx$v znJErcqrl9-AcLeP7|I3;Gw?O@Fo5V=K-|>A!@vflKLPRjc4h_zAgu&6e?vPD0~3&) z1Y~b)X9Byv2*_?+!UR@d1!N7n8c)^!bDTJ%z=p$Z4gu;Odl?qTILZj=Yu^&0Xj`U8%}9- Zb6|X&8VE?h^rO=Rw9!nO!+}~x0|0%hERO&H literal 0 HcmV?d00001 diff --git a/avaframe/tests/test_com4FlowPy.py b/avaframe/tests/test_com4FlowPy.py index 9029fee0b..f4fbf5be9 100644 --- a/avaframe/tests/test_com4FlowPy.py +++ b/avaframe/tests/test_com4FlowPy.py @@ -1,45 +1,49 @@ """ - Pytest for module com4FlowPy +Pytest for module com4FlowPy """ # Load modules import numpy as np +import pathlib import pytest -from numpy.ma.core import ones_like - +import pickle +import os +import rasterio +import geopandas as gpd from avaframe.com4FlowPy import flowClass import avaframe.com4FlowPy.flowCore as flowCore +import avaframe.com4FlowPy.splitAndMerge as SPAM +import avaframe.in2Trans.rasterUtils as IOf + def test_add_os(): - cell = flowClass.Cell(1,1, - np.array([[10,10,10], [10,10,10], [10,10,10]]), 10, - 1,0, None, - 20, 8, 3e-4, 270, - startcell=True) + cell = flowClass.Cell( + 1, + 1, + np.array([[10, 10, 10], [10, 10, 10], [10, 10, 10]]), + 10, + 1, + 0, + None, + 20, + 8, + 3e-4, + 270, + startcell=True, + ) cell.add_os(0.2) assert cell.flux == 1.2 + def test_reverseTopology(): - ''' + """ testing flowCore.reverseTopology() for different examples of dir graphs - ''' - - testGraph = {0: [1,2,3], - 1: [2,4], - 2: [4,5], - 3: [2], - 4: [], - 5: [6], - 6: []} - - testGraphReverse = {0: [], - 1: [0], - 2: [0,1,3], - 3: [0], - 4: [2,1], - 5: [2], - 6: [5]} + """ + + testGraph = {0: [1, 2, 3], 1: [2, 4], 2: [4, 5], 3: [2], 4: [], 5: [6], 6: []} + + testGraphReverse = {0: [], 1: [0], 2: [0, 1, 3], 3: [0], 4: [2, 1], 5: [2], 6: [5]} reverseGraphCalc = flowCore.reverseTopology(testGraph) @@ -49,33 +53,37 @@ def test_reverseTopology(): setCalcChildren = set(testGraphReverse[key]) assert setTestChildren == setCalcChildren - testGraph = {0:[1,2,3], - 1:[4], - 2:[5,6], - 3:[7,8], - 4:[9], - 5:[9,10], - 6:[10], - 7:[11], - 8:[], - 9:[], - 10:[12], - 11:[], - 12:[]} - - testGraphReverse = {0:[], - 1:[0], - 2:[0], - 3:[0], - 4:[1], - 5:[2], - 6:[2], - 7:[3], - 8:[3], - 9:[4,5], - 10:[5,6], - 11:[7], - 12:[10]} + testGraph = { + 0: [1, 2, 3], + 1: [4], + 2: [5, 6], + 3: [7, 8], + 4: [9], + 5: [9, 10], + 6: [10], + 7: [11], + 8: [], + 9: [], + 10: [12], + 11: [], + 12: [], + } + + testGraphReverse = { + 0: [], + 1: [0], + 2: [0], + 3: [0], + 4: [1], + 5: [2], + 6: [2], + 7: [3], + 8: [3], + 9: [4, 5], + 10: [5, 6], + 11: [7], + 12: [10], + } reverseGraphCalc = flowCore.reverseTopology(testGraph) @@ -87,93 +95,60 @@ def test_reverseTopology(): def test_backTracking(): - ''' + """ testing flowCore.backTracking() for different examples of dir graphs - basic graphs are the same as in test_reverseTopology() with added 'infra' values as valueDicts - ''' - - testGraph = {0: [1,2,3], - 1: [2,4], - 2: [4,5], - 3: [2], - 4: [], - 5: [6], - 6: []} - - testValsIn = {0: 0, - 1: 0, - 2: 0, - 3: 0, - 4: 3, - 5: 0, - 6: 2} - - testValsBT = {0: 3, - 1: 3, - 2: 3, - 3: 3, - 4: 3, - 5: 2, - 6: 2} - + """ + + testGraph = {0: [1, 2, 3], 1: [2, 4], 2: [4, 5], 3: [2], 4: [], 5: [6], 6: []} + + testValsIn = {0: 0, 1: 0, 2: 0, 3: 0, 4: 3, 5: 0, 6: 2} + + testValsBT = {0: 3, 1: 3, 2: 3, 3: 3, 4: 3, 5: 2, 6: 2} + calcValsBT = flowCore.backTracking(testGraph, testValsIn) - for key,item in calcValsBT.items(): + for key, item in calcValsBT.items(): assert calcValsBT[key] == testValsBT[key] - - testGraph = {0:[1,2,3], - 1:[4], - 2:[5,6], - 3:[7,8], - 4:[9], - 5:[9,10], - 6:[10], - 7:[11], - 8:[], - 9:[], - 10:[12], - 11:[], - 12:[]} - - testValsIn = {0:0, - 1:0, - 2:0, - 3:0, - 4:0, - 5:0, - 6:0, - 7:0, - 8:0, - 9:1, - 10:0, - 11:3, - 12:2} - - testValsBT = {0:3, - 1:1, - 2:2, - 3:3, - 4:1, - 5:2, - 6:2, - 7:3, - 8:0, - 9:1, - 10:2, - 11:3, - 12:2} - + + testGraph = { + 0: [1, 2, 3], + 1: [4], + 2: [5, 6], + 3: [7, 8], + 4: [9], + 5: [9, 10], + 6: [10], + 7: [11], + 8: [], + 9: [], + 10: [12], + 11: [], + 12: [], + } + + testValsIn = {0: 0, 1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0, 9: 1, 10: 0, 11: 3, 12: 2} + + testValsBT = {0: 3, 1: 1, 2: 2, 3: 3, 4: 1, 5: 2, 6: 2, 7: 3, 8: 0, 9: 1, 10: 2, 11: 3, 12: 2} + calcValsBT = flowCore.backTracking(testGraph, testValsIn) - for key,item in calcValsBT.items(): + for key, item in calcValsBT.items(): assert calcValsBT[key] == testValsBT[key] def test_calculation(): dem = np.array( - [[40, 40, 40, 40, 40], [30, 30, 30, 30, 30], [20, 20, 20, 20, 20], [10, 10, 10, 10, 10], [0, 0, 0, 0, 0]]) + [ + [40, 40, 40, 40, 40], + [30, 30, 30, 30, 30], + [20, 20, 20, 20, 20], + [10, 10, 10, 10, 10], + [0, 0, 0, 0, 0], + ] + ) infra = None pra = np.array([[0, 0, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]]) alpha = 10 @@ -194,7 +169,7 @@ def test_calculation(): } fluxDistOldVersionBool = False previewMode = False - outputs = ['travelLengthMin', 'flux'] + outputs = ["travelLengthMin", "flux"] forestArray = None forestParams = None relOutputParams = { @@ -224,14 +199,14 @@ def test_calculation(): relOutputParams, ] - flux = ones_like(dem) * -9999. + flux = np.ones_like(dem) * -9999.0 flux[[1, 2, 3], [2, 2, 2]] = 1 depFluxSum = np.zeros_like(dem) - routFluxSum = np.where(flux == 1, 1., 0.) - travelLengthMin = ones_like(dem) * -9999. + routFluxSum = np.where(flux == 1, 1.0, 0.0) + travelLengthMin = np.ones_like(dem) * -9999.0 travelLengthMin[1, 2] = 0 - travelLengthMin[2, 2] = np.sqrt(cellsize ** 2) - travelLengthMin[3, 2] = 2 * np.sqrt(cellsize ** 2) + travelLengthMin[2, 2] = np.sqrt(cellsize**2) + travelLengthMin[3, 2] = 2 * np.sqrt(cellsize**2) results = flowCore.calculation(args) assert len(results) == 16 @@ -239,11 +214,309 @@ def test_calculation(): assert np.all(results[10] == routFluxSum) assert np.all(results[11] == depFluxSum) assert np.all(results[8] == travelLengthMin) - assert np.all(results[7] == np.ones_like(flux) * -9999) + assert results[7] == None + + +def createTestRaster(pathTestFolder, rasterName): + + # first create test raster and save in test folder + testRaster = np.zeros((10, 10)) + + testRaster[2:5, 2:5] = 1 + testRaster[6:9, 5:9] = 2 + testRaster[5, 5:8] = 2 + testRaster[0, 0] = 3 + cellsize = 10 + nrows, ncols = testRaster.shape + + header = { + "cellsize": cellsize, + "nrows": nrows, + "ncols": ncols, + "xllcenter": 0, + "yllcenter": 0, + "nodata_value": -9999, + "driver": "GTiff", + "crs": "EPSG:4326", + } + # convert lower-left center to upper-left corner + x_ul = header["xllcenter"] - cellsize / 2 + y_ul = header["yllcenter"] + nrows * cellsize - cellsize / 2 + + transform = rasterio.transform.from_origin(x_ul, y_ul, cellsize, cellsize) + header["transform"] = transform + + # write flipped raster, the read raster function does also flip the raster + IOf.writeResultToRaster(header, testRaster, pathTestFolder / rasterName, useCompression=False, flip=True) + del testRaster + return header + + +def test_tileRaster(): + pathTestFolder = pathlib.Path("avaframe/tests/data/testCom4") + rasterName = "testRaster" + ext = ".tif" + + tileName = "testTile" + pathTempFolder = pathlib.Path("avaframe/tests/data/testCom4/tmp") + xDim = 4 + yDim = 4 + U = 1 + if os.path.exists(pathTempFolder) is False: + os.makedirs(pathTempFolder) + + createTestRaster(pathTestFolder, rasterName) + + SPAM.tileRaster(pathTestFolder / f"{rasterName}{ext}", tileName, pathTempFolder, xDim, yDim, U) + mergedRaster = SPAM.mergeRaster(pathTempFolder, tileName) + + testData = IOf.readRaster(pathTestFolder / f"{rasterName}{ext}", noDataToNan=False) + testRaster = testData["rasterData"] + + nTiles = pickle.load(open(pathTempFolder / "nTiles", "rb")) + ext00 = pickle.load(open(pathTempFolder / "ext_0_0", "rb")) + ext03 = pickle.load(open(pathTempFolder / "ext_0_3", "rb")) + ext10 = pickle.load(open(pathTempFolder / "ext_1_0", "rb")) + ext21 = pickle.load(open(pathTempFolder / "ext_2_1", "rb")) + + tile00 = np.load(pathTempFolder / "testTile_0_0.npy") + tile03 = np.load(pathTempFolder / "testTile_0_3.npy") + tile03Test = testRaster[0:yDim, 2 * xDim - 2 * U : 3 * xDim - 2 * U] + tile00Test = testRaster[0:xDim, 0:yDim] + + assert np.all(testRaster == mergedRaster) + assert nTiles == (3, 3) + assert tile00.shape == (4, 4) + assert np.all(tile00 == tile00Test) + assert np.all(tile03 == tile03Test) + assert ext00 == ((0, xDim), (0, yDim)) + assert ext03 == ((0, xDim), (2 * yDim - 2 * U, 3 * yDim - 2 * U)) + assert ext10 == ((xDim - 2 * U, 2 * xDim - 2 * U), (0, yDim)) + assert ext21 == ((2 * xDim - 4 * U, 3 * xDim - 4 * U), (yDim - 2 * U, 2 * yDim - 2 * U)) + + +def test_mergeDict(): + pathTestFolder = pathlib.Path("avaframe/tests/data/testCom4") + rasterName = "testRaster" + ext = ".tif" + pathRaster = pathTestFolder / (rasterName + ext) + tileName = "testTile" + pathTempFolder = pathlib.Path("avaframe/tests/data/testCom4/tmp") + xDim = 4 + yDim = 4 + U = 1 + if os.path.exists(pathTempFolder) is False: + os.makedirs(pathTempFolder) + + createTestRaster(pathTestFolder, rasterName) + + dictName = "testDict" + + SPAM.tileRaster(pathRaster, tileName, pathTempFolder, xDim, yDim, U) + nTiles = pickle.load(open(pathTempFolder / "nTiles", "rb")) + + for i in range(nTiles[0] + 1): + for j in range(nTiles[1] + 1): + tile = np.load(pathTempFolder / f"{tileName}_{i}_{j}.npy") + rows, cols = np.where(tile > 0) + dictSmallTile = {(r, c): tile[r, c] for r, c in zip(rows, cols)} + saveDict = open(pathTempFolder / ("%s_%s_%s.pickle" % (dictName, i, j)), "wb") + pickle.dump(dictSmallTile, saveDict) + saveDict.close() + + mergedDict = SPAM.mergeDict(pathTempFolder, dictName) + + mergedDictRef = { + (0, 0): 3, + (2, 2): 1, + (2, 3): 1, + (2, 4): 1, + (3, 2): 1, + (3, 3): 1, + (3, 4): 1, + (4, 2): 1, + (4, 3): 1, + (4, 4): 1, + (6, 5): 2, + (6, 6): 2, + (6, 7): 2, + (6, 8): 2, + (7, 5): 2, + (7, 6): 2, + (7, 7): 2, + (7, 8): 2, + (8, 5): 2, + (8, 6): 2, + (8, 7): 2, + (8, 8): 2, + (5, 5): 2, + (5, 6): 2, + (5, 7): 2, + } + + assert mergedDict.keys() == mergedDictRef.keys() + for k in mergedDict: + assert np.all(mergedDict[k] == mergedDictRef[k]) + + # manipulate dictionaries that they overlap + + for i in range(nTiles[0] + 1): + for j in range(nTiles[1] + 1): + tile = np.load(pathTempFolder / f"{tileName}_{i}_{j}.npy") + rows, cols = np.where(tile > 0) + dictSmallTile = {(r, c): tile[r, c] for r, c in zip(rows, cols)} + if i == 2 and j == 1: + dictSmallTile[(1, 2)] = [1, 2] + dictSmallTile[(1, 1)] = [1, 2] + saveDict = open(pathTempFolder / ("%s_%s_%s.pickle" % (dictName, i, j)), "wb") + pickle.dump(dictSmallTile, saveDict) + saveDict.close() + + mergedDict = SPAM.mergeDict(pathTempFolder, dictName) + + mergedDictRef = { + (5, 4): [1, 2], + (5, 3): [1, 2], + (0, 0): 3, + (2, 2): 1, + (2, 3): 1, + (2, 4): 1, + (3, 2): 1, + (3, 3): 1, + (3, 4): 1, + (4, 2): 1, + (4, 3): 1, + (4, 4): 1, + (6, 5): 2, + (6, 6): 2, + (6, 7): 2, + (6, 8): 2, + (7, 5): 2, + (7, 6): 2, + (7, 7): 2, + (7, 8): 2, + (8, 5): 2, + (8, 6): 2, + (8, 7): 2, + (8, 8): 2, + (5, 5): 2, + (5, 6): 2, + (5, 7): 2, + } + assert mergedDict.keys() == mergedDictRef.keys() + for k in mergedDict: + assert np.all(mergedDict[k] == mergedDictRef[k]) + + +def test_mergeDictToRaster(): + pathTempFolder = pathlib.Path("avaframe/tests/data/testCom4/tmp") + tileName = "testTile" + dictName = "testDict" + + nTiles = pickle.load(open(pathTempFolder / "nTiles", "rb")) + + for i in range(nTiles[0] + 1): + for j in range(nTiles[1] + 1): + tile = np.load(pathTempFolder / f"{tileName}_{i}_{j}.npy") + rows, cols = np.where(tile > 0) + dictSmallTile = {(r, c): tile[r, c] for r, c in zip(rows, cols)} + saveDict = open(pathTempFolder / ("%s_%s_%s.pickle" % (dictName, i, j)), "wb") + pickle.dump(dictSmallTile, saveDict) + saveDict.close() + + mergedRaster = SPAM.mergeDictToRaster(pathTempFolder, dictName) + + testData = IOf.readRaster(pathlib.Path("avaframe/tests/data/testCom4/testRaster.tif")) + testRaster = testData["rasterData"] + # due to no overlap: + testRaster[testRaster > 0] = 1 + assert np.all(mergedRaster == testRaster) + + # manipulate dictionaries that they overlap + + for i in range(nTiles[0] + 1): + for j in range(nTiles[1] + 1): + tile = np.load(pathTempFolder / f"{tileName}_{i}_{j}.npy") + rows, cols = np.where(tile > 0) + dictSmallTile = {(r, c): tile[r, c] for r, c in zip(rows, cols)} + if i == 2 and j == 1: + dictSmallTile[(1, 2)] = [1, 2] + dictSmallTile[(1, 1)] = [1, 2] + saveDict = open(pathTempFolder / ("%s_%s_%s.pickle" % (dictName, i, j)), "wb") + pickle.dump(dictSmallTile, saveDict) + saveDict.close() + + mergedRaster = SPAM.mergeDictToRaster(pathTempFolder, dictName) + testRaster[5, 4] = 2 + testRaster[5, 3] = 2 + + assert np.all(mergedRaster == testRaster) + + +def test_mergeDictToPolygon(): + pathTestFolder = pathlib.Path("avaframe/tests/data/testCom4") + rasterName = "testRaster" + ext = ".tif" + pathRaster = pathTestFolder / (rasterName + ext) + tileName = "testTile" + pathTempFolder = pathlib.Path("avaframe/tests/data/testCom4/tmp") + xDim = 4 + yDim = 4 + U = 1 + if os.path.exists(pathTempFolder) is False: + os.makedirs(pathTempFolder) + + rasterHeader = createTestRaster(pathTestFolder, rasterName) + + dictName = "testDict" + + SPAM.tileRaster(pathRaster, tileName, pathTempFolder, xDim, yDim, U) + nTiles = pickle.load(open(pathTempFolder / "nTiles", "rb")) + + for i in range(nTiles[0] + 1): + for j in range(nTiles[1] + 1): + tile = np.load(pathTempFolder / f"{tileName}_{i}_{j}.npy") + rows, cols = np.where(tile > 0) + dictSmallTile = {(r, c): tile[r, c] for r, c in zip(rows, cols)} + saveDict = open(pathTempFolder / ("%s_%s_%s.pickle" % (dictName, i, j)), "wb") + pickle.dump(dictSmallTile, saveDict) + saveDict.close() + + gdfPathPolygons = SPAM.mergeDictToPolygon(pathTempFolder, dictName, rasterHeader) + refPolygons = gpd.read_file(pathTestFolder / "refPolygon.geojson") + + assert len(gdfPathPolygons) == 3 + assert np.all(gdfPathPolygons["PRA_id"] == refPolygons["PRA_id"]) + assert np.all(gdfPathPolygons == refPolygons) + + # manipulate dictionaries that the polygons overlap + + for i in range(nTiles[0] + 1): + for j in range(nTiles[1] + 1): + tile = np.load(pathTempFolder / f"{tileName}_{i}_{j}.npy") + rows, cols = np.where(tile > 0) + dictSmallTile = {(r, c): tile[r, c] for r, c in zip(rows, cols)} + if i == 2 and j == 1: + dictSmallTile[(1, 2)] = [1, 2] + dictSmallTile[(1, 1)] = [1, 2] + saveDict = open(pathTempFolder / ("%s_%s_%s.pickle" % (dictName, i, j)), "wb") + pickle.dump(dictSmallTile, saveDict) + saveDict.close() + + gdfPathPolygons = SPAM.mergeDictToPolygon(pathTempFolder, dictName, rasterHeader) + refPolygons = gpd.read_file(pathTestFolder / "refPolygon_manipulated.geojson") + + assert len(gdfPathPolygons) == 3 + assert np.all(gdfPathPolygons["PRA_id"] == refPolygons["PRA_id"]) + assert np.all(gdfPathPolygons.geometry.geom_equals(refPolygons.geometry)) -if __name__=='__main__': +if __name__ == "__main__": test_add_os() test_reverseTopology() test_backTracking() test_calculation() + test_tileRaster() + test_mergeDict() + test_mergeDictToRaster() + test_mergeDictToPolygon() From 02fe89f4a817c63446c6de624eec7ce42e6c4429 Mon Sep 17 00:00:00 2001 From: Paula Date: Mon, 18 May 2026 16:09:54 +0200 Subject: [PATCH 4/6] delete release volume min and max output option --- avaframe/com4FlowPy/com4FlowPy.py | 184 ++++++++++++++------------ avaframe/com4FlowPy/com4FlowPyCfg.ini | 4 +- avaframe/com4FlowPy/flowClass.py | 51 +++---- avaframe/com4FlowPy/flowCore.py | 62 +-------- avaframe/runCom4FlowPy.py | 52 +++++--- docs/moduleCom4FlowPy.rst | 2 - 6 files changed, 157 insertions(+), 198 deletions(-) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index 935bbc7d1..9930fe626 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -5,6 +5,7 @@ mainly handling input of data, model params and output of model results """ + # Load modules import pathlib import numpy as np @@ -54,7 +55,9 @@ def com4FlowPyMain(cfgPath, cfgSetup): # Flow-Py parameters modelParameters["alpha"] = cfgSetup.getfloat("alpha") # float(cfgSetup["alpha"]) modelParameters["exp"] = cfgSetup.getfloat("exp") # float(cfgSetup["exp"]) - modelParameters["flux_threshold"] = cfgSetup.getfloat("flux_threshold") # float(cfgSetup["flux_threshold"]) + modelParameters["flux_threshold"] = cfgSetup.getfloat( + "flux_threshold" + ) # float(cfgSetup["flux_threshold"]) modelParameters["max_z"] = cfgSetup.getfloat("max_z") # float(cfgSetup["max_z"]) # Flags for use of Forest and/or Infrastructure @@ -91,15 +94,15 @@ def com4FlowPyMain(cfgPath, cfgSetup): modelPaths["tempDir"] = cfgPath["tempDir"] modelPaths["uid"] = cfgPath["uid"] modelPaths["timeString"] = cfgPath["timeString"] - modelPaths["outputFileList"] = cfgPath["outputFiles"].split('|') + modelPaths["outputFileList"] = cfgPath["outputFiles"].split("|") modelPaths["outputNoDataValue"] = cfgPath["outputNoDataValue"] modelPaths["useCompression"] = cfgPath["useCompression"] modelPaths["outputFileFormat"] = cfgPath["outputFileFormat"] if modelPaths["outputFileFormat"] in [".asc", ".ASC"]: - modelPaths["outputFileFormat"] = '.asc' + modelPaths["outputFileFormat"] = ".asc" else: - modelPaths["outputFileFormat"] = '.tif' + modelPaths["outputFileFormat"] = ".tif" # check if 'customDirs' are used - alternative is 'default' AvaFrame Folder Structure modelPaths["useCustomDirs"] = True if cfgPath["customDirs"] == "True" else False @@ -109,9 +112,15 @@ def com4FlowPyMain(cfgPath, cfgSetup): # Multiprocessing Options MPOptions = {} MPOptions["nCPU"] = cfgSetup.getint("cpuCount") # int(cfgSetup["cpuCount"]) #number of CPUs to use - MPOptions["procPerCPU"] = cfgSetup.getint("procPerCPUCore") # int(cfgSetup["procPerCPUCore"]) #processes per core - MPOptions["chunkSize"] = cfgSetup.getint("chunkSize") # int(cfgSetup["chunkSize"]) # default task size for MP - MPOptions["maxChunks"] = cfgSetup.getint("maxChunks") # int(cfgSetup["maxChunks"]) # max number of tasks for MP + MPOptions["procPerCPU"] = cfgSetup.getint( + "procPerCPUCore" + ) # int(cfgSetup["procPerCPUCore"]) #processes per core + MPOptions["chunkSize"] = cfgSetup.getint( + "chunkSize" + ) # int(cfgSetup["chunkSize"]) # default task size for MP + MPOptions["maxChunks"] = cfgSetup.getint( + "maxChunks" + ) # int(cfgSetup["maxChunks"]) # max number of tasks for MP # check if calculation with infrastructure if modelParameters["infraBool"]: @@ -126,12 +135,24 @@ def com4FlowPyMain(cfgPath, cfgSetup): forestParams["forestModule"] = cfgSetup.get("forestModule") modelPaths["forestPath"] = cfgPath["forestPath"] # 'forestFriction' and 'forestDetrainment' parameters - forestParams["maxAddedFriction"] = cfgSetup.getfloat("maxAddedFrictionFor") # float(cfgSetup["maxAddedFr.."]) - forestParams["minAddedFriction"] = cfgSetup.getfloat("minAddedFrictionFor") # float(cfgSetup["minAddedFr.."]) - forestParams["velThForFriction"] = cfgSetup.getfloat("velThForFriction") # float(cfgSetup["velThForFriction"]) - forestParams["maxDetrainment"] = cfgSetup.getfloat("maxDetrainmentFor") # float(cfgSetup["maxDetrainmentFor"]) - forestParams["minDetrainment"] = cfgSetup.getfloat("minDetrainmentFor") # float(cfgSetup["minDetrainmentFor"]) - forestParams["velThForDetrain"] = cfgSetup.getfloat("velThForDetrain") # float(cfgSetup["velThForDetrain"]) + forestParams["maxAddedFriction"] = cfgSetup.getfloat( + "maxAddedFrictionFor" + ) # float(cfgSetup["maxAddedFr.."]) + forestParams["minAddedFriction"] = cfgSetup.getfloat( + "minAddedFrictionFor" + ) # float(cfgSetup["minAddedFr.."]) + forestParams["velThForFriction"] = cfgSetup.getfloat( + "velThForFriction" + ) # float(cfgSetup["velThForFriction"]) + forestParams["maxDetrainment"] = cfgSetup.getfloat( + "maxDetrainmentFor" + ) # float(cfgSetup["maxDetrainmentFor"]) + forestParams["minDetrainment"] = cfgSetup.getfloat( + "minDetrainmentFor" + ) # float(cfgSetup["minDetrainmentFor"]) + forestParams["velThForDetrain"] = cfgSetup.getfloat( + "velThForDetrain" + ) # float(cfgSetup["velThForDetrain"]) # 'forestFrictionLayer' parameter forestParams["fFrLayerType"] = cfgSetup.get("forestFrictionLayerType") # skipForestDist - no forest friciton effect assumed while distance along path <= skipForestDist @@ -165,11 +186,6 @@ def com4FlowPyMain(cfgPath, cfgSetup): 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 ... @@ -319,7 +335,10 @@ def checkInputLayerDimensions(modelParameters, modelPaths): if modelParameters["forestBool"]: _forestHeader = IOf.readRasterHeader(modelPaths["forestPath"]) - if _demHeader["ncols"] == _forestHeader["ncols"] and _demHeader["nrows"] == _forestHeader["nrows"]: + if ( + _demHeader["ncols"] == _forestHeader["ncols"] + and _demHeader["nrows"] == _forestHeader["nrows"] + ): log.info("Forest Layer ok!") else: log.error("Error: Forest Layer doesn't match DEM!") @@ -327,7 +346,10 @@ def checkInputLayerDimensions(modelParameters, modelPaths): if modelParameters["varUmaxBool"]: _varUmaxHeader = IOf.readRasterHeader(modelPaths["varUmaxPath"]) - if _demHeader["ncols"] == _varUmaxHeader["ncols"] and _demHeader["nrows"] == _varUmaxHeader["nrows"]: + if ( + _demHeader["ncols"] == _varUmaxHeader["ncols"] + and _demHeader["nrows"] == _varUmaxHeader["nrows"] + ): log.info("uMax Limit Layer ok!") else: log.error("Error: uMax Limit Layer doesn't match DEM!") @@ -335,7 +357,10 @@ def checkInputLayerDimensions(modelParameters, modelPaths): if modelParameters["varAlphaBool"]: _varAlphaHeader = IOf.readRasterHeader(modelPaths["varAlphaPath"]) - if _demHeader["ncols"] == _varAlphaHeader["ncols"] and _demHeader["nrows"] == _varAlphaHeader["nrows"]: + if ( + _demHeader["ncols"] == _varAlphaHeader["ncols"] + and _demHeader["nrows"] == _varAlphaHeader["nrows"] + ): log.info("variable Alpha Layer ok!") else: log.error("Error: variable Alpha Layer doesn't match DEM!") @@ -343,7 +368,10 @@ def checkInputLayerDimensions(modelParameters, modelPaths): if modelParameters["varExponentBool"]: _varExponentHeader = IOf.readRasterHeader(modelPaths["varExponentPath"]) - if _demHeader["ncols"] == _varExponentHeader["ncols"] and _demHeader["nrows"] == _varExponentHeader["nrows"]: + if ( + _demHeader["ncols"] == _varExponentHeader["ncols"] + and _demHeader["nrows"] == _varExponentHeader["nrows"] + ): log.info("variable exponent Layer ok!") else: log.error("Error: variable exponent Layer doesn't match DEM!") @@ -352,8 +380,10 @@ def checkInputLayerDimensions(modelParameters, modelPaths): log.info("========================") except Exception as ex: - log.error("could not read all required Input Layers, please re-check files and paths provided in .ini files") - log.error('Error occured: %s' % ex) + log.error( + "could not read all required Input Layers, please re-check files and paths provided in .ini files" + ) + log.error("Error occured: %s" % ex) # return sys.exit(1) @@ -369,17 +399,17 @@ def checkInputParameterValues(modelParameters, modelPaths): modelPaths: dict contains paths to input files """ - alpha = modelParameters['alpha'] - if (alpha < 0 or alpha > 90): + alpha = modelParameters["alpha"] + if alpha < 0 or alpha > 90: log.error("Error: Alpha value is not within a physically sensible range ([0,90]).") sys.exit(1) - zDelta = modelParameters['max_z'] - if (zDelta < 0 or zDelta > 8848): + zDelta = modelParameters["max_z"] + if zDelta < 0 or zDelta > 8848: log.error("Error: zDeltaMaxLimit value is not within a physically sensible range ([0,8848]).") sys.exit(1) - exp = modelParameters['exp'] + exp = modelParameters["exp"] if exp < 0: log.error("Error: Exponent value is not within a physically sensible range (> 0).") sys.exit(1) @@ -388,7 +418,7 @@ def checkInputParameterValues(modelParameters, modelPaths): if modelParameters["varAlphaBool"]: data = IOf.readRaster(modelPaths["varAlphaPath"]) - rasterValues = data["rasterData"] + rasterValues = data["rasterData"] rasterValues[rasterValues < 0] = np.nan # handle different noData values if np.any(rasterValues > 90, where=~np.isnan(rasterValues)): log.error("Error: Not all Alpha-raster values are within a physically sensible range ([0,90]),\ @@ -397,15 +427,17 @@ def checkInputParameterValues(modelParameters, modelPaths): if modelParameters["varUmaxBool"]: data = IOf.readRaster(modelPaths["varUmaxPath"]) - rasterValues = data["rasterData"] + rasterValues = data["rasterData"] rasterValues[rasterValues < 0] = np.nan - if modelParameters["varUmaxType"].lower() == 'umax': + if modelParameters["varUmaxType"].lower() == "umax": _maxVal = 1500 # ~sqrt(8848*2*9.81) else: _maxVal = 8848 if np.any(rasterValues > _maxVal, where=~np.isnan(rasterValues)): - log.error("Error: Not all zDeltaMaxLimit-raster values are within a physically sensible range \ - ([0, 8848 m] or [0, 1500 m/s]), in respective startcells the general zDeltaMax value is used.") + log.error( + "Error: Not all zDeltaMaxLimit-raster values are within a physically sensible range \ + ([0, 8848 m] or [0, 1500 m/s]), in respective startcells the general zDeltaMax value is used." + ) _checkVarParams = False if _checkVarParams: @@ -452,16 +484,24 @@ def tileInputLayers(modelParameters, modelPaths, rasterAttributes, tilingParamet log.info("---------------------") SPAM.tileRaster(modelPaths["demPath"], "dem", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) - SPAM.tileRaster(modelPaths["releasePathWork"], "init", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U, isInit=True) + SPAM.tileRaster( + modelPaths["releasePathWork"], "init", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U, isInit=True + ) if modelParameters["infraBool"]: SPAM.tileRaster(modelPaths["infraPath"], "infra", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) if modelParameters["varUmaxBool"]: - SPAM.tileRaster(modelPaths["varUmaxPath"], "varUmax", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) + SPAM.tileRaster( + modelPaths["varUmaxPath"], "varUmax", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U + ) if modelParameters["varAlphaBool"]: - SPAM.tileRaster(modelPaths["varAlphaPath"], "varAlpha", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) + SPAM.tileRaster( + modelPaths["varAlphaPath"], "varAlpha", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U + ) if modelParameters["varExponentBool"]: - SPAM.tileRaster(modelPaths["varExponentPath"], "varExponent", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) + 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"]: @@ -526,7 +566,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): contains model input parameters (from .ini - file) """ _uid = modelPaths["uid"] - _outputs = set(modelPaths['outputFileList']) + _outputs = set(modelPaths["outputFileList"]) _outputNoDataValue = modelPaths["outputNoDataValue"] demHeader = IOf.readRasterHeader(modelPaths["demPath"]) @@ -549,7 +589,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): # 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: + if "cellCounts" in _outputs: cellCounts = defineNotAffectedCells(cellCounts, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -557,11 +597,11 @@ def mergeAndWriteResults(modelPaths, modelOptions): 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: + if "zDelta" in _outputs: zDelta = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta") zDelta = defineNotAffectedCells(zDelta, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( @@ -570,12 +610,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), flip=True, useCompression=useCompression, - ) + ) del zDelta del output log.info("com4_{}_{}_zdelta is written".format(_uid, _ts)) - if 'flux' in _outputs: + if "flux" in _outputs: flux = SPAM.mergeRaster(modelPaths["tempDir"], "res_flux") flux = defineNotAffectedCells(flux, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( @@ -584,12 +624,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): 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: + if "zDeltaSum" in _outputs: zDeltaSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta_sum", method="sum") zDeltaSum = defineNotAffectedCells(zDeltaSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( @@ -603,7 +643,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): del output log.info("com4_{}_{}_zDeltaSum is written".format(_uid, _ts)) - if 'routFluxSum' in _outputs: + if "routFluxSum" in _outputs: routFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_rout_flux_sum", method="sum") routFluxSum = defineNotAffectedCells(routFluxSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( @@ -617,7 +657,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): del output log.info("com4_{}_{}_routFluxSum is written".format(_uid, _ts)) - if 'depFluxSum' in _outputs: + if "depFluxSum" in _outputs: depFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_dep_flux_sum", method="sum") depFluxSum = defineNotAffectedCells(depFluxSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( @@ -659,7 +699,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): del output log.info("com4_{}_{}_fpTravelAngleMin is written".format(_uid, _ts)) - if 'slTravelAngle' in _outputs: + if "slTravelAngle" in _outputs: slTa = SPAM.mergeRaster(modelPaths["tempDir"], "res_sl") slTa = defineNotAffectedCells(slTa, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( @@ -701,7 +741,6 @@ def mergeAndWriteResults(modelPaths, modelOptions): del output log.info("com4_{}_{}_travelLengthMin is written".format(_uid, _ts)) - if modelOptions["infraBool"]: backcalc = SPAM.mergeRaster(modelPaths["tempDir"], "res_backcalc") backcalc = defineNotAffectedCells(backcalc, cellCounts, noDataValue=_outputNoDataValue) @@ -711,14 +750,13 @@ 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 = SPAM.mergeRaster(modelPaths["tempDir"], "res_forestInt", method="min") forestInteraction = defineNotAffectedCells( forestInteraction, cellCounts, noDataValue=_outputNoDataValue ) @@ -728,12 +766,11 @@ 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( @@ -742,7 +779,6 @@ def mergeAndWriteResults(modelPaths, modelOptions): 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) @@ -757,36 +793,6 @@ def mergeAndWriteResults(modelPaths, modelOptions): 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 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) @@ -812,11 +818,13 @@ def checkConvertReleaseShp2Tif(modelPaths): dem = IOf.readRaster(modelPaths["demPath"]) demHeader = dem["header"] - dem['originalHeader'] = demHeader + dem["originalHeader"] = demHeader releaseLine = shpConv.SHP2Array(modelPaths["releasePath"], "releasePolygon") thresholdPointInPoly = 0.01 - releaseLine = gT.prepareArea(releaseLine, dem, thresholdPointInPoly, combine=True, checkOverlap=False) + releaseLine = gT.prepareArea( + releaseLine, dem, thresholdPointInPoly, combine=True, checkOverlap=False + ) # give the same header as the dem releaseArea = np.flipud(releaseLine["rasterData"]) if demHeader["driver"] == "AAIGrid": diff --git a/avaframe/com4FlowPy/com4FlowPyCfg.ini b/avaframe/com4FlowPy/com4FlowPyCfg.ini index a56d56b39..97e248553 100644 --- a/avaframe/com4FlowPy/com4FlowPyCfg.ini +++ b/avaframe/com4FlowPy/com4FlowPyCfg.ini @@ -227,8 +227,6 @@ outputNoDataValue = -9999 # depFluxSum # travelLengthMin # fpTravelAngleMin -# relVolMin -# relVolMax # relIdPolygon # relIdCount # if forestInteraction: forestInteraction is automatically added to outputs @@ -262,7 +260,7 @@ useCustomPathDEM = False workDir = demPath = releasePath = -releaseIdPath = +relIdPath = infraPath = forestPath = varUmaxPath = diff --git a/avaframe/com4FlowPy/flowClass.py b/avaframe/com4FlowPy/flowClass.py index c4e21e3e7..8f6c6f6bd 100644 --- a/avaframe/com4FlowPy/flowClass.py +++ b/avaframe/com4FlowPy/flowClass.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- """ - Functions for calculations at the 'cell level' (vgl. D'Amboise et al., 2022) +Functions for calculations at the 'cell level' (vgl. D'Amboise et al., 2022) """ import numpy as np @@ -28,7 +28,6 @@ def __init__( 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 @@ -37,8 +36,12 @@ def __init__( * bool --> isStart * Cell --> startCell """ - self.rowindex = rowindex # index of the Cell in row-direction (i.e. local y-index in the calculation domain) - self.colindex = colindex # index of the Cell in column-direction (i.e. local x-index in the calculation domain) + self.rowindex = ( + rowindex # index of the Cell in row-direction (i.e. local y-index in the calculation domain) + ) + self.colindex = ( + colindex # index of the Cell in column-direction (i.e. local x-index in the calculation domain) + ) self.dem_ng = dem_ng # elevation values in the 3x3 neigbourhood around the Cell self.altitude = dem_ng[1, 1] # elevation value of the cell (central cell of 3x3 neighbourhood) self.cellsize = cellsize # cellsize in meters @@ -60,22 +63,23 @@ def __init__( self.fluxDistOldVersionBool = fluxDistOldVersionBool - self.tanAlpha = np.tan(np.deg2rad(self.alpha)) # moved to constructor, so this doesn't have to be calculated on + self.tanAlpha = np.tan( + np.deg2rad(self.alpha) + ) # moved to constructor, so this doesn't have to be calculated on # every iteration of calc_z_delta(self) self.min_distance = 0 # minimal distance to start-cell (i.e. along shortest path) min_distance >= self.minDistXYZ = 0 # minimal distance to start-cell (Actual 3D lenght, not projected!!) self.max_distance = 0 # NOTE: self.max_distance is never used - maybe remove!? - self.min_gamma = 0 # NOTE: self.min_gamma (assumingly minimal travel angle to cell) never used - maybe remove!? + self.min_gamma = ( + 0 # NOTE: self.min_gamma (assumingly minimal travel angle to cell) never used - maybe remove!? + ) self.max_gamma = 0 self.sl_gamma = 0 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 @@ -117,7 +121,9 @@ def __init__( elif forestParams["fFrLayerType"] == "relative": self.AlphaFor = self.alpha + FSI - self.AlphaFor = max(self.AlphaFor, self.alpha) # Friction in Forest can't be lower than without forest + self.AlphaFor = max( + self.AlphaFor, self.alpha + ) # Friction in Forest can't be lower than without forest self.tanAlphaFor = np.tan(np.deg2rad(self.AlphaFor)) # NOTE: This is a quick hack to check if all values for Detrainment are set to 0 (as provided in the @@ -191,10 +197,6 @@ 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 @@ -213,7 +215,7 @@ def calcDistMin(self, calc3D=False): _dy = abs(parent.rowindex - self.rowindex) * self.cellsize _dz = abs(parent.altitude - self.altitude) _ldistMin.append(math.sqrt(_dx * _dx + _dy * _dy) + parent.min_distance) - _lDistMinXYZ.append(math.sqrt(_dy*_dy + _dy*_dy + _dz*_dz) + parent.minDistXYZ) + _lDistMinXYZ.append(math.sqrt(_dy * _dy + _dy * _dy + _dz * _dz) + parent.minDistXYZ) self.min_distance = np.amin(_ldistMin) self.minDistXYZ = np.amin(_lDistMinXYZ) else: @@ -255,8 +257,8 @@ def calc_z_delta(self): self.z_gamma = self.altitude - self.dem_ng ds = np.array([[self._SQRT2, 1, self._SQRT2], [1, 0, 1], [self._SQRT2, 1, self._SQRT2]]) - if (not self.is_start): - if (not self.forestBool): + if not self.is_start: + if not self.forestBool: self.calcDistMin() else: self.calcDistMin(calc3D=True) @@ -270,7 +272,7 @@ def calc_z_delta(self): _tanAlpha = self.tanAlpha if self.forestModule in ["forestFriction", "forestDetrainment"]: - if (not self.is_start) and (self.FSI > 0.) and (self.skipForestDist < self.minDistXYZ): + if (not self.is_start) and (self.FSI > 0.0) and (self.skipForestDist < self.minDistXYZ): # if forestBool, we assume that forestFriciton is activated # and if FSI > 0 then we also calculate _tanAlpha with forestEffect # NOTE: We also don't assume a forest Effect on potential Start Zells, since this should @@ -284,7 +286,9 @@ def calc_z_delta(self): _slope = (_rest - self.minAddedFrictionForest) / (0 - self.noFrictionEffectZDelta) # y = mx + b, shere z_delta is the x friction = max(self.minAddedFrictionForest, _slope * self.z_delta + _rest) - _alpha_calc = self.alpha + max(0, friction) # NOTE: not sure what this does, seems redundant! + _alpha_calc = self.alpha + max( + 0, friction + ) # NOTE: not sure what this does, seems redundant! else: _alpha_calc = self.alpha + self.minAddedFrictionForest @@ -302,8 +306,7 @@ def calc_z_delta(self): self.z_delta_neighbour[self.z_delta_neighbour > self.max_z_delta] = self.max_z_delta def calc_tanbeta(self): - """calculates the normalized terrain based routing - """ + """calculates the normalized terrain based routing""" _ds = np.array([[self._SQRT2, 1, self._SQRT2], [1, 1, 1], [self._SQRT2, 1, self._SQRT2]]) _distance = _ds * self.cellsize @@ -330,7 +333,9 @@ def calc_persistence(self): dx = parent.colindex - self.colindex dy = parent.rowindex - self.rowindex - self.no_flow[dy + 1, dx + 1] = 0 # 3x3 Matrix of ones, every parent gets a 0, no flow to a parent field + self.no_flow[dy + 1, dx + 1] = ( + 0 # 3x3 Matrix of ones, every parent gets a 0, no flow to a parent field + ) maxweight = parent.z_delta # Old Calculation @@ -460,7 +465,7 @@ def forest_detrainment(self): _noDetrainmentEffectZdelta = self.noDetrainmentEffectZdelta # detrainment effect scaled to forest, 0 for non-forest - _rest = (self.maxAddedDetrainmentForest * self.FSI) + _rest = self.maxAddedDetrainmentForest * self.FSI # rise over run (should be negative slope) slope = (_rest - self.minAddedDetrainmentForest) / (0 - _noDetrainmentEffectZdelta) # y=mx+b, where zDelta is x diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index 318fd00f3..4652ebac8 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -145,7 +145,6 @@ def run(optTuple): 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) @@ -204,11 +203,6 @@ def run(optTuple): else: relIdArray = None - if relVolBool: - relVolArray = release.copy() - else: - relVolArray = None - varParams = { "varUmaxBool": varUmaxBool, "varUmaxArray": varUmaxArray, @@ -220,8 +214,6 @@ def run(optTuple): relOutputParams = { "relIdBool": relIdBool, "relIdArray": relIdArray, - "relVolBool": relVolBool, - "relVolArray": relVolArray, } # convert release areas to binary (0: no release areas, 1: release areas) @@ -295,8 +287,6 @@ 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 = [] @@ -313,8 +303,6 @@ def run(optTuple): travelLengthMaxList = [] travelLengthMinList = [] processedStartCellIdList = [] - relVolMinList = [] - relVolMaxList = [] if forestInteraction: forestIntList = [] @@ -336,9 +324,7 @@ def run(optTuple): depFluxSumList.append(res[11]) processedStartCellIdList.append(res[12]) if forestInteraction: - forestIntList.append(res[15]) - relVolMinList.append(res[13]) - relVolMaxList.append(res[14]) + forestIntList.append(res[13]) logging.info("Calculation finished, getting results.") for i in range(len(zDeltaList)): @@ -373,14 +359,7 @@ 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: @@ -405,8 +384,6 @@ 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: @@ -509,9 +486,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): 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] @@ -555,16 +530,6 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): else: travelLengthMaxArray = None - if "relVolMin" in outputs: - relVolMinArray = np.ones_like(dem, dtype=np.float32) * -9999 - else: - relVolMinArray = None - - if "relVolMax" in outputs: - relVolMaxArray = np.zeros_like(dem, dtype=np.float32) - else: - relVolMaxArray = None - if infraBool: backcalc = np.ones_like(dem, dtype=np.int32) * -9999 else: @@ -612,10 +577,6 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): startcellId = relIdArray[row_idx, col_idx] else: startcellId = None - if relVolBool: - startcellVol = relVolArray[row_idx, col_idx] - else: - startcellVol = None startcell = Cell( row_idx, @@ -633,7 +594,6 @@ 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 @@ -674,8 +634,6 @@ 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) @@ -726,7 +684,6 @@ 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, ) ) @@ -781,19 +738,6 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): 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 @@ -842,8 +786,6 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): routFluxSumArray, depFluxSumArray, startCellIdDict, - relVolMinArray, - relVolMaxArray, forestIntArray, ) diff --git a/avaframe/runCom4FlowPy.py b/avaframe/runCom4FlowPy.py index aab8ca789..5a3f2b338 100644 --- a/avaframe/runCom4FlowPy.py +++ b/avaframe/runCom4FlowPy.py @@ -26,7 +26,7 @@ import avaframe.in3Utils.geoTrans as gT -def main(avalancheDir=''): +def main(avalancheDir=""): """this is a wrapper around com4FlowPy.py that handles the following tasks: * reading inputs from (local_)avaframeCfg.ini and (local_)com4FlowPyCfg.ini * constructing cfgPath and cfgSetup dictionaries for passing to com4FlowPy.com4FlowPyMain() @@ -53,18 +53,18 @@ def main(avalancheDir=''): # if "useCustomPaths" == False, we also use the AvaDir Info for the # creation of the simulaiton uid - if avalancheDir != '': - cfgMain['MAIN']['avalancheDir'] = avalancheDir + if avalancheDir != "": + cfgMain["MAIN"]["avalancheDir"] = avalancheDir else: # Load avalanche directory from general configuration file avalancheDir = cfgMain["MAIN"]["avalancheDir"] - cfg['GENERAL']['avaDir'] = cfgMain["MAIN"]["avalancheDir"] + cfg["GENERAL"]["avaDir"] = cfgMain["MAIN"]["avalancheDir"] uid = cfgUtils.cfgHash(cfg) # Clean input directory of old work and output files from module initProj.cleanModuleFiles(avalancheDir, com4FlowPy, deleteOutput=False) # Start logging - log = logUtils.initiateLogger(avalancheDir, logName+'_'+uid) + log = logUtils.initiateLogger(avalancheDir, logName + "_" + uid) log.info("==================================") log.info("MAIN SCRIPT") log.info("Current avalanche: %s", avalancheDir) @@ -90,7 +90,11 @@ def main(avalancheDir=''): # check if simulation with same uid already has results folder if os.path.isdir(cfgPath["resDir"]): log.info("folder with same name already exists - aborting") - log.info("simulation results folder with same .ini parameters already exists: simulation {}".format(uid)) + log.info( + "simulation results folder with same .ini parameters already exists: simulation {}".format( + uid + ) + ) sys.exit(1) else: fU.makeADir(cfgPath["resDir"]) @@ -98,12 +102,12 @@ def main(avalancheDir=''): fU.makeADir(cfgPath["tempDir"]) # writing config to .json file - successToJSON = writeCfgJSON(cfg, uid, cfgPath['outDir']) + successToJSON = writeCfgJSON(cfg, uid, cfgPath["outDir"]) if successToJSON is True: - log.info('wrote config to {}/{}.json'.format(cfgPath['outDir'], uid)) + log.info("wrote config to {}/{}.json".format(cfgPath["outDir"], uid)) else: - log.info('could not write config to {}/{}.json'.format(cfgPath['outDir'], uid)) + log.info("could not write config to {}/{}.json".format(cfgPath["outDir"], uid)) log.error("Exception occurred: %s", str(successToJSON), exc_info=True) cfgPath["deleteTemp"] = "False" @@ -132,14 +136,18 @@ def main(avalancheDir=''): except Exception as e: return e - log = logUtils.initiateLogger(workDir, logName+'_'+uid) + log = logUtils.initiateLogger(workDir, logName + "_" + uid) timeString = datetime.now().strftime("%Y%m%d_%H%M%S") try: os.makedirs(workDir / "res_{}".format(uid)) # (time_string)) - res_dir = workDir / "res_{}".format(uid) # (time_string) + res_dir = workDir / "res_{}".format(uid) # (time_string) except FileExistsError: - log.info("simulation results folder with same .ini parameters already exists: simulation {}".format(uid)) + log.info( + "simulation results folder with same .ini parameters already exists: simulation {}".format( + uid + ) + ) sys.exit(1) try: os.makedirs(workDir / res_dir / "temp") @@ -152,9 +160,9 @@ def main(avalancheDir=''): successToJSON = writeCfgJSON(cfg, uid, workDir) if successToJSON is True: - log.info('wrote config to {}/{}.json'.format(workDir, uid)) + log.info("wrote config to {}/{}.json".format(workDir, uid)) else: - log.info('could not write config to {}/{}.json'.format(workDir, uid)) + log.info("could not write config to {}/{}.json".format(workDir, uid)) log.error("Exception occurred: %s", str(successToJSON), exc_info=True) cfgPath["workDir"] = pathlib.Path(workDir) @@ -270,7 +278,9 @@ def readFlowPyinputs(avalancheDir, cfgFlowPy, log): if cfgFlowPy.getboolean("GENERAL", "variableExponent") is True: 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" + message = ( + f"There is no variable EXPONENT file in supported format provided in {avalancheDir}/EXP" + ) log.error(message) raise AssertionError(message) log.info("variable Exponent file is: %s" % varExponentPath) @@ -339,7 +349,7 @@ def checkOutputFilesFormat(strOutputFiles): """ try: - setA = set(strOutputFiles.split('|')) + setA = set(strOutputFiles.split("|")) setB = set( [ "zDelta", @@ -355,24 +365,22 @@ def checkOutputFilesFormat(strOutputFiles): "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): + 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') + raise ValueError("outputFiles defined in .ini have wrong format - using default settings") except ValueError: # else we return the default options - return 'zDelta|cellCounts|travelLength|fpTravelAngle' + return "zDelta|cellCounts|travelLength|fpTravelAngle" def writeCfgJSON(cfg, uid, workDir): @@ -399,7 +407,7 @@ def writeCfgJSON(cfg, uid, workDir): cfgDict = cfgUtils.convertConfigParserToDict(cfg) try: - with open(workDir / "{}.json".format(uid), 'w') as outfile: + with open(workDir / "{}.json".format(uid), "w") as outfile: jsonDict = json.dumps(cfgDict, sort_keys=True, ensure_ascii=True) outfile.write(jsonDict) diff --git a/docs/moduleCom4FlowPy.rst b/docs/moduleCom4FlowPy.rst index ed1b0e745..002b3ac7c 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -245,8 +245,6 @@ 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) From db91a91037f4e40c9642aa9ca0088d2dd201feab Mon Sep 17 00:00:00 2001 From: Paula Date: Mon, 18 May 2026 16:22:17 +0200 Subject: [PATCH 5/6] comments AH Update avaframe/tests/test_com4FlowPy.py Co-authored-by: qltysh[bot] <168846912+qltysh[bot]@users.noreply.github.com> fix bug --- avaframe/com4FlowPy/flowCore.py | 18 +++++++++++------- avaframe/tests/test_com4FlowPy.py | 4 ++-- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index 4652ebac8..8fe55cc21 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -287,7 +287,8 @@ def run(optTuple): travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 if forestInteraction: forestIntArray = np.ones_like(dem, dtype=np.float32) * -9999 - processedStartCellIdDict = {} + if relOutputParams["relIdBool"]: + processedStartCellIdDict = {} zDeltaList = [] fluxList = [] @@ -302,7 +303,8 @@ def run(optTuple): slTravelAngleList = [] travelLengthMaxList = [] travelLengthMinList = [] - processedStartCellIdList = [] + if relOutputParams["relIdBool"]: + processedStartCellIdList = [] if forestInteraction: forestIntList = [] @@ -322,7 +324,8 @@ def run(optTuple): fpTravelAngleMinList.append(res[9]) routFluxSumList.append(res[10]) depFluxSumList.append(res[11]) - processedStartCellIdList.append(res[12]) + if relOutputParams["relIdBool"]: + processedStartCellIdList.append(res[12]) if forestInteraction: forestIntList.append(res[13]) @@ -368,10 +371,11 @@ def run(optTuple): 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 + if relOutputParams["relIdBool"]: + 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) diff --git a/avaframe/tests/test_com4FlowPy.py b/avaframe/tests/test_com4FlowPy.py index f4fbf5be9..709478acd 100644 --- a/avaframe/tests/test_com4FlowPy.py +++ b/avaframe/tests/test_com4FlowPy.py @@ -209,12 +209,12 @@ def test_calculation(): travelLengthMin[3, 2] = 2 * np.sqrt(cellsize**2) results = flowCore.calculation(args) - assert len(results) == 16 + assert len(results) == 14 assert np.all(results[1] == flux) assert np.all(results[10] == routFluxSum) assert np.all(results[11] == depFluxSum) assert np.all(results[8] == travelLengthMin) - assert results[7] == None + assert results[7] is None def createTestRaster(pathTestFolder, rasterName): From 71dd27e17bef61615fd4a8d80ade6d49c84e7e81 Mon Sep 17 00:00:00 2001 From: Paula Date: Mon, 18 May 2026 19:03:53 +0200 Subject: [PATCH 6/6] clean and make polygon geometry valid delete temp files delete unued import --- avaframe/com4FlowPy/splitAndMerge.py | 9 +++++-- .../tests/data/testCom4/refPolygon.geojson | 8 +++--- .../testCom4/refPolygon_manipulated.geojson | 8 +++--- .../data/testCom4/refPolygon_manipulated.qmd | 27 +++++++++++++++++++ avaframe/tests/test_com4FlowPy.py | 2 +- 5 files changed, 43 insertions(+), 11 deletions(-) create mode 100644 avaframe/tests/data/testCom4/refPolygon_manipulated.qmd diff --git a/avaframe/com4FlowPy/splitAndMerge.py b/avaframe/com4FlowPy/splitAndMerge.py index 5a65100a5..17d1bc134 100644 --- a/avaframe/com4FlowPy/splitAndMerge.py +++ b/avaframe/com4FlowPy/splitAndMerge.py @@ -12,7 +12,6 @@ import numpy as np import avaframe.in2Trans.rasterUtils as IOf import shapely -import shapely.ops import geopandas as gpd # create local logger @@ -324,6 +323,7 @@ def mergeDictToPolygon(inDirPath, fName, demHeader): mergedDict = mergeDict(inDirPath, fName) cellsize = demHeader["cellsize"] pathPolygons = {} + buffer = 1e-10 for (row, col), praIds in mergedDict.items(): # get a polygon around every cell contained in mergedDict @@ -342,7 +342,12 @@ def mergeDictToPolygon(inDirPath, fName, demHeader): for pid, polys in pathPolygons.items(): # merge all polygons belonging to a PRA ID - pathPolygons[pid] = shapely.ops.unary_union(polys) + unionPoly = shapely.union_all(polys) + bufferedPoly = unionPoly.buffer(buffer, cap_style="flat").buffer(-buffer, cap_style="square") + bufferedPoly = shapely.make_valid(bufferedPoly) + cleanedPoly = shapely.set_precision(bufferedPoly, 1e-3) + cleanedPoly = shapely.remove_repeated_points(cleanedPoly, tolerance=1e-3) + pathPolygons[pid] = cleanedPoly gdfPathPolygons = gpd.GeoDataFrame( {"PRA_id": list(pathPolygons.keys())}, diff --git a/avaframe/tests/data/testCom4/refPolygon.geojson b/avaframe/tests/data/testCom4/refPolygon.geojson index bbfab2e48..9a6ed08aa 100644 --- a/avaframe/tests/data/testCom4/refPolygon.geojson +++ b/avaframe/tests/data/testCom4/refPolygon.geojson @@ -1,10 +1,10 @@ { "type": "FeatureCollection", -"name": "test", +"name": "refPolygon", "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }, "features": [ -{ "type": "Feature", "properties": { "PRA_id": 3.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 5.0, -5.0 ], [ 5.0, 5.0 ], [ -5.0, 5.0 ], [ -5.0, -5.0 ], [ 5.0, -5.0 ] ] ] } }, -{ "type": "Feature", "properties": { "PRA_id": 1.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 15.0, 35.0 ], [ 15.0, 45.0 ], [ 25.0, 45.0 ], [ 35.0, 45.0 ], [ 45.0, 45.0 ], [ 45.0, 35.0 ], [ 45.0, 25.0 ], [ 45.0, 15.0 ], [ 35.0, 15.0 ], [ 25.0, 15.0 ], [ 15.0, 15.0 ], [ 15.0, 25.0 ], [ 15.0, 35.0 ] ] ] } }, -{ "type": "Feature", "properties": { "PRA_id": 2.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 45.0, 65.0 ], [ 45.0, 75.0 ], [ 45.0, 85.0 ], [ 55.0, 85.0 ], [ 65.0, 85.0 ], [ 75.0, 85.0 ], [ 85.0, 85.0 ], [ 85.0, 75.0 ], [ 85.0, 65.0 ], [ 85.0, 55.0 ], [ 75.0, 55.0 ], [ 75.0, 45.0 ], [ 65.0, 45.0 ], [ 55.0, 45.0 ], [ 45.0, 45.0 ], [ 45.0, 55.0 ], [ 45.0, 65.0 ] ] ] } } +{ "type": "Feature", "properties": { "PRA_id": 3.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -5.0, 5.0 ], [ 5.0, 5.0 ], [ 5.0, -5.0 ], [ -5.0, -5.0 ], [ -5.0, 5.0 ] ] ] } }, +{ "type": "Feature", "properties": { "PRA_id": 1.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 45.0, 45.0 ], [ 45.0, 15.0 ], [ 15.0, 15.0 ], [ 15.0, 45.0 ], [ 45.0, 45.0 ] ] ] } }, +{ "type": "Feature", "properties": { "PRA_id": 2.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 85.0, 85.0 ], [ 85.0, 55.0 ], [ 75.0, 55.0 ], [ 75.0, 45.0 ], [ 45.0, 45.0 ], [ 45.0, 85.0 ], [ 85.0, 85.0 ] ] ] } } ] } diff --git a/avaframe/tests/data/testCom4/refPolygon_manipulated.geojson b/avaframe/tests/data/testCom4/refPolygon_manipulated.geojson index caca8e148..f116d62ff 100644 --- a/avaframe/tests/data/testCom4/refPolygon_manipulated.geojson +++ b/avaframe/tests/data/testCom4/refPolygon_manipulated.geojson @@ -1,10 +1,10 @@ { "type": "FeatureCollection", -"name": "test", +"name": "refPolygon_manipulated", "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }, "features": [ -{ "type": "Feature", "properties": { "PRA_id": 3.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 5.0, -5.0 ], [ 5.0, 5.0 ], [ -5.0, 5.0 ], [ -5.0, -5.0 ], [ 5.0, -5.0 ] ] ] } }, -{ "type": "Feature", "properties": { "PRA_id": 1.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 15.0, 35.0 ], [ 15.0, 45.0 ], [ 25.0, 45.0 ], [ 25.0, 55.0 ], [ 35.0, 55.0 ], [ 45.0, 55.0 ], [ 45.0, 45.0 ], [ 45.0, 35.0 ], [ 45.0, 25.0 ], [ 45.0, 15.0 ], [ 35.0, 15.0 ], [ 25.0, 15.0 ], [ 15.0, 15.0 ], [ 15.0, 25.0 ], [ 15.0, 35.0 ] ] ] } }, -{ "type": "Feature", "properties": { "PRA_id": 2.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 45.0, 65.0 ], [ 45.0, 75.0 ], [ 45.0, 85.0 ], [ 55.0, 85.0 ], [ 65.0, 85.0 ], [ 75.0, 85.0 ], [ 85.0, 85.0 ], [ 85.0, 75.0 ], [ 85.0, 65.0 ], [ 85.0, 55.0 ], [ 75.0, 55.0 ], [ 75.0, 45.0 ], [ 65.0, 45.0 ], [ 55.0, 45.0 ], [ 45.0, 45.0 ], [ 35.0, 45.0 ], [ 25.0, 45.0 ], [ 25.0, 55.0 ], [ 35.0, 55.0 ], [ 45.0, 55.0 ], [ 45.0, 65.0 ] ] ] } } +{ "type": "Feature", "properties": { "PRA_id": 3.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -5.0, 5.0 ], [ 5.0, 5.0 ], [ 5.0, -5.0 ], [ -5.0, -5.0 ], [ -5.0, 5.0 ] ] ] } }, +{ "type": "Feature", "properties": { "PRA_id": 1.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 25.0, 45.0 ], [ 25.0, 55.0 ], [ 45.0, 55.0 ], [ 45.0, 15.0 ], [ 15.0, 15.0 ], [ 15.0, 45.0 ], [ 25.0, 45.0 ] ] ] } }, +{ "type": "Feature", "properties": { "PRA_id": 2.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 45.0, 85.0 ], [ 85.0, 85.0 ], [ 85.0, 55.0 ], [ 75.0, 55.0 ], [ 75.0, 45.0 ], [ 25.0, 45.0 ], [ 25.0, 55.0 ], [ 45.0, 55.0 ], [ 45.0, 85.0 ] ] ] } } ] } diff --git a/avaframe/tests/data/testCom4/refPolygon_manipulated.qmd b/avaframe/tests/data/testCom4/refPolygon_manipulated.qmd new file mode 100644 index 000000000..2f2471ab8 --- /dev/null +++ b/avaframe/tests/data/testCom4/refPolygon_manipulated.qmd @@ -0,0 +1,27 @@ + + + + + + dataset + + + + + + + + + GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]] + +proj=longlat +datum=WGS84 +no_defs + 3452 + 4326 + EPSG:4326 + WGS 84 + longlat + EPSG:7030 + true + + + + diff --git a/avaframe/tests/test_com4FlowPy.py b/avaframe/tests/test_com4FlowPy.py index 709478acd..61e36fb71 100644 --- a/avaframe/tests/test_com4FlowPy.py +++ b/avaframe/tests/test_com4FlowPy.py @@ -487,7 +487,7 @@ def test_mergeDictToPolygon(): assert len(gdfPathPolygons) == 3 assert np.all(gdfPathPolygons["PRA_id"] == refPolygons["PRA_id"]) - assert np.all(gdfPathPolygons == refPolygons) + assert np.all(gdfPathPolygons.geometry.geom_equals(refPolygons.geometry)) # manipulate dictionaries that the polygons overlap