diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py old mode 100755 new mode 100644 index a572b9f24..d47f67f4d --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -77,6 +77,27 @@ def com4FlowPyMain(cfgPath, cfgSetup): # Flag for use of old flux distribution version modelParameters["fluxDistOldVersionBool"] = cfgSetup.getboolean("fluxDistOldVersion") + modelParameters["calcGeneration"] = cfgSetup.getboolean("calcGeneration") + modelParameters["calcThalweg"] = cfgSetup.getboolean("calcThalweg") + if modelParameters["calcThalweg"]: + modelParameters["thalwegReleaseArea"] = cfgSetup.getboolean("thalwegReleaseArea") + modelParameters["thalwegSaveRam"] = cfgSetup.getboolean("thalwegSaveRam") + modelParameters["videoRelId"] = cfgSetup.get("videoRelId", fallback="") + modelParameters["videoDataVariable"] = cfgSetup.get("videoDataVariable") + else: + modelParameters["thalwegReleaseArea"] = False + modelParameters["thalwegSaveRam"] = False + modelParameters["thalwegVariables"] = cfgSetup.get("thalwegVariables") + modelParameters["thalwegCenterOf"] = cfgSetup.get("thalwegCenterOf") + + if modelParameters["thalwegSaveRam"]: + if modelParameters["videoRelId"] != "": + message = f"If thalwegSaveRam is True, no video data is stored, please check the configuration settings!" + log.error(message) + raise ValueError(message) + + # modelParameters["infra"] = cfgSetup["infra"] + # modelParameters["forest"] = cfgSetup["forest"] # compute engine: "python" (default, Cell-based) or "numba" (JIT kernel) modelParameters["engine"] = cfgSetup.get("engine", "python").strip().lower() @@ -106,6 +127,7 @@ def com4FlowPyMain(cfgPath, cfgSetup): modelPaths["outputFileFormat"] = ".asc" else: modelPaths["outputFileFormat"] = ".tif" + modelPaths["thalwegDir"] = cfgPath["thalwegDir"] # check if 'customDirs' are used - alternative is 'default' AvaFrame Folder Structure modelPaths["useCustomDirs"] = True if cfgPath["customDirs"] == "True" else False @@ -134,7 +156,6 @@ def com4FlowPyMain(cfgPath, cfgSetup): forestParams = {} # check if calculation with forest if modelParameters["forestBool"]: - forestParams["forestModule"] = cfgSetup.get("forestModule") modelPaths["forestPath"] = cfgPath["forestPath"] # 'forestFriction' and 'forestDetrainment' parameters @@ -182,7 +203,10 @@ def com4FlowPyMain(cfgPath, cfgSetup): else: modelPaths["varExponentPath"] = "" - if "relIdPolygon" in modelPaths["outputFileList"] or "relIdCount" in modelPaths["outputFileList"]: + # conditions if relId is used + _outputPolygon = "relIdPolygon" in modelPaths["outputFileList"] + _outputCount = "relIdCount" in modelPaths["outputFileList"] + if _outputPolygon or _outputCount or modelParameters["thalwegReleaseArea"]: modelPaths["relIdPath"] = cfgPath["relIdPath"] modelParameters["outputRelIdBool"] = True else: @@ -213,6 +237,9 @@ def com4FlowPyMain(cfgPath, cfgSetup): demHeader = IOf.readRasterHeader(modelPaths["demPath"]) rasterAttributes["nodata"] = demHeader["nodata_value"] rasterAttributes["cellsize"] = demHeader["cellsize"] + rasterAttributes["xllcenter"] = demHeader["xllcenter"] + rasterAttributes["yllcenter"] = demHeader["yllcenter"] + rasterAttributes["nrows"] = demHeader["nrows"] # tile input layers and write tiles (pickled np.arrays) to temp Folder nTiles = tileInputLayers(modelParameters, modelPaths, rasterAttributes, tilingParameters) @@ -295,7 +322,7 @@ def startLogging(modelParameters, forestParams, modelPaths, MPOptions): for param, value in MPOptions.items(): log.info(f"{'%s:'%param : <20}{value : <5}") # log.info("{}:\t{}".format(param,value)) - log.info("------------------------") + log.info("------------------------") log.info(f"{'WorkDir:' : <12}{'%s'%modelPaths['workDir'] : <5}") log.info(f"{'ResultsDir:' : <12}{'%s'%modelPaths['resDir'] : <5}") # log.info("WorkDir: {}".format(modelPaths["workDir"])) @@ -326,16 +353,27 @@ def checkInputLayerDimensions(modelParameters, modelPaths): if _demHeader["ncols"] == _relHeader["ncols"] and _demHeader["nrows"] == _relHeader["nrows"]: log.info("DEM and Release Layer ok!") else: - log.error("Error: Release Layer doesn't match DEM!") - sys.exit(1) + message = "Error: Release Layer doesn't match DEM!" + log.error(message) + raise ValueError(message) + + if modelParameters["outputRelIdBool"]: + _relIdHeader = IOf.readRasterHeader(modelPaths["relIdPath"]) + if _demHeader["ncols"] == _relIdHeader["ncols"] and _demHeader["nrows"] == _relIdHeader["nrows"]: + log.info("Release ID Layer ok!") + else: + message = "Error: Release ID Layer doesn't match DEM!" + log.error(message) + raise ValueError(message) if modelParameters["infraBool"]: _infraHeader = IOf.readRasterHeader(modelPaths["infraPath"]) if _demHeader["ncols"] == _infraHeader["ncols"] and _demHeader["nrows"] == _infraHeader["nrows"]: log.info("Infra Layer ok!") else: - log.error("Error: Infra Layer doesn't match DEM!") - sys.exit(1) + message = "Error: Infra Layer doesn't match DEM!" + log.error(message) + raise ValueError(message) if modelParameters["forestBool"]: _forestHeader = IOf.readRasterHeader(modelPaths["forestPath"]) @@ -345,8 +383,9 @@ def checkInputLayerDimensions(modelParameters, modelPaths): ): log.info("Forest Layer ok!") else: - log.error("Error: Forest Layer doesn't match DEM!") - sys.exit(1) + message = "Error: Forest Layer doesn't match DEM!" + log.error(message) + raise ValueError(message) if modelParameters["varUmaxBool"]: _varUmaxHeader = IOf.readRasterHeader(modelPaths["varUmaxPath"]) @@ -356,8 +395,9 @@ def checkInputLayerDimensions(modelParameters, modelPaths): ): log.info("uMax Limit Layer ok!") else: - log.error("Error: uMax Limit Layer doesn't match DEM!") - sys.exit(1) + message = "Error: uMax Limit Layer doesn't match DEM!" + log.error(message) + raise ValueError(message) if modelParameters["varAlphaBool"]: _varAlphaHeader = IOf.readRasterHeader(modelPaths["varAlphaPath"]) @@ -367,8 +407,9 @@ def checkInputLayerDimensions(modelParameters, modelPaths): ): log.info("variable Alpha Layer ok!") else: - log.error("Error: variable Alpha Layer doesn't match DEM!") - sys.exit(1) + message = "Error: variable Alpha Layer doesn't match DEM!" + log.error(message) + raise ValueError(message) if modelParameters["varExponentBool"]: _varExponentHeader = IOf.readRasterHeader(modelPaths["varExponentPath"]) @@ -378,8 +419,9 @@ def checkInputLayerDimensions(modelParameters, modelPaths): ): log.info("variable exponent Layer ok!") else: - log.error("Error: variable exponent Layer doesn't match DEM!") - sys.exit(1) + message = "Error: variable exponent Layer doesn't match DEM!" + log.error(message) + raise ValueError(message) log.info("========================") @@ -389,7 +431,9 @@ def checkInputLayerDimensions(modelParameters, modelPaths): ) log.error("Error occured: %s" % ex) # return - sys.exit(1) + raise ValueError( + "could not read all required Input Layers, please re-check files and paths provided in .ini files" + ) def checkInputParameterValues(modelParameters, modelPaths, validParamRanges): """check if the input parameters are valid @@ -479,30 +523,151 @@ def tileInputLayers(modelParameters, modelPaths, rasterAttributes, tilingParamet log.info("Start Tiling...") log.info("---------------------") + if modelParameters["thalwegReleaseArea"]: + _relIdRasterDict = IOf.readRaster(modelPaths["relIdPath"]) + _relIdRaster = _relIdRasterDict["rasterData"] + exList, eyList = SPAM.getTileEnds(modelPaths["tempDir"], _tileCOLS, _tileROWS, _U, _relIdRaster) + + SPAM.tileRasterWithIndices(modelPaths["demPath"], "dem", modelPaths["tempDir"], exList, eyList, _U) + SPAM.tileRasterWithIndices( + modelPaths["releasePathWork"], + "init", + modelPaths["tempDir"], + exList, + eyList, + _U, + isInit=True, + ) - SPAM.tileRaster(modelPaths["demPath"], "dem", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U) - SPAM.tileRaster( - modelPaths["releasePathWork"], "init", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U, isInit=True - ) + if modelParameters["infraBool"]: + SPAM.tileRasterWithIndices( + modelPaths["infraPath"], + "infra", + modelPaths["tempDir"], + exList, + eyList, + _U, + ) + if modelParameters["varUmaxBool"]: + SPAM.tileRasterWithIndices( + modelPaths["varUmaxPath"], + "varUmax", + modelPaths["tempDir"], + exList, + eyList, + _U, + ) + if modelParameters["varAlphaBool"]: + SPAM.tileRasterWithIndices( + modelPaths["varAlphaPath"], + "varAlpha", + modelPaths["tempDir"], + exList, + eyList, + _U, + ) + if modelParameters["varExponentBool"]: + SPAM.tileRasterWithIndices( + modelPaths["varExponentPath"], + "varExponent", + modelPaths["tempDir"], + exList, + eyList, + _U, + ) + if modelParameters["forestBool"]: + SPAM.tileRasterWithIndices( + modelPaths["forestPath"], + "forest", + modelPaths["tempDir"], + exList, + eyList, + _U, + ) + if modelParameters["outputRelIdBool"]: + SPAM.tileRasterWithIndices( + modelPaths["relIdPath"], + "relId", + modelPaths["tempDir"], + exList, + eyList, + _U, + ) - 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 - ) - if modelParameters["varAlphaBool"]: + else: SPAM.tileRaster( - modelPaths["varAlphaPath"], "varAlpha", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U + modelPaths["demPath"], + "dem", + modelPaths["tempDir"], + _tileCOLS, + _tileROWS, + _U, ) - if modelParameters["varExponentBool"]: SPAM.tileRaster( - modelPaths["varExponentPath"], "varExponent", modelPaths["tempDir"], _tileCOLS, _tileROWS, _U + modelPaths["releasePathWork"], + "init", + modelPaths["tempDir"], + _tileCOLS, + _tileROWS, + _U, + isInit=True, ) - 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) + + 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, + ) + if modelParameters["varAlphaBool"]: + SPAM.tileRaster( + modelPaths["varAlphaPath"], + "varAlpha", + modelPaths["tempDir"], + _tileCOLS, + _tileROWS, + _U, + ) + if modelParameters["varExponentBool"]: + 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("==================================") @@ -536,7 +701,17 @@ def performModelCalculation(nTiles, modelParameters, modelPaths, rasterAttribute for i in range(nTiles[0] + 1): for j in range(nTiles[1] + 1): - optList.append((i, j, modelParameters, modelPaths, rasterAttributes, forestParams, MPOptions)) + optList.append( + ( + i, + j, + modelParameters, + modelPaths, + rasterAttributes, + forestParams, + MPOptions, + ) + ) log.info(" >> Start Calculation << ") log.info("-------------------------") @@ -771,7 +946,8 @@ def mergeAndWriteResults(modelPaths, modelOptions): if "relIdPolygon" in _outputs: pathPolygons = SPAM.mergeDictToPolygon(modelPaths["tempDir"], "res_startCellIdDict", outputHeader) pathPolygons.to_file( - modelPaths["resDir"] / "com4_{}_{}_relIdPolygon.geojson".format(_uid, _ts), driver="GeoJSON" + modelPaths["resDir"] / "com4_{}_{}_pathPolygons.geojson".format(_uid, _ts), + driver="GeoJSON", ) del pathPolygons log.info("com4_{}_{}_relIdPolygon is written".format(_uid, _ts)) @@ -812,7 +988,6 @@ def checkConvertReleaseShp2Tif(modelPaths): # the release is a shp polygon, we need to convert it to a raster # releaseLine = shpConv.readLine(releasePath, 'releasePolygon', demDict) if modelPaths["releasePath"].suffix == ".shp": - dem = IOf.readRaster(modelPaths["demPath"]) demHeader = dem["header"] dem["originalHeader"] = demHeader diff --git a/avaframe/com4FlowPy/com4FlowPyCfg.ini b/avaframe/com4FlowPy/com4FlowPyCfg.ini index 6c3e86e9b..64ec456e0 100644 --- a/avaframe/com4FlowPy/com4FlowPyCfg.ini +++ b/avaframe/com4FlowPy/com4FlowPyCfg.ini @@ -179,6 +179,40 @@ fluxDistOldVersion = False # numba : JIT-compiled kernel (double precision), bit-for-bit vs 'python' engine = python +#++++++++++++ Calculate with generations +# If calcGenerations = True, a different order of cells in a path are calculated. +# The results can vary when computing with generations. +# Additionally, the generation (iteration step) can be derived, which is required +# to get thalweg information. +calcGeneration = False +# set to True to compute thalwegs and save them as pickle files +# To compute thalwegs, calcGeneration needs to be set to True. +calcThalweg = False +# You can choose how the thalweg is computed: with the center of energy, center of flux and/or +# center of zdelta. The format should be: thalwegCenterOf = ['zdelta','energy', 'flux'] +thalwegCenterOf = ['zdelta','energy', 'flux'] +# The saved variables can be chosen in thalwegVariables. Possible variables are. +# ['col', 'row', 'x', 'y', 'flux', 'fluxSum' 'flowEnergy', 'altitude', 'travelLength', 'zDelta', 'gamma', 'flowEnergyArray', 'zDeltaArray', 'fluxArray'] +# the arrays contain the respective values of the path and need much memory +thalwegVariables = ['x', 'y', 'z', 's', 'zDelta'] +# if thalwegReleaseArea is True, one thalweg is computed for each continuous release area +# (in Inputs/RELID a raster file needs to be provided containing the release area ID) +# if thalwegReleaseArea is False, a thalweg is computed for each startcell +thalwegReleaseArea = True +# only compute thalweg variable x and y to safe RAM (also no video data) +thalwegSaveRam = True +# release ids for that data is stored to make a gif +# (multiple ids can be separated by a |, e.g. 1|3|7; if empty, no video data are saved) +# video data is only stored when thalwegSaveRam is False +videoRelId = +# define one output variable for the video data +# possible options: +# z_delta +# flux +# min_distance (max travel length) +# max_gamma (max travelAngle) +videoDataVariable = z_delta + #++++++++++++ Parameters for Tiling # tileSize: size of tiles in x and y direction in meters (if total size of) x # or y of input DEM is larger than tileSize, then the input raster diff --git a/avaframe/com4FlowPy/flowClass.py b/avaframe/com4FlowPy/flowClass.py index 8f6c6f6bd..b969e67a7 100644 --- a/avaframe/com4FlowPy/flowClass.py +++ b/avaframe/com4FlowPy/flowClass.py @@ -10,7 +10,6 @@ class Cell: - def __init__( self, rowindex, @@ -57,7 +56,7 @@ def __init__( self.z_delta = z_delta self.alpha = float(alpha) - self.exp = int(exp) + self.exp = float(exp) self.max_z_delta = float(max_z_delta) self.flux_threshold = float(flux_threshold) @@ -84,7 +83,6 @@ def __init__( # if FSI != None AND forestParams != None - then self.ForestBool = True and forestParams and # FSI are accordingly initialized if (FSI is not None) and (forestParams is not None): - self.forestBool = True self.forestModule = forestParams["forestModule"] self.skipForestDist = forestParams["skipForestDist"] @@ -115,7 +113,6 @@ def __init__( self.noDetrainmentEffectZdelta = (_vThDe * _vThDe) / _sqrt2xG elif self.forestModule == "forestFrictionLayer": - if forestParams["fFrLayerType"] == "absolute": self.AlphaFor = FSI elif forestParams["fFrLayerType"] == "relative": @@ -264,7 +261,6 @@ def calc_z_delta(self): self.calcDistMin(calc3D=True) if self.forestBool: - if self.forestModule == "forestFrictionLayer": if (not self.is_start) and (self.skipForestDist < self.minDistXYZ): _tanAlpha = self.tanAlphaFor @@ -319,6 +315,11 @@ def calc_tanbeta(self): if abs(np.sum(self.tan_beta)) > 0: self.r_t = self.tan_beta**self.exp / np.sum(self.tan_beta**self.exp) + def calcFlowEnergy(self): + # calculate flow energy (corresponding to kinetic energy) + # analog to: kin_energy = mass * velocity² / 2 + self.flowEnergy = self.flux * self.z_delta * 9.81 + def calc_persistence(self): """ calculates persistence-based routing @@ -407,6 +408,7 @@ def calc_distribution(self): # substituted by self.flux_threshold???? self.flux = max(0.0003, self.flux - self.detrainment) + self.calcFlowEnergy() threshold = self.flux_threshold if np.sum(self.r_t) > 0: self.dist = (self.persistence * self.r_t) / np.sum(self.persistence * self.r_t) * self.flux @@ -449,10 +451,10 @@ def calc_distribution(self): row_local, col_local = np.where(self.dist >= threshold) return ( - self.rowindex - 1 + row_local, - self.colindex - 1 + col_local, - self.dist[row_local, col_local], - self.z_delta_neighbour[row_local, col_local], + list(self.rowindex - 1 + row_local), + list(self.colindex - 1 + col_local), + list(self.dist[row_local, col_local]), + list(self.z_delta_neighbour[row_local, col_local]), ) def forest_detrainment(self): diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index 99d0c9ae8..a5ad1fbea 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -8,21 +8,25 @@ import sys import numpy as np import logging -import os -import platform +import pathlib import gc import psutil import time import pickle +from itertools import zip_longest from multiprocessing import Pool from avaframe.com4FlowPy.flowClass import Cell +from avaframe.com4FlowPy.flowPath import Path +log = logging.getLogger(__name__) -def get_start_idx(dem, release): +def get_start_idx(dem, release, relIdArray=None, calcThalweg=False): """Sort Release Pixels by altitude and return the result as lists for the Rows and Columns, starting with the highest altitude + If releaseIds are provided, sort by release Id to ensure that all cells belonging to a + segmented PRA are computed after each other. Parameters ----------- @@ -30,6 +34,10 @@ def get_start_idx(dem, release): Digital Elevation Model to gain information about altitude release: numpy array The release layer, release pixels need int value > 0 + relIdArray: numpy array + release Ids + calcThalweg: bool + flag if thalweg is computed Returns ----------- @@ -41,16 +49,27 @@ def get_start_idx(dem, release): row_list, col_list = np.where(release > 0) # Gives back the indices of the release areas if len(row_list) > 0: altitude_list = [] + relIdList = [] for i in range(len(row_list)): altitude_list.append(dem[row_list[i], col_list[i]]) - altitude_list, row_list, col_list = list( - zip(*sorted(zip(altitude_list, row_list, col_list), reverse=True)) - ) - # Sort this lists by altitude + if relIdArray is not None and calcThalweg: + relIdList.append(relIdArray[row_list[i], col_list[i]]) + + # sort this list by releaseId + if relIdArray is not None and calcThalweg: + relIdList, altitude_list, row_list, col_list = list( + zip(*sorted(zip(relIdList, altitude_list, row_list, col_list), reverse=True)) + ) + else: + # Sort this lists by altitude + altitude_list, row_list, col_list = list( + zip(*sorted(zip(altitude_list, row_list, col_list), reverse=True)) + ) + return row_list, col_list -def split_release(release, pieces): +def split_release(release, pieces, relIdArray, calcThalweg): """Split the release layer in several tiles. The area is determined by the number of release pixels in it, so that every tile has the same amount of release pixels in it. @@ -66,12 +85,20 @@ def split_release(release, pieces): The release tiles have still the size of the original layer, so no split for the DEM is needed. + If thalweg is computed, and release Ids are provided, cells belonging to one + release Ids are not divided into separate release_lists (for different chunks), + so the thalweg can be computed for one segmented PRA. + Parameters ----------- release: np.array a binary 0|1 array with release pixels designated by '1' pieces: int number of chunck in which the release layer should be split + relIdArray: numpy array + release Ids of segmented relase areas + calcThalweg: bool + flag if thalweg is calculated Returns ----------- @@ -79,33 +106,63 @@ def split_release(release, pieces): contains the tiles(arrays) [array0, array1, ..] """ - # Flatten the array and compute the cumulative sum - flat_release = release.flatten() - cumulative_sum = np.cumsum(flat_release) + if calcThalweg and relIdArray is not None: + # release split for thalweg computation - total_sum = cumulative_sum[-1] - sum_per_split = total_sum / pieces + uniqueIds, counts = np.unique(relIdArray[release == 1], return_counts=True) - release_list = [] - start_index = 0 + pieces = np.minimum(pieces, len(uniqueIds)) - for i in range(1, pieces): - # Find the split point in the flattened array - split_index = np.searchsorted(cumulative_sum, sum_per_split * i) + idCount = list(zip(uniqueIds, counts)) + idCount.sort(key=lambda x: x[1], reverse=True) - # Create a new array for this split - split_flat = np.zeros_like(flat_release) - split_flat[start_index:split_index] = flat_release[start_index:split_index] + # prepare lists for ids and number of cells + _numberCells = np.zeros(pieces, dtype=int) + _ids = [[] for _ in range(pieces)] - # Reshape the flat array back to 2D and add to the list - release_list.append(split_flat.reshape(release.shape)) + # add relId to this chunk that has less cells yet + for id, count in idCount: + idx = np.argmin(_numberCells) + _ids[idx].append(id) + _numberCells[idx] += count + + # write release cells + release_list = [] + for idsChunk in _ids: + release_piece = np.zeros_like(release, dtype=release.dtype) + if idsChunk: + id_piece = np.isin(relIdArray, idsChunk) + release_piece[id_piece] = release[id_piece] + release_list.append(release_piece) + + else: + # Flatten the array and compute the cumulative sum + flat_release = release.flatten() + cumulative_sum = np.cumsum(flat_release) + + total_sum = cumulative_sum[-1] + sum_per_split = total_sum / pieces - start_index = split_index + release_list = [] + start_index = 0 - # Handle the last piece - split_flat = np.zeros_like(flat_release) - split_flat[start_index:] = flat_release[start_index:] - release_list.append(split_flat.reshape(release.shape)) + for i in range(1, pieces): + # Find the split point in the flattened array + split_index = np.searchsorted(cumulative_sum, sum_per_split * i) + + # Create a new array for this split + split_flat = np.zeros_like(flat_release) + split_flat[start_index:split_index] = flat_release[start_index:split_index] + + # Reshape the flat array back to 2D and add to the list + release_list.append(split_flat.reshape(release.shape)) + + start_index = split_index + + # Handle the last piece + split_flat = np.zeros_like(flat_release) + split_flat[start_index:] = flat_release[start_index:] + release_list.append(split_flat.reshape(release.shape)) return release_list @@ -146,6 +203,22 @@ def run(optTuple): fluxDistOldVersionBool = optTuple[2]["fluxDistOldVersionBool"] relIdBool = optTuple[2]["outputRelIdBool"] previewMode = optTuple[2]["previewMode"] + calcGeneration = optTuple[2]["calcGeneration"] + calcThalweg = optTuple[2]["calcThalweg"] + if calcThalweg: + thalwegParameters = { + "thalwegDir": optTuple[3]["thalwegDir"], + "thalwegCenterOf": optTuple[2]["thalwegCenterOf"], + "thalwegVariables": optTuple[2]["thalwegVariables"], + "calcRelID": optTuple[2]["thalwegReleaseArea"], + "thalwegSaveRam": optTuple[2]["thalwegSaveRam"], + "videoRelId": optTuple[2]["videoRelId"], + "videoDataVariable": optTuple[2]["videoDataVariable"], + } + else: + thalwegParameters = { + "thalwegSaveRam": False, + } # Temp-Dir (all input files are located here and results are written back in here) tempDir = optTuple[3]["tempDir"] @@ -154,6 +227,7 @@ def run(optTuple): outputs = optTuple[3]["outputFileList"] # raster-layer Attributes + rasterAttributes = optTuple[4] cellsize = float(optTuple[4]["cellsize"]) nodata = float(optTuple[4]["nodata"]) @@ -161,6 +235,8 @@ def run(optTuple): dem = np.load(tempDir / ("dem_%s_%s.npy" % (optTuple[0], optTuple[1]))) release = np.load(tempDir / ("init_%s_%s.npy" % (optTuple[0], optTuple[1]))) + extentTile = np.load(tempDir / ("ext_%s_%s" % (optTuple[0], optTuple[1])), allow_pickle=True) + rasterAttributes["extentTile"] = extentTile if infraBool: infra = np.load(tempDir / ("infra_%s_%s.npy" % (optTuple[0], optTuple[1]))) else: @@ -220,6 +296,7 @@ def run(optTuple): # every positive value >0 is interpreted as release area release[release < 0] = 0 release[release == nodata] = 0 # added in case nodata is non-negative + release[np.isnan(release)] = 0 release[release > 0] = 1 nRel = np.sum(release) @@ -233,7 +310,7 @@ def run(optTuple): chunkSize=MPOptions["chunkSize"], ) - release_list = split_release(release, nChunks) + release_list = split_release(release, nChunks, relIdArray, calcThalweg) # select compute engine: "numba" JIT kernel, else the default Python (Cell) path. # numba does not (yet) implement infra/back-calculation, previewMode or the @@ -273,8 +350,7 @@ def run(optTuple): exp, flux_threshold, max_z_delta, - nodata, - cellsize, + rasterAttributes, infraBool, forestBool, varParams, @@ -284,6 +360,9 @@ def run(optTuple): forestParams, outputs, relOutputParams, + calcGeneration, + calcThalweg, + thalwegParameters, ] for release_sub in release_list ], @@ -384,7 +463,6 @@ def run(optTuple): np.minimum(forestIntArray, forestIntList[i]), np.maximum(forestIntArray, forestIntList[i]), ) - if "relIdPolygon" in outputs or "relIdCount" in outputs: for key in processedStartCellIdList[i]: if key in processedStartCellIdDict: @@ -394,26 +472,51 @@ def run(optTuple): processedStartCellIdDict[key] = processedStartCellIdList[i][key] if relOutputParams["relIdBool"]: - saveDict = open(tempDir / ("res_startCellIdDict_%s_%s.pickle" % (optTuple[0], optTuple[1])), "wb") + 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) - np.save(tempDir / ("res_rout_flux_sum_%s_%s" % (optTuple[0], optTuple[1])), routFluxSumArray) - np.save(tempDir / ("res_dep_flux_sum_%s_%s" % (optTuple[0], optTuple[1])), depFluxSumArray) + np.save( + tempDir / ("res_rout_flux_sum_%s_%s" % (optTuple[0], optTuple[1])), + routFluxSumArray, + ) + np.save( + tempDir / ("res_dep_flux_sum_%s_%s" % (optTuple[0], optTuple[1])), + depFluxSumArray, + ) np.save(tempDir / ("res_flux_%s_%s" % (optTuple[0], optTuple[1])), fluxArray) np.save(tempDir / ("res_count_%s_%s" % (optTuple[0], optTuple[1])), countArray) - np.save(tempDir / ("res_fp_max_%s_%s" % (optTuple[0], optTuple[1])), fpTravelAngleMaxArray) - np.save(tempDir / ("res_fp_min_%s_%s" % (optTuple[0], optTuple[1])), fpTravelAngleMinArray) + np.save( + tempDir / ("res_fp_max_%s_%s" % (optTuple[0], optTuple[1])), + fpTravelAngleMaxArray, + ) + np.save( + tempDir / ("res_fp_min_%s_%s" % (optTuple[0], optTuple[1])), + fpTravelAngleMinArray, + ) 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_travel_length_max_%s_%s" % (optTuple[0], optTuple[1])), + travelLengthMaxArray, + ) + np.save( + tempDir / ("res_travel_length_min_%s_%s" % (optTuple[0], optTuple[1])), + travelLengthMinArray, + ) if infraBool: np.save(tempDir / ("res_backcalc_%s_%s" % (optTuple[0], optTuple[1])), backcalc) if forestInteraction: - np.save(tempDir / ("res_forestInt_%s_%s" % (optTuple[0], optTuple[1])), forestIntArray) + np.save( + tempDir / ("res_forestInt_%s_%s" % (optTuple[0], optTuple[1])), + forestIntArray, + ) def calculation(args): @@ -432,18 +535,20 @@ def calculation(args): - args[4] (float) - exponent - args[5] (float) - threshold of minimum flux - args[6] (float) - maximum of zDelta - - args[7] (float) - nodata values of rasters - - args[8] (float) - cellsize of rasters - - args[9] (bool) - flag for calculation with/without infrastructure - - args[10] (bool) - flag for calculation with/without forest - - args[11] (dict) - contains flags and numpy arrays for variable input parameters (Alpha, exp, uMax) - - args[12] (bool) - flag for computing flux distribution with old version - - args[13] (bool) - flag for previewMode / fast Calculation - - - 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 + - args[7] (float) - raster attributes + - args[8] (bool) - flag for calculation with/without infrastructure + - args[9] (bool) - flag for calculation with/without forest + - args[10] (dict) - contains flags and numpy arrays for variable input parameters (Alpha, exp, uMax) + - args[11] (bool) - flag for computing flux distribution with old version + - args[12] (bool) - flag for previewMode / fast Calculation + + - args[13] (numpy array) - contains forest information (None if forestBool=False) + - args[14] (dict) - contains parameters for forest interaction models (None if forestBool=False) + - args[15] (list) - output names + - args[16] (dict) - contains flags and rasters for release - information outputs + - args[17] (bool) - flag for computing each generation + - args[18] (bool) - flag for computing thalweg + - args[19] (dict) - thalweg parameters Returns ----------- @@ -471,6 +576,7 @@ def calculation(args): minimum of the count a forested cell is hit (only returned if args[18]["forestInteraction"]==True) """ + log = logging.getLogger(__name__) # helper function for backTracking, a bit slower than inline but improves # readability by avoiding repetitions @@ -498,25 +604,29 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): exp = args[4] flux_threshold = args[5] max_z_delta = args[6] - nodata = args[7] - cellsize = args[8] - infraBool = args[9] - forestBool = args[10] - varUmaxBool = args[11]["varUmaxBool"] - varUmaxArray = args[11]["varUmaxArray"] - varAlphaBool = args[11]["varAlphaBool"] - varAlphaArray = args[11]["varAlphaArray"] - varExponentBool = args[11]["varExponentBool"] - varExponentArray = args[11]["varExponentArray"] - fluxDistOldVersionBool = args[12] - previewMode = args[13] - outputs = args[16] - relIdArray = args[17]["relIdArray"] - relIdBool = args[17]["relIdBool"] + rasterAttributes = args[7] + cellsize = rasterAttributes["cellsize"] + nodata = rasterAttributes["nodata"] + infraBool = args[8] + forestBool = args[9] + varUmaxBool = args[10]["varUmaxBool"] + varUmaxArray = args[10]["varUmaxArray"] + varAlphaBool = args[10]["varAlphaBool"] + varAlphaArray = args[10]["varAlphaArray"] + varExponentBool = args[10]["varExponentBool"] + varExponentArray = args[10]["varExponentArray"] + fluxDistOldVersionBool = args[11] + previewMode = args[12] + outputs = args[15] + relIdArray = args[16]["relIdArray"] + relIdBool = args[16]["relIdBool"] + calcGeneration = args[17] + calcThalweg = args[18] + thalwegParameters = args[19] if forestBool: - forestArray = args[14] - forestParams = args[15] + forestArray = args[13] + forestParams = args[14] forestInteraction = forestParams["forestInteraction"] else: forestInteraction = False @@ -572,11 +682,20 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): # Core # NOTE-TODO: row_list, col_list are tuples - rethink variable naming - row_list, col_list = get_start_idx(dem, release) - + row_list, col_list = get_start_idx(dem, release, relIdArray, calcThalweg) + + generationListRelId = [] + colListRelId = [] + rowListRelId = [] + fluxListRelId = [] + zdeltaListRelId = [] + travelLengthMaxListRelId = [] startcell_idx = 0 + if calcThalweg: + timeThalweg = 0.0 + nextRowIdx = row_list[0] + nextColIdx = col_list[0] while startcell_idx < len(row_list): - if infraBool: # if infraBool - here we initialize a directed graph structure pathTopology = {} # topology of path as directed graph @@ -584,7 +703,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): processedCells = {} # dictionary of cells that have been processed already zDeltaPathArray = np.zeros_like(dem, dtype=np.float32) - cell_list = [] + row_idx = row_list[startcell_idx] col_idx = col_list[startcell_idx] dem_ng = dem[row_idx - 1 : row_idx + 2, col_idx - 1 : col_idx + 2] # neighbourhood DEM @@ -637,144 +756,515 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): processedCells[(startcell.rowindex, startcell.colindex)] = 1 # list of flowClass.Cell() Objects that is contains the "path" for each release-cell - cell_list.append(startcell) + if calcGeneration: + cellList = [startcell] # list of parents for current iteration + genList = [cellList] # list of all cells (which are calculated), organised in generations + childList = [] # list of childs of the current iteration + childIndex = {} + if thalwegParameters["thalwegSaveRam"]: + colThalwegLists = [] + rowThalwegLists = [] + fluxThalwegLists = [] + zdeltaThalwegLists = [] + travelLengthMaxThalwegLists = [] + + for gen, cellList in enumerate(genList): + if thalwegParameters["thalwegSaveRam"]: + colThalwegGen = [] + rowThalwegGen = [] + fluxThalwegGen = [] + zdeltaThalwegGen = [] + travelLengthMaxThalwegGen = [] + for idx, cell in enumerate(cellList): + 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 + # number of child cells (successors) to currently processed cell + row, col, flux, z_delta = cell.calc_distribution() + + if len(row) > 0: + # mass, row, col = list(zip(*sorted(zip( mass, row, col), reverse=False))) + z_delta, flux, row, col = list( + zip(*sorted(zip(z_delta, flux, row, col), reverse=False)) + ) + # Sort this lists by elh, to start with the highest cell + + if infraBool: + # if the current cell is not already in the dir-graph, then we add it here + updateInfraDirGraph(cell.rowindex, cell.colindex) + + newRow = [] + newCol = [] + newFlux = [] + newZDelta = [] + + for r, c, f, zd in zip(row, col, flux, z_delta): + key = (r, c) + + if key in childIndex: + child = childList[childIndex[key]] + child.add_os(f) + child.add_parent(cell) + + if infraBool: + updateInfraDirGraph(r, c, cell.rowindex, cell.colindex) + + if zd > child.z_delta: + child.z_delta = zd + + else: + newRow.append(r) + newCol.append(c) + newFlux.append(f) + newZDelta.append(zd) + + row = newRow + col = newCol + flux = newFlux + z_delta = newZDelta + + + # TODO: we could put this checking part in an extra function, if we can move updateInfraDirGraph + """ + # I substitute this with checking for the cell in the dict (part above) + for i in range(len(childList)): # Check if Cell already exists in childList + k = 0 + while k < len(row): + if row[k] == childList[i].rowindex and col[k] == childList[i].colindex: + childList[i].add_os(flux[k]) + childList[i].add_parent(cell) + + if infraBool: + updateInfraDirGraph(row[k], col[k], cell.rowindex, cell.colindex) + + if z_delta[k] > childList[i].z_delta: + childList[i].z_delta = z_delta[k] + row = np.delete(row, k) + col = np.delete(col, k) + flux = np.delete(flux, k) + z_delta = np.delete(z_delta, k) + else: + k += 1 + """ + + for k in range(len(row)): + dem_ng = dem[row[k] - 1 : row[k] + 2, col[k] - 1 : col[k] + 2] # neighbourhood DEM + + # This bit handles edge cases and noData-values in the DEM!! this is an important piece of code, since + # no-data handling is expected (by some users/applications) to behave like here: + # i.e. if nodata in the 3x3 neighbourhood --> no calculation + if (nodata in dem_ng) or np.size(dem_ng) < 9: + continue + if infraBool: + updateInfraDirGraph(row[k], col[k], cell.rowindex, cell.colindex) - if infraBool: - # adding start-cell as "root-node" to directed graph of the modeled process path - 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) + # if the current child cell is already in processedCells + # just add +1 to the visit-counter, else add it to the + # processedCells dictionary with visit-count = 1 + if (row[k], col[k]) in processedCells: + processedCells[(row[k], col[k])] += 1 + else: + processedCells[(row[k], col[k])] = 1 + + childIndex[(row[k], col[k])] = len( + childList) # important that it's before childlist.append(...) + childList.append( + Cell( + row[k], + col[k], + dem_ng, + cellsize, + flux[k], + z_delta[k], + cell, + alpha, + exp, + flux_threshold, + max_z_delta, + startcell, + FSI=( + forestArray[row[k], col[k]] + if isinstance(forestArray, np.ndarray) + else None + ), + forestParams=forestParams, + ) + ) + + # TODO: writing arrays in a separate function? + routFluxSumArray[cell.rowindex, cell.colindex] += cell.flux + depFluxSumArray[cell.rowindex, cell.colindex] += cell.fluxDep + 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 + ) + zDeltaPathArray[cell.rowindex, cell.colindex] = max( + zDeltaPathArray[cell.rowindex, cell.colindex], cell.z_delta + ) + if "fpTravelAngleMax" in outputs or "fpTravelAngle" in outputs: + fpTravelAngleMaxArray[cell.rowindex, cell.colindex] = max( + fpTravelAngleMaxArray[cell.rowindex, cell.colindex], + cell.max_gamma, + ) + if "fpTravelAngleMin" in outputs: + if fpTravelAngleMinArray[cell.rowindex, cell.colindex] >= 0 and cell.max_gamma >= 0: + fpTravelAngleMinArray[cell.rowindex, cell.colindex] = min( + fpTravelAngleMinArray[cell.rowindex, cell.colindex], + cell.max_gamma, + ) + else: + fpTravelAngleMinArray[cell.rowindex, cell.colindex] = max( + fpTravelAngleMinArray[cell.rowindex, cell.colindex], + cell.max_gamma, + ) + slTravelAngleArray[cell.rowindex, cell.colindex] = max( + slTravelAngleArray[cell.rowindex, cell.colindex], cell.sl_gamma + ) + if "travelLengthMax" in outputs or "travelLength" in outputs: + travelLengthMaxArray[cell.rowindex, cell.colindex] = max( + travelLengthMaxArray[cell.rowindex, cell.colindex], + cell.min_distance, + ) + if "travelLengthMin" in outputs: + if ( + travelLengthMinArray[cell.rowindex, cell.colindex] >= 0 + and cell.min_distance >= 0 + ): + travelLengthMinArray[cell.rowindex, cell.colindex] = min( + travelLengthMinArray[cell.rowindex, cell.colindex], + cell.min_distance, + ) + else: + travelLengthMinArray[cell.rowindex, cell.colindex] = max( + travelLengthMinArray[cell.rowindex, cell.colindex], + cell.min_distance, + ) + + # TODO: why does the cell count not work as without generation-computation? + if processedCells[(cell.rowindex, cell.colindex)] == 1: + countArray[cell.rowindex, cell.colindex] += int(1) + elif ( + processedCells[(cell.rowindex, cell.colindex)] > 1 + and countArray[cell.rowindex, cell.colindex] <= 0 + ): + countArray[cell.rowindex, cell.colindex] += int(1) + + if forestInteraction: + if forestIntArray[cell.rowindex, cell.colindex] >= 0 and cell.forestIntCount >= 0: + forestIntArray[cell.rowindex, cell.colindex] = min( + forestIntArray[cell.rowindex, cell.colindex], + cell.forestIntCount, + ) + else: + forestIntArray[cell.rowindex, cell.colindex] = max( + forestIntArray[cell.rowindex, cell.colindex], + cell.forestIntCount, + ) + if thalwegParameters["thalwegSaveRam"]: + colThalwegGen.append(cell.colindex) + rowThalwegGen.append(cell.rowindex) + fluxThalwegGen.append(cell.flux) + zdeltaThalwegGen.append(cell.z_delta) + travelLengthMaxThalwegGen.append(cell.min_distance) + + if len(childList) > 0: + cellList = childList + genList.append(cellList) + childList = [] + childIndex = {} + + if thalwegParameters["thalwegSaveRam"]: + colThalwegLists.append(colThalwegGen) + rowThalwegLists.append(rowThalwegGen) + fluxThalwegLists.append(fluxThalwegGen) + zdeltaThalwegLists.append(zdeltaThalwegGen) + travelLengthMaxThalwegLists.append(travelLengthMaxThalwegGen) + + # empty last generation in genList to save RAM + if gen > 1: + genList[gen - 1] = [] + + if calcThalweg and thalwegParameters["calcRelID"]: + if thalwegParameters["thalwegSaveRam"]: + colListRelId = [ + (colThisCell or []) + (colBefore or []) + for colThisCell, colBefore in zip_longest(colListRelId, colThalwegLists) + ] + rowListRelId = [ + (rowThisCell or []) + (rowBefore or []) + for rowThisCell, rowBefore in zip_longest(rowListRelId, rowThalwegLists) + ] + fluxListRelId = [ + (fluxThisCell or []) + (fluxBefore or []) + for fluxThisCell, fluxBefore in zip_longest(fluxListRelId, fluxThalwegLists) + ] + zdeltaListRelId = [ + (zdeltaThisCell or []) + (zdeltaBefore or []) + for zdeltaThisCell, zdeltaBefore in zip_longest(zdeltaListRelId, zdeltaThalwegLists) + ] + travelLengthMaxListRelId = [ + (travelLengthMaxThisCell or []) + (travelLengthMaxBefore or []) + for travelLengthMaxThisCell, travelLengthMaxBefore in zip_longest( + travelLengthMaxListRelId, travelLengthMaxThalwegLists + ) + ] + else: + # zip the generationLists within one release Id + generationListRelId = [ + (generationThisCell or []) + (generationBefore or []) + for generationThisCell, generationBefore in zip_longest(generationListRelId, genList) + ] + + # check if the next startcell has the same startcellId + if startcell_idx + 1 < len(row_list): + nextRowIdx = row_list[startcell_idx + 1] + nextColIdx = col_list[startcell_idx + 1] + lastStartcell = False else: - startCellIdDict[(cell.rowindex, cell.colindex)] = np.array([startcellId]) + # if this was the last startcell, we also want to compute the thalweg! + lastStartcell = True + if startcellId != relIdArray[nextRowIdx, nextColIdx] or lastStartcell: + log.info(f"Finished computing PRA with ID {startcellId}. Start computing its thalweg!") + timeThawlegStart = time.time() + if thalwegParameters["thalwegSaveRam"]: + listsRelId = { + "row": rowListRelId, + "col": colListRelId, + "flux": fluxListRelId, + "zdelta": zdeltaListRelId, + "travelLengthMax": travelLengthMaxListRelId, + } + path = Path( + dem, + row_list[startcell_idx], + col_list[startcell_idx], + None, + rasterAttributes, + countArray, + startcellId, + listsRelId, + cellList[0], + ) + path.calcAndSaveThalwegData(thalwegParameters) + del path + colListRelId = [] + rowListRelId = [] + fluxListRelId = [] + zdeltaListRelId = [] + travelLengthMaxListRelId = [] + listsRelId = {} - # calculate flux, z_delta from current cell (cell) to child-cells - # lenght of row, col, flux, and z_delta vectors correspond to - # number of child cells (successors) to currently processed cell - row, col, flux, z_delta = cell.calc_distribution() + else: + if str(int(startcellId)) in list(thalwegParameters["videoRelId"].split("|")): + saveGenerationVideoData(generationListRelId, + startcellId, + dem, + rasterAttributes, + outDir=thalwegParameters["thalwegDir"] / "videoData", + variable=thalwegParameters["videoDataVariable"]) + path = Path( + dem, + row_list[startcell_idx], + col_list[startcell_idx], + generationListRelId, + rasterAttributes, + countArray, + startcellId, + ) + path.calcAndSaveThalwegData(thalwegParameters) + del path + generationListRelId = [] + timeThalwegEnd = time.time() + timeThalweg += timeThalwegEnd - timeThawlegStart + log.info( + f"Finished computing thalweg of PRA with ID {startcellId}, it took {np.round(timeThalwegEnd - timeThawlegStart, 1)} s." + ) + + elif calcThalweg: + path = Path( + dem, + row_list[startcell_idx], + col_list[startcell_idx], + genList, + rasterAttributes, + countArray, + ) + path.calcAndSaveThalwegData(thalwegParameters) - if len(flux) > 0: # i.e. if there are child cells - # Sort this lists by z_delta, to start with the highest cell - z_delta, flux, row, col = list(zip(*sorted(zip(z_delta, flux, row, col), reverse=False))) + else: + cellList = [] + cellIndex = {} + cellList.append(startcell) if infraBool: - # if the current cell is not already in the dir-graph, then we add it here - updateInfraDirGraph(cell.rowindex, cell.colindex) + # adding start-cell as "root-node" to directed graph of the modeled process path + updateInfraDirGraph(startcell.rowindex, startcell.colindex) + + for idx, cell in enumerate(cellList): + 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]) - # check if child cells already exist - for i in range(idx, len(cell_list)): - k = 0 - while k < len(row): - 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) + row, col, flux, z_delta = cell.calc_distribution() - if infraBool: - updateInfraDirGraph(row[k], col[k], cell.rowindex, cell.colindex) + if len(flux) > 0: + # mass, row, col = list(zip(*sorted(zip( mass, row, col), reverse=False))) + z_delta, flux, row, col = list(zip(*sorted(zip(z_delta, flux, row, col), reverse=False))) + # Sort this lists by elh, to start with the highest cell - if z_delta[k] > cell_list[i].z_delta: - cell_list[i].z_delta = z_delta[k] - row = np.delete(row, k) - col = np.delete(col, k) - flux = np.delete(flux, k) - z_delta = np.delete(z_delta, k) - else: - k += 1 + if infraBool: + # if the current cell is not already in the dir-graph, then we add it here + updateInfraDirGraph(cell.rowindex, cell.colindex) - for k in range(len(row)): - dem_ng = dem[row[k] - 1 : row[k] + 2, col[k] - 1 : col[k] + 2] # neighbourhood DEM + # check if cell already exists + newRow = [] + newCol = [] + newFlux = [] + newZDelta = [] - # This bit handles edge cases and noData-values in the DEM!! this is an important piece of code, since - # no-data handling is expected (by some users/applications) to behave like here: - # i.e. if nodata in the 3x3 neighbourhood --> no calculation - if (nodata in dem_ng) or np.size(dem_ng) < 9: - continue + for r, c, f, zd in zip(row, col, flux, z_delta): + key = (r, c) - if infraBool: - updateInfraDirGraph(row[k], col[k], cell.rowindex, cell.colindex) + if key in cellIndex: + cellExist = cellList[cellIndex[key]] + cellExist.add_os(f) + cellExist.add_parent(cell) - # if the current child cell is already in processedCells - # just add +1 to the visit-counter, else add it to the - # processedCells dictionary with visit-count = 1 - if (row[k], col[k]) in processedCells: - processedCells[(row[k], col[k])] += 1 - 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, + if infraBool: + updateInfraDirGraph(r, c, cellExist.rowindex, cellExist.colindex) + + if zd > cellExist.z_delta: + cellExist.z_delta = zd + + else: + newRow.append(r) + newCol.append(c) + newFlux.append(f) + newZDelta.append(zd) + + row = newRow + col = newCol + flux = newFlux + z_delta = newZDelta + + + for k in range(len(row)): + dem_ng = dem[row[k] - 1 : row[k] + 2, col[k] - 1 : col[k] + 2] # neighbourhood DEM + + # This bit handles edge cases and noData-values in the DEM!! this is an important piece of + # code, sinceno-data handling is expected (by some users/applications) to behave like here: + # i.e. if nodata in the 3x3 neighbourhood --> no calculation + if (nodata in dem_ng) or np.size(dem_ng) < 9: + continue + + if infraBool: + updateInfraDirGraph(row[k], col[k], cell.rowindex, cell.colindex) + + # if the current child cell is already in processedCells + # just add +1 to the visit-counter, else add it to the + # processedCells dictionary with visit-count = 1 + if (row[k], col[k]) in processedCells: + processedCells[(row[k], col[k])] += 1 + else: + processedCells[(row[k], col[k])] = 1 + cellIndex[(row[k], col[k])] = len( + cellList) # important that it's before childlist.append(...) + cellList.append( + Cell( + row[k], + col[k], + dem_ng, + cellsize, + flux[k], + z_delta[k], + cell, + alpha, + exp, + flux_threshold, + max_z_delta, + startcell, + FSI=forestArray[row[k], col[k]] if isinstance(forestArray, np.ndarray) else None, + forestParams=forestParams, + ) ) + zDeltaArray[cell.rowindex, cell.colindex] = max( + zDeltaArray[cell.rowindex, cell.colindex], cell.z_delta ) - - 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) - routFluxSumArray[cell.rowindex, cell.colindex] += cell.flux - depFluxSumArray[cell.rowindex, cell.colindex] += cell.fluxDep - zDeltaPathArray[cell.rowindex, cell.colindex] = max( - zDeltaPathArray[cell.rowindex, cell.colindex], cell.z_delta - ) - if "fpTravelAngleMax" in outputs or "fpTravelAngle" in outputs: - fpTravelAngleMaxArray[cell.rowindex, cell.colindex] = max( - fpTravelAngleMaxArray[cell.rowindex, cell.colindex], cell.max_gamma + fluxArray[cell.rowindex, cell.colindex] = max( + fluxArray[cell.rowindex, cell.colindex], cell.flux ) - if "fpTravelAngleMin" in outputs: - if fpTravelAngleMinArray[cell.rowindex, cell.colindex] >= 0 and cell.max_gamma >= 0: - fpTravelAngleMinArray[cell.rowindex, cell.colindex] = min( - fpTravelAngleMinArray[cell.rowindex, cell.colindex], cell.max_gamma - ) - else: - fpTravelAngleMinArray[cell.rowindex, cell.colindex] = max( - fpTravelAngleMinArray[cell.rowindex, cell.colindex], cell.max_gamma - ) - slTravelAngleArray[cell.rowindex, cell.colindex] = max( - slTravelAngleArray[cell.rowindex, cell.colindex], cell.sl_gamma - ) - if "travelLengthMax" in outputs or "travelLength" in outputs: - travelLengthMaxArray[cell.rowindex, cell.colindex] = max( - travelLengthMaxArray[cell.rowindex, cell.colindex], cell.min_distance + routFluxSumArray[cell.rowindex, cell.colindex] += cell.flux + depFluxSumArray[cell.rowindex, cell.colindex] += cell.fluxDep + zDeltaPathArray[cell.rowindex, cell.colindex] = max( + zDeltaPathArray[cell.rowindex, cell.colindex], cell.z_delta ) - if "travelLengthMin" in outputs: - if travelLengthMinArray[cell.rowindex, cell.colindex] >= 0 and cell.min_distance >= 0: - travelLengthMinArray[cell.rowindex, cell.colindex] = min( - travelLengthMinArray[cell.rowindex, cell.colindex], cell.min_distance + if "fpTravelAngleMax" in outputs or "fpTravelAngle" in outputs: + fpTravelAngleMaxArray[cell.rowindex, cell.colindex] = max( + fpTravelAngleMaxArray[cell.rowindex, cell.colindex], + cell.max_gamma, ) - else: - 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) + if "fpTravelAngleMin" in outputs: + if fpTravelAngleMinArray[cell.rowindex, cell.colindex] >= 0 and cell.max_gamma >= 0: + fpTravelAngleMinArray[cell.rowindex, cell.colindex] = min( + fpTravelAngleMinArray[cell.rowindex, cell.colindex], + cell.max_gamma, + ) + else: + fpTravelAngleMinArray[cell.rowindex, cell.colindex] = max( + fpTravelAngleMinArray[cell.rowindex, cell.colindex], + cell.max_gamma, + ) + slTravelAngleArray[cell.rowindex, cell.colindex] = max( + slTravelAngleArray[cell.rowindex, cell.colindex], cell.sl_gamma + ) - if forestInteraction: - if forestIntArray[cell.rowindex, cell.colindex] >= 0 and cell.forestIntCount >= 0: - forestIntArray[cell.rowindex, cell.colindex] = min( - forestIntArray[cell.rowindex, cell.colindex], cell.forestIntCount - ) - else: - forestIntArray[cell.rowindex, cell.colindex] = max( - forestIntArray[cell.rowindex, cell.colindex], cell.forestIntCount + if "travelLengthMax" in outputs or "travelLength" in outputs: + travelLengthMaxArray[cell.rowindex, cell.colindex] = max( + travelLengthMaxArray[cell.rowindex, cell.colindex], + cell.min_distance, ) + if "travelLengthMin" in outputs: + if travelLengthMinArray[cell.rowindex, cell.colindex] >= 0 and cell.min_distance >= 0: + travelLengthMinArray[cell.rowindex, cell.colindex] = min( + travelLengthMinArray[cell.rowindex, cell.colindex], + cell.min_distance, + ) + else: + 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) + + if forestInteraction: + if forestIntArray[cell.rowindex, cell.colindex] >= 0 and cell.forestIntCount >= 0: + forestIntArray[cell.rowindex, cell.colindex] = min( + forestIntArray[cell.rowindex, cell.colindex], + cell.forestIntCount, + ) + else: + forestIntArray[cell.rowindex, cell.colindex] = max( + forestIntArray[cell.rowindex, cell.colindex], + cell.forestIntCount, + ) if infraBool: # if 'infraBool' is True - i.e. calculation is performed with infrastructure information @@ -798,15 +1288,17 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): # if this is the case, then we exclude the affected release cell(s) from further processing and update # the row_list, col_list variables containing the release cells that should be processed release[zDeltaArray > 0] = 0 - row_list, col_list = get_start_idx(dem, release) + row_list, col_list = get_start_idx(dem, release, relIdArray, calcThalweg) - zDeltaPathList.append(zDeltaPathArray) - del cell_list, processedCells, zDeltaPathArray + if "zDeltaSum" in outputs: + zDeltaPathList.append(zDeltaPathArray) + del processedCells, zDeltaPathArray startcell_idx += 1 - for zDeltaPathArray in zDeltaPathList: - zDeltaSumArray += zDeltaPathArray + if "zDeltaSum" in outputs: + for zDeltaPathArray in zDeltaPathList: + zDeltaSumArray += zDeltaPathArray gc.collect() return ( @@ -1012,3 +1504,129 @@ def reverseTopology(topologyDict): reverseGraph[child].append(parentNode) return reverseGraph + + +def saveGenerationVideoData(genList, relId, dem, rasterAttributes, + outDir, variable="z_delta"): + """ + Saves, per generation, both a cumulative raster snapshot (history) and + a generation-only raster snapshot (just the cells of that generation), + plus the center-of-flux (thalweg) position, for a given relId. + Before saving, both the raster stacks and the row + coordinates are flipped along the row axis (upside down), so the saved + ``.npz`` data is already in the same up/down orientation as the + other com4FlowPy output rasters. + + + Parameters + ----------- + genList: list + list of cell lists (one list per generation), as built in calculation() + relId: int + release ID for which the data is saved + dem: np.array + DEM of the tile (used for the array shape) + rasterAttributes: dict + contains, among other things, "extentTile" for converting tile + coordinates to full-DEM coordinates + outDir: pathlib.Path + target directory + variable: str + "z_delta" or "flux" (attribute name of the Cell class) + + Returns + --------- + outFile: pathlib.Path + Path to the saved `video data file. + """ + ((startY, _), (startX, _)) = rasterAttributes["extentTile"] + nGen = len(genList) + nRows = dem.shape[0] + + frameStackHistory = np.zeros((nGen, *dem.shape), dtype=np.float32) # cumulative + frameStackCurrent = np.zeros((nGen, *dem.shape), dtype=np.float32) # this generation only + rowCoF = np.zeros(nGen, dtype=np.float32) + colCoF = np.zeros(nGen, dtype=np.float32) + rowCoE = np.zeros(nGen, dtype=np.float32) + colCoE = np.zeros(nGen, dtype=np.float32) + + runningMax = np.zeros_like(dem, dtype=np.float32) + + for gen, cellList in enumerate(genList): + if len(cellList) == 0: + frameStackHistory[gen] = runningMax + # the thalweg location is not moved compared to the generation before + rowCoF[gen] = rowCoF[gen - 1] if gen > 0 else np.nan + colCoF[gen] = colCoF[gen - 1] if gen > 0 else np.nan + rowCoE[gen] = rowCoE[gen - 1] if gen > 0 else np.nan + colCoE[gen] = colCoE[gen - 1] if gen > 0 else np.nan + continue + + rowArr = [] + colArr = [] + valArr = [] + fluxWeights = [] + energyWeights = [] + + for cell in cellList: + rowArr.append(cell.rowindex) + colArr.append(cell.colindex) + valArr.append(getattr(cell, variable)) + fluxWeights.append(cell.flux) + energyWeights.append(cell.flowEnergy) + + rowArr = np.asarray(rowArr) + colArr = np.asarray(colArr) + valArr = np.asarray(valArr) + fluxWeights = np.asarray(fluxWeights) + energyWeights = np.asarray(energyWeights) + + # generation-only snapshot (not cumulative) + currentArr = np.zeros_like(dem, dtype=np.float32) + currentArr[rowArr, colArr] = valArr + frameStackCurrent[gen] = currentArr + + # cumulative snapshot (history, including current generation) + runningMax[rowArr, colArr] = np.maximum(runningMax[rowArr, colArr], valArr) + frameStackHistory[gen] = runningMax + + # center of flux (weighted average position, weighted by flux) + if fluxWeights.sum() > 0: + rowCoF[gen] = np.average(rowArr, weights=fluxWeights) + colCoF[gen] = np.average(colArr, weights=fluxWeights) + else: + rowCoF[gen] = np.average(rowArr) + colCoF[gen] = np.average(colArr) + + # center of energy + if energyWeights.sum() > 0: + rowCoE[gen] = np.average(rowArr, weights=energyWeights) + colCoE[gen] = np.average(colArr, weights=energyWeights) + else: + rowCoE[gen] = np.average(rowArr) + colCoE[gen] = np.average(colArr) + + # --- flip everything from the simulation's internal (upside-down) + # row orientation to the standard output orientation, BEFORE applying + # the tile-to-full-DEM offset --------------------------------------- + frameStackHistory = np.flip(frameStackHistory, axis=1) + frameStackCurrent = np.flip(frameStackCurrent, axis=1) + + rowCoF = (nRows - 1) - rowCoF + rowCoE = (nRows - 1) - rowCoE + + rowCoFFull = rowCoF + startY + colCoFFull = colCoF + startX + rowCoEFull = rowCoE + startY + colCoEFull = colCoE + startX + + outDir = pathlib.Path(outDir) + outDir.mkdir(parents=True, exist_ok=True) + outFile = outDir / f"videoData_{variable}_{int(relId)}.npz" + np.savez_compressed(outFile, + framesHistory=frameStackHistory, + framesCurrent=frameStackCurrent, + rowCoF=rowCoFFull, colCoF=colCoFFull, + rowCoE=rowCoEFull, colCoE=colCoEFull, ) + log.info(f"Video data for relId {int(relId)} saved: {outFile}") + return outFile diff --git a/avaframe/com4FlowPy/flowPath.py b/avaframe/com4FlowPy/flowPath.py new file mode 100644 index 000000000..973618e4b --- /dev/null +++ b/avaframe/com4FlowPy/flowPath.py @@ -0,0 +1,388 @@ +import numpy as np +import pickle +import logging + +import avaframe.in3Utils.geoTrans as gT + +log = logging.getLogger(__name__) + + +class Path: + """Class contains a path, containing one startcell and corresponding child cells""" + + def __init__( + self, + dem, + startcellRow, + startcellCol, + genList, + rasterAttributes, + countArray, + relId=None, + listsRelId=None, + exampleCell=None, + ): + """initializes a GMF path, that belongs to a startcell + + Parameters + ---------- + dem: numpy array + Digital elevation model + startcellRow: int + Row index of startcell + startcellCol: int + Column index of startcell + genList: list + contains all cells that belong to the path (per generation an extra list) + rasterAttributes: dict + contains information about the input rasters + """ + self.dem = dem + self.cellsize = rasterAttributes["cellsize"] + self.nrows = rasterAttributes["nrows"] + self.rasterAttributes = rasterAttributes + + self.genList = genList + self.startcellRow = startcellRow + self.startcellCol = startcellCol + self.relId = int(relId) + self.pathRaster = np.where(countArray > 0, countArray, np.nan) + + if self.genList is None: + self.rowList = listsRelId["row"] + self.colList = listsRelId["col"] + self.fluxList = listsRelId["flux"] + self.zdeltaList = listsRelId["zdelta"] + self.travelLengthList = listsRelId["travelLengthMax"] + self.alpha = exampleCell.alpha + self.exp = exampleCell.exp + self.maxZDelta = exampleCell.max_z_delta + self.numberGen = len(self.rowList) + else: + self.alpha = genList[0][0].alpha + self.exp = genList[0][0].exp + self.maxZDelta = genList[0][0].max_z_delta + self.numberGen = len(genList) + + self.zDeltaGeneration = [] + self.fluxGeneration = [] + self.depFluxGeneration = [] + self.travelLengthGeneration = [] + self.flowEnergyGeneration = [] + self.rowGeneration = [] + self.colGeneration = [] + self.altitudeGeneration = [] + self.gammaGeneration = [] + self.flux_gen = [] + + self.zDeltaArray = np.zeros_like(self.dem, dtype=np.float32) + self.flowEnergyArray = np.zeros_like(self.dem, dtype=np.float32) + self.fluxArray = np.zeros_like(self.dem, dtype=np.float32) + self.routFluxSumArray = np.zeros_like(self.dem, dtype=np.float32) + self.depFluxSumArray = np.zeros_like(self.dem, dtype=np.float32) + + """ + self.travel_length_array = np.zeros_like(self.dem, dtype=np.float32) + self.generation_array = np.full_like(self.dem, np.nan, dtype=np.float32) + """ + + def getGenerationList(self, variable, generation=None): + """write lists with size and format of genList containing specific parameters + (the main list contains lists for every generation) + + Parameters + ----------- + variable: string + for the variable is the generation list created + generation: int + generation that is extracted (if None, all generations are added) + + Returns + ----------- + variableGeneration: list + contains all parameter values of a path (in generation structure) + """ + + variableGeneration = [] + if generation is None: + for cellList in self.genList: + listVariable = self.getListFromCellList(cellList, variable) + variableGeneration.append(listVariable) + else: + cellList = self.genList[generation] + variableGeneration = self.getListFromCellList(cellList, variable) + return variableGeneration + + def getListFromCellList(self, cellList, variable): + listVariable = [] + + for cell in cellList: + if variable == "zDelta": + listVariable.append(cell.z_delta) + elif variable == "flux": + listVariable.append(cell.flux) + elif variable in ["travelLength", "s"]: + listVariable.append(cell.min_distance) + elif variable in ["altitude", "z"]: + listVariable.append(cell.altitude) + elif variable == "row": + listVariable.append(cell.rowindex) + elif variable == "col": + listVariable.append(cell.colindex) + elif variable == "gamma": + listVariable.append(cell.max_gamma) + elif variable == "flowEnergy": + listVariable.append(cell.flowEnergy) + else: + log.error(f"variable {variable} can not be computed to a generation list") + return listVariable + + def getPathArrays(self): + """write arrays with size of DEM, containing the maximum of the variable values of every path + value 0 means, the path does not hit the cell + TODO: only calculate 'important'/output arrays + """ + for gen, cellList in enumerate(self.genList): + for cell in cellList: + self.zDeltaArray[cell.rowindex, cell.colindex] = max( + self.zDeltaArray[cell.rowindex, cell.colindex], cell.z_delta + ) + self.flowEnergyArray[cell.rowindex, cell.colindex] = max( + self.flowEnergyArray[cell.rowindex, cell.colindex], cell.flowEnergy + ) + self.fluxArray[cell.rowindex, cell.colindex] = max( + self.fluxArray[cell.rowindex, cell.colindex], cell.flux + ) + self.routFluxSumArray[cell.rowindex, cell.colindex] += cell.flux + self.depFluxSumArray[cell.rowindex, cell.colindex] += cell.fluxDep + + """ + self.travel_length_array[cell.rowindex, cell.colindex] = max(self.travel_length_array[cell.rowindex, cell.colindex], cell.min_distance) + self.generation_array[cell.rowindex, cell.colindex] = gen + """ + + def calcThalwegCenterof(self, variable, variableCo): + """calculates for a specific variable the center of a specific variable (thalweg) + + Parameters + ---------- + variable: list + variable, which is centered (in format genList) + variableCo: list + center of variableCo is calculated (variable is weighted) (in format genList) + + Returns + ---------- + variableSum: numpy array + sum of variable per generation + coVar: numpy array + centered variable (per generation) + """ + + coVar = np.zeros(self.numberGen) + variableSum = np.zeros(self.numberGen) + for gen in range(0, self.numberGen): + var = np.array(variable[gen]) + co = np.array(variableCo[gen]) + variableSum[gen] = np.sum(var) + variableCoSum = np.sum(co) + if variableCoSum > 0: # flow_energy and zdelta are 0 in generation 0 + # coVar[gen] = 1 / variableCoSum * np.sum(var * co) + coVar[gen] = np.average(var, weights=co) + else: + # TODO: does this makes sense?? + coVar[gen] = np.average(var) + return variableSum, coVar + + def getCenterofs(self, variables, centerOfs): + """ + calculate sum of variable for every iteration step/ generation and + center of energy, flux and zDelta for the following variables: + + Parameters + ---------- + variables: list + List of variables that should be weighted (with center of energy and flux) + """ + + # self.getVariablesGeneration() + + for varName in variables: + if varName in [ + "s", + "z", + "x", + "y", + "flowEnergyArray", + "zDeltaArray", + "fluxArray", + "routFluxSumArray", + "depFluxSumArray", + ]: + continue + if varName == "depFluxSum": + variables.append("depFlux") + continue + if varName == "fluxSum": + variables.append("flux") + continue + + values = self.getGenerationList(varName) + + if "CoE" in centerOfs: + self.energyGenList = self.getGenerationList("flowEnergy") + sumE, coE = self.calcThalwegCenterof(values, self.energyGenList) + # TODO: zdelta is 0 in generation 1, so the first value does not make sense / + # -> now last indices are deleted in postprocessing + setattr(self, f"{varName}CoE", coE) + if "CoF" in centerOfs: + self.fluxGenList = self.getGenerationList("flux") + sumF, coF = self.calcThalwegCenterof(values, self.fluxGenList) + setattr(self, f"{varName}CoF", coF) + if "CoZd" in centerOfs: + self.zDeltaGenList = self.getGenerationList("zDelta") + sumZd, coZd = self.calcThalwegCenterof(values, self.zDeltaGenList) + setattr(self, f"{varName}CoZd", coZd) + + def correctIndicesTile(self, row, col): + """ + correct row and col from the tile to the whole DEM extent + + Parameters + -------------- + row: numpy array + row in the tile + col: numpy array + col in the tile + + Returns + ------------- + rowLarge: numpy array + row in the whole DEM extent + colLarge: numpy array + col in the whole DEM extent + """ + ((sY, _), (sX, _)) = self.rasterAttributes["extentTile"] + + rowLarge = row + sY + colLarge = col + sX + + return (rowLarge, colLarge) + + def saveDict(self, saveDir, centerOfs, variables): + """ + save thalweg data. (One file per thalweg) + + Parameters + ------------ + saveDir: pathlib.PosixPath + directory, in which the thalweg data is saved + centerOfs: list + contains the center-of-variable names that are saved + variables: list + contains the variable names that are saved + """ + + thalwegData = { + "alpha": round(self.alpha, 1), + "exponent": self.exp, + "zDeltaMax": round(self.maxZDelta, 1), + # 'crs': self.crs, + "numberGen": self.numberGen, + } + variables = variables + centerOfs = centerOfs + + for co in centerOfs: + for varName in variables: + if varName in [ + "flowEnergyArray", + "zDeltaArray", + "fluxArray", + "routFluxSumArray", + "depFluxSumArray", + ]: + if np.any(getattr(self, f"{varName}")) is False: + self.getPathArrays() + value = getattr(self, f"{varName}") + elif varName == "z": + value = getattr(self, f"altitude{co}") + elif varName == "s": + value = getattr(self, f"travelLength{co}") + else: + value = getattr(self, f"{varName}{co}") + thalwegData[f"{varName}"] = value + + # output file name and save teh pickle file + if self.relId is None: + outName = f"thalwegData_{co}_{self.startcellRow}_{self.startcellCol}.pickle" + else: + outName = f"thalwegData_{co}_{self.relId}.pickle" + with open(saveDir / (outName), "wb") as handle: + pickle.dump(thalwegData, handle, protocol=pickle.HIGHEST_PROTOCOL) + + def calcAndSaveThalwegData(self, thalwegParameters): + """main function for paths & thalwegs: calculates the thalweg and saves the data + + Parameters: + ------------ + thalwegParameters: dict + contains information to calculate and save the thalweg data (from .ini file) + """ + saveDir = thalwegParameters["thalwegDir"] + if thalwegParameters["thalwegSaveRam"]: + # only compute thalweg location for coF + # TODO: do we only want to compute coF or also coE and coZd? + variables = ["x", "y", "travelLength", "zdelta"] + cos = ["cof"] + else: + cos = eval(thalwegParameters["thalwegCenterOf"]) + variables = eval(thalwegParameters["thalwegVariables"]) + centerOfs = [] + + for co in cos: + co.lower() + if co in ["energy", "coe"]: + centerOf = "CoE" + elif co in ["flux", "cof"]: + centerOf = "CoF" + elif co in ["zdelta", "cozd"]: + centerOf = "CoZd" + else: + message = f"{co} is a not valid thalweg parameter" + log.error(message) + raise ValueError(message) + centerOfs.append(centerOf) + + if "s" in variables: + variables.append("travelLength") + if "z" in variables: + variables.append("altitude") + if "x" in variables or "y" in variables: + variables.append("col") + variables.append("row") + + if thalwegParameters["thalwegSaveRam"]: + _, self.colCoF = self.calcThalwegCenterof(self.colList, self.fluxList) + _, self.rowCoF = self.calcThalwegCenterof(self.rowList, self.fluxList) + _, self.zdeltaCoF = self.calcThalwegCenterof(self.zdeltaList, self.fluxList) + _, self.travelLengthCoF = self.calcThalwegCenterof(self.travelLengthList, self.fluxList) + else: + self.getCenterofs(variables, centerOfs) + # empty generation list to safe RAM + self.genList = [] + for co in centerOfs: + # convert column and row to coordinates s, y + # TODO: when there is more than one tile, think if the other outputs need to be corrected?? + colCentered = getattr(self, f"col{co}") + rowCentered = getattr(self, f"row{co}") + rowLarge, colLarge = self.correctIndicesTile(rowCentered, colCentered) + + x, y = gT.indicesToCoords(colLarge, rowLarge, self.rasterAttributes) + setattr(self, f"x{co}", x) + setattr(self, f"y{co}", y) + setattr(self, f"col{co}", colLarge) + setattr(self, f"row{co}", rowLarge) + # update y coordinate + self.saveDict(saveDir, centerOfs, variables) + log.debug(f"thalweg data saved in {saveDir}") diff --git a/avaframe/com4FlowPy/splitAndMerge.py b/avaframe/com4FlowPy/splitAndMerge.py index 3158e75ef..08212737f 100644 --- a/avaframe/com4FlowPy/splitAndMerge.py +++ b/avaframe/com4FlowPy/splitAndMerge.py @@ -178,6 +178,252 @@ def tileRaster(fNameIn, fNameOut, dirName, xDim, yDim, U, isInit=False): # return largeRaster +def tileRasterWithIndices(fNameIn, fNameOut, dirName, exList, eyList, U, isInit=False): + """ + divides a raster into tiles and saves the tiles + the tile size is determined by the indices of the tile ends and can vary within the tiles. + + Parameters + ----------- + fNameIn : str + path to raster that is tiled + fNameOut: str + name of saved raster file + dirName: str + path to folder, where tiled raster is saved (temp - folder) + exList: list + contains end indices of tiles in x dimension (number of raster columns) + eyList: list + contains end indices of tiles in y dimension (number of raster columns) + U: int + size of tile overlapping (number of raster cells) + isInit: bool + if isInit is True, edges are assigned to -9999 (default: False) + """ + + log.info("tile rasters using ends!") + + # largeRaster, largeHeader = iof.f_readASC(fNameIn, dType='float') + largeData = IOf.readRaster(fNameIn, noDataToNan=False) + largeRaster = largeData["rasterData"] + + i, j, imax, jmax = 0, 0, 0, 0 + sX, sY, eX, eY = 0, 0, 0, 0 + + JMAX = len(exList) - 1 + IMAX = len(eyList) - 1 + + if isInit is False: + for eY in eyList: + for eX in exList: + rangeRowsCols = ((sY, eY), (sX, eX)) + pickle.dump(rangeRowsCols, open(dirName / ("ext_%i_%i" % (i, j)), "wb")) + + np.save(dirName / ("%s_%i_%i" % (fNameOut, i, j)), largeRaster[sY:eY, sX:eX]) + log.info("saved %s - TileNr.: %i_%i", fNameOut, i, j) + + sX = eX - 2 * U + jmax = max(j, jmax) + j += 1 + sX, j, eX = 0, 0, 0 + sY = eY - 2 * U + imax = max(i, imax) + i += 1 + else: + for eY in eyList: + for eX in exList: + rangeRowsCols = ((sY, eY), (sX, eX)) + pickle.dump(rangeRowsCols, open(dirName / ("ext_%i_%i" % (i, j)), "wb")) + + initRas = largeRaster[sY:eY, sX:eX].copy() + if j != JMAX: + initRas[:, -U:] = -9999 # Rand im Osten + if i != 0: + initRas[0:U, :] = -9999 # Rand im Norden + if j != 0: + initRas[:, 0:U] = -9999 # Rand im Westen + if i != IMAX: + initRas[-U:, :] = -9999 # Rand im Sueden + + np.save(dirName / ("%s_%i_%i" % (fNameOut, i, j)), initRas) + del initRas + log.info("saved %s - TileNr.: %i_%i", fNameOut, i, j) + + sX = eX - 2 * U + jmax = max(j, jmax) + j += 1 + sX, j, eX = 0, 0, 0 + sY = eY - 2 * U + imax = max(i, imax) + i += 1 + + pickle.dump((imax, jmax), open(dirName / "nTiles", "wb")) + log.info("finished tiling %s: nTiles=%s" % (fNameOut, (imax + 1) * (jmax + 1))) + log.info("----------------------------") + del largeRaster + gc.collect() + + +def getTileEnds(dirName, xDim, yDim, U, relIdRaster): + """ + computes the indices of tiling ends for a raster. + The tiles can not split a release cells with the same release ID. + + Parameters + ----------- + dirName: str + path to folder, where tiled raster is saved (temp - folder) + xDim: int + minimum size of one tile in x dimension (number of raster columns) + yDim: int + minimum size of one tile in y dimension (number of raster rows) + U: int + size of tile overlapping (number of raster cells) + + Returns + ------------- + exList: list + contains end x - indices of the tiles + eyList: list + contains end y - indices of the tiles + """ + log.info("get tile ends") + i, j, imax, jmax = 0, 0, 0, 0 + sX, sY, eX, eY = 0, 0, 0, 0 + + nrows, ncols = relIdRaster.shape[0], relIdRaster.shape[1] + pickle.dump((nrows, ncols), open(dirName / "extentLarge", "wb")) + + I, J, IMAX, JMAX = 0, 0, 0, 0 + + # compute the maximum or tiles in x and y direction are possible + while eY < nrows: + eY = sY + yDim + while eX < ncols: + eX = sX + xDim + + sX = eX - 2 * U + JMAX = max(J, JMAX) + J += 1 + sX, J, eX = 0, 0, 0 + sY = eY - 2 * U + IMAX = max(I, IMAX) + I += 1 + + sX, sY, eX, eY = 0, 0, 0, 0 + exList = [] + eyList = [] + i, j, imax, jmax = 0, 0, 0, 0 + + while eY < nrows: + eY = sY + yDim + shiftCountY = 0 + # iterate as long as necessary that the y - end of the tile does not cut a PRA + while True: + # TODO: loop stops if ey = nrows?? + if i == IMAX: + # the border tile + eyList.append(eY) + break + + relIdSY = 0 + relIdEY = eY - U + if j != 0: + relIdSY = sY + U + + mask = np.zeros_like(relIdRaster) + mask[relIdSY:relIdEY, :] = 1 + idsIn, idsOut = getMaskedRasters(mask, relIdRaster) + + if np.any(np.isin(idsIn, idsOut)): + shiftCountY += 1 + eY += 1 + + # relIdEY is relIdSY in the next iteration, so we only need to search for the end indices + else: + eyList.append(eY) + break + if shiftCountY > 0: + log.info( + f"tiling would divide a release area, tiling size (in y direction) changed by {shiftCountY} cells." + ) + + sY = eY - 2 * U + imax = max(i, imax) + i += 1 + + while eX < ncols: + eX = sX + xDim + + shiftCountX = 0 + + while True: + # iterate as long as necessary that the x - end of the tile does not cut a PRA + if j == JMAX: + # the border tile + exList.append(eX) + break + relIdSX = 0 + relIdEX = eX - U + if j != 0: + relIdSX = sX + U + + mask = np.zeros_like(relIdRaster) + mask[:, relIdSX:relIdEX] = 1 + idsIn, idsOut = getMaskedRasters(mask, relIdRaster) + + if np.any(np.isin(idsIn, idsOut)): + shiftCountX += 1 + eX += 1 + else: + exList.append(eX) + break + + if shiftCountX > 0: + log.info( + f"tiling would divide a release area, tiling size (in x direction) changed by {shiftCountX} cells." + ) + sX = eX - 2 * U + jmax = max(j, jmax) + j += 1 + + pickle.dump((imax, jmax), open(dirName / "nTiles", "wb")) + log.info("nTiles=%s" % ((imax + 1) * (jmax + 1))) + log.info("----------------------------") + gc.collect() + return exList, eyList + + +def getMaskedRasters(mask, raster): + """ + get unique values of a raster within and outside a mask + + Parameters + ----------- + mask : numpy array + mask that contains 1 (inside) and 0 (outside) + raster : numpy array + raster values + + Returns + ------------- + idsIn: numpy array + unique values of raster inside mask + idsOut: numpy array + unique values of raster outside mask + + """ + + relIdInside = np.where(mask == 1, raster, -9999) + relIdOutside = np.where(mask == 0, raster, -9999) + + idsIn = np.unique(relIdInside) + idsOut = np.unique(relIdOutside) + idsIn = idsIn[idsIn > 0] + idsOut = idsOut[idsOut > 0] + return idsIn, idsOut + + def mergeRaster(inDirPath, fName, method="max"): """ Merges the results for each tile to one array using the @@ -230,8 +476,14 @@ def mergeRaster(inDirPath, fName, method="max"): 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), + 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( @@ -268,7 +520,10 @@ def mergeDict(inDirPath, fName): smallDict = pickle.load(file) if bool(smallDict): for cellindSmall in smallDict: - cellind = (cellindSmall[0] + pos[0][0], cellindSmall[1] + pos[1][0]) + cellind = ( + cellindSmall[0] + pos[0][0], + cellindSmall[1] + pos[1][0], + ) if cellind in mergedDict: mergedDict[cellind] = np.append(smallDict[cellindSmall], mergedDict[cellind]) else: diff --git a/avaframe/in3Utils/geoTrans.py b/avaframe/in3Utils/geoTrans.py index f3c0428fe..521c4b789 100644 --- a/avaframe/in3Utils/geoTrans.py +++ b/avaframe/in3Utils/geoTrans.py @@ -291,7 +291,11 @@ def remeshData(rasterDict, cellSizeNew, remeshOption="griddata", interpMethod="c yGrid = yGrid[mask] z = zCopy[mask] zNew = sp.interpolate.griddata( - (xGrid, yGrid), z, (xGridNew, yGridNew), method=interpMethod, fill_value=header["nodata_value"] + (xGrid, yGrid), + z, + (xGridNew, yGridNew), + method=interpMethod, + fill_value=header["nodata_value"], ) elif remeshOption == "RectBivariateSpline": if np.isnan(z).any(): @@ -482,7 +486,10 @@ def remeshRaster(rasterFile, cfgSim, typeIndicator="DEM", onlySearch=False, lega else: log.info("Using rasterio resampling") remeshedRaster = remeshDataRio( - rasterFile, cszRasterNew, cfgSim["GENERAL"]["remeshInterpMethod"], larger=False + rasterFile, + cszRasterNew, + cfgSim["GENERAL"]["remeshInterpMethod"], + larger=False, ) flipArg = False @@ -1591,7 +1598,13 @@ def rotateRaster(rasterDict, theta, deg=True): # project data on this new grid rotatedZ, _ = projectOnGrid( - xTheta, yTheta, rasterDict["rasterData"], csz=csz, xllc=xllc, yllc=yllc, interp="bilinear" + xTheta, + yTheta, + rasterDict["rasterData"], + csz=csz, + xllc=xllc, + yllc=yllc, + interp="bilinear", ) rotatedRaster = {"header": header, "rasterData": rotatedZ} @@ -2224,3 +2237,33 @@ def interpolateLineLinear(lineDict, distance): lineDict["y"] = Y1 return lineDict + + +def indicesToCoords(col, row, header): + """ + transform indeces (row and col) to coordinates (x, y) + considering the flipped (upside down) rasters + TODO: is there already a function for this calculation?? + + Parameters + ---------- + col: numpy array + column indices of cells belonging to path + row: numpy array + row indices of cells belonging to path + header: dict + header of a raster + + Returns + ---------- + x: numpy array + x coordinates (in m) of cells belonging to path + y: numpy array + y coordinates (in m) of cells belonging to path + """ + cellsize = header["cellsize"] + xllcorner = header["xllcenter"] - cellsize / 2 + yllcorner = header["yllcenter"] - cellsize / 2 + x = xllcorner + col * cellsize + y = yllcorner + row * cellsize + return x, y diff --git a/avaframe/out3Plot/outCom4Gif.py b/avaframe/out3Plot/outCom4Gif.py new file mode 100644 index 000000000..baa9b5ea7 --- /dev/null +++ b/avaframe/out3Plot/outCom4Gif.py @@ -0,0 +1,744 @@ +""" functions to create a GIF for com4FlowPys generation data""" + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.animation as animation +import matplotlib.colors as mcolors +import matplotlib.patheffects as pe +from avaframe.runScripts.runComputeDist import outFile +from matplotlib.gridspec import GridSpec +from matplotlib.collections import LineCollection +import logging +import os +import pathlib +from matplotlib.colors import LightSource +from PIL import Image +from cmcrameri import cm as cmapCrameri + +import avaframe.in2Trans.rasterUtils as rasterUtils +from avaframe.in3Utils import cfgUtils +import avaframe.in1Data.getInput as getInput +import avaframe.in3Utils.fileHandlerUtils as fU +import avaframe.in3Utils.geoTrans as gT +import pickle + +import avaframe.out3Plot.outCom4Gif as outCom4Gif + +log = logging.getLogger(__name__) + + + +def addHillshadeSimple(dem, cellSize, cfgHS): + """ + Computes a hillshade array from the DEM + + Parameters + ----------- + dem: numpy array + DEM raster data + cellSize: float + cellsize of DEM + cfgHS: configparser object + configuration for hillshade + """ + lightSource = LightSource(azdeg=cfgHS.getfloat("azimuth"), altdeg=cfgHS.getfloat("altitude")) + hillshade = lightSource.hillshade(dem, vert_exag=cfgHS.getfloat("vertExag"), dx=cellSize, dy=cellSize) + return hillshade + + +def searchResFolder(avalanchedir, module="com4FlowPy"): + """ + search for a com4Flowpy result folder, return the respective simhash if only one exists + + Parameters + ---------- + avalanchedir: str or pathlib.Path + path to avalanche project + module: str + module name (default: com4FlowPy) + + Returns + ------- + oneResFolder: bool + True if exactly one result folder exists, otherwise False + simHash: str + the simhash of the simulation, if one result folder exists + """ + dir = avalanchedir / "Outputs" / module / "peakFiles" + + resFolders = [] + for filename in os.listdir(dir): + if filename.startswith("res_"): + resFolders.append(filename) + if len(resFolders) == 0: + message = f"No results {module} folder found in {dir}" + log.error(message) + raise FileNotFoundError(message) + elif len(resFolders) > 1: + oneResFolder = False + simHash = "" + elif len(resFolders) == 1: + oneResFolder = True + simHash = str(resFolders[0]).split("res_", 1)[1] + return oneResFolder, simHash + + +def removeContour(contourHolder): + """ + Removes the previously drawn contour artist (if any), so a new one + can be drawn for the current frame without accumulating old contours. + + Parameters + ----------- + contourHolder: dict + mutable dict with key "artist", holding the current contour artist + (or None). Passed by reference so state persists across calls. + """ + if contourHolder["artist"] is not None: + try: + contourHolder["artist"].remove() + except AttributeError: + for coll in contourHolder["artist"].collections: + coll.remove() + contourHolder["artist"] = None + return contourHolder + + +def figureToRgbArray(fig): + """ + Render the current state of a matplotlib figure to an RGB numpy array. + + Parameters + ---------- + fig : matplotlib.figure.Figure + Figure to render (must already have been drawn/updated). + + Returns + ------- + numpy.ndarray + RGB image array of shape (height, width, 3), dtype uint8. + """ + fig.canvas.draw() + buf = np.asarray(fig.canvas.buffer_rgba()) + return buf[..., :3].copy() + + +def buildGlobalPalette(frameArrays, colors=256, sampleEveryFrame=5, sampleEveryPixel=8): + """ + Build a single shared color palette from a subsample of frames and pixels. + + Using one global palette (instead of letting PIL quantize each GIF + frame independently) prevents visible color flicker between frames, + since every frame is then mapped onto the exact same set of colors. + + Both frames and pixels within each frame are subsampled before + building the palette, since using all pixels of many full-resolution + frames can require more contiguous memory than PIL can allocate for + a single quantization image. + + Parameters + ---------- + frameArrays : list of numpy.ndarray + List of RGB frame arrays, each of shape (height, width, 3). + colors : int, optional + Number of palette colors (GIF hard limit is 256). Default 256. + sampleEveryFrame : int, optional + Use only every Nth frame to build the palette. Default 5. + sampleEveryPixel : int, optional + Use only every Nth pixel (in both height and width) within each + sampled frame. Default 8. + + Returns + ------- + paletteImg: PIL.Image.Image + A palette-mode image whose palette should be reused for all + frames via ``Image.quantize(palette=...)``. + """ + sampledFrames = frameArrays[::max(1, sampleEveryFrame)] + + # subsample pixels within each frame too, so the total number of + # pixels stays small regardless of frame resolution or frame count + pixelSamples = [ + frame[::sampleEveryPixel, ::sampleEveryPixel, :].reshape(-1, 3) + for frame in sampledFrames + ] + stacked = np.concatenate(pixelSamples, axis=0) + + # reshape into a reasonably square-ish image instead of a single row, + # which avoids extremely wide (1, N) images that some PIL/C backends + # fail to allocate contiguous memory for + nPixels = stacked.shape[0] + width = int(np.ceil(np.sqrt(nPixels))) + height = int(np.ceil(nPixels / width)) + padded = np.zeros((height * width, 3), dtype=np.uint8) + padded[:nPixels] = stacked + + paletteSourceImg = Image.fromarray(padded.reshape(height, width, 3), mode="RGB") + paletteImg = paletteSourceImg.quantize(colors=colors, method=Image.MEDIANCUT) + return paletteImg + + +def loadExtendedProfile(extendedProfilePicklePath): + """ + Load the extended thalweg profile (avaProfileMass) from a pickle file. + + Parameters + ---------- + extendedProfilePicklePath : str or pathlib.Path + Path to the pickle file containing the extended profile dict, + with keys "x", "y", "indStartMassAverage", "indEndMassAverage" + (same structure as used in avalancheThalwegPlot). + + Returns + ------- + dict + The loaded avaProfileMass dictionary. + """ + with open(extendedProfilePicklePath, "rb") as handle: + avaProfileMass = pickle.load(handle) + return avaProfileMass + + +def coordsToRowCol(x, y, rasterHeader): + """ + Convert real-world x/y coordinates to raster row/col array indices. + + Assumes the standard AvaFrame raster convention: row 0 corresponds to + the northernmost (top) row, col 0 to the westernmost (left) column, + matching how ``imshow`` displays the raster arrays elsewhere in this + module (no ``extent`` set, i.e. plotted in pixel/array index space). + + Parameters + ---------- + x, y : numpy.ndarray or float + Real-world coordinates. + rasterHeader : dict + Raster header as returned by ``rasterUtils.readRaster``, must + contain "xllcenter", "yllcenter", "cellsize" and "nrows". + + Returns + ------- + col, row : numpy.ndarray or float + Corresponding array indices (float, not rounded, so subpixel + positions along the profile are preserved for smooth plotting). + """ + cellSize = rasterHeader["cellsize"] + col = (x - rasterHeader["xllcenter"]) / cellSize + row = rasterHeader["nrows"] - 1 - (y - rasterHeader["yllcenter"]) / cellSize + return col, row + + +def addThalwegExtension(ax, avaProfileMass, rasterHeader): + """ + Plot the top and bottom thalweg extensions plus the center-of-mass + path onto an existing axes, in the same row/col pixel coordinate + system used by the generation video's imshow layers. + + Mirrors the styling used in avaframe.out3Plot.outCom1DFA's + avalancheThalwegPlot (colored outline via path_effects), so the + final GIF frame looks consistent with AvaFrame's standard thalweg + plots. + + Parameters + ---------- + ax : matplotlib.axes.Axes + Axes to draw on (same axes as the animated raster layers). + avaProfileMass : dict + Extended profile data with keys "x", "y", + "indStartMassAverage", "indEndMassAverage". + rasterHeader : dict + Raster header used to convert avaProfileMass["x"]/["y"] into + row/col indices via coordsToRowCol. + + Returns + ------- + None + """ + indStart = avaProfileMass["indStartMassAverage"] + indEnd = avaProfileMass["indEndMassAverage"] + + colProfile, rowProfile = coordsToRowCol( + avaProfileMass["x"], avaProfileMass["y"], rasterHeader + ) + + topLine, = ax.plot( + colProfile[: indStart + 1], + rowProfile[: indStart + 1], + "-b.", + # color="blue", + zorder=20, + label="top extension", + lw=2.5, + path_effects=[pe.Stroke(linewidth=2, foreground="b"), pe.Normal()], + ) + bottomLine, = ax.plot( + colProfile[indEnd:], + rowProfile[indEnd:], + "-g.", + # color="green", + zorder=20, + label="bottom extension", + lw=2.5, + path_effects=[pe.Stroke(linewidth=2, foreground="g"), pe.Normal()], + ) + centerLine, = ax.plot( + colProfile[indStart: indEnd + 1], + rowProfile[indStart: indEnd + 1], + "-k.", + # color="black", + zorder=20, + label="center of flux path", + lw=2.5, + path_effects=[pe.Stroke(linewidth=2, foreground="k"), pe.Normal()], + ) + return topLine, bottomLine, centerLine + + +def getProfileData(dem, profilePickleData): + """ + Compute travel distance (s), elevation (z) and get energy-height + (zDelta) values per generation. + + Parameters + ---------- + dem : dict + DEM dictionary containing header and the raster data + profilePickleData : dict + Thalweg data loaded from the pickle file, must contain the x, y and zDelta values + per generation/iteration, (in real-world coordinates.) + + Returns + ------- + profile: dict + contains numpy arrays with: + x: x-coordinates of thalweg + y: y-coordinates of thalweg + zDelta: zDelta values for each thalweg + s: Cumulative travel distance along the (x, y) trajectory, one + value per generation, starting at 0. + z: DEM elevation at each (x, y) position, one value per generation, + obtained via bilinear interpolation on the DEM. + """ + x = np.asarray(profilePickleData["x"], dtype=np.float64) + y = np.asarray(profilePickleData["y"], dtype=np.float64) + zDelta = np.asarray(profilePickleData["zDelta"], dtype=np.float32) + + # cumulative travel distance: s[0] = 0, then running sum of + # step-wise Euclidean distances between consecutive (x, y) points + dx = np.diff(x) + dy = np.diff(y) + stepDist = np.sqrt(dx ** 2 + dy ** 2) + s = np.concatenate([[0.0], np.cumsum(stepDist)]).astype(np.float32) + + # project (x, y) onto the DEM to get elevation z, via bilinear interpolation + points = {"x": x, "y": y} + points, ioob = gT.projectOnRaster(dem, points, interp="bilinear") + if ioob > 0: + log.warning(f"{ioob} thalweg points were out of bounds of the DEM during elevation projection") + z = points["z"].astype(np.float32) + + profile = {"x": x, "y": y, "z": z, "s": s, "zDelta": zDelta} + + return profile + + +def setup2DProfileAxis(ax2, sFull, zFull, zDeltaFull, sharedCmap, norm): + """ + Set up a 2D profile axis (elevation and z + zDelta vs. + travel distance). + + Parameters + ---------- + ax2 : matplotlib.axes.Axes + Axes to set up. + sFull, zFull, zDeltaFull : numpy.ndarray + travel distance, elevation and zDelta arrays that are plotted. + + Returns + ------- + zLine, zVelLine : matplotlib.lines.Line2D + Line artists for the elevation profile and the z + zDelta + ("velocity altitude") profile, to be updated per frame. + currentProfilePoint : matplotlib.lines.Line2D + Marker artist for the current generation's position on the + profile. + """ + # set general axis limits + validMask = ~np.isnan(zFull) + sMargin = 0.05 * (np.nanmax(sFull) - np.nanmin(sFull) + 1e-6) + zMin = np.nanmin(zFull[validMask]) + zMax = np.nanmax((zFull + zDeltaFull)[validMask]) + zMargin = 0.05 * (zMax - zMin + 1e-6) + + ax2.set_xlim(np.nanmin(sFull) - sMargin, np.nanmax(sFull) + sMargin) + ax2.set_ylim(zMin - zMargin, zMax + zMargin) + ax2.set_xlabel("Travel distance (horizontally projected) [m]") + ax2.set_ylabel("Elevation [m]") + # ax2.set_title("Thalweg elevation and energy line") + + ax2.hlines(zFull[validMask][-1], 0, sFull[-1], + colors="grey", linestyles="dotted", linewidths=1, + label="Runout length") + ax2.vlines(0, zFull[validMask][-1], zFull[validMask][0], + colors="grey", linestyles="dashed", linewidths=1, + label="Elevation drop") + (zLine,) = ax2.plot([], [], color="black", lw=2, label="Elevation z") + zVelDummy = ax2.plot([], [], color="blue", lw=2, label="Energy line height (+ z) \n(color indicates velocity)") + + segments = np.empty((0, 2, 2)) + + (currentProfilePoint,) = ax2.plot([], [], "o", color="k", markersize=8, + markeredgecolor="white", markeredgewidth=1.5, + zorder=20) + ax2.legend(loc="upper right") + zVelLine = LineCollection( + segments, + cmap=sharedCmap, + norm=norm, + linewidth=3, + ) + ax2.add_collection(zVelLine) + + return zLine, zVelLine, currentProfilePoint + + +def makeGenerationVideo(cfg=None, avalancheDir="", ax=None): + """ + Create an animated GIF of com4FlowPy generation data for one release ID. + The animation is saved as a GIF file to disk + + For a single release ID, renders one frame per generation + (or per group of generations, see ``gensPerFrame``) showing: + + - a static hillshade of the DEM as background, + - the cumulative raster values of all previous generations + - the raster values of the current generation only and outlined with a contour line, + - the center-of-flux or center-of-energy trajectory ("thalweg") + accumulated up to the current generation, plus a marker at the + current position. + + Parameters + ---------- + cfg : configparser.ConfigParser, optional + Full configuration object for the GIF module. + avalancheDir : str or pathlib.Path, optional + Path to the avalanche project directory. + """ + if cfg is None: + cfg = cfgUtils.getModuleConfig(outCom4Gif) + cfgHS = cfg["HILLSHADE"] + cfgGen = cfg["GENERAL"] + cfgPath = cfg["PATH"] + + showProfilePanel = cfgGen.getboolean("show2DProfilePanel", fallback=True) + + if avalancheDir == "": + # Load avalanche directory from general configuration file + cfgMain = cfgUtils.getGeneralConfig() + avalancheDir = cfgMain["MAIN"]["avalancheDir"] + avalancheDir = pathlib.Path(avalancheDir) + + if cfgPath.get("demPath", fallback="") == "" or cfgPath.get("demPath", fallback="") is None: + demFile = getInput.getDEMPath(avalancheDir) + else: + demFile = cfgPath["demPath"] + + resFolderUnique, simhash = searchResFolder(avalancheDir, module="com4FlowPy") + + if resFolderUnique: + resHash = simhash + else: + resHash = cfgGen.get("simhash") + if resHash == "" or resHash is None: + message = "Please provide a valid simhash to simulation results." + log.error(message) + raise ValueError(message) + + outDir = avalancheDir / "Outputs" / "com4FlowPy" + thalwegDataDir = outDir / "peakFiles" / f"res_{resHash}" / "thalwegData" + videoDataDir = thalwegDataDir / "videoData" + + videoDataVariable = cfgGen.get("videoDataVariable") + relId = cfgGen.get("relId") + centerOf = cfgGen.get("centerOf") + + if videoDataVariable == "velocity": + videoDataVariableFile = "z_delta" + else: + videoDataVariableFile = videoDataVariable + + npzFile = videoDataDir / f"videoData_{videoDataVariableFile}_{str(relId)}.npz" + + if showProfilePanel: + outFile2D = "2D" + else: + outFile2D = "" + + if cfgPath.get("outVideoPath", fallback="") == "" or cfgPath.get("outVideoPath", fallback="") is None: + videoOutputDir = outDir / "reports" + outFile = videoOutputDir / f"videoData{outFile2D}_{resHash}_{videoDataVariable}_{str(relId)}_{centerOf}.gif" + else: + outFile = cfgPath["outVideoPath"] + + fU.makeADir(videoOutputDir) + + clabelDict = {"flux": "flux", "z_delta": "energy line height zDelta [m]", "min_distance": "travel length [m]", + "max_gamma": "tavel angle [°]", "velocity": "velocity [m/s]", } + + fps = cfgGen.getint("fps") + gensPerFrame = cfgGen.getint("gensPerFrame") + + data = np.load(npzFile) + frameStackHistory = data["framesHistory"] + frameStackCurrent = data["framesCurrent"] + + if cfgGen.get("videoDataVariable") == "velocity": + # compute Zdelta to velocity + for gen, (histArr, currArr) in enumerate(zip(frameStackHistory, frameStackCurrent)): + frameStackCurrent[gen] = (currArr * 2 * 9.81) ** 0.5 + frameStackHistory[gen] = (histArr * 2 * 9.81) ** 0.5 + + demDict = rasterUtils.readRaster(demFile) + + if centerOf.lower() == "coe": + coRow = data["rowCoE"] + coCol = data["colCoE"] + centerOfPickle = "CoE" + else: + coRow = data["rowCoF"] + coCol = data["colCoF"] + centerOfPickle = "CoF" + + if gensPerFrame > 1: + nFramesOut = int(np.ceil(len(frameStackHistory) / gensPerFrame)) + frameStackHistory = frameStackHistory[::gensPerFrame][:nFramesOut] + frameStackCurrent = frameStackCurrent[::gensPerFrame][:nFramesOut] + coRow = coRow[::gensPerFrame][:nFramesOut] + coCol = coCol[::gensPerFrame][:nFramesOut] + + if showProfilePanel: + # read thalweg data and also use subset (generation per frame) + profilePicklePath = thalwegDataDir / f"thalwegData_{centerOfPickle}_{relId}.pickle" + thalwegData = np.load(profilePicklePath, allow_pickle=True) + + profile = getProfileData( + demDict, thalwegData + ) + + if gensPerFrame > 1: + sProfile = profile["s"][::gensPerFrame][:nFramesOut] + zProfile = profile["z"][::gensPerFrame][:nFramesOut] + zDeltaProfile = profile["zDelta"][::gensPerFrame][:nFramesOut] + + dem = np.flipud(demDict["rasterData"]) + cellSize = demDict["header"]["cellsize"] + hillshade = addHillshadeSimple(dem, cellSize, cfgHS) + + # fixed color range, computed once from the GLOBAL max across all frames + vmin = 0 + vmax = np.nanmax(frameStackHistory) + sharedCmap = cmapCrameri.batlow.reversed() + + norm = mcolors.Normalize( + vmin=0, + vmax=vmax, + clip=True + ) + + thalwegColor = cfgGen.get("thalwegColor") + thalwegWidth = cfgGen.getfloat("thalwegWidth") + pointColor = cfgGen.get("pointColor") + pointSize = cfgGen.getfloat("pointSize") + outlineColor = cfgGen.get("outlineColor") + outlineWidth = cfgGen.getfloat("outlineWidth") + historyAlpha = cfgGen.getfloat("historyAlpha") + currentAlpha = cfgGen.getfloat("currentAlpha") + historyPointsColor = cfgGen.get("historyPointsColor") + historyPointsAlpha = cfgGen.getfloat("historyPointsAlpha") + historyPointsSize = cfgGen.getfloat("historyPointsSize") + + if showProfilePanel: + fig = plt.figure(figsize=(14, 10)) + + gs = GridSpec( + 1, 2, + width_ratios=[1.0, 1.8], + wspace=0.5 + ) + + ax1 = fig.add_subplot(gs[0, 0]) + ax2 = fig.add_subplot(gs[0, 1]) + + else: + fig, ax1 = plt.subplots(figsize=(8, 10)) + + # mapping panel + ax1.imshow(hillshade, cmap="gray", vmin=0, vmax=1) + imHistory = ax1.imshow(np.where(frameStackHistory[0] > 0, frameStackHistory[0], np.nan), + cmap=sharedCmap, norm=norm, alpha=historyAlpha) + imCurrent = ax1.imshow(np.where(frameStackCurrent[0] > 0, frameStackCurrent[0], np.nan), + cmap=sharedCmap, norm=norm, alpha=currentAlpha) + + imHistory.set_clim(vmin, vmax) + imCurrent.set_clim(vmin, vmax) + + thalwegLine, = ax1.plot([], [], "-", color=thalwegColor, linewidth=thalwegWidth, + alpha=0.8, label="Thalweg (path so far)") + historyPoints = ax1.scatter([], [], s=historyPointsSize, color=historyPointsColor, + alpha=historyPointsAlpha, zorder=15, + label="Center of flux (history)") + currentPoint, = ax1.plot([], [], "o", color=pointColor, markersize=pointSize, + markeredgecolor="white", markeredgewidth=1.5, + label="Center of flux (current)") + + titleText = ax1.set_title("") + ax1.legend(loc="lower left") + ax1.set_axis_off() + + fig.canvas.draw() # triggers layout computation + + # get proper position for colorbar + p1 = ax1.get_position() + if showProfilePanel: + ax1.set_aspect("auto") + p2 = ax2.get_position() + + ax1.set_position([ + p1.x0, + p2.y0, + p1.width, + p2.height + ]) + p1 = ax1.get_position() + gap = p2.x0 - p1.x1 + else: + gap = 0.1 + + cax_width = 0.015 + cax = fig.add_axes([ + p1.x1 + 0.08 * gap, + p1.y0, + cax_width, + p1.height * 0.95, + ]) + + gradient = np.linspace(vmin, vmax, 256).reshape(-1, 1) + cax.imshow(gradient, aspect="auto", cmap=sharedCmap, origin="lower", + extent=[0, 1, vmin, vmax]) + cax.set_xticks([]) + cax.yaxis.tick_right() + cax.yaxis.set_label_position("right") + cax.set_ylabel(clabelDict[cfgGen.get("videoDataVariable")]) + + if showProfilePanel: + # --- right panel: growing z / zDelta profile ------------------------ + zLine, zVelLine, currentProfilePoint = setup2DProfileAxis(ax2, sProfile, zProfile, zDeltaProfile, sharedCmap, + norm) + + contourHolder = {"artist": None} + + + def update(gen): + history = frameStackHistory[gen] + current = frameStackCurrent[gen] + + imHistory.set_data(np.where(history > 0, history, np.nan)) + imCurrent.set_data(np.where(current > 0, current, np.nan)) + imHistory.set_clim(vmin, vmax) + imCurrent.set_clim(vmin, vmax) + + removeContour(contourHolder) + mask = (current > 0).astype(np.float32) + if mask.max() > 0: + contourHolder["artist"] = ax1.contour(mask, levels=[0.5], + colors=outlineColor, + linewidths=outlineWidth) + + validMask = ~np.isnan(coRow[: gen + 1]) + thalwegLine.set_data(coCol[: gen + 1][validMask], coRow[: gen + 1][validMask]) + + historyMask = ~np.isnan(coRow[:gen]) + historyOffsets = np.column_stack([coCol[:gen][historyMask], coRow[:gen][historyMask]]) + historyPoints.set_offsets(historyOffsets) + + currentPoint.set_data([coCol[gen]], [coRow[gen]]) + titleText.set_text(f"Generation {gen * gensPerFrame}") + + artists = (imHistory, imCurrent, thalwegLine, historyPoints, currentPoint, titleText) + + if showProfilePanel: + zLine.set_data(sProfile[: gen + 1], zProfile[: gen + 1]) + + x = sProfile[:gen + 1] + y = zProfile[:gen + 1] + zDeltaProfile[:gen + 1] + + points = np.column_stack([x, y]).reshape(-1, 1, 2) + + segments = np.concatenate( + [points[:-1], points[1:]], + axis=1, + ) + + zVelLine.set_segments(segments) + + if videoDataVariable == "velocity": + values = np.sqrt( + 2 * 9.81 * zDeltaProfile[:gen] + ) + else: + values = zDeltaProfile[:gen] + + zVelLine.set_array(values) + + currentProfilePoint.set_data([sProfile[gen]], [zProfile[gen] + zDeltaProfile[gen]]) + artists += (zLine, zVelLine, currentProfilePoint) + + return artists + + frameArrays = [] + for gen in range(len(frameStackHistory)): + update(gen) + frameArrays.append(figureToRgbArray(fig)) + + # --- final extra frame: last generation state + thalweg extension --- + if cfgGen.getboolean("showExtendedThalweg"): + if cfgPath.get("extendedProfilePicklePath", "") == "": + extendedProfilePicklePath = thalwegDataDir / f"extended_thalwegData_{centerOfPickle}_{relId}.pickle" + else: + extendedProfilePicklePath = cfgPath.get("extendedProfilePicklePath") + avaProfileMass = loadExtendedProfile(extendedProfilePicklePath) + topLine, bottomLine, centerLine = addThalwegExtension(ax1, avaProfileMass, demDict["header"]) + ax1.legend( + handles=[topLine, bottomLine, centerLine], + loc="lower left") + fig.canvas.draw() + frameArrays.append(figureToRgbArray(fig)) # just ONE extra frame now + # ----------------------------------------------------------------- + + plt.close(fig) + + basePalette = buildGlobalPalette(frameArrays) + pilFrames = [] + for arr in frameArrays: + img = Image.fromarray(arr, mode="RGB") + imgQuantized = img.quantize(palette=basePalette, dither=Image.NONE) + pilFrames.append(imgQuantized) + + normalDuration = int(1000 / fps) + finalFrameDuration = cfgGen.getint("finalFrameDurationMs", fallback=3000) # ms, e.g. 3 seconds + + if cfgGen.getboolean("showExtendedThalweg"): + # last frame (the extension frame) gets a longer duration + durations = [normalDuration] * (len(pilFrames) - 1) + [finalFrameDuration] + else: + durations = [normalDuration] * len(pilFrames) + + pilFrames[0].save( + outFile, + save_all=True, + append_images=pilFrames[1:], + duration=durations, + loop=0, # play once, do not loop; use loop=0 for infinite looping + ) + + print(f"Video saved: {outFile}") diff --git a/avaframe/out3Plot/outCom4GifCfg.ini b/avaframe/out3Plot/outCom4GifCfg.ini new file mode 100644 index 000000000..46f66196f --- /dev/null +++ b/avaframe/out3Plot/outCom4GifCfg.ini @@ -0,0 +1,59 @@ +# configuration for creating a GIF out from com4FlowPy's results videoData + +[GENERAL] +# if show2DProfilePanel is True, the profile of the thalweg is plotted additionally in a separate panel +# data is loaded from the respective thalweg data in the com4FlowPy results +show2DProfilePanel = True +# release id for that the video is created (only one is possible) +relId = +# simulation hash id of the simulation for that +simhash = +# provide the variable that should be plotted (it needs to exist in the videoData folder) +# possible options: +# z_delta +# flux +# min_distance (max travel length) +# max_gamma (max travelAngle)variable = z_delta +# velocity (computed from zDelta) +videoDataVariable = velocity +# choose if the thalweg defined as center of energy (coE) or center of flux (coF) is plotted +centerOf = coF +# frames per second of the output video (playback speed) +fps = 5 +# number of generation steps merged into a single video frame +gensPerFrame = 3 + +# color, linewidth, point size and transparency for the thalweg data +thalwegColor = red +thalwegWidth = 2 + +pointColor = red +pointSize = 8 + +outlineColor = white +outlineWidth = 1. + +historyAlpha = 0.5 +currentAlpha = 0.95 + +historyPointsColor = k +historyPointsAlpha = 0.3 +historyPointsSize = 7 + +# if True, the last frame shows the extended (and resampled) path (that is read from the extendedProfilePicklePath file +# or from the respective thalwegData in the result folder if extendedProfilePicklePath is empty) +showExtendedPath = True +# the duraion of the last frame if the extended path is shown +finalFrameDurationMs = 5000 + +[HILLSHADE] +azimuth = 315 +altitude = 45 +vertExag = 1.0 + +[PATH] +# if paths are empty the files are expected in the Avaframe structure +demPath = +videoDataPath = +outVideoPath = +extendedProfilePicklePath = \ No newline at end of file diff --git a/avaframe/runCom4FlowPy.py b/avaframe/runCom4FlowPy.py index bfb000b5a..a6882d75c 100644 --- a/avaframe/runCom4FlowPy.py +++ b/avaframe/runCom4FlowPy.py @@ -154,9 +154,14 @@ def main(avalancheDir="", cfg=None): if successToJSON is True: 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) + if cfgSetup["calcThalweg"] == "True": + cfgPath["thalwegDir"] = cfgPath["resDir"] / "thalwegData" + fU.makeADir(cfgPath["thalwegDir"]) + else: + cfgPath["thalwegDir"] = "" cfgPath["deleteTemp"] = "False" cfgPath["uid"] = uid @@ -250,6 +255,12 @@ def main(avalancheDir="", cfg=None): shutil.rmtree(temp_dir) fU.makeADir(temp_dir) + if cfgSetup["calcThalweg"] is True: + thalwegDir = workDir / res_dir / "thalwegData" + fU.makeADir(thalwegDir) + else: + thalwegDir = "" + # writing config to .json file successToJSON = writeCfgJSON(cfg, uid, workDir) @@ -260,6 +271,10 @@ def main(avalancheDir="", cfg=None): log.error("Exception occurred: %s", str(successToJSON), exc_info=True) cfgPath["workDir"] = pathlib.Path(workDir) + if cfgSetup["calcThalweg"] is True: + cfgPath["thalwegDir"] = pathlib.Path(thalwegDir) + else: + cfgPath["thalwegDir"] = None cfgPath["outDir"] = pathlib.Path(res_dir) cfgPath["resDir"] = cfgPath["outDir"] cfgPath["tempDir"] = pathlib.Path(temp_dir) @@ -408,9 +423,14 @@ def readFlowPyinputs(avalancheDir, cfgFlowPy, log): cfgPath["forestPath"] = forestPath # read release ID raster - if "relIdPolygon" in cfgFlowPy["PATHS"]["outputFiles"].split("|") or "relIdCount" in cfgFlowPy["PATHS"][ - "outputFiles" - ].split("|"): + if ( + "relIdPolygon" in cfgFlowPy["PATHS"]["outputFiles"].split("|") + or "relIdCount" in cfgFlowPy["PATHS"]["outputFiles"].split("|") + or ( + cfgFlowPy.getboolean("GENERAL", "thalwegReleaseArea") + and cfgFlowPy.getboolean("GENERAL", "calcThalweg") + ) + ): 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" diff --git a/avaframe/runScripts/runCreateGIFCom4.py b/avaframe/runScripts/runCreateGIFCom4.py new file mode 100644 index 000000000..c737b6298 --- /dev/null +++ b/avaframe/runScripts/runCreateGIFCom4.py @@ -0,0 +1,32 @@ +""" +Create a GIF out of com4FlowPy video data results +""" +import time + +# Local imports + +from avaframe.in3Utils import cfgUtils +from avaframe.in3Utils import logUtils +import avaframe.out3Plot.outCom4Gif as outCom4Gif + +# Time the whole routine +startTime = time.time() + +# log file name; leave empty to use default runLog.log +logName = 'runCreateGIFCom4' + +# Load avalanche directory from general configuration file +cfgMain = cfgUtils.getGeneralConfig() +avalancheDir = cfgMain['MAIN']['avalancheDir'] + +# Start logging +log = logUtils.initiateLogger(avalancheDir, logName) +log.info('MAIN SCRIPT') +log.info('Current avalanche: %s', avalancheDir) + +# TODO: add an option to run FlowPy before (and modify the GIF config to create a video for this simualtion) + +# Load configuration for hybrid model +cfg = cfgUtils.getModuleConfig(outCom4Gif) + +outCom4Gif.makeGenerationVideo(avalancheDir=avalancheDir, cfg=cfg) diff --git a/avaframe/tests/test_com4FlowPy.py b/avaframe/tests/test_com4FlowPy.py index a154db42c..cb086dbd6 100644 --- a/avaframe/tests/test_com4FlowPy.py +++ b/avaframe/tests/test_com4FlowPy.py @@ -136,9 +136,37 @@ def test_backTracking(): 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} + 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} + 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) @@ -157,7 +185,15 @@ def test_calculation(): ] ) 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]]) + 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 exp = 99 fluxTh = 0.001 @@ -185,6 +221,7 @@ def test_calculation(): "relVolBool": False, "relVolArray": None, } + rasterAttributes = {"cellsize": cellsize, "nodata": nodata} args = [ dem, infra, @@ -193,8 +230,7 @@ def test_calculation(): exp, fluxTh, zDeltaMax, - nodata, - cellsize, + rasterAttributes, infraBool, forestBool, variableParameters, @@ -204,6 +240,9 @@ def test_calculation(): forestParams, outputs, relOutputParams, + False, + False, + None, ] flux = np.ones_like(dem) * -9999.0 @@ -225,7 +264,6 @@ def test_calculation(): def createTestRaster(pathTestFolder, rasterName): - # first create test raster and save in test folder testRaster = np.zeros((10, 10)) @@ -299,7 +337,10 @@ def test_tileRaster(tmp_path): 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)) + assert ext21 == ( + (2 * xDim - 4 * U, 3 * xDim - 4 * U), + (yDim - 2 * U, 2 * yDim - 2 * U), + ) def test_mergeDict(tmp_path): @@ -563,6 +604,8 @@ def test_runCom4FlowPy(tmp_path): "cpuCount": "1", "tileSize": "15000", "tileOverlap": "5000", + "thalwegreleasearea": "False", + "calcThalweg": "False", } cfg["PATHS"] = { "outputFiles": "zDelta", @@ -752,6 +795,231 @@ def test_runCom4FlowPy(tmp_path): for key in resDictTest6: assert resDictTest6[key] == resDict[key] +def test_getMaskedRasters(): + raster = np.array( + [ + [1, 2, 3], + [4, 5, 0], + [0, 7, 8], + ] + ) + mask = np.array( + [ + [1, 1, 0], + [1, 0, 0], + [0, 0, 1], + ] + ) + + idsIn, idsOut = SPAM.getMaskedRasters(mask, raster) + + # inside mask: (0,0)=1, (0,1)=2, (1,0)=4, (2,2)=8 + assert list(idsIn) == [1, 2, 4, 8] + # outside mask: (0,2)=3, (1,1)=5, (1,2)=0, (2,0)=0, (2,1)=7 -> zeros/negatives dropped + assert list(idsOut) == [3, 5, 7] + + raster = np.array([[1, 2], [3, 4]]) + mask = np.ones_like(raster) + + idsIn, idsOut = SPAM.getMaskedRasters(mask, raster) + + assert list(idsIn) == [1, 2, 3, 4] + assert list(idsOut) == [] + + +def test_getTileEnds(tmp_path): + # same raster/geometry as used implicitly by test_tileRaster, so the + # expected exList/eyList can be cross-checked against the ext_i_j values + # asserted there. + pathTestFolder = tmp_path / "data" / "testCom4" + rasterName = "testRaster" + pathTempFolder = pathTestFolder / "tmp" + os.makedirs(pathTempFolder, exist_ok=True) + + createTestRaster(pathTestFolder, rasterName) + testData = IOf.readRaster(pathTestFolder / f"{rasterName}.tif", noDataToNan=False) + relIdRaster = testData["rasterData"] + + xDim, yDim, U = 4, 4, 1 + + exList, eyList = SPAM.getTileEnds(pathTempFolder, xDim, yDim, U, relIdRaster) + + assert exList == [6, 10] + assert eyList == [6, 10] + + nTiles = pickle.load(open(pathTempFolder / "nTiles", "rb")) + assert nTiles == (1, 1) + + # relId over entire raster + + relIdRaster = np.ones((10, 10)) + exList, eyList = SPAM.getTileEnds(pathTempFolder, xDim, yDim, U, relIdRaster) + + assert exList == [11] + assert eyList == [11] + nTiles = pickle.load(open(pathTempFolder / "nTiles", "rb")) + assert nTiles == (0, 0) + + extentLarge = pickle.load(open(pathTempFolder / "extentLarge", "rb")) + assert extentLarge == (10, 10) + + pathTempFolder = tmp_path / "tmp" + os.makedirs(pathTempFolder, exist_ok=True) + + relIdRaster = np.zeros((10, 10)) + # release area (ID 5) straddles the naive x-tile boundary (column 3/4) + relIdRaster[2:5, 2:6] = 5 + + xDim, yDim, U = 4, 4, 1 + + exList, eyList = SPAM.getTileEnds(pathTempFolder, xDim, yDim, U, relIdRaster) + + assert exList == [7, 9, 11] + assert eyList == [6, 8, 10] + + relIdRaster = np.zeros((10, 10)) + # release area fully inside what would become the first naive tile + relIdRaster[0:2, 0:2] = 7 + + xDim, yDim, U = 4, 4, 1 + + exList, eyList = SPAM.getTileEnds(pathTempFolder, xDim, yDim, U, relIdRaster) + + assert exList == [4, 6, 8, 10] + assert eyList == [4, 6, 8, 10] + + +def test_tileRasterWithIndices(tmp_path): + pathTestFolder = tmp_path / "data" / "testCom4" + rasterName = "testRaster" + ext = ".tif" + pathTempFolder = pathTestFolder / "tmp" + os.makedirs(pathTempFolder, exist_ok=True) + + createTestRaster(pathTestFolder, rasterName) + + fNameIn = pathTestFolder / f"{rasterName}{ext}" + fNameOut = "testTile" + U = 1 + # same end-indices SPAM.getTileEnds would produce for this raster/config + exList = [4, 6, 8, 10] + eyList = [4, 6, 8, 10] + + SPAM.tileRasterWithIndices(fNameIn, fNameOut, pathTempFolder, exList, eyList, U, isInit=False) + + testData = IOf.readRaster(fNameIn, noDataToNan=False) + testRaster = testData["rasterData"] + + nTiles = pickle.load(open(pathTempFolder / "nTiles", "rb")) + assert nTiles == (3, 3) + + # corner tile (0,0): rows 0:4, cols 0:4 + ext00 = pickle.load(open(pathTempFolder / "ext_0_0", "rb")) + assert ext00 == ((0, 4), (0, 4)) + tile00 = np.load(pathTempFolder / "testTile_0_0.npy") + assert tile00.shape == (4, 4) + assert np.all(tile00 == testRaster[0:4, 0:4]) + + # interior tile (1,2): rows 2:6, cols 4:8 + ext12 = pickle.load(open(pathTempFolder / "ext_1_2", "rb")) + assert ext12 == ((2, 6), (4, 8)) + tile12 = np.load(pathTempFolder / "testTile_1_2.npy") + assert np.all(tile12 == testRaster[2:6, 4:8]) + + # bottom-right corner tile (3,3): rows 6:10, cols 6:10 + ext33 = pickle.load(open(pathTempFolder / "ext_3_3", "rb")) + assert ext33 == ((6, 10), (6, 10)) + tile33 = np.load(pathTempFolder / "testTile_3_3.npy") + assert np.all(tile33 == testRaster[6:10, 6:10]) + + # with tiling init + pathTestFolder = tmp_path / "data" / "testCom4" + rasterName = "testRaster" + ext = ".tif" + pathTempFolder = pathTestFolder / "tmp" + os.makedirs(pathTempFolder, exist_ok=True) + + createTestRaster(pathTestFolder, rasterName) + + fNameIn = pathTestFolder / f"{rasterName}{ext}" + fNameOut = "testTileInit" + U = 1 + exList = [4, 6, 8, 10] + eyList = [4, 6, 8, 10] + + SPAM.tileRasterWithIndices(fNameIn, fNameOut, pathTempFolder, exList, eyList, U, isInit=True) + + testData = IOf.readRaster(fNameIn, noDataToNan=False) + testRaster = testData["rasterData"] + # test raster only contains values >= 0, so any -9999 found below must + # come from the edge-nulling logic, not from the source data + assert not np.any(testRaster == -9999) + + # --- corner tile (0,0): east edge nulled (j != JMAX) and south edge + # nulled (i != IMAX); north/west untouched (i == 0, j == 0) + tile00 = np.load(pathTempFolder / "testTileInit_0_0.npy") + assert np.all(tile00[:, -U:] == -9999) # east + assert np.all(tile00[-U:, :] == -9999) # south + assert np.all(tile00[0, :-U] == testRaster[0, 0:4][:-U]) # north untouched + assert np.all(tile00[:-U, 0] == testRaster[0:4, 0][:-U]) # west untouched + + # --- interior tile (1,2): all four edges nulled + tile12 = np.load(pathTempFolder / "testTileInit_1_2.npy") + assert np.all(tile12[:, -U:] == -9999) # east + assert np.all(tile12[0:U, :] == -9999) # north + assert np.all(tile12[:, 0:U] == -9999) # west + assert np.all(tile12[-U:, :] == -9999) # south + + # --- bottom-right corner tile (3,3): north edge nulled (i != 0) and + # west edge nulled (j != 0); east/south untouched (j == JMAX, i == IMAX) + tile33 = np.load(pathTempFolder / "testTileInit_3_3.npy") + assert np.all(tile33[0:U, :] == -9999) # north + assert np.all(tile33[:, 0:U] == -9999) # west + assert np.all(tile33[-1, U:] == testRaster[9, 6:10][U:]) # south untouched + assert np.all(tile33[U:, -1] == testRaster[6:10, 9][U:]) # east untouched + + # workflow: first get tiles, then make tiling + createTestRaster(pathTestFolder, rasterName) + testData = IOf.readRaster(pathTestFolder / f"{rasterName}.tif", noDataToNan=False) + relIdRaster = testData["rasterData"] + + xDim, yDim, U = 4, 4, 1 + fNameIn = pathTestFolder / f"{rasterName}.tif" + fNameOut = "testTile" + + exList, eyList = SPAM.getTileEnds(pathTempFolder, xDim, yDim, U, relIdRaster) + SPAM.tileRasterWithIndices(fNameIn, fNameOut, pathTempFolder, exList, eyList, U, isInit=True) + + tile00 = np.load(pathTempFolder / "testTile_0_0.npy") + tile01 = np.load(pathTempFolder / "testTile_0_1.npy") + tile10 = np.load(pathTempFolder / "testTile_1_0.npy") + tile11 = np.load(pathTempFolder / "testTile_1_1.npy") + + ids00 = np.unique(tile00) + ids01 = np.unique(tile01) + ids10 = np.unique(tile10) + ids11 = np.unique(tile11) + + for id in ids00: + if id <= 0: + continue + assert id not in ids01 + assert id not in ids10 + assert id not in ids11 + + for id in ids10: + if id <= 0: + continue + assert id not in ids01 + assert id not in ids00 + assert id not in ids11 + + for id in ids01: + if id <= 0: + continue + assert id not in ids00 + assert id not in ids10 + assert id not in ids11 def testCompareRasters(monkeypatch): """Test the comparison of two raster arrays. @@ -1227,6 +1495,241 @@ def test_checkVariableInputParameters_varExponent(monkeypatch): com4FlowPy.checkVariableInputParameters(modelParameters, makeModelPaths(), validParamRanges) +def test_get_start_idx_sortedByAltitudeDescending(): + """ test that release pixels are returned sorted by altitude, highest first """ + dem = np.array([ + [100.0, 200.0, 150.0], + [300.0, 50.0, 400.0], + [10.0, 500.0, 20.0], + ]) + release = np.array([ + [1, 1, 0], + [1, 0, 1], + [0, 1, 0], + ]) + + row_list, col_list = flowCore.get_start_idx(dem, release) + + # release pixels are at (0,0)=100, (0,1)=200, (1,0)=300, (1,2)=400, (2,1)=500 + # sorted descending by altitude: 500, 400, 300, 200, 100 + expectedOrder = [(2, 1), (1, 2), (1, 0), (0, 1), (0, 0)] + actualOrder = list(zip(row_list, col_list)) + + assert actualOrder == expectedOrder + + +def test_get_start_idx_noReleasePixels(): + """ test that empty row/col lists are returned when there are no release pixels """ + dem = np.array([[100.0, 200.0], [300.0, 400.0]]) + release = np.zeros((2, 2), dtype=int) + + row_list, col_list = flowCore.get_start_idx(dem, release) + + assert len(row_list) == 0 + assert len(col_list) == 0 + + +def test_get_start_idx_singleReleasePixel(): + """ test behavior with exactly one release pixel """ + dem = np.array([[10.0, 20.0], [30.0, 40.0]]) + release = np.array([[0, 0], [0, 1]]) + + row_list, col_list = flowCore.get_start_idx(dem, release) + + assert list(row_list) == [1] + assert list(col_list) == [1] + + +def test_get_start_idx_sortedByRelIdWhenThalwegRequested(): + """ test that with relIdArray and calcThalweg=True, pixels are grouped/sorted by release Id + (descending), rather than purely by altitude + """ + dem = np.array([ + [100.0, 200.0, 150.0], + [300.0, 50.0, 400.0], + ]) + release = np.array([ + [1, 1, 1], + [1, 0, 1], + ]) + relIdArray = np.array([ + [1, 2, 1], + [2, 0, 3], + ]) + + row_list, col_list = flowCore.get_start_idx(dem, release, relIdArray=relIdArray, calcThalweg=True) + + # cells and their relId: (0,0)->1, (0,1)->2, (0,2)->1, (1,0)->2, (1,2)->3 + # sorted primarily by relId descending: relId 3 first, then relId 2s, then relId 1s + resultRelIds = [relIdArray[r, c] for r, c in zip(row_list, col_list)] + assert resultRelIds == sorted(resultRelIds, reverse=True) + + # relId 3 has exactly one cell -> must be (1,2) + assert (row_list[0], col_list[0]) == (1, 2) + + # cells belonging to the same relId are contiguous in the output + idxRelId1 = [i for i, rid in enumerate(resultRelIds) if rid == 1] + assert idxRelId1 == list(range(min(idxRelId1), max(idxRelId1) + 1)) + idxRelId2 = [i for i, rid in enumerate(resultRelIds) if rid == 2] + assert idxRelId2 == list(range(min(idxRelId2), max(idxRelId2) + 1)) + + +def test_get_start_idx_relIdArrayProvidedButCalcThalwegFalse(): + """ test that relIdArray is ignored (falls back to altitude sort) when calcThalweg is False """ + dem = np.array([ + [100.0, 200.0], + [300.0, 400.0], + ]) + release = np.array([ + [1, 1], + [1, 1], + ]) + relIdArray = np.array([ + [5, 5], + [1, 1], + ]) + + row_list, col_list = flowCore.get_start_idx(dem, release, relIdArray=relIdArray, calcThalweg=False) + + # falls back to pure altitude-descending sort, ignoring relIdArray + expectedOrder = [(1, 1), (1, 0), (0, 1), (0, 0)] + actualOrder = list(zip(row_list, col_list)) + assert actualOrder == expectedOrder + + +def test_split_release_evenSplit_byPixelCount(): + """ test the default (non-thalweg) branch: release cells split roughly evenly + by cumulative pixel count into `pieces` chunks + """ + release = np.zeros((2, 10), dtype=int) + release[0, :] = 1 # 10 release pixels along the first row + + release_list = flowCore.split_release(release, pieces=2, relIdArray=None, calcThalweg=False) + + assert len(release_list) == 2 + for piece in release_list: + assert piece.shape == release.shape + + # every original release pixel appears in exactly one piece + totalReconstructed = sum(piece.sum() for piece in release_list) + assert totalReconstructed == release.sum() + + # no overlap between pieces + combined = np.zeros_like(release) + for piece in release_list: + assert np.all((combined & piece) == 0) # no pixel assigned twice + combined = combined | piece + np.testing.assert_array_equal(combined, release) + + +def test_split_release_evenSplit_singlePiece(): + """ test that pieces=1 returns the whole release layer unchanged (in a single piece) """ + release = np.array([ + [1, 0, 1], + [0, 1, 0], + ]) + + release_list = flowCore.split_release(release, pieces=1, relIdArray=None, calcThalweg=False) + + assert len(release_list) == 1 + np.testing.assert_array_equal(release_list[0], release) + + +def test_split_release_evenSplit_unevenPixelCount(): + """ test even splitting when the number of release pixels doesn't divide evenly into pieces """ + release = np.zeros((1, 7), dtype=int) + release[0, :] = 1 # 7 release pixels, split into 3 pieces + + release_list = flowCore.split_release(release, pieces=3, relIdArray=None, calcThalweg=False) + + assert len(release_list) == 3 + totalReconstructed = sum(piece.sum() for piece in release_list) + assert totalReconstructed == 7 + + # no overlaps, full reconstruction + combined = np.zeros_like(release) + for piece in release_list: + combined = combined | piece + np.testing.assert_array_equal(combined, release) + + +def test_split_release_thalweg_groupsByReleaseId(): + """ test the thalweg branch: cells belonging to the same relId stay together in one chunk, + and chunks are balanced by total cell count + """ + release = np.ones((2, 6), dtype=int) + # 3 release areas of sizes 6, 4, 2 + relIdArray = np.array([ + [1, 1, 1, 2, 2, 2], + [1, 1, 1, 2, 3, 3], + ]) + + release_list = flowCore.split_release(release, pieces=2, relIdArray=relIdArray, calcThalweg=True) + + assert len(release_list) == 2 + + # every cell belonging to a given relId must end up entirely within a single chunk + for relId in np.unique(relIdArray): + # find the mask of all cells belonging to this release Id + cellsWithThisId = (relIdArray == relId) + + # count in how many of the output chunks these cells actually show up + numChunksContainingId = 0 + for piece in release_list: + pixelsInThisChunk = piece[cellsWithThisId] + if np.any(pixelsInThisChunk > 0): + numChunksContainingId += 1 + + # a release area must not be split across multiple chunks + assert numChunksContainingId == 1 + + # every release pixel appears in exactly one piece, none lost/duplicated + combined = np.zeros_like(release) + for piece in release_list: + assert np.all((combined & piece) == 0) + combined = combined | piece + np.testing.assert_array_equal(combined, release) + + +def test_split_release_thalweg_piecesClampedToUniqueIdCount(): + """ test that pieces is clamped down to the number of unique release Ids when + there are fewer distinct release areas than requested pieces + """ + release = np.array([ + [1, 1], + [0, 0], + ]) + relIdArray = np.array([ + [1, 1], + [0, 0], + ]) + # only 1 unique release id, but pieces=5 requested + release_list = flowCore.split_release(release, pieces=5, relIdArray=relIdArray, calcThalweg=True) + + # should be clamped down to 1 chunk (np.minimum(pieces, len(uniqueIds))) + assert len(release_list) == 1 + np.testing.assert_array_equal(release_list[0], release) + + +def test_split_release_thalweg_balancesChunkSizes(): + """ test that with several release areas of different sizes, the two most balanced + combinations of areas end up in different chunks (greedy balancing by count) + """ + release = np.ones((1, 15), dtype=int) + # release areas of very different sizes: 10, 3, 2 + relIdArray = np.array([[1] * 10 + [2] * 3 + [3] * 2]) + + release_list = flowCore.split_release(release, pieces=2, relIdArray=relIdArray, calcThalweg=True) + + counts = [piece.sum() for piece in release_list] + assert sum(counts) == 15 + # the big area (10 cells) should end up alone in one chunk since combining it with + # anything else would make that chunk more unbalanced (greedy assigns smallest chunk first) + assert 10 in counts + # other chunk should contain the two smaller areas (3+2=5) + assert 5 in counts + + if __name__ == "__main__": test_add_os() test_reverseTopology() @@ -1243,3 +1746,17 @@ def test_checkVariableInputParameters_varExponent(monkeypatch): test_checkGlobalParameters_exp() test_checkGlobalParameters_fluxThreshold() test_checkGlobalParameters_errorMessageMentionsValidRange() + test_getMaskedRasters() + test_getTileEnds(tmpDir) + test_tileRasterWithIndices(tmpDir) + test_get_start_idx_sortedByAltitudeDescending() + test_get_start_idx_noReleasePixels() + test_get_start_idx_singleReleasePixel() + test_get_start_idx_sortedByRelIdWhenThalwegRequested() + test_get_start_idx_relIdArrayProvidedButCalcThalwegFalse() + test_split_release_evenSplit_byPixelCount() + test_split_release_evenSplit_singlePiece() + test_split_release_evenSplit_unevenPixelCount() + test_split_release_thalweg_groupsByReleaseId() + test_split_release_thalweg_piecesClampedToUniqueIdCount() + test_split_release_thalweg_balancesChunkSizes() diff --git a/avaframe/tests/test_flowPath.py b/avaframe/tests/test_flowPath.py new file mode 100644 index 000000000..596fa5b02 --- /dev/null +++ b/avaframe/tests/test_flowPath.py @@ -0,0 +1,556 @@ +import numpy as np +import pytest +import pickle +import logging +from types import SimpleNamespace + +# adjust this import to match your actual module path +from avaframe.com4FlowPy.flowPath import Path # noqa: F401 (placeholder alias, see note below) + + +# --------------------------------------------------------------------------- +# helpers +# --------------------------------------------------------------------------- + +def makeCell(z_delta=0.0, flux=0.0, min_distance=0.0, altitude=0.0, + rowindex=0, colindex=0, max_gamma=0.0, flowEnergy=0.0, + fluxDep=0.0, alpha=25.0, exp=8, max_z_delta=100.0): + """ build a minimal fake cell object with the attributes Path expects """ + return SimpleNamespace( + z_delta=z_delta, flux=flux, min_distance=min_distance, altitude=altitude, + rowindex=rowindex, colindex=colindex, max_gamma=max_gamma, + flowEnergy=flowEnergy, fluxDep=fluxDep, alpha=alpha, exp=exp, + max_z_delta=max_z_delta, + ) + + +def makeRasterAttributes(cellsize=10.0, nrows=5, ncols=5, extentTile=((0, 5), (0, 5))): + return {"cellsize": cellsize, "nrows": nrows, "ncols": ncols, "extentTile": extentTile} + + +def test_Path_init_withGenList(): + """ test Path init when genList is provided (normal, non-RAM-saving branch) """ + dem = np.zeros((5, 5)) + countArray = np.zeros((5, 5)) + countArray[1, 1] = 3 + countArray[2, 2] = 5 + + cell1 = makeCell(alpha=27.0, exp=8, max_z_delta=150.0) + cell2 = makeCell(alpha=27.0, exp=8, max_z_delta=150.0) + genList = [[cell1, cell2], [cell1]] + + rasterAttributes = makeRasterAttributes() + + p = Path(dem, 2, 3, genList, rasterAttributes, countArray, relId=7) + + assert p.cellsize == 10.0 + assert p.nrows == 5 + assert p.startcellRow == 2 + assert p.startcellCol == 3 + assert p.relId == 7 + assert p.alpha == 27.0 + assert p.exp == 8 + assert p.maxZDelta == 150.0 + assert p.numberGen == 2 + + # pathRaster: countArray values kept where >0, else nan + assert p.pathRaster[1, 1] == 3 + assert p.pathRaster[2, 2] == 5 + assert np.isnan(p.pathRaster[0, 0]) + + # output arrays initialized to zero, correct shape/dtype + for arrName in ["zDeltaArray", "flowEnergyArray", "fluxArray", "routFluxSumArray", "depFluxSumArray"]: + arr = getattr(p, arrName) + assert arr.shape == dem.shape + assert arr.dtype == np.float32 + np.testing.assert_array_equal(arr, np.zeros_like(dem, dtype=np.float32)) + + +def test_Path_init_withoutGenList(): + """ test Path init when genList is None (RAM-saving branch, uses listsRelId/exampleCell) """ + dem = np.zeros((5, 5)) + countArray = np.zeros((5, 5)) + rasterAttributes = makeRasterAttributes() + + listsRelId = { + "row": [1, 2, 3], + "col": [1, 2, 3], + "flux": [0.1, 0.2, 0.3], + "zdelta": [10.0, 20.0, 30.0], + "travelLengthMax": [5.0, 10.0, 15.0], + } + exampleCell = makeCell(alpha=30.0, exp=9, max_z_delta=200.0) + + p = Path(dem, 0, 0, None, rasterAttributes, countArray, relId=1, + listsRelId=listsRelId, exampleCell=exampleCell) + + assert p.rowList == listsRelId["row"] + assert p.colList == listsRelId["col"] + assert p.fluxList == listsRelId["flux"] + assert p.zdeltaList == listsRelId["zdelta"] + assert p.travelLengthList == listsRelId["travelLengthMax"] + assert p.alpha == 30.0 + assert p.exp == 9 + assert p.maxZDelta == 200.0 + assert p.numberGen == 3 + + +# --------------------------------------------------------------------------- +# getListFromCellList +# --------------------------------------------------------------------------- + +def _minimalPath(): + """ build a minimal, cheaply-constructed Path instance for method-level tests """ + dem = np.zeros((4, 4)) + countArray = np.zeros((4, 4)) + cell = makeCell() + genList = [[cell]] + rasterAttributes = makeRasterAttributes(nrows=4, ncols=4, extentTile=((2, 6), (3, 7))) + return Path(dem, 0, 0, genList, rasterAttributes, countArray, relId=1) + + +def test_getListFromCellList_zDelta(): + """ test that variable 'zDelta' maps to cell.z_delta """ + dummyPath = _minimalPath() + cells = [makeCell(z_delta=1.0, flux=2.0, min_distance=3.0, altitude=4.0, + rowindex=5, colindex=6, max_gamma=7.0, flowEnergy=8.0), + makeCell(z_delta=10.0, flux=20.0, min_distance=30.0, altitude=40.0, + rowindex=50, colindex=60, max_gamma=70.0, flowEnergy=80.0)] + + result = dummyPath.getListFromCellList(cells, "zDelta") + assert result == [1.0, 10.0] + + result = dummyPath.getListFromCellList(cells, "flux") + assert result == [2.0, 20.0] + + result = dummyPath.getListFromCellList(cells, "travelLength") + assert result == [3.0, 30.0] + + result = dummyPath.getListFromCellList(cells, "row") + assert result == [5, 50] + + result = dummyPath.getListFromCellList(cells, "col") + assert result == [6, 60] + + +# --------------------------------------------------------------------------- +# getGenerationList +# --------------------------------------------------------------------------- + +def test_getGenerationList_allGenerations(): + """ test that getGenerationList (no generation arg) returns nested per-generation lists """ + dem = np.zeros((4, 4)) + countArray = np.zeros((4, 4)) + gen0 = [makeCell(flux=1.0), makeCell(flux=2.0)] + gen1 = [makeCell(flux=3.0)] + genList = [gen0, gen1] + rasterAttributes = makeRasterAttributes(nrows=4, ncols=4) + + p = Path(dem, 0, 0, genList, rasterAttributes, countArray, relId=1) + + result = p.getGenerationList("flux") + + assert result == [[1.0, 2.0], [3.0]] + + +def test_getGenerationList_singleGeneration(): + """ test that getGenerationList with a specific generation index returns a flat list """ + dem = np.zeros((4, 4)) + countArray = np.zeros((4, 4)) + gen0 = [makeCell(flux=1.0), makeCell(flux=2.0)] + gen1 = [makeCell(flux=3.0)] + genList = [gen0, gen1] + rasterAttributes = makeRasterAttributes(nrows=4, ncols=4) + + p = Path(dem, 0, 0, genList, rasterAttributes, countArray, relId=1) + + result = p.getGenerationList("flux", generation=1) + + assert result == [3.0] + + +# --------------------------------------------------------------------------- +# getPathArrays +# --------------------------------------------------------------------------- + +def test_getPathArrays(): + """ test that getPathArrays fills max/sum arrays correctly from the genList cells """ + dem = np.zeros((4, 4)) + countArray = np.zeros((4, 4)) + rasterAttributes = makeRasterAttributes(nrows=4, ncols=4) + + # same (row, col) hit twice, across two generations -> max/sum behavior should differ + cellA = makeCell(rowindex=1, colindex=1, z_delta=5.0, flowEnergy=50.0, flux=0.5, fluxDep=0.1) + cellB = makeCell(rowindex=1, colindex=1, z_delta=8.0, flowEnergy=20.0, flux=0.3, fluxDep=0.2) + cellC = makeCell(rowindex=2, colindex=3, z_delta=1.0, flowEnergy=1.0, flux=1.0, fluxDep=1.0) + genList = [[cellA], [cellB, cellC]] + + p = Path(dem, 0, 0, genList, rasterAttributes, countArray, relId=1) + p.getPathArrays() + + zDeltaArrayRef = np.zeros((4, 4)) + zDeltaArrayRef[1, 1] = 8.0 + zDeltaArrayRef[2, 3] = 1.0 + + # max of z_delta at (1,1) across both hits + assert p.zDeltaArray[1, 1] == 8.0 + # max of flowEnergy at (1,1) + assert p.flowEnergyArray[1, 1] == 50.0 + # max of flux at (1,1) + assert p.fluxArray[1, 1] == pytest.approx(0.5) + # sums accumulate across hits + assert p.routFluxSumArray[1, 1] == pytest.approx(0.5 + 0.3) + assert p.depFluxSumArray[1, 1] == pytest.approx(0.1 + 0.2) + + # cell hit only once + assert p.zDeltaArray[2, 3] == 1.0 + assert p.routFluxSumArray[2, 3] == pytest.approx(1.0) + + # untouched cell stays zero + assert p.zDeltaArray[0, 0] == 0.0 + assert p.routFluxSumArray[0, 0] == 0.0 + + assert np.all(p.zDeltaArray == zDeltaArrayRef) + + +# --------------------------------------------------------------------------- +# calcThalwegCenterof +# --------------------------------------------------------------------------- + +def test_calcThalwegCenterof_weightedAverage(): + """ test weighted average computation when weights are non-zero """ + p = _minimalPath() + p.numberGen = 2 + + variable = [[1.0, 3.0], [10.0]] + variableCo = [[1.0, 1.0], [5.0]] # generation 0: equal weights -> avg 2.0; generation 1: single value + + variableSum, coVar = p.calcThalwegCenterof(variable, variableCo) + + np.testing.assert_allclose(variableSum, [4.0, 10.0]) + np.testing.assert_allclose(coVar, [2.0, 10.0]) + + +def test_calcThalwegCenterof_zeroWeightFallback(): + """ test that a generation with all-zero weights falls back to a plain average """ + p = _minimalPath() + p.numberGen = 1 + + variable = [[2.0, 4.0, 6.0]] + variableCo = [[0.0, 0.0, 0.0]] # sum of weights is 0 -> fallback branch + + variableSum, coVar = p.calcThalwegCenterof(variable, variableCo) + + assert variableSum[0] == pytest.approx(12.0) + assert coVar[0] == pytest.approx(4.0) # plain average of [2,4,6] + + +# --------------------------------------------------------------------------- +# getCenterofs +# --------------------------------------------------------------------------- + +def test_getCenterofs_setsWeightedAttributes_col(): + """ test that getCenterofs computes and sets e.g. colCoF from the generation data """ + dem = np.zeros((4, 4)) + countArray = np.zeros((4, 4)) + rasterAttributes = makeRasterAttributes(nrows=4, ncols=4) + gen0 = [makeCell(rowindex=1, colindex=4, flux=1.0), makeCell(rowindex=3, colindex=12, flux=3.0)] + gen1 = [makeCell(rowindex=5, colindex=5, flux=5.0)] + genList = [gen0, gen1] + p = Path(dem, 0, 0, genList, rasterAttributes, countArray, relId=1) + variables = ["col"] + p.getCenterofs(variables, ["CoF"]) + # col weighted by flux -> weighted avg of [1,3] w/ weights [1,3] = (1*1+3*3)/4 = 2.5; gen1 single value 5 + assert hasattr(p, "colCoF") + np.testing.assert_allclose(p.colCoF, [10.0, 5.0]) + + variables = ["flux"] + p.getCenterofs(variables, ["CoF"]) + + # flux weighted by itself -> weighted avg of [1,3] w/ weights [1,3] = (1*1+3*3)/4 = 2.5; gen1 single value 5 + assert hasattr(p, "fluxCoF") + np.testing.assert_allclose(p.fluxCoF, [2.5, 5.0]) + + variables = ["row"] + p.getCenterofs(variables, ["CoF"]) + # row weighted by flux -> weighted avg of [1,3] w/ weights [1,3] = (1*1+3*3)/4 = 2.5; gen1 single value 5 + assert hasattr(p, "rowCoF") + np.testing.assert_allclose(p.rowCoF, [2.5, 5.0]) + + +def test_getCenterofs_skipsExcludedVariables(): + """ test that variables in the exclusion list are skipped (no attribute is set) """ + dem = np.zeros((4, 4)) + countArray = np.zeros((4, 4)) + rasterAttributes = makeRasterAttributes(nrows=4, ncols=4) + genList = [[makeCell(flux=1.0)]] + + p = Path(dem, 0, 0, genList, rasterAttributes, countArray, relId=1) + + variables = ["zDeltaArray", "x"] + p.getCenterofs(variables, ["CoF"]) + + assert not hasattr(p, "xCoF") + assert not hasattr(p, "zDeltaArrayCoF") + + +def test_getCenterofs_expandsDepFluxSumAndFluxSum(): + """ test that 'depFluxSum'/'fluxSum' trigger appending 'depFlux'/'flux' to the variables list """ + dem = np.zeros((4, 4)) + countArray = np.zeros((4, 4)) + rasterAttributes = makeRasterAttributes(nrows=4, ncols=4) + genList = [[makeCell(flux=2.0, fluxDep=1.0)]] + + p = Path(dem, 0, 0, genList, rasterAttributes, countArray, relId=1) + + variables = ["fluxSum"] + p.getCenterofs(variables, ["CoF"]) + + # "flux" should have been appended and processed + assert "flux" in variables + assert hasattr(p, "fluxCoF") + + +# --------------------------------------------------------------------------- +# correctIndicesTile +# --------------------------------------------------------------------------- + +def test_correctIndicesTile(): + """ test that row/col indices are correctly offset by the tile's extent """ + p = _minimalPath() + p.rasterAttributes["extentTile"] = ((100, 200), (50, 150)) + + row = np.array([0, 1, 2]) + col = np.array([0, 5, 10]) + + rowLarge, colLarge = p.correctIndicesTile(row, col) + + np.testing.assert_array_equal(rowLarge, row + 100) + np.testing.assert_array_equal(colLarge, col + 50) + + +# --------------------------------------------------------------------------- +# saveDict +# --------------------------------------------------------------------------- + +def test_saveDict_withRelId(tmp_path): + """ test that saveDict writes a pickle file named by relId and with correct content """ + p = _minimalPath() + p.alpha = 27.456 + p.exp = 8 + p.maxZDelta = 123.456 + p.numberGen = 4 + p.relId = 42 + + # attributes needed for variable "x" and "s" with centerOf "CoF" + p.xCoF = np.array([1.0, 2.0, 3.0]) + p.travelLengthCoF = np.array([0.0, 10.0, 20.0]) + + p.saveDict(tmp_path, ["CoF"], ["x", "s"]) + + outFile = tmp_path / "thalwegData_CoF_42.pickle" + assert outFile.is_file() + + with open(outFile, "rb") as f: + data = pickle.load(f) + + assert data["alpha"] == pytest.approx(27.5, abs=0.05) # rounded to 1 decimal + assert data["exponent"] == 8 + assert data["zDeltaMax"] == pytest.approx(123.5, abs=0.05) + assert data["numberGen"] == 4 + np.testing.assert_allclose(data["x"], [1.0, 2.0, 3.0]) + np.testing.assert_allclose(data["s"], [0.0, 10.0, 20.0]) + + +def test_saveDict_withoutRelId(tmp_path): + """ test that saveDict falls back to startcellRow/startcellCol naming when relId is None """ + p = _minimalPath() + p.relId = None + p.startcellRow = 7 + p.startcellCol = 9 + p.xCoE = np.array([1.0]) + + p.saveDict(tmp_path, ["CoE"], ["x"]) + + outFile = tmp_path / "thalwegData_CoE_7_9.pickle" + assert outFile.is_file() + + +# --------------------------------------------------------------------------- +# calcAndSaveThalwegData +# --------------------------------------------------------------------------- + +def test_calcAndSaveThalwegData_invalidCoRaises(tmp_path): + """ test that an invalid thalweg 'centerOf' parameter raises ValueError """ + p = _minimalPath() + + thalwegParameters = { + "thalwegDir": tmp_path, + "thalwegSaveRam": False, + "thalwegCenterOf": "['notValid']", + "thalwegVariables": "['x']", + } + + with pytest.raises(ValueError): + p.calcAndSaveThalwegData(thalwegParameters) + + +def test_calcAndSaveThalwegData_saveRamBranch(tmp_path): + """ test the thalwegSaveRam=True branch: computes CoF attributes from the RAM-saving lists + and writes output via saveDict, using the real gT.indicesToCoords for coordinate conversion. + """ + dem = np.zeros((10, 10)) + countArray = np.zeros((10, 10)) + header = {"cellsize": 10.0, "xllcenter": 0.0, "yllcenter": 0.0} + rasterAttributes = { + "cellsize": header["cellsize"], + "xllcenter": header["xllcenter"], + "yllcenter": header["yllcenter"], + "nrows": 10, + "ncols": 10, + "extentTile": ((0, 10), (0, 10)), + } + + listsRelId = { + "row": [[1, 2, 3]], + "col": [[1, 2, 3]], + "flux": [[1.0, 1.0, 1.0]], + "zdelta": [[10.0, 20.0, 30.0]], + "travelLengthMax": [[0.0, 5.0, 10.0]], + } + exampleCell = makeCell(alpha=25.0, exp=8, max_z_delta=100.0) + + p = Path(dem, 0, 0, None, rasterAttributes, countArray, relId=3, + listsRelId=listsRelId, exampleCell=exampleCell) + + thalwegParameters = { + "thalwegDir": tmp_path, + "thalwegSaveRam": True, + } + + p.calcAndSaveThalwegData(thalwegParameters) + + assert hasattr(p, "colCoF") + assert hasattr(p, "rowCoF") + assert hasattr(p, "zdeltaCoF") + assert hasattr(p, "travelLengthCoF") + assert hasattr(p, "xCoF") + assert hasattr(p, "yCoF") + + # since row/col/flux are all equal-weighted with equal flux ([1,2,3] weighted equally by [1,1,1]), + # the weighted center-of-flux is just the mean: row = col = 2 + assert p.colCoF == pytest.approx(2.0) + assert p.rowCoF == pytest.approx(2.0) + + # tile offset is (0,0) here, so rowLarge/colLarge == rowCoF/colCoF + # x = xllcorner + col*cellsize = -5.0 + 2*10.0 = 15.0 (xllcenter=0, cellsize=10 -> xllcorner=-5) + assert p.xCoF == pytest.approx(15.0) + assert p.yCoF == pytest.approx(15.0) + + # genList should have been emptied to save RAM + assert p.genList == [] + + outFile = tmp_path / "thalwegData_CoF_3.pickle" + assert outFile.is_file() + + +def test_getGenerationList_threeGenerations_allGenerations(): + """ test getGenerationList (no generation arg) with a 3-generation genList, + including a generation that reuses cell1 from generation 0 + """ + dem = np.zeros((5, 5)) + countArray = np.zeros((5, 5)) + + cell1 = makeCell(z_delta=1.0, flux=1.0) + cell2 = makeCell(z_delta=2.0, flux=2.0) + cell3 = makeCell(z_delta=3.0, flux=3.0) + cell4 = makeCell(z_delta=4.0, flux=4.0) + cell5 = makeCell(z_delta=5.0, flux=5.0) + + genList = [[cell1, cell2], [cell3, cell4, cell5], [cell1]] + rasterAttributes = makeRasterAttributes(nrows=5, ncols=5) + + p = Path(dem, 0, 0, genList, rasterAttributes, countArray, relId=1) + + result = p.getGenerationList("zDelta") + + assert result == [[1.0, 2.0], [3.0, 4.0, 5.0], [1.0]] + + +def test_getGenerationList_threeGenerations_singleGeneration(): + """ test getGenerationList with an explicit generation index into a 3-generation genList """ + dem = np.zeros((5, 5)) + countArray = np.zeros((5, 5)) + + cell1 = makeCell(flux=1.0) + cell2 = makeCell(flux=2.0) + cell3 = makeCell(flux=3.0) + cell4 = makeCell(flux=4.0) + cell5 = makeCell(flux=5.0) + + genList = [[cell1, cell2], [cell3, cell4, cell5], [cell1]] + rasterAttributes = makeRasterAttributes(nrows=5, ncols=5) + + p = Path(dem, 0, 0, genList, rasterAttributes, countArray, relId=1) + + # generation 0: 2 cells + assert p.getGenerationList("flux", generation=0) == [1.0, 2.0] + # generation 1: 3 cells + assert p.getGenerationList("flux", generation=1) == [3.0, 4.0, 5.0] + # generation 2: reuses cell1 -> 1 cell + assert p.getGenerationList("flux", generation=2) == [1.0] + + +def test_getPathArrays_threeGenerations_repeatedCell(): + """ test getPathArrays with a 3-generation genList where cell1 (row=0, col=0) appears + in both generation 0 and generation 2, verifying max/sum behavior across the repeat + """ + dem = np.zeros((5, 5)) + countArray = np.zeros((5, 5)) + + cell1 = makeCell(rowindex=0, colindex=0, z_delta=5.0, flowEnergy=50.0, flux=1.0, fluxDep=0.5) + cell2 = makeCell(rowindex=1, colindex=1, z_delta=2.0, flowEnergy=20.0, flux=2.0, fluxDep=0.2) + cell3 = makeCell(rowindex=2, colindex=2, z_delta=3.0, flowEnergy=30.0, flux=3.0, fluxDep=0.3) + cell4 = makeCell(rowindex=3, colindex=3, z_delta=4.0, flowEnergy=40.0, flux=4.0, fluxDep=0.4) + cell5 = makeCell(rowindex=4, colindex=4, z_delta=1.0, flowEnergy=10.0, flux=0.5, fluxDep=0.1) + + genList = [[cell1, cell2], [cell3, cell4, cell5], [cell1]] + rasterAttributes = makeRasterAttributes(nrows=5, ncols=5) + + p = Path(dem, 0, 0, genList, rasterAttributes, countArray, relId=1) + p.getPathArrays() + + # cell1 is hit twice (gen 0 and gen 2) at (0,0), with identical values both times + assert p.zDeltaArray[0, 0] == 5.0 + assert p.flowEnergyArray[0, 0] == 50.0 + assert p.fluxArray[0, 0] == pytest.approx(1.0) + # sums accumulate over both hits + assert p.routFluxSumArray[0, 0] == pytest.approx(1.0 + 1.0) + assert p.depFluxSumArray[0, 0] == pytest.approx(0.5 + 0.5) + + # cells hit only once, spread across generations + assert p.zDeltaArray[1, 1] == 2.0 + assert p.zDeltaArray[2, 2] == 3.0 + assert p.zDeltaArray[3, 3] == 4.0 + assert p.zDeltaArray[4, 4] == 1.0 + assert p.routFluxSumArray[4, 4] == pytest.approx(0.5) + + +def test_calcThalwegCenterof_threeGenerations(): + """ test calcThalwegCenterof directly using values shaped like a 3-generation genList + (generation sizes 2, 3, 1) with flux as the weighting variable + """ + p = _minimalPath() + p.numberGen = 3 + + # zDelta values per generation, matching [[cell1,cell2],[cell3,cell4,cell5],[cell1]] + zDeltaGen = [[1.0, 2.0], [3.0, 4.0, 5.0], [1.0]] + fluxGen = [[1.0, 2.0], [3.0, 4.0, 5.0], [1.0]] + + variableSum, coVar = p.calcThalwegCenterof(zDeltaGen, fluxGen) + + # gen 0: sum=3.0, weighted avg = (1*1+2*2)/(1+2) = 5/3 + # gen 1: sum=12.0, weighted avg = (3*3+4*4+5*5)/(3+4+5) = 50/12 + # gen 2: sum=1.0, weighted avg = 1.0 (single value) + np.testing.assert_allclose(variableSum, [3.0, 12.0, 1.0]) + np.testing.assert_allclose(coVar, [5.0 / 3.0, 50.0 / 12.0, 1.0]) diff --git a/avaframe/tests/test_geoTrans.py b/avaframe/tests/test_geoTrans.py index ea81ced81..ac4fc2e97 100644 --- a/avaframe/tests/test_geoTrans.py +++ b/avaframe/tests/test_geoTrans.py @@ -1266,3 +1266,50 @@ def test_getNormalMesh(capfd): atol=atol, ) assert TestNZ + + +def test_indicesToCoords_basic(): + """ test indicesToCoords with simple round-number values """ + header = {"cellsize": 10.0, "xllcenter": 0.0, "yllcenter": 0.0} + + col = np.array([0.0, 1.0, 2.5]) + row = np.array([0.0, 1.0, 2.5]) + + x, y = geoTrans.indicesToCoords(col, row, header) + + # xllcorner = 0 - 5 = -5.0, yllcorner = -5.0 + # x = xllcorner + col * cellsize = -5.0 + col*10.0 + expectedX = np.array([-5.0, 5.0, 20.0]) + expectedY = np.array([-5.0, 5.0, 20.0]) + + np.testing.assert_allclose(x, expectedX) + np.testing.assert_allclose(y, expectedY) + + +def test_indicesToCoords_nonZeroLowerLeftCenter(): + """ test indicesToCoords with a non-zero xllcenter/yllcenter """ + header = {"cellsize": 5.0, "xllcenter": 100.0, "yllcenter": 200.0} + + col = np.array([0.0, 2.0]) + row = np.array([0.0, 3.0]) + + x, y = geoTrans.indicesToCoords(col, row, header) + + # xllcorner = 100 - 2.5 = 97.5, yllcorner = 200 - 2.5 = 197.5 + expectedX = 97.5 + col * 5.0 + expectedY = 197.5 + row * 5.0 + + np.testing.assert_allclose(x, expectedX) + np.testing.assert_allclose(y, expectedY) + + +def test_indicesToCoords_scalarInputs(): + """ test indicesToCoords also works with plain Python floats (not just arrays) """ + header = {"cellsize": 2.0, "xllcenter": 0.0, "yllcenter": 0.0} + + x, y = geoTrans.indicesToCoords(5.0, 2.0, header) + + # xllcorner = -1.0, yllcorner = -1.0 + # x = -1.0 + 5.0*2.0 = 9.0, y = -1.0 + 2.0*2.0 = 3.0 + assert x == pytest.approx(9.0) + assert y == pytest.approx(3.0) diff --git a/docs/moduleCom4FlowPy.rst b/docs/moduleCom4FlowPy.rst index 02b494f69..112a59d8f 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -123,6 +123,7 @@ ii) additional modules (forest, infrastructure) - ``forest``: if set to ``True`` the runout calculation is performed with the *forest module* (a forest layer has to be provided) - ``infra``: if set to ``True`` the calculation is performend with the *backcalculation module* (an infrastructure layer has to be provided) +- ``calcGeneration``: if set to ``True`` the calculation (iteration per cell) is done per generation (iteration). The results can vary. This computation is required for deriving thalwegs. if ``infra`` is set to ``True`` the infrastructure layer has to be provided either in ``avalancheDir/INPUTS/INFRA`` (if ``useCustomPaths=False``) or at the defined ``infraPath`` (if ``useCustomPaths=True``). The layer has to be of the same resolution and extent as the other input layers; infrastructure cells have to be coded with values > 0, while @@ -187,6 +188,7 @@ If ``forestInteraction = True``, an additional output Layer is computed, which r iv) variable parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + There are options to set for each path variable parameters: - alpha (``variableAlpha = True``), @@ -198,7 +200,59 @@ If the value of the variable layer in the cell that is assigned to a release cel When ``variableUmaxLim = True``, the type of the provided parameter is required: ``varUmaxParameter = uMax`` (in m/s) or ``varUmaxParameter = zDeltaMax`` (in m). (A layer containing release cells is still required). -v) tiling and multiprocessing parameters + +v) thalweg output +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +To generate thalweg data, set both ``calcGeneration`` and ``calcThalweg`` to ``True``. + +.. note:: + + Enabling ``calcGeneration = True`` may lead to slight differences in the results because the raster cells are processed in a different iteration order. + +The ``thalwegReleaseArea`` parameter determines how thalwegs are generated: + +* If ``thalwegReleaseArea is True``, one thalweg is computed for each continuous release area. In this case, the ``Inputs/RELID`` directory must contain a raster defining the release area IDs, where each continuous release area has a unique ID. +* If ``thalwegReleaseArea is False``, one thalweg is computed for each start cell. + +The generated thalweg data are stored in a separate ``thalwegData`` directory within the output folder. +Each computed thalweg is saved as a dictionary in a separate pickle file. + +The thalweg can be calculated using one of three center definitions at each iteration step that is selected using the ``thalwegCenterOf`` parameter: + +* center of energy +* center of flux +* center of velocity altitude (``zDelta``) + +Each output dictionary always contains the Flow-Py input parameters ``alpha``, ``zDeltaMax``, +and ``exponent`` used for the path simulation. Additional variables can be specified using the +``thalwegVariables`` parameter. +The following variables are available: + +* ``col`` - column index of the thalweg +* ``row`` - row index of the thalweg +* ``x`` - x coordinate of the thalweg (in the coordinate system of the PRA-raster) +* ``y`` - y coordinate of the thalweg (in the coordinate system of the PRA-raster) +* ``z`` or ``altitude`` - elevation along thalweg +* ``flux`` - flux along thalweg +* ``zDelta`` - velocity altitude along thalweg +* ``flowEnergy`` - flow energy (:math:`= flux * zdelta * g`, with gravitational acceleration :math:`g`) along thalweg +* ``gamma`` - travel angle along thalweg +* ``s`` or ``travelLength`` - horizontally projected travel length along thalweg +* ``flowEnergyArray`` - flow energy of the path (2 dimensional array) +* ``zDeltaArray`` - velocity altitude of the path (2 dimensional array) +* ``fluxArray`` - flux of the path (2 dimensional array) + +When a thalweg is computed and ``thalwegReleaseArea is True``, a GIF of the evolution of the flow path and the thalweg can be created. +The data for the GIF are stored, when a respective ``videoRelId`` is selected. +To create the GIF, modify ``out3Plot/(local_)outCom4GifCfg.ini`` respectively and run: + + :: + + python runScripts/runCreateGIFCom4.py + + +vi) tiling and multiprocessing parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ If the model extent (i.e. number of cells and/or rows in the input layers) is larger than ``tileSize``, then :py:mod:`com4FlowPy`