From 81d5952992af03360e079aaf876f50ea65eebe75 Mon Sep 17 00:00:00 2001 From: paulasp Date: Tue, 8 Jul 2025 17:17:36 +0200 Subject: [PATCH] add output minimum of travelangle and travellength fix bug fix bug doc for additional output files small adaptation null data values are 0 Revert "null data values are 0" This reverts commit aa46b4b38f3a200368d877b67402ce05bb490ae0. travelLength, fpTravelAngle, slTravelAngle, flux have value -9999 in cells that are not hit by a process Update docs/moduleCom4FlowPy.rst Co-authored-by: ahuber-bfw <112931921+ahuber-bfw@users.noreply.github.com> Update docs/moduleCom4FlowPy.rst Co-authored-by: ahuber-bfw <112931921+ahuber-bfw@users.noreply.github.com> minor minor mod to runCom4FlowPy.py to include new outputs as valid output options only comupte rasters when specified in outputs add test test calculation --- avaframe/com4FlowPy/com4FlowPy.py | 56 ++++++--- avaframe/com4FlowPy/com4FlowPyCfg.ini | 6 +- avaframe/com4FlowPy/flowCore.py | 171 +++++++++++++++++++------- avaframe/runCom4FlowPy.py | 3 +- avaframe/tests/test_com4FlowPy.py | 51 +++++++- docs/moduleCom4FlowPy.rst | 6 +- 6 files changed, 227 insertions(+), 66 deletions(-) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index 99403b7fc..146eaa63d 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -518,13 +518,15 @@ def mergeAndWriteResults(modelPaths, modelOptions): # Merge calculated tiles zDelta = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta") flux = SPAM.mergeRaster(modelPaths["tempDir"], "res_flux") - cellCounts = SPAM.mergeRaster(modelPaths["tempDir"], "res_count", method='sum') - zDeltaSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta_sum", method='sum') - routFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_rout_flux_sum", method='sum') - depFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_dep_flux_sum", method='sum') - fpTa = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp") + cellCounts = SPAM.mergeRaster(modelPaths["tempDir"], "res_count", method="sum") + zDeltaSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta_sum", method="sum") + routFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_rout_flux_sum", method="sum") + depFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_dep_flux_sum", method="sum") + fpTaMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp_max") + fpTaMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp_min", method="min") slTa = SPAM.mergeRaster(modelPaths["tempDir"], "res_sl") - travelLength = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length") + travelLengthMax = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length_max") + travelLengthMin = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length_min", method="min") if modelOptions["infraBool"]: backcalc = SPAM.mergeRaster(modelPaths["tempDir"], "res_backcalc") @@ -542,7 +544,7 @@ def mergeAndWriteResults(modelPaths, modelOptions): demHeader = IOf.readRasterHeader(modelPaths["demPath"]) outputHeader = demHeader.copy() outputHeader["nodata_value"]=-9999 - + if _oF == ".asc": outputHeader["driver"] = "AAIGrid" elif _oF == ".tif": @@ -566,20 +568,42 @@ def mergeAndWriteResults(modelPaths, modelOptions): if 'depFluxSum' in _outputs: output = IOf.writeResultToRaster(outputHeader, depFluxSum, modelPaths["resDir"] / "com4_{}_{}_depFluxSum".format(_uid, _ts), flip=True) - if 'fpTravelAngle' in _outputs: - output = IOf.writeResultToRaster(outputHeader, fpTa, - modelPaths["resDir"] / "com4_{}_{}_fpTravelAngle".format(_uid, _ts), flip=True) + if "fpTravelAngle" in _outputs or "fpTravelAngleMax" in _outputs: + output = IOf.writeResultToRaster( + outputHeader, + fpTaMax, + modelPaths["resDir"] / "com4_{}_{}_fpTravelAngleMax".format(_uid, _ts), + flip=True, + ) + if "fpTravelAngleMin" in _outputs: + output = IOf.writeResultToRaster( + outputHeader, + fpTaMin, + modelPaths["resDir"] / "com4_{}_{}_fpTravelAngleMin".format(_uid, _ts), + flip=True, + ) if 'slTravelAngle' in _outputs: output = IOf.writeResultToRaster(outputHeader, slTa, modelPaths["resDir"] / "com4_{}_{}_slTravelAngle".format(_uid, _ts), flip=True) - if 'travelLength' in _outputs: - output = IOf.writeResultToRaster(outputHeader, travelLength, - modelPaths["resDir"] / "com4_{}_{}_travelLength".format(_uid, _ts), flip=True) + if "travelLength" in _outputs or "travelLengthMax" in _outputs: + output = IOf.writeResultToRaster( + outputHeader, + travelLengthMax, + modelPaths["resDir"] / "com4_{}_{}_travelLengthMax".format(_uid, _ts), + flip=True, + ) + if "travelLengthMin" in _outputs: + output = IOf.writeResultToRaster( + outputHeader, + travelLengthMin, + modelPaths["resDir"] / "com4_{}_{}_travelLengthMin".format(_uid, _ts), + flip=True, + ) # NOTE: # if not modelOptions["infraBool"]: # if no infra - # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("cell_counts%s" %(output_format)),cell_counts) - # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("z_delta_sum%s" %(output_format)),z_delta_sum) + # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("cell_counts%s" %(output_format)),cell_counts) + # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("z_delta_sum%s" %(output_format)),z_delta_sum) if modelOptions["infraBool"]: # if infra output = IOf.writeResultToRaster(outputHeader, backcalc, modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts), flip=True) @@ -676,4 +700,4 @@ def deleteTempFolder(tempFolderPath): log.info("deletion of temp folder {} failed".format(tempFolderPath)) log.info(" isDir:{} isTemp:{}}".format(isDir, validTemp)) - log.info("+++++++++++++++++++++++") \ No newline at end of file + log.info("+++++++++++++++++++++++") diff --git a/avaframe/com4FlowPy/com4FlowPyCfg.ini b/avaframe/com4FlowPy/com4FlowPyCfg.ini index 96f7b9e36..19e17d887 100644 --- a/avaframe/com4FlowPy/com4FlowPyCfg.ini +++ b/avaframe/com4FlowPy/com4FlowPyCfg.ini @@ -205,16 +205,18 @@ maxChunks = 500 outputFileFormat = .tif # define the different output files that are written to disk -# default = 'zDelta|cellCounts|travelLength|fpTravelAngle' +# default = 'zDelta|cellCounts|travelLengthMax|fpTravelAngleMax' # additional options: # slTravelAngle # flux # zDeltaSum # routFluxSum # depFluxSum +# travelLengthMin +# fpTravelAngleMin # if forestInteraction: forestInteraction is automatically added to outputs # if infra: backCalculation is automatically added to output -outputFiles = zDelta|cellCounts|travelLength|fpTravelAngle +outputFiles = zDelta|cellCounts|travelLengthMax|fpTravelAngleMax #++++++++++++ Custom paths True/False # default: False diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index 9876e315b..731195178 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -151,6 +151,9 @@ def run(optTuple): # Temp-Dir (all input files are located here and results are written back in here) tempDir = optTuple[3]["tempDir"] + # List of output layers + outputs = optTuple[3]["outputFileList"] + # raster-layer Attributes cellsize = float(optTuple[4]["cellsize"]) nodata = float(optTuple[4]["nodata"]) @@ -229,14 +232,24 @@ def run(optTuple): results = pool.map( calculation, [ - [ # TODO: write in dicts: - dem, infra, release_sub, - alpha, exp, flux_threshold, max_z_delta, - nodata, cellsize, - infraBool, forestBool, - varParams, fluxDistOldVersionBool, + [ # TODO: write in dicts: + dem, + infra, + release_sub, + alpha, + exp, + flux_threshold, + max_z_delta, + nodata, + cellsize, + infraBool, + forestBool, + varParams, + fluxDistOldVersionBool, previewMode, - forestArray, forestParams, + forestArray, + forestParams, + outputs, ] for release_sub in release_list ], @@ -248,16 +261,18 @@ def run(optTuple): # Move this part into a separate function # initializing arrays for storing the results from the multiprocessing step zDeltaArray = np.zeros_like(dem, dtype=np.float32) - fluxArray = np.zeros_like(dem, dtype=np.float32) + fluxArray = np.ones_like(dem, dtype=np.float32) * -9999 countArray = np.zeros_like(dem, dtype=np.int32) zDeltaSumArray = np.zeros_like(dem, dtype=np.float32) routFluxSumArray = np.zeros_like(dem, dtype=np.float32) depFluxSumArray = np.zeros_like(dem, dtype=np.float32) if infraBool: backcalc = np.ones_like(dem, dtype=np.int32) * -9999 - fpTravelAngleArray = np.zeros_like(dem, dtype=np.float32) - slTravelAngleArray = np.zeros_like(dem, dtype=np.float32) - travelLengthArray = np.zeros_like(dem, dtype=np.float32) + fpTravelAngleMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 + fpTravelAngleMinArray = np.ones_like(dem, dtype=np.float32) * -9999 + slTravelAngleArray = np.ones_like(dem, dtype=np.float32) * -9999 + travelLengthMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 + travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 if forestInteraction: forestIntArray = np.ones_like(dem, dtype=np.float32) * -9999 @@ -269,9 +284,11 @@ def run(optTuple): depFluxSumList = [] if infraBool: backcalcList = [] - fpTravelAngleList = [] + fpTravelAngleMaxList = [] + fpTravelAngleMinList = [] slTravelAngleList = [] - travelLengthList = [] + travelLengthMaxList = [] + travelLengthMinList = [] if forestInteraction: forestIntList = [] @@ -284,13 +301,15 @@ def run(optTuple): zDeltaSumList.append(res[3]) if infraBool: backcalcList.append(res[4]) - fpTravelAngleList.append(res[5]) + fpTravelAngleMaxList.append(res[5]) slTravelAngleList.append(res[6]) - travelLengthList.append(res[7]) - routFluxSumList.append(res[8]) - depFluxSumList.append(res[9]) + travelLengthMaxList.append(res[7]) + travelLengthMinList.append(res[8]) + fpTravelAngleMinList.append(res[9]) + routFluxSumList.append(res[10]) + depFluxSumList.append(res[11]) if forestInteraction: - forestIntList.append(res[10]) + forestIntList.append(res[12]) logging.info("Calculation finished, getting results.") for i in range(len(zDeltaList)): @@ -302,9 +321,23 @@ def run(optTuple): depFluxSumArray += depFluxSumList[i] if infraBool: backcalc = np.maximum(backcalc, backcalcList[i]) - fpTravelAngleArray = np.maximum(fpTravelAngleArray, fpTravelAngleList[i]) + if "fpTravelAngleMax" in outputs or "fpTravelAngle" in outputs: + fpTravelAngleMaxArray = np.maximum(fpTravelAngleMaxArray, fpTravelAngleMaxList[i]) + if "fpTravelAngleMin" in outputs: + fpTravelAngleMinArray = np.where( + (travelLengthMinArray >= 0) & (fpTravelAngleMinList[i] >= 0), + np.minimum(fpTravelAngleMinArray, fpTravelAngleMinList[i]), + np.maximum(fpTravelAngleMinArray, fpTravelAngleMinList[i]), + ) slTravelAngleArray = np.maximum(slTravelAngleArray, slTravelAngleList[i]) - travelLengthArray = np.maximum(travelLengthArray, travelLengthList[i]) + if "travelLengthMax" in outputs or "travelLength" in outputs: + travelLengthMaxArray = np.maximum(travelLengthMaxArray, travelLengthMaxList[i]) + if "travelLengthMin" in outputs: + travelLengthMinArray = np.where( + (travelLengthMinArray >= 0) & (travelLengthMinList[i] >= 0), + np.minimum(travelLengthMinArray, travelLengthMinList[i]), + np.maximum(travelLengthMinArray, travelLengthMinList[i]), + ) if forestInteraction: forestIntArray = np.where((forestIntArray >= 0) & (forestIntList[i] >= 0), np.minimum(forestIntArray, forestIntList[i]), @@ -317,9 +350,11 @@ def run(optTuple): 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_%s_%s" % (optTuple[0], optTuple[1])), fpTravelAngleArray) + 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_%s_%s" % (optTuple[0], optTuple[1])), travelLengthArray) + 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: @@ -352,6 +387,7 @@ def calculation(args): - args[14] (numpy array) - contains forest information (None if forestBool=False) - args[15] (dict) - contains parameters for forest interaction models (None if forestBool=False) + - args[16] (list) - output names Returns ----------- @@ -418,6 +454,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): varExponentArray = args[11]['varExponentArray'] fluxDistOldVersionBool = args[12] previewMode = args[13] + outputs = args[16] if forestBool: forestArray = args[14] @@ -433,13 +470,15 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): zDeltaPathList = [] routFluxSumArray = np.zeros_like(dem, dtype=np.float32) depFluxSumArray = np.zeros_like(dem, dtype=np.float32) - fluxArray = np.zeros_like(dem, dtype=np.float32) + fluxArray = np.ones_like(dem, dtype=np.float32) * -9999 countArray = np.zeros_like(dem, dtype=np.int32) - fpTravelAngleArray = np.zeros_like(dem, dtype=np.float32) # fp = Flow Path - slTravelAngleArray = np.zeros_like(dem, dtype=np.float32) # sl = Straight Line + fpTravelAngleMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 # fp = Flow Path + fpTravelAngleMinArray = np.ones_like(dem, dtype=np.float32) * -9999 # fp = Flow Path + slTravelAngleArray = np.ones_like(dem, dtype=np.float32) * -9999 # sl = Straight Line - travelLengthArray = np.zeros_like(dem, dtype=np.float32) + travelLengthMinArray = np.ones_like(dem, dtype=np.float32) * -9999 + travelLengthMaxArray = np.ones_like(dem, dtype=np.float32) * -9999 if infraBool: backcalc = np.ones_like(dem, dtype=np.int32) * -9999 @@ -515,7 +554,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): 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))) - + if infraBool: # if the current cell is not already in the dir-graph, then we add it here updateInfraDirGraph(cell.rowindex, cell.colindex) @@ -530,7 +569,7 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): if infraBool: updateInfraDirGraph(row[k], col[k], cell.rowindex, cell.colindex) - + if z_delta[k] > cell_list[i].z_delta: cell_list[i].z_delta = z_delta[k] row = np.delete(row, k) @@ -576,12 +615,34 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): 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) - fpTravelAngleArray[cell.rowindex, cell.colindex] = max(fpTravelAngleArray[cell.rowindex, cell.colindex], - cell.max_gamma) + 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) - travelLengthArray[cell.rowindex, cell.colindex] = max(travelLengthArray[cell.rowindex, cell.colindex], - cell.min_distance) + 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) @@ -601,10 +662,10 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): for key, val in updatedInfraValues.items(): backcalc[key[0], key[1]] = max(backcalc[key[0], key[1]], val) # writing max-values to back-tracking array - + del pathTopology, infraValues, updatedInfraValues gc.collect() - + if previewMode: # if the 'previewMode' is On/'True', then we check here if the current modeled process zones already # includes other release Cells (i.e. if release cells are "hit from above") @@ -624,11 +685,36 @@ def updateInfraDirGraph(row, col, parentRow=None, parentCol=None): gc.collect() if forestInteraction: - return zDeltaArray, fluxArray, countArray, zDeltaSumArray, backcalc, fpTravelAngleArray, slTravelAngleArray, \ - travelLengthArray, routFluxSumArray, depFluxSumArray, forestIntArray + return ( + zDeltaArray, + fluxArray, + countArray, + zDeltaSumArray, + backcalc, + fpTravelAngleMaxArray, + slTravelAngleArray, + travelLengthMaxArray, + travelLengthMinArray, + fpTravelAngleMinArray, + routFluxSumArray, + depFluxSumArray, + forestIntArray, + ) else: - return zDeltaArray, fluxArray, countArray, zDeltaSumArray, backcalc, fpTravelAngleArray, slTravelAngleArray, \ - travelLengthArray, routFluxSumArray, depFluxSumArray + return ( + zDeltaArray, + fluxArray, + countArray, + zDeltaSumArray, + backcalc, + fpTravelAngleMaxArray, + slTravelAngleArray, + travelLengthMaxArray, + travelLengthMinArray, + fpTravelAngleMinArray, + routFluxSumArray, + depFluxSumArray, + ) def enoughMemoryAvailable(limit=0.05): @@ -762,7 +848,7 @@ def backTracking(topologyDict, infraValDict): valDictSorted = {k: v for k, v in sorted(infraValDict.items(), key= lambda item: item[1], reverse=True)} # reverse the graph topology, so "parents" become "children" reverseGraph = reverseTopology(topologyDict) - + # helper function to recursively traverse the reverseGraph and # propagate the infraValues "upslope" def propagateInfraVal(node, infraValToPropagate, visited): @@ -781,13 +867,12 @@ def propagateInfraVal(node, infraValToPropagate, visited): for parentNode in reverseGraph.get(node, []): propagateInfraVal(parentNode, infraValToPropagate, visited) - for node, val in valDictSorted.items(): if val > 0: propagateInfraVal(node, val, set()) - + return infraValDict - + def reverseTopology(topologyDict): ''' reverse graph topology @@ -815,5 +900,3 @@ def reverseTopology(topologyDict): reverseGraph[child].append(parentNode) return reverseGraph - - diff --git a/avaframe/runCom4FlowPy.py b/avaframe/runCom4FlowPy.py index cdb157af1..3fbde5a2d 100644 --- a/avaframe/runCom4FlowPy.py +++ b/avaframe/runCom4FlowPy.py @@ -321,7 +321,8 @@ def checkOutputFilesFormat(strOutputFiles): try: setA = set(strOutputFiles.split('|')) - setB = set(['zDelta', 'cellCounts', 'fpTravelAngle', 'travelLength', + setB = set(['zDelta', 'cellCounts', 'fpTravelAngle', 'travelLength', 'fpTravelAngleMax', + 'fpTravelAngleMin', 'travelLengthMax', 'travelLengthMin', 'slTravelAngle', 'flux', 'zDeltaSum', 'routFluxSum', 'depFluxSum']) # if there is at least 1 correct outputfile defined, we use the string provided in the .ini file if (setA & setB): diff --git a/avaframe/tests/test_com4FlowPy.py b/avaframe/tests/test_com4FlowPy.py index 9ed33e94d..aabd8c4d1 100644 --- a/avaframe/tests/test_com4FlowPy.py +++ b/avaframe/tests/test_com4FlowPy.py @@ -5,6 +5,7 @@ # Load modules import numpy as np import pytest +from numpy.ma.core import ones_like from avaframe.com4FlowPy import flowClass import avaframe.com4FlowPy.flowCore as flowCore @@ -170,7 +171,55 @@ def test_backTracking(): assert calcValsBT[key] == testValsBT[key] +def test_calculation(): + dem = np.array( + [[40, 40, 40, 40, 40], [30, 30, 30, 30, 30], [20, 20, 20, 20, 20], [10, 10, 10, 10, 10], [0, 0, 0, 0, 0]]) + infra = None + pra = np.array([[0, 0, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]]) + alpha = 10 + exp = 99 + fluxTh = 0.001 + zDeltaMax = 8000 + nodata = -9999 + cellsize = 10 + infraBool = False + forestBool = False + variableParameters = { + "varUmaxBool": False, + "varUmaxArray": None, + "varAlphaBool": False, + "varAlphaArray": None, + "varExponentBool": False, + "varExponentArray": None, + } + fluxDistOldVersionBool = False + previewMode = False + outputs = ['travelLengthMin', 'flux'] + forestArray = None + forestParams = None + args = [dem, infra, pra, alpha, exp, fluxTh, zDeltaMax, nodata, cellsize, infraBool, forestBool, variableParameters, + fluxDistOldVersionBool, previewMode, forestArray, forestParams, outputs] + + flux = ones_like(dem) * -9999. + flux[[1, 2, 3], [2, 2, 2]] = 1 + depFluxSum = np.zeros_like(dem) + routFluxSum = np.where(flux == 1, 1., 0.) + travelLengthMin = ones_like(dem) * -9999. + travelLengthMin[1, 2] = 0 + travelLengthMin[2, 2] = np.sqrt(cellsize ** 2) + travelLengthMin[3, 2] = 2 * np.sqrt(cellsize ** 2) + results = flowCore.calculation(args) + + assert len(results) == 12 + assert np.all(results[1] == flux) + assert np.all(results[10] == routFluxSum) + assert np.all(results[11] == depFluxSum) + assert np.all(results[8] == travelLengthMin) + assert np.all(results[7] == np.ones_like(flux) * -9999) + + if __name__=='__main__': test_add_os() test_reverseTopology() - test_backTracking() \ No newline at end of file + test_backTracking() + test_calculation() diff --git a/docs/moduleCom4FlowPy.rst b/docs/moduleCom4FlowPy.rst index 25bf419f9..28514dcff 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -229,8 +229,8 @@ By default the following four output layers are written to disk at the end of th - ``zdelta``: the maximum z_delta of all paths for every raster cell (geometric measure of process magnitude, can be associated to kinetic energy/velocity) - ``cellCounts``: number of paths/release cells that route flux through a raster cell -- ``travelLength``: the travel length along the flow path -- ``fpTravelAngle``: the gamma angle along the flow path +- ``travelLengthMax``: the travel length along the flow path (the maximum value of all paths for every raster cell) - **this output is equivalent to the ``travelLength`` output from previous code versions!** +- ``fpTravelAngleMax``: the gamma angle along the flow path (the maximum value of all paths for every raster cell) - **this output is equivalent to the ``fpTravelAngle``output used in previous versions of the code!!** In addition these output layers are also available: @@ -239,6 +239,8 @@ In addition these output layers are also available: - ``slTravelAngle``: gamma angle calculated along a straight-line between release cell and current cell - ``routFluxSum``: routing flux summed up over all paths - ``depFluxSum``: deposited flux summed up over all paths +- ``travelLengthMin``: the travel length along the flow path (the minimum value of all paths for every raster cell) +- ``fpTravelAngleMin``: the gamma angle along the flow path (the minimum value of all paths for every raster cell) If ``forestInteraction = True`` this layer will be written automatically (no need to separately define in ``outputFiles``):