From 4d19b1bbe6354a2549fe9d69360461c55c6a8096 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Latorre=20V=2E?= Date: Thu, 11 Jun 2026 07:25:26 -0400 Subject: [PATCH 1/3] feat(ana5Utils): add extBottomOption to extend the DFA path to the deposit front The mass-averaged path ends short of the deposit once the front decelerates, and the existing bottom extension extrapolates a straight line of fixed relative length (factBottomExt * sMax) clipped at the DEM border, independently of where the avalanche actually stopped. This adds extBottomOption = 1: the front of the deposit is located in the peak flow thickness field (flow-thickness-weighted centroid of the flow cells in the lowest part of the flow elevation range, with a lowest-cell fallback for flat deposits) and the path is extended to it along a Dijkstra least-cost path over the DEM, penalizing uphill steps and cells away from the flow footprint. The distance from the flow is computed in meters (distance_transform_edt with sampling=cellsize) so the penalties are resolution independent. Default behaviour is unchanged (extBottomOption = 0); the option is dispatched in extendDFAPath, which gains an optional fieldFT argument, and the pft field is read in generatePathAndSplitpoint only when the option is active. Used for the corridor-scale thalweg extraction on Route 115-CH (Chile), validated there against 208 manually delineated reference thalwegs (median endpoint distance 64 m). --- avaframe/ana5Utils/DFAPathGeneration.py | 259 +++++++++++++++++++- avaframe/ana5Utils/DFAPathGenerationCfg.ini | 23 +- avaframe/tests/test_DFAPathGeneration.py | 123 ++++++++++ docs/moduleAna5Utils.rst | 22 +- 4 files changed, 416 insertions(+), 11 deletions(-) diff --git a/avaframe/ana5Utils/DFAPathGeneration.py b/avaframe/ana5Utils/DFAPathGeneration.py index 625ac29cb..dcf555d8e 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 @@ -84,6 +87,10 @@ def generatePathAndSplitpoint(avalancheDir, cfgDFAPath, cfgMain, runDFAModule): 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 + fieldFT = None + if cfgDFAPath['PATH'].getint('extBottomOption', fallback=0) == 1: + fieldFT = readPeakFT(avalancheDir, simName) # get the mass average path avaProfileMass, particlesIni = generateMassAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInfo=cfgDFAPath['PATH'].getboolean('addVelocityInfo')) @@ -93,7 +100,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, + fieldFT=fieldFT) # resample path and keep track of start and end of mass averaged part avaProfileMass = resamplePath(cfgDFAPath['PATH'], dem, avaProfileMass) # get split point @@ -306,7 +314,37 @@ def getMassAvgPathFromFields(fieldsList, fieldHeader, dem): return avaProfileMass -def extendDFAPath(cfg, avaProfile, dem, particlesIni): +def readPeakFT(avalancheDir, simName, comModule='com1DFA'): + """ read the peak flow thickness field (pft) of one simulation + + Parameters + ----------- + avalancheDir: str or pathlib path + avalanche directory + simName: str + simulation name + comModule: str + computational module name (subdirectory of Outputs) + + Returns + -------- + fieldFT: numpy array + peak flow thickness raster (same grid as the simulation dem), + None if no pft peak field is available for this simulation + """ + inputDir = pathlib.Path(avalancheDir, 'Outputs', comModule, 'peakFiles') + peakFilesDF = fU.makeSimDF(inputDir, avaDir=avalancheDir) + # 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 in %s' % (simName, inputDir)) + return None + return IOf.readRaster(peakFilesDF.loc[index[0], 'files'])['rasterData'] + + +def extendDFAPath(cfg, avaProfile, dem, particlesIni, fieldFT=None): """ extend the DFA path at the top and bottom avaProfile with x, y, z, s information @@ -316,6 +354,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 fixed-length extrapolation + method, 1 for the least-cost extension to the deposit front (requires fieldFT) - 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 +369,8 @@ def extendDFAPath(cfg, avaProfile, dem, particlesIni): dem dict particlesIni: dict initial particles dict + fieldFT: numpy array, optional + peak flow thickness field on the dem grid, only used if extBottomOption = 1 Returns -------- @@ -339,7 +381,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 fieldFT is None: + log.warning('extBottomOption is 1 but no flow thickness field was provided, ' + 'falling back to the fixed-length bottom extension') + avaProfile = extendProfileBottom(cfg, dem, avaProfile) + else: + avaProfile = extendProfileToFront(cfg, dem, avaProfile, fieldFT) + else: + avaProfile = extendProfileBottom(cfg, dem, avaProfile) return avaProfile @@ -523,6 +573,203 @@ def extendProfileBottom(cfg, dem, profile): return profile +def extendProfileToFront(cfg, dem, profile, fieldFT): + """ extend the DFA path at the bottom to the front of the deposit + + Locate the front of the deposit in the 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 fixed-length 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 + fieldFT: numpy array + flow thickness field (peak or last time step) 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 fieldFT.shape != zRaster.shape: + message = '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(fieldFT, zRaster, cfg.getfloat('ftThreshold'), + cfg.getfloat('lowFrontFraction')) + if frontRow is None: + log.warning('No flow cell above ftThreshold was found, ' + 'falling back to the fixed-length 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), fieldFT, 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 fixed-length 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(fieldFT, demRaster, ftThreshold, lowFrontFraction): + """ locate the front of the deposit in a 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 + ----------- + fieldFT: numpy array + 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 = (fieldFT > 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 = fieldFT[rows, cols] * lowBand + 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, fieldFT, 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 + fieldFT: numpy array + 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(~(fieldFT > 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 +899,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..ad9df8781 100644 --- a/avaframe/tests/test_DFAPathGeneration.py +++ b/avaframe/tests/test_DFAPathGeneration.py @@ -136,6 +136,129 @@ 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)) + fieldFT = np.zeros((10, 11)) + # flow tongue along row 5, thicker towards the front + fieldFT[5, 1:9] = np.array([0.5, 0.5, 0.5, 0.5, 1., 2., 3., 5.]) + frontRow, frontCol = DFAPathGeneration.findFlowFront(fieldFT, 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(fieldFT, 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)) + fieldFT = np.zeros((7, 7)) + fieldFT[0, :] = 1. + fieldFT[:, 0] = 1. + fieldFT[:, 6] = 1. + cellPath = DFAPathGeneration.leastCostPath((3, 0), (3, 6), fieldFT, demRaster, csz, 0.01, 10., 1.) + assert cellPath[0] == (3, 0) + assert cellPath[-1] == (3, 6) + assert all(fieldFT[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), fieldFT, 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)) + fieldFT = np.zeros((10, 11)) + fieldFT[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, + fieldFT=fieldFT) + # 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 fixed-length 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, + fieldFT=fieldFT) + assert avaProfileExt['x'][-1] > 14. + assert np.all(np.isfinite(avaProfileExt['z'])) + assert np.all(np.diff(avaProfileExt['s']) > 0) + + # without a flow thickness field, option 1 falls back to the fixed-length 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, + fieldFT=fieldFT) + 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 simulation is found by its hash (the index of the configuration dataframe) + fieldFT = DFAPathGeneration.readPeakFT(tmp_path, '0123456789') + assert fieldFT.shape == (2, 3) + # and by its full simulation name + fieldFT = DFAPathGeneration.readPeakFT(tmp_path, 'relA_0123456789_C_M_null_dfa') + assert fieldFT.shape == (2, 3) + # no pft available for the simulation: returns None + assert DFAPathGeneration.readPeakFT(tmp_path, 'someOtherSim') is None + + def test_resamplePath(): """""" # setup required inputs diff --git a/docs/moduleAna5Utils.rst b/docs/moduleAna5Utils.rst index 464ff823c..93d7b7068 100644 --- a/docs/moduleAna5Utils.rst +++ b/docs/moduleAna5Utils.rst @@ -119,10 +119,24 @@ 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. Since the mass-averaged path + ends short of the deposit once the front decelerates, 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), + and then connects the end of the mass-averaged path 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 fixed-length extension of option 0. Resampling ========== From 3b56dc1e57e1d4f64827ce31c348b86214857eb7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Latorre=20V=2E?= Date: Sat, 20 Jun 2026 09:53:18 -0400 Subject: [PATCH 2/3] perf(ana5Utils): build peak file dataframe once for extBottomOption Move the makeSimDF call out of the per-simulation loop in generatePathAndSplitpoint: the peak files are parsed once and the dataframe is passed to readPeakFT, instead of re-reading and re-parsing all peak files for every simulation. Document why findFlowFront needs no divide-by-zero guard on the flow-weight sum. --- avaframe/ana5Utils/DFAPathGeneration.py | 28 ++++++++++++++---------- avaframe/tests/test_DFAPathGeneration.py | 9 +++++--- 2 files changed, 23 insertions(+), 14 deletions(-) diff --git a/avaframe/ana5Utils/DFAPathGeneration.py b/avaframe/ana5Utils/DFAPathGeneration.py index dcf555d8e..a047d4b96 100644 --- a/avaframe/ana5Utils/DFAPathGeneration.py +++ b/avaframe/ana5Utils/DFAPathGeneration.py @@ -83,14 +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 fieldFT = None - if cfgDFAPath['PATH'].getint('extBottomOption', fallback=0) == 1: - fieldFT = readPeakFT(avalancheDir, simName) + if extendToFront: + fieldFT = readPeakFT(peakFilesDF, simName) # get the mass average path avaProfileMass, particlesIni = generateMassAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInfo=cfgDFAPath['PATH'].getboolean('addVelocityInfo')) @@ -314,17 +322,15 @@ def getMassAvgPathFromFields(fieldsList, fieldHeader, dem): return avaProfileMass -def readPeakFT(avalancheDir, simName, comModule='com1DFA'): - """ read the peak flow thickness field (pft) of one simulation +def readPeakFT(peakFilesDF, simName): + """ get the peak flow thickness field (pft) of one simulation from the peak file dataframe Parameters ----------- - avalancheDir: str or pathlib path - avalanche directory + peakFilesDF: pandas DataFrame + peak files of all simulations (from fU.makeSimDF), parsed once for the whole run simName: str simulation name - comModule: str - computational module name (subdirectory of Outputs) Returns -------- @@ -332,14 +338,12 @@ def readPeakFT(avalancheDir, simName, comModule='com1DFA'): peak flow thickness raster (same grid as the simulation dem), None if no pft peak field is available for this simulation """ - inputDir = pathlib.Path(avalancheDir, 'Outputs', comModule, 'peakFiles') - peakFilesDF = fU.makeSimDF(inputDir, avaDir=avalancheDir) # 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 in %s' % (simName, inputDir)) + log.warning('No pft peak field found for simulation %s' % simName) return None return IOf.readRaster(peakFilesDF.loc[index[0], 'files'])['rasterData'] @@ -686,6 +690,8 @@ def findFlowFront(fieldFT, demRaster, ftThreshold, lowFrontFraction): elev = demRaster[rows, cols] lowBand = elev <= elev.min() + lowFrontFraction * (elev.max() - elev.min()) weight = fieldFT[rows, cols] * lowBand + # weight.sum() is always positive: the lowest flow cell is in lowBand and, being a flow + # cell, carries fieldFT > 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]: diff --git a/avaframe/tests/test_DFAPathGeneration.py b/avaframe/tests/test_DFAPathGeneration.py index ad9df8781..d528a3993 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(): @@ -249,14 +250,16 @@ def test_readPeakFT(tmp_path): 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) - fieldFT = DFAPathGeneration.readPeakFT(tmp_path, '0123456789') + fieldFT = DFAPathGeneration.readPeakFT(peakFilesDF, '0123456789') assert fieldFT.shape == (2, 3) # and by its full simulation name - fieldFT = DFAPathGeneration.readPeakFT(tmp_path, 'relA_0123456789_C_M_null_dfa') + fieldFT = DFAPathGeneration.readPeakFT(peakFilesDF, 'relA_0123456789_C_M_null_dfa') assert fieldFT.shape == (2, 3) # no pft available for the simulation: returns None - assert DFAPathGeneration.readPeakFT(tmp_path, 'someOtherSim') is None + assert DFAPathGeneration.readPeakFT(peakFilesDF, 'someOtherSim') is None def test_resamplePath(): From 8a5fd52d8caa568334cfc8b54457421cead2192d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Latorre=20V=2E?= Date: Mon, 22 Jun 2026 21:46:40 -0400 Subject: [PATCH 3/3] refactor(ana5Utils): name the bottom-extension field peak flow thickness The least-cost bottom extension uses the peak flow thickness field (pft), so rename fieldFT to fieldPFT and make the docstrings and messages say peak flow thickness consistently. Relabel the option-0 extension as straight-line instead of fixed-length, rename nrows/ncols to nRows/nCols in leastCostPath, and update the docs (the mass-averaged path ends at the last center of mass point; document lowFrontFraction and its default). --- avaframe/ana5Utils/DFAPathGeneration.py | 76 ++++++++++++------------ avaframe/tests/test_DFAPathGeneration.py | 44 +++++++------- docs/moduleAna5Utils.rst | 21 +++---- 3 files changed, 71 insertions(+), 70 deletions(-) diff --git a/avaframe/ana5Utils/DFAPathGeneration.py b/avaframe/ana5Utils/DFAPathGeneration.py index a047d4b96..1157190fb 100644 --- a/avaframe/ana5Utils/DFAPathGeneration.py +++ b/avaframe/ana5Utils/DFAPathGeneration.py @@ -96,9 +96,9 @@ def generatePathAndSplitpoint(avalancheDir, cfgDFAPath, cfgMain, runDFAModule): 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 - fieldFT = None + fieldPFT = None if extendToFront: - fieldFT = readPeakFT(peakFilesDF, simName) + fieldPFT = readPeakFT(peakFilesDF, simName) # get the mass average path avaProfileMass, particlesIni = generateMassAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInfo=cfgDFAPath['PATH'].getboolean('addVelocityInfo')) @@ -109,7 +109,7 @@ def generatePathAndSplitpoint(avalancheDir, cfgDFAPath, cfgMain, runDFAModule): 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, - fieldFT=fieldFT) + fieldPFT=fieldPFT) # resample path and keep track of start and end of mass averaged part avaProfileMass = resamplePath(cfgDFAPath['PATH'], dem, avaProfileMass) # get split point @@ -334,7 +334,7 @@ def readPeakFT(peakFilesDF, simName): Returns -------- - fieldFT: numpy array + fieldPFT: numpy array peak flow thickness raster (same grid as the simulation dem), None if no pft peak field is available for this simulation """ @@ -348,7 +348,7 @@ def readPeakFT(peakFilesDF, simName): return IOf.readRaster(peakFilesDF.loc[index[0], 'files'])['rasterData'] -def extendDFAPath(cfg, avaProfile, dem, particlesIni, fieldFT=None): +def extendDFAPath(cfg, avaProfile, dem, particlesIni, fieldPFT=None): """ extend the DFA path at the top and bottom avaProfile with x, y, z, s information @@ -358,8 +358,8 @@ def extendDFAPath(cfg, avaProfile, dem, particlesIni, fieldFT=None): 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 fixed-length extrapolation - method, 1 for the least-cost extension to the deposit front (requires fieldFT) + - 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 @@ -373,7 +373,7 @@ def extendDFAPath(cfg, avaProfile, dem, particlesIni, fieldFT=None): dem dict particlesIni: dict initial particles dict - fieldFT: numpy array, optional + fieldPFT: numpy array, optional peak flow thickness field on the dem grid, only used if extBottomOption = 1 Returns @@ -386,12 +386,12 @@ def extendDFAPath(cfg, avaProfile, dem, particlesIni, fieldFT=None): avaProfile, _ = gT.prepareLine(dem, avaProfile, distance=resampleDistance, Point=None) avaProfile = extendProfileTop(cfg.getint('extTopOption'), particlesIni, avaProfile) if cfg.getint('extBottomOption', fallback=0) == 1: - if fieldFT is None: - log.warning('extBottomOption is 1 but no flow thickness field was provided, ' - 'falling back to the fixed-length bottom extension') + 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, fieldFT) + avaProfile = extendProfileToFront(cfg, dem, avaProfile, fieldPFT) else: avaProfile = extendProfileBottom(cfg, dem, avaProfile) return avaProfile @@ -577,15 +577,15 @@ def extendProfileBottom(cfg, dem, profile): return profile -def extendProfileToFront(cfg, dem, profile, fieldFT): +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 flow thickness field (findFlowFront) + 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 fixed-length extension + reached, the profile falls back to the straight-line extension (extendProfileBottom). Parameters @@ -601,8 +601,8 @@ def extendProfileToFront(cfg, dem, profile, fieldFT): dem dict profile: dict profile to extend - fieldFT: numpy array - flow thickness field (peak or last time step) on the same grid as the dem + fieldPFT: numpy array + peak flow thickness field (pft) on the same grid as the dem Returns -------- @@ -612,8 +612,8 @@ def extendProfileToFront(cfg, dem, profile, fieldFT): header = dem['header'] csz = header['cellsize'] zRaster = dem['rasterData'] - if fieldFT.shape != zRaster.shape: - message = 'Flow thickness field and dem do not have the same shape' + 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 @@ -621,21 +621,21 @@ def extendProfileToFront(cfg, dem, profile, fieldFT): yLast = profile['y'][-1] sLast = profile['s'][-1] # locate the deposit front - frontRow, frontCol = findFlowFront(fieldFT, zRaster, cfg.getfloat('ftThreshold'), + 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 fixed-length bottom extension') + '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), fieldFT, zRaster, csz, + 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 fixed-length bottom extension') + '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 @@ -657,8 +657,8 @@ def extendProfileToFront(cfg, dem, profile, fieldFT): return profile -def findFlowFront(fieldFT, demRaster, ftThreshold, lowFrontFraction): - """ locate the front of the deposit in a flow thickness field +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 @@ -669,8 +669,8 @@ def findFlowFront(fieldFT, demRaster, ftThreshold, lowFrontFraction): Parameters ----------- - fieldFT: numpy array - flow thickness field + fieldPFT: numpy array + peak flow thickness field demRaster: numpy array dem raster of the same shape ftThreshold: float @@ -683,15 +683,15 @@ def findFlowFront(fieldFT, demRaster, ftThreshold, lowFrontFraction): frontRow, frontCol: int cell of the front, (None, None) if there is no flow above ftThreshold """ - flow = (fieldFT > ftThreshold) & np.isfinite(demRaster) + 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 = fieldFT[rows, cols] * lowBand + weight = fieldPFT[rows, cols] * lowBand # weight.sum() is always positive: the lowest flow cell is in lowBand and, being a flow - # cell, carries fieldFT > ftThreshold >= 0, so no divide-by-zero guard is needed + # 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]: @@ -702,7 +702,7 @@ def findFlowFront(fieldFT, demRaster, ftThreshold, lowFrontFraction): return frontRow, frontCol -def leastCostPath(startCell, goalCell, fieldFT, demRaster, csz, ftThreshold, upSlopePenalty, +def leastCostPath(startCell, goalCell, fieldPFT, demRaster, csz, ftThreshold, upSlopePenalty, flowDistPenalty): """ Dijkstra least-cost path between two cells of the dem grid @@ -714,8 +714,8 @@ def leastCostPath(startCell, goalCell, fieldFT, demRaster, csz, ftThreshold, upS ----------- startCell, goalCell: tuple (row, col) of the start and goal cells - fieldFT: numpy array - flow thickness field used to build the flow footprint + fieldPFT: numpy array + peak flow thickness field used to build the flow footprint demRaster: numpy array dem raster of the same shape csz: float @@ -737,14 +737,14 @@ def leastCostPath(startCell, goalCell, fieldFT, demRaster, csz, ftThreshold, upS log.error(message) raise AssertionError(message) # distance (m) from the flow footprint, used to keep the path on the deposit - distToFlow = distance_transform_edt(~(fieldFT > ftThreshold), sampling=csz) - nrows, ncols = demRaster.shape + 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 = np.full((nRows, nCols), np.inf) costGrid[startCell] = 0. - previousCell = np.full((nrows, ncols, 2), -1, dtype=np.int64) + 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) @@ -754,7 +754,7 @@ def leastCostPath(startCell, goalCell, fieldFT, demRaster, csz, ftThreshold, upS continue for dRow, dCol in neighbours: nRow, nCol = row + dRow, col + dCol - if 0 <= nRow < nrows and 0 <= nCol < ncols and demValid[nRow, nCol]: + 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 diff --git a/avaframe/tests/test_DFAPathGeneration.py b/avaframe/tests/test_DFAPathGeneration.py index d528a3993..4460dd1be 100644 --- a/avaframe/tests/test_DFAPathGeneration.py +++ b/avaframe/tests/test_DFAPathGeneration.py @@ -141,16 +141,16 @@ 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)) - fieldFT = np.zeros((10, 11)) + fieldPFT = np.zeros((10, 11)) # flow tongue along row 5, thicker towards the front - fieldFT[5, 1:9] = np.array([0.5, 0.5, 0.5, 0.5, 1., 2., 3., 5.]) - frontRow, frontCol = DFAPathGeneration.findFlowFront(fieldFT, demRaster, 0.01, 0.05) + 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(fieldFT, np.ones((10, 11)), 0.01, 0.05) + 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 @@ -174,19 +174,19 @@ def test_leastCostPath(): # the channel instead of cutting straight through the no-flow area csz = 5. demRaster = np.zeros((7, 7)) - fieldFT = np.zeros((7, 7)) - fieldFT[0, :] = 1. - fieldFT[:, 0] = 1. - fieldFT[:, 6] = 1. - cellPath = DFAPathGeneration.leastCostPath((3, 0), (3, 6), fieldFT, demRaster, csz, 0.01, 10., 1.) + 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(fieldFT[row, col] > 0 for row, col in cellPath) + 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), fieldFT, demBarrier, csz, 0.01, 10., 1.) + cellPath = DFAPathGeneration.leastCostPath((3, 0), (3, 6), fieldPFT, demBarrier, csz, 0.01, 10., 1.) assert cellPath == [] @@ -203,8 +203,8 @@ def test_extendProfileToFront(): 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)) - fieldFT = np.zeros((10, 11)) - fieldFT[5, 1:9] = np.array([0.5, 0.5, 0.5, 0.5, 1., 2., 3., 5.]) + 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])} @@ -212,24 +212,24 @@ def test_extendProfileToFront(): particlesIni, _ = gT.projectOnRaster(dem, particlesIni, interp='bilinear') avaProfileExt = DFAPathGeneration.extendDFAPath(cfg['PATH'], avaProfile, dem, particlesIni, - fieldFT=fieldFT) + 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 fixed-length extension takes over so + # 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, - fieldFT=fieldFT) + fieldPFT=fieldPFT) assert avaProfileExt['x'][-1] > 14. assert np.all(np.isfinite(avaProfileExt['z'])) assert np.all(np.diff(avaProfileExt['s']) > 0) - # without a flow thickness field, option 1 falls back to the fixed-length extension + # 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, @@ -238,7 +238,7 @@ def test_extendProfileToFront(): 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, - fieldFT=fieldFT) + fieldPFT=fieldPFT) assert np.allclose(avaProfileFallback['x'], avaProfileOpt0['x']) assert np.allclose(avaProfileFallback['y'], avaProfileOpt0['y']) @@ -253,11 +253,11 @@ def test_readPeakFT(tmp_path): # 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) - fieldFT = DFAPathGeneration.readPeakFT(peakFilesDF, '0123456789') - assert fieldFT.shape == (2, 3) + fieldPFT = DFAPathGeneration.readPeakFT(peakFilesDF, '0123456789') + assert fieldPFT.shape == (2, 3) # and by its full simulation name - fieldFT = DFAPathGeneration.readPeakFT(peakFilesDF, 'relA_0123456789_C_M_null_dfa') - assert fieldFT.shape == (2, 3) + 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 diff --git a/docs/moduleAna5Utils.rst b/docs/moduleAna5Utils.rst index 93d7b7068..697ccfb08 100644 --- a/docs/moduleAna5Utils.rst +++ b/docs/moduleAna5Utils.rst @@ -127,16 +127,17 @@ We also extend the path at the bottom, to have some buffer in the runout area. T 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. Since the mass-averaged path - ends short of the deposit once the front decelerates, 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), - and then connects the end of the mass-averaged path 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 fixed-length extension of option 0. +* ``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 ==========