diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index b448049ed..6170631ba 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -158,6 +158,18 @@ def com4FlowPyMain(cfgPath, cfgSetup): else: modelPaths["varExponentPath"] = "" + if "relIdPolygon" in modelPaths["outputFileList"] or "relIdCount" in modelPaths["outputFileList"]: + modelPaths["relIdPath"] = cfgPath["relIdPath"] + modelParameters["outputRelIdBool"] = True + else: + modelPaths["relIdPath"] = "" + modelParameters["outputRelIdBool"] = False + + if "relVolMin" in modelPaths["outputFileList"] or "relVolMax" in modelPaths["outputFileList"]: + modelParameters["outputRelVolBool"] = True + else: + modelParameters["outputRelVolBool"] = False + # TODO: provide some kind of check for the model Parameters # i.e. * sensible value ranges # * contradicting options ... @@ -452,6 +464,8 @@ def tileInputLayers(modelParameters, modelPaths, rasterAttributes, tilingParamet SPAM.tileRaster(modelPaths["varExponentPath"], "varExponent", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) if modelParameters["forestBool"]: SPAM.tileRaster(modelPaths["forestPath"], "forest", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) + if modelParameters["outputRelIdBool"]: + SPAM.tileRaster(modelPaths["relIdPath"], "relId", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) log.info("Finished Tiling All Input Rasters.") log.info("==================================") @@ -515,56 +529,40 @@ def mergeAndWriteResults(modelPaths, modelOptions): _outputs = set(modelPaths['outputFileList']) _outputNoDataValue = modelPaths["outputNoDataValue"] - log.info(" merging results ...") - log.info("-------------------------") - - # Merge calculated tiles - zDelta = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta") - flux = SPAM.mergeRaster(modelPaths["tempDir"], "res_flux") - cellCounts = SPAM.mergeRaster(modelPaths["tempDir"], "res_count", method="sum") - zDeltaSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta_sum", method="sum") - routFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_rout_flux_sum", method="sum") - depFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_dep_flux_sum", method="sum") - fpTaMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp_max") - fpTaMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp_min", method="min") - slTa = SPAM.mergeRaster(modelPaths["tempDir"], "res_sl") - travelLengthMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length_max") - travelLengthMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length_min", method="min") - - if modelOptions["infraBool"]: - backcalc = SPAM.mergeRaster(modelPaths["tempDir"], "res_backcalc") - - if modelOptions["forestInteraction"]: - forestInteraction = SPAM.mergeRaster(modelPaths["tempDir"], "res_forestInt", method='min') - - # Write Output Files to Disk - log.info("-------------------------") - log.info(" writing output files ...") - log.info("-------------------------") - _oF = modelPaths["outputFileFormat"] - _ts = modelPaths["timeString"] - demHeader = IOf.readRasterHeader(modelPaths["demPath"]) outputHeader = demHeader.copy() outputHeader["nodata_value"] = _outputNoDataValue + _oF = modelPaths["outputFileFormat"] + _ts = modelPaths["timeString"] + + useCompression = modelPaths["useCompression"] if _oF == ".asc": outputHeader["driver"] = "AAIGrid" elif _oF == ".tif": outputHeader["driver"] = "GTiff" - useCompression = modelPaths["useCompression"] + log.info(" merging and writing results ...") + log.info("-------------------------") - if 'flux' in _outputs: - flux = defineNotAffectedCells(flux, cellCounts, noDataValue=_outputNoDataValue) + # Merge calculated tiles + # compute cellCounts and don't delete because it is used for defining not affected cells + # other rasters (and polygons) are deleted after writing (to reduce computation time and used RAM) + cellCounts = SPAM.mergeRaster(modelPaths["tempDir"], "res_count", method="sum") + if 'cellCounts' in _outputs: + cellCounts = defineNotAffectedCells(cellCounts, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, - flux, - modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), + cellCounts, + modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), flip=True, useCompression=useCompression, - ) + ) + del output + log.info("com4_{}_{}_cellCounts is written".format(_uid, _ts)) + if 'zDelta' in _outputs: + zDelta = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta") zDelta = defineNotAffectedCells(zDelta, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -572,17 +570,27 @@ def mergeAndWriteResults(modelPaths, modelOptions): modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), flip=True, useCompression=useCompression, - ) - if 'cellCounts' in _outputs: - cellCounts = defineNotAffectedCells(cellCounts, cellCounts, noDataValue=_outputNoDataValue) + ) + del zDelta + del output + log.info("com4_{}_{}_zdelta is written".format(_uid, _ts)) + + if 'flux' in _outputs: + flux = SPAM.mergeRaster(modelPaths["tempDir"], "res_flux") + flux = defineNotAffectedCells(flux, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, - cellCounts, - modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), + flux, + modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), flip=True, useCompression=useCompression, - ) + ) + del flux + del output + log.info("com4_{}_{}_flux is written".format(_uid, _ts)) + if 'zDeltaSum' in _outputs: + zDeltaSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta_sum", method="sum") zDeltaSum = defineNotAffectedCells(zDeltaSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -591,7 +599,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del zDeltaSum + del output + log.info("com4_{}_{}_zDeltaSum is written".format(_uid, _ts)) + if 'routFluxSum' in _outputs: + routFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_rout_flux_sum", method="sum") routFluxSum = defineNotAffectedCells(routFluxSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -600,7 +613,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del routFluxSum + del output + log.info("com4_{}_{}_routFluxSum is written".format(_uid, _ts)) + if 'depFluxSum' in _outputs: + depFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_dep_flux_sum", method="sum") depFluxSum = defineNotAffectedCells(depFluxSum, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -609,7 +627,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del depFluxSum + del output + log.info("com4_{}_{}_depFluxSum is written".format(_uid, _ts)) + if "fpTravelAngle" in _outputs or "fpTravelAngleMax" in _outputs: + fpTaMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp_max") fpTaMax = defineNotAffectedCells(fpTaMax, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -618,7 +641,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del fpTaMax + del output + log.info("com4_{}_{}_fpTravelAngleMax is written".format(_uid, _ts)) + if "fpTravelAngleMin" in _outputs: + fpTaMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp_min", method="min") fpTaMin = defineNotAffectedCells(fpTaMin, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -627,7 +655,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del fpTaMin + del output + log.info("com4_{}_{}_fpTravelAngleMin is written".format(_uid, _ts)) + if 'slTravelAngle' in _outputs: + slTa = SPAM.mergeRaster(modelPaths["tempDir"], "res_sl") slTa = defineNotAffectedCells(slTa, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -636,7 +669,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del slTa + del output + log.info("com4_{}_{}_slTravelAngle is written".format(_uid, _ts)) + if "travelLength" in _outputs or "travelLengthMax" in _outputs: + travelLengthMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length_max") travelLengthMax = defineNotAffectedCells(travelLengthMax, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -645,7 +683,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del travelLengthMax + del output + log.info("com4_{}_{}_travelLengthMax is written".format(_uid, _ts)) + if "travelLengthMin" in _outputs: + travelLengthMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length_min", method="min") travelLengthMin = defineNotAffectedCells(travelLengthMin, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -654,12 +697,13 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) + del travelLengthMin + del output + log.info("com4_{}_{}_travelLengthMin is written".format(_uid, _ts)) - # NOTE: - # if not modelOptions["infraBool"]: # if no infra - # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("cell_counts%s" %(output_format)),cell_counts) - # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("z_delta_sum%s" %(output_format)),z_delta_sum) - if modelOptions["infraBool"]: # if infra + + if modelOptions["infraBool"]: + backcalc = SPAM.mergeRaster(modelPaths["tempDir"], "res_backcalc") backcalc = defineNotAffectedCells(backcalc, cellCounts, noDataValue=_outputNoDataValue) output = IOf.writeResultToRaster( outputHeader, @@ -667,8 +711,14 @@ def mergeAndWriteResults(modelPaths, modelOptions): modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts), flip=True, useCompression=useCompression, - ) + ) + del backcalc + del output + log.info("com4_{}_{}_backcalculation is written".format(_uid, _ts)) + + if modelOptions["forestInteraction"]: + forestInteraction = SPAM.mergeRaster(modelPaths["tempDir"], "res_forestInt", method='min') forestInteraction = defineNotAffectedCells( forestInteraction, cellCounts, noDataValue=_outputNoDataValue ) @@ -678,8 +728,69 @@ def mergeAndWriteResults(modelPaths, modelOptions): modelPaths["resDir"] / "com4_{}_{}_forestInteraction".format(_uid, _ts), flip=True, useCompression=useCompression, + ) + del forestInteraction + del output + log.info("com4_{}_{}_forestInteraction is written".format(_uid, _ts)) + + + if "relIdPolygon" in _outputs: + pathPolygons = SPAM.mergeDictToPolygon(modelPaths["tempDir"], "res_startCellIdDict", outputHeader) + pathPolygons.to_file( + modelPaths["resDir"] / "com4_{}_{}_pathPolygons.geojson".format(_uid, _ts), driver="GeoJSON" + ) + del pathPolygons + log.info("com4_{}_{}_pathPolygons is written".format(_uid, _ts)) + + + if "relIdCount" in _outputs: + countRelId = SPAM.mergeDictToRaster(modelPaths["tempDir"], "res_startCellIdDict") + countRelId = defineNotAffectedCells(countRelId, cellCounts, noDataValue=_outputNoDataValue) + output = IOf.writeResultToRaster( + outputHeader, + countRelId, + modelPaths["resDir"] / "com4_{}_{}_countRelId".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) + del countRelId + del output + log.info("com4_{}_{}_countRelId is written".format(_uid, _ts)) + + + if "relVolMin" in _outputs: + relVolMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_relVol_min", method="min") + relVolMin = defineNotAffectedCells(relVolMin, cellCounts, noDataValue=_outputNoDataValue) + output = IOf.writeResultToRaster( + outputHeader, + relVolMin, + modelPaths["resDir"] / "com4_{}_{}_relVolMin".format(_uid, _ts), + flip=True, + useCompression=useCompression, ) - del output + del relVolMin + del output + log.info("com4_{}_{}_relVolMin is written".format(_uid, _ts)) + + + if "relVolMax" in _outputs: + relVolMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_relVol_max") + relVolMax = defineNotAffectedCells(relVolMax, cellCounts, noDataValue=_outputNoDataValue) + output = IOf.writeResultToRaster( + outputHeader, + relVolMax, + modelPaths["resDir"] / "com4_{}_{}_relVolMax".format(_uid, _ts), + flip=True, + useCompression=useCompression, + ) + del relVolMax + del output + log.info("com4_{}_{}_relVolMax is written".format(_uid, _ts)) + + # NOTE: + # if not modelOptions["infraBool"]: # if no infra + # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("cell_counts%s" %(output_format)),cell_counts) + # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("z_delta_sum%s" %(output_format)),z_delta_sum) def checkConvertReleaseShp2Tif(modelPaths): diff --git a/avaframe/com4FlowPy/com4FlowPyCfg.ini b/avaframe/com4FlowPy/com4FlowPyCfg.ini index c4920e7d8..a56d56b39 100644 --- a/avaframe/com4FlowPy/com4FlowPyCfg.ini +++ b/avaframe/com4FlowPy/com4FlowPyCfg.ini @@ -227,8 +227,14 @@ outputNoDataValue = -9999 # depFluxSum # travelLengthMin # fpTravelAngleMin +# relVolMin +# relVolMax +# relIdPolygon +# relIdCount # if forestInteraction: forestInteraction is automatically added to outputs # if infra: backCalculation is automatically added to output +# if relVolMin or relVolMax is in outputFiles, the Volume of the PRA should be provided in the raster file in the REL folder +# if relIdCount or relIdPolygon is in outputFiles, the ids of the PRAs should be provided in the raster file in the RELID folder outputFiles = zDelta|cellCounts|travelLengthMax|fpTravelAngleMax #++++++++++++ Custom paths True/False @@ -256,6 +262,7 @@ useCustomPathDEM = False workDir = demPath = releasePath = +releaseIdPath = infraPath = forestPath = varUmaxPath = diff --git a/avaframe/com4FlowPy/flowClass.py b/avaframe/com4FlowPy/flowClass.py index 94cfa2dd6..c4e21e3e7 100644 --- a/avaframe/com4FlowPy/flowClass.py +++ b/avaframe/com4FlowPy/flowClass.py @@ -10,14 +10,25 @@ class Cell: + def __init__( self, - rowindex, colindex, - dem_ng, cellsize, - flux, z_delta, parent, - alpha, exp, flux_threshold, max_z_delta, - startcell, fluxDistOldVersionBool=False, - FSI=None, forestParams=None, + rowindex, + colindex, + dem_ng, + cellsize, + flux, + z_delta, + parent, + alpha, + exp, + flux_threshold, + max_z_delta, + startcell, + fluxDistOldVersionBool=False, + FSI=None, + forestParams=None, + startcellVol=None, ): """constructor for the Cell class that describes a raster cell that is hit by the GMF. the constructor function is called every time a new instance of type 'Cell' is @@ -62,6 +73,9 @@ def __init__( self._SQRT2 = np.sqrt(2.0) self._RAD90 = np.deg2rad(90.0) + self.startcellVolMin = startcellVol + self.startcellVolMax = startcellVol + # NOTE: Forest Interaction included here # if FSI != None AND forestParams != None - then self.ForestBool = True and forestParams and # FSI are accordingly initialized @@ -177,6 +191,10 @@ def add_parent(self, parent): if parent.forestIntCount < (self.forestIntCount - self.isForest): self.forestIntCount = parent.forestIntCount + self.isForest + def calc_startCellVol(self, startcellVolNew): + self.startcellVolMin = min(self.startcellVolMin, startcellVolNew) + self.startcellVolMax = max(self.startcellVolMax, startcellVolNew) + def calcDistMin(self, calc3D=False): """ function calculates the projected horizontal (self.min_distance) and 3D (self.minDistXYZ) length @@ -416,7 +434,7 @@ def calc_distribution(self): self.dist[self.dist < threshold] = 0 if np.sum(self.dist) != self.flux and count > 0: # correction/flux conservation for potential rounding losses or gains - # (self.flux - np.sum(self.dist)) will either be negative or positive + # (self.flux - np.sum(self.dist)) will either be negative or positive # depending on the direction of the rounding error self.dist[self.dist >= threshold] += (self.flux - np.sum(self.dist)) / count if count == 0: @@ -446,4 +464,4 @@ def forest_detrainment(self): # rise over run (should be negative slope) slope = (_rest - self.minAddedDetrainmentForest) / (0 - _noDetrainmentEffectZdelta) # y=mx+b, where zDelta is x - self.detrainment = max(self.minAddedDetrainmentForest, slope * self.z_delta + _rest) \ No newline at end of file + self.detrainment = max(self.minAddedDetrainmentForest, slope * self.z_delta + _rest) diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index a9a2315f7..822d9bc23 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -13,6 +13,7 @@ import gc import psutil import time +import pickle if os.name == "nt": from multiprocessing.pool import Pool as Pool @@ -146,6 +147,8 @@ def run(optTuple): varAlphaBool = optTuple[2]["varAlphaBool"] varExponentBool = optTuple[2]["varExponentBool"] fluxDistOldVersionBool = optTuple[2]["fluxDistOldVersionBool"] + relIdBool = optTuple[2]["outputRelIdBool"] + relVolBool = optTuple[2]["outputRelVolBool"] previewMode = optTuple[2]["previewMode"] # Temp-Dir (all input files are located here and results are written back in here) @@ -199,6 +202,16 @@ def run(optTuple): else: varExponentArray = None + if relIdBool: + relIdArray = np.load(tempDir / ("relId_%s_%s.npy" % (optTuple[0], optTuple[1]))) + else: + relIdArray = None + + if relVolBool: + relVolArray = release.copy() + else: + relVolArray = None + varParams = { 'varUmaxBool': varUmaxBool, 'varUmaxArray': varUmaxArray, @@ -207,6 +220,12 @@ def run(optTuple): 'varExponentBool': varExponentBool, 'varExponentArray': varExponentArray, } + relOutputParams = { + "relIdBool": relIdBool, + "relIdArray": relIdArray, + "relVolBool": relVolBool, + "relVolArray": relVolArray, + } # convert release areas to binary (0: no release areas, 1: release areas) # every positive value >0 is interpreted as release area @@ -250,6 +269,7 @@ def run(optTuple): forestArray, forestParams, outputs, + relOutputParams ] for release_sub in release_list ], @@ -275,6 +295,9 @@ def run(optTuple): travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 if forestInteraction: forestIntArray = np.ones_like(dem, dtype=np.float32) * -9999 + relVolMinArray = np.ones_like(dem, dtype=np.float32) * -9999 + relVolMaxArray = np.zeros_like(dem, dtype=np.float32) + processedStartCellIdDict = {} zDeltaList = [] fluxList = [] @@ -289,6 +312,9 @@ def run(optTuple): slTravelAngleList = [] travelLengthMaxList = [] travelLengthMinList = [] + processedStartCellIdList = [] + relVolMinList = [] + relVolMaxList = [] if forestInteraction: forestIntList = [] @@ -308,8 +334,11 @@ def run(optTuple): fpTravelAngleMinList.append(res[9]) routFluxSumList.append(res[10]) depFluxSumList.append(res[11]) + processedStartCellIdList.append(res[12]) if forestInteraction: - forestIntList.append(res[12]) + forestIntList.append(res[15]) + relVolMinList.append(res[13]) + relVolMaxList.append(res[14]) logging.info("Calculation finished, getting results.") for i in range(len(zDeltaList)): @@ -342,7 +371,26 @@ def run(optTuple): forestIntArray = np.where((forestIntArray >= 0) & (forestIntList[i] >= 0), np.minimum(forestIntArray, forestIntList[i]), np.maximum(forestIntArray, forestIntList[i])) + if "relVolMin" in outputs: + relVolMinArray = np.where( + (relVolMinArray >= 0) & (relVolMinList[i] >= 0), + np.minimum(relVolMinArray, relVolMinList[i]), + np.maximum(relVolMinArray, relVolMinList[i]), + ) + if "relVolMax" in outputs: + relVolMaxArray = np.maximum(relVolMaxArray, relVolMaxList[i]) + if "relIdPolygon" in outputs or "relIdCount" in outputs: + for key in processedStartCellIdList[i]: + if key in processedStartCellIdDict: + ids = np.append(processedStartCellIdList[i][key], processedStartCellIdDict[key]) + processedStartCellIdDict[key] = np.unique(ids) + else: + processedStartCellIdDict[key] = processedStartCellIdList[i][key] + saveDict = open(tempDir / ("res_startCellIdDict_%s_%s.pickle" % (optTuple[0], optTuple[1])), "wb") + pickle.dump(processedStartCellIdDict, saveDict) + saveDict.close() + del processedStartCellIdDict # Save Calculated tiles np.save(tempDir / ("res_z_delta_%s_%s" % (optTuple[0], optTuple[1])), zDeltaArray) np.save(tempDir / ("res_z_delta_sum_%s_%s" % (optTuple[0], optTuple[1])), zDeltaSumArray) @@ -355,6 +403,8 @@ def run(optTuple): np.save(tempDir / ("res_sl_%s_%s" % (optTuple[0], optTuple[1])), slTravelAngleArray) np.save(tempDir / ("res_travel_length_max_%s_%s" % (optTuple[0], optTuple[1])), travelLengthMaxArray) np.save(tempDir / ("res_travel_length_min_%s_%s" % (optTuple[0], optTuple[1])), travelLengthMinArray) + np.save(tempDir / ("res_relVol_max_%s_%s" % (optTuple[0], optTuple[1])), relVolMaxArray) + np.save(tempDir / ("res_relVol_min_%s_%s" % (optTuple[0], optTuple[1])), relVolMinArray) if infraBool: np.save(tempDir / ("res_backcalc_%s_%s" % (optTuple[0], optTuple[1])), backcalc) if forestInteraction: @@ -388,6 +438,7 @@ def calculation(args): - args[14] (numpy array) - contains forest information (None if forestBool=False) - args[15] (dict) - contains parameters for forest interaction models (None if forestBool=False) - args[16] (list) - output names + - args[17] (dict) - contains flags and rasters for release - information outputs Returns ----------- @@ -455,6 +506,10 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): fluxDistOldVersionBool = args[12] previewMode = args[13] outputs = args[16] + relIdArray = args[17]["relIdArray"] + relVolArray = args[17]["relVolArray"] + relIdBool = args[17]["relIdBool"] + relVolBool = args[17]["relVolBool"] if forestBool: forestArray = args[14] @@ -480,6 +535,9 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 travelLengthMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 + relVolMinArray = np.ones_like(dem, dtype=np.float32) * -9999 + relVolMaxArray = np.zeros_like(dem, dtype=np.float32) + if infraBool: backcalc = np.ones_like(dem, dtype=np.int32) * -9999 else: @@ -492,12 +550,15 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): if forestInteraction: forestIntArray = np.ones_like(dem, dtype=np.float32) * -9999 + else: + forestIntArray = None # Core # NOTE-TODO: row_list, col_list are tuples - rethink variable naming row_list, col_list = get_start_idx(dem, release) startcell_idx = 0 + startCellIdDict = {} while startcell_idx < len(row_list): if infraBool: @@ -525,18 +586,37 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): startcell_idx += 1 continue + if relIdBool: + startcellId = relIdArray[row_idx, col_idx] + else: + startcellId = None + if relVolBool: + startcellVol = relVolArray[row_idx, col_idx] + else: + startcellVol = None + startcell = Cell( - row_idx, col_idx, - dem_ng, cellsize, - 1, 0, None, - alpha, exp, flux_threshold, max_z_delta, - startcell=True, fluxDistOldVersionBool=fluxDistOldVersionBool, + row_idx, + col_idx, + dem_ng, + cellsize, + 1, + 0, + None, + alpha, + exp, + flux_threshold, + max_z_delta, + startcell=True, + fluxDistOldVersionBool=fluxDistOldVersionBool, FSI=forestArray[row_idx, col_idx] if isinstance(forestArray, np.ndarray) else None, forestParams=forestParams, + startcellVol=startcellVol, ) # dictionary of all the cells that have been processed and the number of times the cell has been visited processedCells[(startcell.rowindex, startcell.colindex)] = 1 + # list of flowClass.Cell() Objects that is contains the "path" for each release-cell cell_list.append(startcell) @@ -545,6 +625,14 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): updateInfraDirGraph(startcell.rowindex, startcell.colindex) for idx, cell in enumerate(cell_list): + if relIdBool: + if (cell.rowindex, cell.colindex) in startCellIdDict: + startcellIdList = np.append( + startCellIdDict[(cell.rowindex, cell.colindex)], startcellId) + startCellIdDict[(cell.rowindex, cell.colindex)] = np.unique( + startcellIdList) + else: + startCellIdDict[(cell.rowindex, cell.colindex)] = np.array([startcellId]) # calculate flux, z_delta from current cell (cell) to child-cells # lenght of row, col, flux, and z_delta vectors correspond to @@ -566,6 +654,8 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): if row[k] == cell_list[i].rowindex and col[k] == cell_list[i].colindex: cell_list[i].add_os(flux[k]) cell_list[i].add_parent(cell) + if relVolBool: + cell_list[i].calc_startCellVol(startcellVol) if infraBool: updateInfraDirGraph(row[k], col[k], cell.rowindex, cell.colindex) @@ -599,16 +689,26 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): else: processedCells[(row[k], col[k])] = 1 - cell_list.append(Cell( - row[k], col[k], - dem_ng, cellsize, - flux[k], z_delta[k], - cell, - alpha, exp, flux_threshold, max_z_delta, - startcell, fluxDistOldVersionBool=fluxDistOldVersionBool, - FSI=forestArray[row[k], col[k]] if isinstance(forestArray, np.ndarray) else None, - forestParams=forestParams, - )) + cell_list.append( + Cell( + row[k], + col[k], + dem_ng, + cellsize, + flux[k], + z_delta[k], + cell, + alpha, + exp, + flux_threshold, + max_z_delta, + startcell, + fluxDistOldVersionBool=fluxDistOldVersionBool, + FSI=forestArray[row[k], col[k]] if isinstance(forestArray, np.ndarray) else None, + forestParams=forestParams, + startcellVol=startcellVol, + ) + ) zDeltaArray[cell.rowindex, cell.colindex] = max(zDeltaArray[cell.rowindex, cell.colindex], cell.z_delta) fluxArray[cell.rowindex, cell.colindex] = max(fluxArray[cell.rowindex, cell.colindex], cell.flux) @@ -643,6 +743,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): travelLengthMinArray[cell.rowindex, cell.colindex] = max( travelLengthMinArray[cell.rowindex, cell.colindex], cell.min_distance ) + if processedCells[(cell.rowindex, cell.colindex)] == 1: countArray[cell.rowindex, cell.colindex] += int(1) @@ -653,6 +754,19 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): else: forestIntArray[cell.rowindex, cell.colindex] = max(forestIntArray[cell.rowindex, cell.colindex], cell.forestIntCount) + if "relVolMax" in outputs: + relVolMaxArray[cell.rowindex, cell.colindex] = max( + relVolMaxArray[cell.rowindex, cell.colindex], cell.startcellVolMax + ) + if "relVolMin" in outputs: + if relVolMinArray[cell.rowindex, cell.colindex] >= 0 and cell.startcellVolMin >= 0: + relVolMinArray[cell.rowindex, cell.colindex] = min( + relVolMinArray[cell.rowindex, cell.colindex], cell.startcellVolMin + ) + else: + relVolMinArray[cell.rowindex, cell.colindex] = max( + relVolMinArray[cell.rowindex, cell.colindex], cell.startcellVolMin + ) if infraBool: # if 'infraBool' is True - i.e. calculation is performed with infrastructure information @@ -683,38 +797,24 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): zDeltaSumArray += zDeltaPathArray gc.collect() - - if forestInteraction: - return ( - zDeltaArray, - fluxArray, - countArray, - zDeltaSumArray, - backcalc, - fpTravelAngleMaxArray, - slTravelAngleArray, - travelLengthMaxArray, - travelLengthMinArray, - fpTravelAngleMinArray, - routFluxSumArray, - depFluxSumArray, - forestIntArray, - ) - else: - return ( - zDeltaArray, - fluxArray, - countArray, - zDeltaSumArray, - backcalc, - fpTravelAngleMaxArray, - slTravelAngleArray, - travelLengthMaxArray, - travelLengthMinArray, - fpTravelAngleMinArray, - routFluxSumArray, - depFluxSumArray, - ) + return ( + zDeltaArray, + fluxArray, + countArray, + zDeltaSumArray, + backcalc, + fpTravelAngleMaxArray, + slTravelAngleArray, + travelLengthMaxArray, + travelLengthMinArray, + fpTravelAngleMinArray, + routFluxSumArray, + depFluxSumArray, + startCellIdDict, + relVolMinArray, + relVolMaxArray, + forestIntArray, + ) def enoughMemoryAvailable(limit=0.05): diff --git a/avaframe/com4FlowPy/splitAndMerge.py b/avaframe/com4FlowPy/splitAndMerge.py index a4a82c246..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 c9f45570e..52a35bc49 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -191,6 +191,7 @@ the following folder structure inside the ``avalancheDir`` directory inside whic Inputs/ ElevationModel - digital elevation model (.asc) REL/ - release area file (can be either .asc, .tif, or .shp) + RELID/ - release ID file (can be either .asc, .tif) RES/ - forest structure information (FSI) (.asc or .tif) INFRA/ - infrastructure layer (.asc or .tif) ALPHA/ - variable alpha angle layer (.tif) @@ -206,6 +207,7 @@ inputs and working directories/model outputs in different places, which might be - ``workDir`` :math:`\ldots` working directory (a ``temp/`` folder, model log and model results will be written here) - ``demPath`` :math:`\ldots` path to input DEM (must be ``.asc`` currently) - ``releasePath`` :math:`\ldots` path to release area raster (``.asc, .tif``) +- ``releaseIdPath`` :math:`\ldots` path to release ID raster (``.asc, .tif``) - ``infraPath`` :math:`\ldots` path to infrastructure raster (``.asc, .tif``) (required if ``infra = True``) - ``forestPath`` :math:`\ldots` path to forest (FSI) raster (``.asc, .tif``) (required if ``forest = True``) - ``varAlphaPath`` :math:`\ldots` path to variable alpha angle raster (``.tif``) (required if ``variableAlpha = True``) @@ -243,6 +245,10 @@ In addition these output layers are also available: - ``depFluxSum``: deposited flux summed up over all paths - ``travelLengthMin``: the travel length along the flow path (the minimum value of all paths for every raster cell) - ``fpTravelAngleMin``: the gamma angle along the flow path (the minimum value of all paths for every raster cell) +- ``relVolMin``: the minimum of the tracked release volume that route through the raster cell (the volume is provided in the release file in the REL folder) +- ``relVolMax``: the maximum of the tracked release volume that route through the raster cell (the volume is provided in the release file in the REL folder) +- ``relIdPolygon``: polygons (*.geojson) that cover the affected process belonging to one release ID (can include multiple release cells; release IDs are provided in the RELID folder) +- ``relIdCount``: number of paths belonging to different release IDs that route flux through a raster cell (release IDs are provided in the RELID folder) If ``forestInteraction = True`` this layer will be written automatically (no need to separately define in ``outputFiles``):