diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index b448049ed..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 @@ -158,6 +179,13 @@ 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 + # TODO: provide some kind of check for the model Parameters # i.e. * sensible value ranges # * contradicting options ... @@ -307,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!") @@ -315,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!") @@ -323,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!") @@ -331,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!") @@ -340,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) @@ -357,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) @@ -376,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]),\ @@ -385,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: @@ -440,18 +484,28 @@ 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"]: + SPAM.tileRaster(modelPaths["relIdPath"], "relId", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) log.info("Finished Tiling All Input Rasters.") log.info("==================================") @@ -512,59 +566,43 @@ 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"] - 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, ) - if 'zDelta' in _outputs: + 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, @@ -573,16 +611,26 @@ def mergeAndWriteResults(modelPaths, modelOptions): 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, ) - if 'zDeltaSum' in _outputs: + 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 +639,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) - if 'routFluxSum' in _outputs: + 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 +653,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) - if 'depFluxSum' in _outputs: + 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 +667,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 +681,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 +695,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) - if 'slTravelAngle' in _outputs: + 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 +709,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 +723,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 +737,12 @@ 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, @@ -668,7 +751,12 @@ def mergeAndWriteResults(modelPaths, modelOptions): 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 ) @@ -679,7 +767,36 @@ def mergeAndWriteResults(modelPaths, modelOptions): flip=True, useCompression=useCompression, ) - del output + 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)) + + # 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): @@ -701,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": @@ -730,11 +849,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 ----------- @@ -755,9 +874,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 @@ -772,7 +891,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("+++++++++++++++++++++++") diff --git a/avaframe/com4FlowPy/com4FlowPyCfg.ini b/avaframe/com4FlowPy/com4FlowPyCfg.ini index c4920e7d8..97e248553 100644 --- a/avaframe/com4FlowPy/com4FlowPyCfg.ini +++ b/avaframe/com4FlowPy/com4FlowPyCfg.ini @@ -227,8 +227,12 @@ outputNoDataValue = -9999 # depFluxSum # travelLengthMin # fpTravelAngleMin +# 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 +260,7 @@ useCustomPathDEM = False workDir = demPath = releasePath = +relIdPath = infraPath = forestPath = varUmaxPath = diff --git a/avaframe/com4FlowPy/flowClass.py b/avaframe/com4FlowPy/flowClass.py index 94cfa2dd6..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 @@ -10,14 +10,24 @@ 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, ): """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 @@ -26,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 @@ -49,13 +63,17 @@ 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 @@ -103,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 @@ -195,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: @@ -237,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) @@ -252,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 @@ -266,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 @@ -284,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 @@ -312,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 @@ -416,7 +439,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: @@ -442,8 +465,8 @@ 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 - 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..8fe55cc21 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,7 @@ def run(optTuple): varAlphaBool = optTuple[2]["varAlphaBool"] varExponentBool = optTuple[2]["varExponentBool"] fluxDistOldVersionBool = optTuple[2]["fluxDistOldVersionBool"] + relIdBool = optTuple[2]["outputRelIdBool"] previewMode = optTuple[2]["previewMode"] # Temp-Dir (all input files are located here and results are written back in here) @@ -196,6 +198,11 @@ def run(optTuple): else: varExponentArray = None + if relIdBool: + relIdArray = np.load(tempDir / ("relId_%s_%s.npy" % (optTuple[0], optTuple[1]))) + else: + relIdArray = None + varParams = { "varUmaxBool": varUmaxBool, "varUmaxArray": varUmaxArray, @@ -204,6 +211,10 @@ def run(optTuple): "varExponentBool": varExponentBool, "varExponentArray": varExponentArray, } + relOutputParams = { + "relIdBool": relIdBool, + "relIdArray": relIdArray, + } # convert release areas to binary (0: no release areas, 1: release areas) # every positive value >0 is interpreted as release area @@ -250,6 +261,7 @@ def run(optTuple): forestArray, forestParams, outputs, + relOutputParams, ] for release_sub in release_list ], @@ -275,6 +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 + if relOutputParams["relIdBool"]: + processedStartCellIdDict = {} zDeltaList = [] fluxList = [] @@ -289,6 +303,8 @@ def run(optTuple): slTravelAngleList = [] travelLengthMaxList = [] travelLengthMinList = [] + if relOutputParams["relIdBool"]: + processedStartCellIdList = [] if forestInteraction: forestIntList = [] @@ -308,8 +324,10 @@ def run(optTuple): fpTravelAngleMinList.append(res[9]) routFluxSumList.append(res[10]) depFluxSumList.append(res[11]) + if relOutputParams["relIdBool"]: + processedStartCellIdList.append(res[12]) if forestInteraction: - forestIntList.append(res[12]) + forestIntList.append(res[13]) logging.info("Calculation finished, getting results.") for i in range(len(zDeltaList)): @@ -345,6 +363,19 @@ def run(optTuple): np.maximum(forestIntArray, forestIntList[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] + + 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) @@ -390,6 +421,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 +489,8 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): fluxDistOldVersionBool = args[12] previewMode = args[13] outputs = args[16] + relIdArray = args[17]["relIdArray"] + relIdBool = args[17]["relIdBool"] if forestBool: forestArray = args[14] @@ -466,6 +500,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) @@ -474,32 +512,44 @@ 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 + + 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 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: + 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 +577,11 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): startcell_idx += 1 continue + if relIdBool: + startcellId = relIdArray[row_idx, col_idx] + else: + startcellId = None + startcell = Cell( row_idx, col_idx, @@ -547,6 +602,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): # 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 +611,12 @@ 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 @@ -667,6 +729,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) @@ -713,38 +776,22 @@ 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, + forestIntArray, + ) def enoughMemoryAvailable(limit=0.05): diff --git a/avaframe/com4FlowPy/splitAndMerge.py b/avaframe/com4FlowPy/splitAndMerge.py index a4a82c246..17d1bc134 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 @@ -11,6 +11,8 @@ import gc import numpy as np import avaframe.in2Trans.rasterUtils as IOf +import shapely +import geopandas as gpd # create local logger log = logging.getLogger(__name__) @@ -22,11 +24,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) @@ -176,14 +178,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 @@ -207,7 +209,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): @@ -217,19 +219,140 @@ 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) return mergedRas + + +def mergeDict(inDirPath, fName): + """ + Merges the dictionary-results for each tile to one dictionary + + Parameters + ---------- + 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 + + 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 = {} + buffer = 1e-10 + + 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 + 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())}, + geometry=list(pathPolygons.values()), + crs=demHeader["crs"], + ) + del pathPolygons + return gdfPathPolygons 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/runCom4FlowPy.py b/avaframe/runCom4FlowPy.py index c8af4737a..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) @@ -163,6 +171,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 +238,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) @@ -269,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) @@ -291,6 +302,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) @@ -324,18 +349,38 @@ def checkOutputFilesFormat(strOutputFiles): """ try: - setA = set(strOutputFiles.split('|')) - setB = set(['zDelta', 'cellCounts', 'fpTravelAngle', 'travelLength', 'fpTravelAngleMax', - 'fpTravelAngleMin', 'travelLengthMax', 'travelLengthMin', - 'slTravelAngle', 'flux', 'zDeltaSum', 'routFluxSum', 'depFluxSum']) + setA = set(strOutputFiles.split("|")) + setB = set( + [ + "zDelta", + "cellCounts", + "fpTravelAngle", + "travelLength", + "fpTravelAngleMax", + "fpTravelAngleMin", + "travelLengthMax", + "travelLengthMin", + "slTravelAngle", + "flux", + "zDeltaSum", + "routFluxSum", + "depFluxSum", + "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): @@ -362,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/avaframe/tests/data/testCom4/refPolygon.geojson b/avaframe/tests/data/testCom4/refPolygon.geojson new file mode 100644 index 000000000..9a6ed08aa --- /dev/null +++ b/avaframe/tests/data/testCom4/refPolygon.geojson @@ -0,0 +1,10 @@ +{ +"type": "FeatureCollection", +"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": [ [ [ 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 new file mode 100644 index 000000000..f116d62ff --- /dev/null +++ b/avaframe/tests/data/testCom4/refPolygon_manipulated.geojson @@ -0,0 +1,10 @@ +{ +"type": "FeatureCollection", +"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": [ [ [ 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/data/testCom4/testRaster.tif b/avaframe/tests/data/testCom4/testRaster.tif new file mode 100644 index 000000000..a47562f83 Binary files /dev/null and b/avaframe/tests/data/testCom4/testRaster.tif differ diff --git a/avaframe/tests/test_com4FlowPy.py b/avaframe/tests/test_com4FlowPy.py index 8c4321149..61e36fb71 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,32 +169,354 @@ def test_calculation(): } fluxDistOldVersionBool = False previewMode = False - outputs = ['travelLengthMin', 'flux'] + 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 = 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) == 12 + 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 np.all(results[7] == np.ones_like(flux) * -9999) + assert results[7] is 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.geometry.geom_equals(refPolygons.geometry)) + + # 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() diff --git a/docs/moduleCom4FlowPy.rst b/docs/moduleCom4FlowPy.rst index dfc8e18a0..002b3ac7c 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,8 @@ 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) +- ``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``):