diff --git a/avaframe/ana5Utils/DFAPathGeneration.py b/avaframe/ana5Utils/DFAPathGeneration.py index 625ac29cb..1157190fb 100644 --- a/avaframe/ana5Utils/DFAPathGeneration.py +++ b/avaframe/ana5Utils/DFAPathGeneration.py @@ -4,13 +4,16 @@ # Load modules import math +import heapq import numpy as np import logging import pathlib import shutil +from scipy.ndimage import distance_transform_edt # Local imports from avaframe.in1Data import getInput as gI +import avaframe.in2Trans.rasterUtils as IOf import avaframe.in2Trans.shpConversion as shpConv from avaframe.in3Utils import cfgUtils from avaframe.in3Utils import fileHandlerUtils as fU @@ -80,10 +83,22 @@ def generatePathAndSplitpoint(avalancheDir, cfgDFAPath, cfgMain, runDFAModule): log.error(message) raise FileExistsError(message) + # the peak flow thickness fields are only needed when the path is extended to the deposit + # front; parse the peak files once here instead of re-reading them for every simulation + extendToFront = cfgDFAPath['PATH'].getint('extBottomOption', fallback=0) == 1 + peakFilesDF = None + if extendToFront: + peakFilesDir = pathlib.Path(avalancheDir, 'Outputs', 'com1DFA', 'peakFiles') + peakFilesDF = fU.makeSimDF(peakFilesDir, avaDir=avalancheDir) + for simName, simDFrow in simDF.iterrows(): log.info('Computing avalanche path from simulation: %s', simName) pathFromPart = cfgDFAPath['PATH'].getboolean('pathFromPart') resampleDistance = cfgDFAPath['PATH'].getfloat('nCellsResample') * dem['header']['cellsize'] + # peak flow thickness field, only needed when extending the path to the deposit front + fieldPFT = None + if extendToFront: + fieldPFT = readPeakFT(peakFilesDF, simName) # get the mass average path avaProfileMass, particlesIni = generateMassAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInfo=cfgDFAPath['PATH'].getboolean('addVelocityInfo')) @@ -93,7 +108,8 @@ def generatePathAndSplitpoint(avalancheDir, cfgDFAPath, cfgMain, runDFAModule): # make the parabolic fit parabolicFit = getParabolicFit(cfgDFAPath['PATH'], avaProfileMass, dem) # here the avaProfileMass given in input is overwritten and returns only an x, y, z extended profile - avaProfileMass = extendDFAPath(cfgDFAPath['PATH'], avaProfileMass, dem, particlesIni) + avaProfileMass = extendDFAPath(cfgDFAPath['PATH'], avaProfileMass, dem, particlesIni, + fieldPFT=fieldPFT) # resample path and keep track of start and end of mass averaged part avaProfileMass = resamplePath(cfgDFAPath['PATH'], dem, avaProfileMass) # get split point @@ -306,7 +322,33 @@ def getMassAvgPathFromFields(fieldsList, fieldHeader, dem): return avaProfileMass -def extendDFAPath(cfg, avaProfile, dem, particlesIni): +def readPeakFT(peakFilesDF, simName): + """ get the peak flow thickness field (pft) of one simulation from the peak file dataframe + + Parameters + ----------- + peakFilesDF: pandas DataFrame + peak files of all simulations (from fU.makeSimDF), parsed once for the whole run + simName: str + simulation name + + Returns + -------- + fieldPFT: numpy array + peak flow thickness raster (same grid as the simulation dem), + None if no pft peak field is available for this simulation + """ + # the simulation can be identified by its full name or by its hash (the index of the + # configuration dataframe iterated in generatePathAndSplitpoint is the hash) + isSim = (peakFilesDF['simName'] == simName) | (peakFilesDF['simID'] == simName) + index = peakFilesDF.index[isSim & (peakFilesDF['resType'] == 'pft')] + if len(index) == 0: + log.warning('No pft peak field found for simulation %s' % simName) + return None + return IOf.readRaster(peakFilesDF.loc[index[0], 'files'])['rasterData'] + + +def extendDFAPath(cfg, avaProfile, dem, particlesIni, fieldPFT=None): """ extend the DFA path at the top and bottom avaProfile with x, y, z, s information @@ -316,6 +358,8 @@ def extendDFAPath(cfg, avaProfile, dem, particlesIni): configuration object with: - extTopOption: int, how to extend towards the top? 0 for heighst point method, a for largest runout method + - extBottomOption: int, how to extend towards the bottom? 0 for the straight-line extrapolation + method, 1 for the least-cost extension to the deposit front (requires fieldPFT) - nCellsResample: int, resampling length is given by nCellsResample*demCellSize - nCellsMinExtend: int, when extending towards the bottom, take points at more than nCellsMinExtend*demCellSize from last point to get the direction @@ -329,6 +373,8 @@ def extendDFAPath(cfg, avaProfile, dem, particlesIni): dem dict particlesIni: dict initial particles dict + fieldPFT: numpy array, optional + peak flow thickness field on the dem grid, only used if extBottomOption = 1 Returns -------- @@ -339,7 +385,15 @@ def extendDFAPath(cfg, avaProfile, dem, particlesIni): resampleDistance = cfg.getfloat('nCellsResample') * dem['header']['cellsize'] avaProfile, _ = gT.prepareLine(dem, avaProfile, distance=resampleDistance, Point=None) avaProfile = extendProfileTop(cfg.getint('extTopOption'), particlesIni, avaProfile) - avaProfile = extendProfileBottom(cfg, dem, avaProfile) + if cfg.getint('extBottomOption', fallback=0) == 1: + if fieldPFT is None: + log.warning('extBottomOption is 1 but no peak flow thickness field was provided, ' + 'falling back to the straight-line bottom extension') + avaProfile = extendProfileBottom(cfg, dem, avaProfile) + else: + avaProfile = extendProfileToFront(cfg, dem, avaProfile, fieldPFT) + else: + avaProfile = extendProfileBottom(cfg, dem, avaProfile) return avaProfile @@ -523,6 +577,205 @@ def extendProfileBottom(cfg, dem, profile): return profile +def extendProfileToFront(cfg, dem, profile, fieldPFT): + """ extend the DFA path at the bottom to the front of the deposit + + Locate the front of the deposit in the peak flow thickness field (findFlowFront) + and extend the profile to it along a least-cost path over the dem + (leastCostPath). In contrast to extendProfileBottom, the extension follows + the terrain and the deposit and stops at the front instead of extrapolating + a straight line of fixed relative length. If the front cannot be located or + reached, the profile falls back to the straight-line extension + (extendProfileBottom). + + Parameters + ----------- + cfg: configParser + configuration object with: + + - ftThreshold: float, minimum flow thickness (m) for a cell to belong to the flow footprint + - lowFrontFraction: float, fraction of the flow elevation range defining the front band + - upSlopePenalty: float, cost multiplier on the positive elevation gain of an edge + - flowDistPenalty: float, cost multiplier on the distance (m) of a cell from the flow footprint + dem: dict + dem dict + profile: dict + profile to extend + fieldPFT: numpy array + peak flow thickness field (pft) on the same grid as the dem + + Returns + -------- + profile: dict + extended profile (x, y, z, s) + """ + header = dem['header'] + csz = header['cellsize'] + zRaster = dem['rasterData'] + if fieldPFT.shape != zRaster.shape: + message = 'Peak flow thickness field and dem do not have the same shape' + log.error(message) + raise AssertionError(message) + # get last point + xLast = profile['x'][-1] + yLast = profile['y'][-1] + sLast = profile['s'][-1] + # locate the deposit front + frontRow, frontCol = findFlowFront(fieldPFT, zRaster, cfg.getfloat('ftThreshold'), + cfg.getfloat('lowFrontFraction')) + if frontRow is None: + log.warning('No flow cell above ftThreshold was found, ' + 'falling back to the straight-line bottom extension') + return extendProfileBottom(cfg, dem, profile) + # cell of the last profile point (profile and field share the dem grid) + startRow = min(max(int(round((yLast - header['yllcenter']) / csz)), 0), zRaster.shape[0] - 1) + startCol = min(max(int(round((xLast - header['xllcenter']) / csz)), 0), zRaster.shape[1] - 1) + cellPath = leastCostPath((startRow, startCol), (frontRow, frontCol), fieldPFT, zRaster, csz, + cfg.getfloat('ftThreshold'), cfg.getfloat('upSlopePenalty'), + cfg.getfloat('flowDistPenalty')) + if len(cellPath) < 2: + log.warning('No least-cost path to the deposit front was found, ' + 'falling back to the straight-line bottom extension') + return extendProfileBottom(cfg, dem, profile) + # drop the first cell (the path end itself) and convert to coordinates; the extension + # points are cell centers of valid dem cells, so z is read directly from the raster + rows, cols = np.array(cellPath[1:]).T + xExtBottom = header['xllcenter'] + cols * csz + yExtBottom = header['yllcenter'] + rows * csz + zExtBottom = zRaster[rows, cols] + dx = np.diff(np.append(xLast, xExtBottom)) + dy = np.diff(np.append(yLast, yExtBottom)) + sExtBottom = sLast + np.cumsum(np.sqrt(dx**2 + dy**2)) + log.info('Path extended to the deposit front (%.0f m beyond the mass averaged path end)' + % (sExtBottom[-1] - sLast)) + + # extend profile + profile['x'] = np.append(profile['x'], xExtBottom) + profile['y'] = np.append(profile['y'], yExtBottom) + profile['z'] = np.append(profile['z'], zExtBottom) + profile['s'] = np.append(profile['s'], sExtBottom) + return profile + + +def findFlowFront(fieldPFT, demRaster, ftThreshold, lowFrontFraction): + """ locate the front of the deposit in a peak flow thickness field + + The front is the flow-thickness-weighted centroid of the flow cells lying + in the lowest lowFrontFraction of the flow elevation range. The band always + contains the lowest flow cell; for a flat deposit it covers the whole + footprint, so the front falls back to the centroid of the deposit. If the + centroid falls outside the flow footprint (e.g. between two deposit lobes), + the front is snapped to the nearest cell of the band. + + Parameters + ----------- + fieldPFT: numpy array + peak flow thickness field + demRaster: numpy array + dem raster of the same shape + ftThreshold: float + minimum flow thickness (m) for a cell to belong to the flow footprint + lowFrontFraction: float + fraction of the flow elevation range defining the front band + + Returns + -------- + frontRow, frontCol: int + cell of the front, (None, None) if there is no flow above ftThreshold + """ + flow = (fieldPFT > ftThreshold) & np.isfinite(demRaster) + if not flow.any(): + return None, None + rows, cols = np.nonzero(flow) + elev = demRaster[rows, cols] + lowBand = elev <= elev.min() + lowFrontFraction * (elev.max() - elev.min()) + weight = fieldPFT[rows, cols] * lowBand + # weight.sum() is always positive: the lowest flow cell is in lowBand and, being a flow + # cell, carries fieldPFT > ftThreshold >= 0, so no divide-by-zero guard is needed + frontRow = int(round(np.sum(rows * weight) / weight.sum())) + frontCol = int(round(np.sum(cols * weight) / weight.sum())) + if not flow[frontRow, frontCol]: + # centroid outside the footprint (e.g. two deposit lobes): snap to the band + bandInd = np.flatnonzero(lowBand) + iNear = bandInd[np.argmin((rows[bandInd] - frontRow)**2 + (cols[bandInd] - frontCol)**2)] + frontRow, frontCol = int(rows[iNear]), int(cols[iNear]) + return frontRow, frontCol + + +def leastCostPath(startCell, goalCell, fieldPFT, demRaster, csz, ftThreshold, upSlopePenalty, + flowDistPenalty): + """ Dijkstra least-cost path between two cells of the dem grid + + The cost of an edge is its horizontal length plus a penalty on the positive + elevation gain and a penalty on the distance (m) of the target cell from + the flow footprint, so that the path descends along the deposit. + + Parameters + ----------- + startCell, goalCell: tuple + (row, col) of the start and goal cells + fieldPFT: numpy array + peak flow thickness field used to build the flow footprint + demRaster: numpy array + dem raster of the same shape + csz: float + cell size (m) + ftThreshold: float + minimum flow thickness (m) for a cell to belong to the flow footprint + upSlopePenalty: float + cost multiplier on the positive elevation gain of an edge + flowDistPenalty: float + cost multiplier on the distance (m) of a cell from the flow footprint + + Returns + -------- + cellPath: list + (row, col) cells from start to goal, empty if the goal is unreachable + """ + if upSlopePenalty < 0 or flowDistPenalty < 0: + message = 'The penalty parameters of the least-cost path must not be negative' + log.error(message) + raise AssertionError(message) + # distance (m) from the flow footprint, used to keep the path on the deposit + distToFlow = distance_transform_edt(~(fieldPFT > ftThreshold), sampling=csz) + nRows, nCols = demRaster.shape + demValid = np.isfinite(demRaster) + neighbours = ((-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)) + diagDist = csz * math.sqrt(2.) + costGrid = np.full((nRows, nCols), np.inf) + costGrid[startCell] = 0. + previousCell = np.full((nRows, nCols, 2), -1, dtype=np.int64) + queue = [(0., int(startCell[0]), int(startCell[1]))] + while queue: + cost, row, col = heapq.heappop(queue) + if (row, col) == (goalCell[0], goalCell[1]): + break + if cost > costGrid[row, col]: + continue + for dRow, dCol in neighbours: + nRow, nCol = row + dRow, col + dCol + if 0 <= nRow < nRows and 0 <= nCol < nCols and demValid[nRow, nCol]: + # no diagonal moves across the corner of two nodata cells + if dRow and dCol and not (demValid[row, nCol] and demValid[nRow, col]): + continue + dz = demRaster[nRow, nCol] - demRaster[row, col] + edge = ((diagDist if (dRow and dCol) else csz) + max(dz, 0.) * upSlopePenalty + + distToFlow[nRow, nCol] * flowDistPenalty) + if cost + edge < costGrid[nRow, nCol]: + costGrid[nRow, nCol] = cost + edge + previousCell[nRow, nCol] = (row, col) + heapq.heappush(queue, (cost + edge, nRow, nCol)) + if previousCell[goalCell[0], goalCell[1], 0] < 0 and tuple(goalCell) != tuple(startCell): + return [] + cellPath = [] + row, col = int(goalCell[0]), int(goalCell[1]) + while row >= 0: + cellPath.append((row, col)) + row, col = int(previousCell[row, col, 0]), int(previousCell[row, col, 1]) + cellPath.reverse() + return cellPath + + def getParabolicFit(cfg, avaProfile, dem): """fit a parabola on a set of (s, z) points @@ -652,8 +905,10 @@ def resamplePath(cfg, dem, avaProfile): avaProfile, _ = gT.prepareLine(dem, avaProfile, distance=resampleDistance, Point=None) # make sure we get the good start and end point... prepareLine might make a small error on the s coord indFirst = np.argwhere(avaProfile['s'] >= s0 - resampleDistance/3)[0][0] - # look for the first point in the extension and take the one before - indEnd = np.argwhere(avaProfile['s'] >= sEnd + resampleDistance/3)[0][0]-1 + # look for the first point in the extension and take the one before; if the extension is + # shorter than a resample step, the mass averaged part reaches the last point + indEndCandidates = np.argwhere(avaProfile['s'] >= sEnd + resampleDistance/3) + indEnd = indEndCandidates[0][0]-1 if len(indEndCandidates) > 0 else np.size(avaProfile['s'])-1 avaProfile['indStartMassAverage'] = indFirst avaProfile['indEndMassAverage'] = indEnd return avaProfile diff --git a/avaframe/ana5Utils/DFAPathGenerationCfg.ini b/avaframe/ana5Utils/DFAPathGenerationCfg.ini index 35ef554b4..3e1b91f54 100644 --- a/avaframe/ana5Utils/DFAPathGenerationCfg.ini +++ b/avaframe/ana5Utils/DFAPathGenerationCfg.ini @@ -21,6 +21,12 @@ nCellsResample = 10 # option 1: find the point that will lead to the longest runout extTopOption = 1 +# extension method at the bottom +# option 0: extend the path in the direction of its last points by factBottomExt x sMax +# option 1: extend the path to the front of the deposit with a least-cost path (Dijkstra) over the dem, +# based on the peak flow thickness field (pft) of the simulation +extBottomOption = 0 + # when extending the path at the bottom, extend path in # the direction extracted form the last points of the path # (all points at a distance nCellsMinExtend x cellSize < distance < nCellsMaxExtend x cellSize @@ -30,7 +36,8 @@ nCellsMaxExtend = 20 # this value needs to be chosen in accordance with nCellsResample (if nCellsMaxExtend < nCellsResample no points # might be found for the path extension process) -# for the extrapolation at the bottom, add factBottomExt * sMax to the path (and then check if the extension +# for the extrapolation at the bottom (extBottomOption = 0), add factBottomExt * sMax to the path +# (and then check if the extension # is on the DEM, if not iterate by dichotomy until we find a point on the topography or reach a maximum # iteration number (integer, in this case we pick the last point found inside the dem) # or precision (nBottomExtPrecision * cellSize)) @@ -38,6 +45,17 @@ factBottomExt = 0.3 maxIterationExtBot = 10 nBottomExtPrecision = 10 +# parameters of the bottom extension to the deposit front (extBottomOption = 1) +# minimum flow thickness [m] for a cell to belong to the flow footprint +ftThreshold = 0.01 +# the front is the flow-thickness-weighted centroid of the flow cells lying in the lowest +# lowFrontFraction of the flow elevation range +lowFrontFraction = 0.05 +# penalties of the least-cost path: upSlopePenalty applies to the positive elevation gain [m] +# of an edge, flowDistPenalty to the distance [m] of a cell from the flow footprint +upSlopePenalty = 10. +flowDistPenalty = 5. + # split point finding # first fit a parabola on the non extended path. Start and end point match the profile # the 3rd constraint is given by: @@ -66,7 +84,8 @@ tSteps = 0:5 # resType (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, particles). pta|FT|FM is the minimum for generating the # path from fields, pta|particles is the minimum for generating the path from particles (see pathFromPart parameter) -resType = pta|FT|FM +# pft is additionally required for the bottom extension to the deposit front (extBottomOption = 1) +resType = pta|FT|FM|pft # list of simulations that shall be performed (null, ent, res, entres, available (use all available input data)) simTypeList = null diff --git a/avaframe/tests/test_DFAPathGeneration.py b/avaframe/tests/test_DFAPathGeneration.py index c47b2452a..4460dd1be 100644 --- a/avaframe/tests/test_DFAPathGeneration.py +++ b/avaframe/tests/test_DFAPathGeneration.py @@ -7,6 +7,7 @@ # Local imports import avaframe.ana5Utils.DFAPathGeneration as DFAPathGeneration import avaframe.in3Utils.geoTrans as gT +import avaframe.in3Utils.fileHandlerUtils as fU def test_appendAverageStd(): @@ -136,6 +137,131 @@ def test_extendDFAPath(): assert avaProfileExt['z'][-1] == pytest.approx(0, abs=1e-6) +def test_findFlowFront(): + """""" + # plane sloping towards increasing column index, flat runout zone from column 5 on + demRaster = np.tile(np.array([50., 40., 30., 20., 10., 0., 0., 0., 0., 0., 0.]), (10, 1)) + fieldPFT = np.zeros((10, 11)) + # flow tongue along row 5, thicker towards the front + fieldPFT[5, 1:9] = np.array([0.5, 0.5, 0.5, 0.5, 1., 2., 3., 5.]) + frontRow, frontCol = DFAPathGeneration.findFlowFront(fieldPFT, demRaster, 0.01, 0.05) + # front band is the flat zone (columns 5-8), ft-weighted centroid at column 78/11 + assert frontRow == 5 + assert frontCol == 7 + + # flat deposit: the band covers the whole footprint, ft-weighted centroid at column 83/13 + frontRow, frontCol = DFAPathGeneration.findFlowFront(fieldPFT, np.ones((10, 11)), 0.01, 0.05) + assert (frontRow, frontCol) == (5, 6) + + # two disjoint lobes at the same elevation: the centroid falls between them and is + # snapped to the nearest cell of the front band + twoLobes = np.zeros((10, 11)) + twoLobes[5, 1] = 1. + twoLobes[5, 9] = 1. + frontRow, frontCol = DFAPathGeneration.findFlowFront(twoLobes, np.ones((10, 11)), 0.01, 0.05) + assert (frontRow, frontCol) in [(5, 1), (5, 9)] + + # no flow above threshold + frontRow, frontCol = DFAPathGeneration.findFlowFront(np.zeros((10, 11)), demRaster, 0.01, 0.05) + assert frontRow is None + assert frontCol is None + + +def test_leastCostPath(): + """""" + # flat dem with a flow channel along the top row and the two outer columns; + # the penalty on the distance (in meters) from the flow keeps the path inside + # the channel instead of cutting straight through the no-flow area + csz = 5. + demRaster = np.zeros((7, 7)) + fieldPFT = np.zeros((7, 7)) + fieldPFT[0, :] = 1. + fieldPFT[:, 0] = 1. + fieldPFT[:, 6] = 1. + cellPath = DFAPathGeneration.leastCostPath((3, 0), (3, 6), fieldPFT, demRaster, csz, 0.01, 10., 1.) + assert cellPath[0] == (3, 0) + assert cellPath[-1] == (3, 6) + assert all(fieldPFT[row, col] > 0 for row, col in cellPath) + + # a goal behind a nodata barrier is unreachable + demBarrier = np.zeros((7, 7)) + demBarrier[:, 3] = np.nan + cellPath = DFAPathGeneration.leastCostPath((3, 0), (3, 6), fieldPFT, demBarrier, csz, 0.01, 10., 1.) + assert cellPath == [] + + +def test_extendProfileToFront(): + """""" + # setup required inputs + cfg = configparser.ConfigParser() + cfg['PATH'] = {'nCellsResample': '1', 'extTopOption': '0', 'extBottomOption': '1', + 'nCellsMinExtend': '1', 'nCellsMaxExtend': '20', 'factBottomExt': 0.2, + 'maxIterationExtBot': 10, 'nBottomExtPrecision': 10, + 'ftThreshold': 0.01, 'lowFrontFraction': 0.05, + 'upSlopePenalty': 10., 'flowDistPenalty': 5.} + + dem = {'header': {'xllcenter': 0, 'yllcenter': 0, 'cellsize': 2, 'nrows': 10, 'ncols': 11}, + 'rasterData': np.tile(np.array([50., 40., 30., 20., 10., 0., 0., 0., 0., 0., 0.]), (10, 1))} + # flow tongue along the profile, reaching the flat runout zone (front cell (5, 7)) + fieldPFT = np.zeros((10, 11)) + fieldPFT[5, 1:9] = np.array([0.5, 0.5, 0.5, 0.5, 1., 2., 3., 5.]) + + avaProfile = {'x': np.array([2, 4, 5, 6]), 'y': np.array([10, 10, 10, 10]), + 'z': np.array([40, 30, 25, 20])} + particlesIni = {'x': np.array([1., 0.9]), 'y': np.array([10., 10.])} + particlesIni, _ = gT.projectOnRaster(dem, particlesIni, interp='bilinear') + + avaProfileExt = DFAPathGeneration.extendDFAPath(cfg['PATH'], avaProfile, dem, particlesIni, + fieldPFT=fieldPFT) + # the extension descends along the tongue and ends on the deposit front + assert avaProfileExt['x'][-1] == pytest.approx(14., abs=1e-6) + assert avaProfileExt['y'][-1] == pytest.approx(10., abs=1e-6) + assert avaProfileExt['z'][-1] == pytest.approx(0., abs=1e-6) + assert np.all(np.diff(avaProfileExt['s']) > 0) + + # if the path already ends on the front cell, the straight-line extension takes over so + # that the profile is always extended at the bottom (resamplePath relies on it) + avaProfileEnd = {'x': np.array([8, 10, 12, 14]), 'y': np.array([10, 10, 10, 10]), + 'z': np.array([10., 0., 0., 0.])} + avaProfileExt = DFAPathGeneration.extendDFAPath(cfg['PATH'], avaProfileEnd, dem, particlesIni, + fieldPFT=fieldPFT) + assert avaProfileExt['x'][-1] > 14. + assert np.all(np.isfinite(avaProfileExt['z'])) + assert np.all(np.diff(avaProfileExt['s']) > 0) + + # without a peak flow thickness field, option 1 falls back to the straight-line extension + avaProfileNoField = {'x': np.array([2, 4, 5, 6]), 'y': np.array([10, 10, 10, 10]), + 'z': np.array([40, 30, 25, 20])} + avaProfileFallback = DFAPathGeneration.extendDFAPath(cfg['PATH'], avaProfileNoField, dem, + particlesIni) + cfg['PATH']['extBottomOption'] = '0' + avaProfileOpt0 = {'x': np.array([2, 4, 5, 6]), 'y': np.array([10, 10, 10, 10]), + 'z': np.array([40, 30, 25, 20])} + avaProfileOpt0 = DFAPathGeneration.extendDFAPath(cfg['PATH'], avaProfileOpt0, dem, particlesIni, + fieldPFT=fieldPFT) + assert np.allclose(avaProfileFallback['x'], avaProfileOpt0['x']) + assert np.allclose(avaProfileFallback['y'], avaProfileOpt0['y']) + + +def test_readPeakFT(tmp_path): + """""" + peakDir = tmp_path / 'Outputs' / 'com1DFA' / 'peakFiles' + peakDir.mkdir(parents=True) + content = ('ncols 3\nnrows 2\nxllcenter 0.\nyllcenter 0.\ncellsize 5.\nNODATA_value -9999\n' + '1. 2. 3.\n4. 5. 6.\n') + (peakDir / 'relA_0123456789_C_M_null_dfa_pft.asc').write_text(content) + # the peak files are parsed once and the dataframe is passed to readPeakFT + peakFilesDF = fU.makeSimDF(peakDir, avaDir=tmp_path) + # the simulation is found by its hash (the index of the configuration dataframe) + fieldPFT = DFAPathGeneration.readPeakFT(peakFilesDF, '0123456789') + assert fieldPFT.shape == (2, 3) + # and by its full simulation name + fieldPFT = DFAPathGeneration.readPeakFT(peakFilesDF, 'relA_0123456789_C_M_null_dfa') + assert fieldPFT.shape == (2, 3) + # no pft available for the simulation: returns None + assert DFAPathGeneration.readPeakFT(peakFilesDF, 'someOtherSim') is None + + def test_resamplePath(): """""" # setup required inputs diff --git a/docs/moduleAna5Utils.rst b/docs/moduleAna5Utils.rst index 464ff823c..697ccfb08 100644 --- a/docs/moduleAna5Utils.rst +++ b/docs/moduleAna5Utils.rst @@ -119,10 +119,25 @@ There are two options available to extend the mass-averaged path profile in the distance between a point in the release and the first point of the mass-averaged path profile. -We also extend the path at the bottom, to have some buffer in the runout area. This is done by finding the direction of -the path given by the last few points within the path in the x,y domain (all points at a distance ``nCellsMinExtend`` * -cellSize < distance < ``nCellsMaxExtend`` * cellSize)) and extending in this direction by a given factor -(``factBottomExt``) of the total length of the path :math:`s`. +We also extend the path at the bottom, to have some buffer in the runout area. Two options exist +(``extBottomOption``): + +* ``extBottomOption = 0`` (default): find the direction of + the path given by the last few points within the path in the x,y domain (all points at a distance ``nCellsMinExtend`` * + cellSize < distance < ``nCellsMaxExtend`` * cellSize)) and extend in this direction by a given factor + (``factBottomExt``) of the total length of the path :math:`s`. + +* ``extBottomOption = 1``: extend the path to the front of the deposit. The mass-averaged path ends + at the last center of mass position, which lies behind the deposit front, so this option first + locates the front as the flow-thickness-weighted centroid of the flow cells lying in the lowest + ``lowFrontFraction`` of the flow elevation range (flow cells are those above ``ftThreshold`` in the + peak flow thickness field). A smaller ``lowFrontFraction`` places the front closer to the furthest + reaching deposit cells; the default is 0.05. The end of the mass-averaged path is then connected to + the front with a least-cost path (Dijkstra) over the DEM. The cost of a step is its horizontal + length, penalized by the positive elevation gain (``upSlopePenalty``) and by the distance from the + flow footprint (``flowDistPenalty``), so the extension descends along the deposit. With this option + the path ends at the simulated deposit front instead of at a fixed relative distance. If the front + cannot be located or reached, the path falls back to the straight-line extension of option 0. Resampling ==========