From 4a931c4f1dc3a5dd2635dc9e85814b5bfb492865 Mon Sep 17 00:00:00 2001 From: PaulaSp3 Date: Fri, 29 Aug 2025 13:18:09 +0200 Subject: [PATCH] extend thalweg to top in flow direction resample path comments FSO correct comment adapt test comments AW adapt dou fix typo Co-authored-by: Anna Wirbel <68687423+awirb@users.noreply.github.com> move extension with direction into separate function can compute extension for 2 points use parameter from ini file adapt tests; config ini --- avaframe/ana5Utils/DFAPathGeneration.py | 663 ++++++++++++-------- avaframe/ana5Utils/DFAPathGenerationCfg.ini | 3 + avaframe/ana5Utils/preparePathGeneral.py | 140 +++++ avaframe/com1DFA/DFAtools.py | 41 ++ avaframe/out3Plot/outCom1DFA.py | 10 +- avaframe/out3Plot/outCom3Plots.py | 370 +++++++---- avaframe/out3Plot/plotUtils.py | 1 + avaframe/tests/test_DFAPathGeneration.py | 237 ++++++- avaframe/tests/test_preparePathGeneral.py | 169 +++++ docs/moduleAna5Utils.rst | 43 ++ 10 files changed, 1303 insertions(+), 374 deletions(-) create mode 100644 avaframe/ana5Utils/preparePathGeneral.py create mode 100644 avaframe/tests/test_preparePathGeneral.py diff --git a/avaframe/ana5Utils/DFAPathGeneration.py b/avaframe/ana5Utils/DFAPathGeneration.py index 1157190fb..ad2e74d89 100644 --- a/avaframe/ana5Utils/DFAPathGeneration.py +++ b/avaframe/ana5Utils/DFAPathGeneration.py @@ -1,5 +1,5 @@ """ - Tools for generating an avalanche path from a DFA simulation +Tools for generating an avalanche path from a DFA simulation """ # Load modules @@ -30,7 +30,8 @@ # change log level in calling module to DEBUG to see log messages log = logging.getLogger(__name__) cfgAVA = cfgUtils.getGeneralConfig() -debugPlot = cfgAVA['FLAGS'].getboolean('debugPlot') +debugPlot = cfgAVA["FLAGS"].getboolean("debugPlot") + def generatePathAndSplitpoint(avalancheDir, cfgDFAPath, cfgMain, runDFAModule): """ @@ -50,29 +51,33 @@ def generatePathAndSplitpoint(avalancheDir, cfgDFAPath, cfgMain, runDFAModule): splitPoint: pathlib file path to the split point result saved as a shapefile """ - if runDFAModule: # call DFA module to perform simulations with overrides from DFAPath config - # Clean avalanche directory of old work and output files from module - initProj.cleanModuleFiles(avalancheDir, com1DFA, deleteOutput=True) - # create and read the default com1DFA config (no local is read) - com1DFACfg = cfgUtils.getModuleConfig(com1DFA, avalancheDir, toPrint=False, - onlyDefault=cfgDFAPath['com1DFA_com1DFA_override'].getboolean( - 'defaultConfig')) - # and override with settings from DFAPath config - com1DFACfg, cfgDFAPath = cfgHandling.applyCfgOverride(com1DFACfg, cfgDFAPath, com1DFA, - addModValues=False) - outDir = pathlib.Path(avalancheDir, 'Outputs', 'ana5Utils', 'DFAPath') - fU.makeADir(outDir) - # write configuration to file for documentation - com1DFACfgFile = outDir / 'com1DFAPathGenerationCfg.ini' - with open(com1DFACfgFile, 'w') as configfile: - com1DFACfg.write(configfile) - # call com1DFA and perform simulations - dem, plotDict, reportDictList, simDF = com1DFA.com1DFAMain(cfgMain, cfgInfo=com1DFACfg) - else: # read existing simulation results + if runDFAModule: # call DFA module to perform simulations with overrides from DFAPath config + # Clean avalanche directory of old work and output files from module + initProj.cleanModuleFiles(avalancheDir, com1DFA, deleteOutput=True) + # create and read the default com1DFA config (no local is read) + com1DFACfg = cfgUtils.getModuleConfig( + com1DFA, + avalancheDir, + toPrint=False, + onlyDefault=cfgDFAPath["com1DFA_com1DFA_override"].getboolean("defaultConfig"), + ) + # and override with settings from DFAPath config + com1DFACfg, cfgDFAPath = cfgHandling.applyCfgOverride( + com1DFACfg, cfgDFAPath, com1DFA, addModValues=False + ) + outDir = pathlib.Path(avalancheDir, "Outputs", "ana5Utils", "DFAPath") + fU.makeADir(outDir) + # write configuration to file for documentation + com1DFACfgFile = outDir / "com1DFAPathGenerationCfg.ini" + with open(com1DFACfgFile, "w") as configfile: + com1DFACfg.write(configfile) + # call com1DFA and perform simulations + dem, plotDict, reportDictList, simDF = com1DFA.com1DFAMain(cfgMain, cfgInfo=com1DFACfg) + else: # read existing simulation results # read simulation dem demOri = gI.readDEM(avalancheDir) dem = com1DFA.setDEMoriginToZero(demOri) - dem['originalHeader'] = demOri['header'].copy() + dem["originalHeader"] = demOri["header"].copy() # load DFA results (use runCom1DFA to generate these results for example) # here is an example with com1DFA but another DFA computational module can be used # as long as it produces some pta, particles or FT, FM and FV results @@ -100,31 +105,40 @@ def generatePathAndSplitpoint(avalancheDir, cfgDFAPath, cfgMain, runDFAModule): if extendToFront: fieldPFT = readPeakFT(peakFilesDF, simName) # get the mass average path - avaProfileMass, particlesIni = generateMassAveragePath(avalancheDir, pathFromPart, simName, dem, - addVelocityInfo=cfgDFAPath['PATH'].getboolean('addVelocityInfo')) - avaProfileMass, _ = gT.prepareLine(dem, avaProfileMass, distance=resampleDistance, Point=None) - avaProfileMass['indStartMassAverage'] = 1 - avaProfileMass['indEndMassAverage'] = np.size(avaProfileMass['x']) + avaProfileMass, particlesIni = generateMassAveragePath( + avalancheDir, + pathFromPart, + simName, + dem, + addVelocityInfo=cfgDFAPath["PATH"].getboolean("addVelocityInfo"), + ) + kResample = cfgDFAPath.getint("kResample") + avaProfileMass, _ = gT.prepareLine(dem, avaProfileMass, distance=resampleDistance, Point=None, k=kResample) + avaProfileMass["indStartMassAverage"] = 1 + avaProfileMass["indEndMassAverage"] = np.size(avaProfileMass["x"]) # make the parabolic fit - parabolicFit = getParabolicFit(cfgDFAPath['PATH'], avaProfileMass, dem) + 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, fieldPFT=fieldPFT) # resample path and keep track of start and end of mass averaged part - avaProfileMass = resamplePath(cfgDFAPath['PATH'], dem, avaProfileMass) + avaProfileMass = resamplePath(cfgDFAPath["PATH"], dem, avaProfileMass) # get split point - splitPoint = getSplitPoint(cfgDFAPath['PATH'], avaProfileMass, parabolicFit) + splitPoint = getSplitPoint(cfgDFAPath["PATH"], avaProfileMass, parabolicFit) # make analysis and generate plots - _ = outCom3Plots.generateCom1DFAPathPlot(avalancheDir, cfgDFAPath['PATH'], avaProfileMass, dem, - parabolicFit, splitPoint, simName) + _ = outCom3Plots.generateCom1DFAPathPlot( + avalancheDir, cfgDFAPath["PATH"], avaProfileMass, dem, parabolicFit, splitPoint, simName + ) # now save the path and split point as shapefiles - avaPath,splitPoint = saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem) + avaPath, splitPoint = saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem) return avaPath, splitPoint -def generateMassAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInfo=False, flagAvaDir=True, - comModule='com1DFA'): - """ extract path from fields or particles + +def generateMassAveragePath( + avalancheDir, pathFromPart, simName, dem, addVelocityInfo=False, flagAvaDir=True, comModule="com1DFA" +): + """extract path from fields or particles Parameters ----------- @@ -156,31 +170,33 @@ def generateMassAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocit x, y coord of the initial particles or flow thickness field """ if pathFromPart: - particlesList, timeStepInfo = particleTools.readPartFromPickle(avalancheDir, simName=simName, flagAvaDir=True, - comModule='com1DFA') + particlesList, timeStepInfo = particleTools.readPartFromPickle( + avalancheDir, simName=simName, flagAvaDir=True, comModule="com1DFA" + ) particlesIni = particlesList[0] - log.info('Using particles to generate avalanche path profile') + log.info("Using particles to generate avalanche path profile") # postprocess to extract path and energy line avaProfileMass = getMassAvgPathFromPart(particlesList, addVelocityInfo=addVelocityInfo) else: - particlesList = '' + particlesList = "" # read field - fieldName = ['FT', 'FM'] + fieldName = ["FT", "FM"] if addVelocityInfo: - fieldName.append('FV') - fieldsList, fieldHeader, timeList = com1DFA.readFields(avalancheDir, fieldName, simName=simName, - flagAvaDir=True, comModule='com1DFA') + fieldName.append("FV") + fieldsList, fieldHeader, timeList = com1DFA.readFields( + avalancheDir, fieldName, simName=simName, flagAvaDir=True, comModule="com1DFA" + ) # get fields header - ncols = fieldHeader['ncols'] - nrows = fieldHeader['nrows'] - csz = fieldHeader['cellsize'] + ncols = fieldHeader["ncols"] + nrows = fieldHeader["nrows"] + csz = fieldHeader["cellsize"] # we want the origin to be in (0, 0) as it is in the avaProfile that comes in X, Y = gT.makeCoordinateGrid(0, 0, csz, ncols, nrows) - indNonZero = np.where(fieldsList[0]['FT'] > 0) + indNonZero = np.where(fieldsList[0]["FT"] > 0) # convert this data in a particles style (dict with x, y, z info) - particlesIni = {'x': X[indNonZero], 'y': Y[indNonZero]} + particlesIni = {"x": X[indNonZero], "y": Y[indNonZero]} particlesIni, _ = gT.projectOnRaster(dem, particlesIni) - log.info('Using fields to generate avalanche path profile') + log.info("Using fields to generate avalanche path profile") # postprocess to extract path and energy line avaProfileMass = getMassAvgPathFromFields(fieldsList, fieldHeader, dem) @@ -188,7 +204,7 @@ def generateMassAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocit def getMassAvgPathFromPart(particlesList, addVelocityInfo=False): - """ compute mass averaged path from particles + """compute mass averaged path from particles Also returns the averaged velocity and kinetic energy associated If addVelocityInfo is True, information about velocity and kinetic energy is computed @@ -208,46 +224,46 @@ def getMassAvgPathFromPart(particlesList, addVelocityInfo=False): the avaProfileMass dict (u2, ekin, totEKin) """ - propList = ['x', 'y', 'z', 's', 'sCor'] - propListPart = ['x', 'y', 'z', 'trajectoryLengthXY', 'trajectoryLengthXYCor'] + propList = ["x", "y", "z", "s", "sCor"] + propListPart = ["x", "y", "z", "trajectoryLengthXY", "trajectoryLengthXYCor"] avaProfileMass = {} # do we have velocity info? if addVelocityInfo: - propList.append('u2') - propList.append('ekin') - propListPart.append('u2') - propListPart.append('ekin') - avaProfileMass['totEKin'] = np.empty((0, 1)) + propList.append("u2") + propList.append("ekin") + propListPart.append("u2") + propListPart.append("ekin") + avaProfileMass["totEKin"] = np.empty((0, 1)) # initialize other properties for prop in propList: avaProfileMass[prop] = np.empty((0, 1)) - avaProfileMass[prop + 'std'] = np.empty((0, 1)) + avaProfileMass[prop + "std"] = np.empty((0, 1)) # loop on each particle dictionary (ie each time step saved) for particles in particlesList: - if particles['nPart'] > 0: - m = particles['m'] + if particles["nPart"] > 0: + m = particles["m"] if addVelocityInfo: - ux = particles['ux'] - uy = particles['uy'] - uz = particles['uz'] + ux = particles["ux"] + uy = particles["uy"] + uz = particles["uz"] u = DFAtls.norm(ux, uy, uz) - u2Array = u*u - kineticEneArray = 0.5*m*u2Array - particles['u2'] = u2Array - particles['ekin'] = kineticEneArray + u2Array = u * u + kineticEneArray = 0.5 * m * u2Array + particles["u2"] = u2Array + particles["ekin"] = kineticEneArray # mass-averaged path avaProfileMass = appendAverageStd(propList, avaProfileMass, particles, m, naming=propListPart) if addVelocityInfo: - avaProfileMass['totEKin'] = np.append(avaProfileMass['totEKin'], np.nansum(kineticEneArray)) + avaProfileMass["totEKin"] = np.append(avaProfileMass["totEKin"], np.nansum(kineticEneArray)) return avaProfileMass def getMassAvgPathFromFields(fieldsList, fieldHeader, dem): - """ compute mass averaged path from fields + """compute mass averaged path from fields Also returns the averaged velocity and kinetic energy associated The dem and fieldsList (FT, FM and FV) need to have identical dimensions and cell size. @@ -270,52 +286,52 @@ def getMassAvgPathFromFields(fieldsList, fieldHeader, dem): the avaProfileMass dict (u2, ekin, totEKin) """ # get DEM - demRaster = dem['rasterData'] + demRaster = dem["rasterData"] # get fields header - ncols = fieldHeader['ncols'] - nrows = fieldHeader['nrows'] - xllc = fieldHeader['xllcenter'] - yllc = fieldHeader['yllcenter'] - csz = fieldHeader['cellsize'] + ncols = fieldHeader["ncols"] + nrows = fieldHeader["nrows"] + xllc = fieldHeader["xllcenter"] + yllc = fieldHeader["yllcenter"] + csz = fieldHeader["cellsize"] X, Y = gT.makeCoordinateGrid(xllc, yllc, csz, ncols, nrows) - propList = ['x', 'y', 'z'] + propList = ["x", "y", "z"] avaProfileMass = {} # do we have velocity info? addVelocityInfo = False - if 'FV' in fieldsList[0]: - propList.append('u2') - propList.append('ekin') - avaProfileMass['totEKin'] = np.empty((0, 1)) + if "FV" in fieldsList[0]: + propList.append("u2") + propList.append("ekin") + avaProfileMass["totEKin"] = np.empty((0, 1)) addVelocityInfo = True # initialize other properties for prop in propList: avaProfileMass[prop] = np.empty((0, 1)) - avaProfileMass[prop + 'std'] = np.empty((0, 1)) + avaProfileMass[prop + "std"] = np.empty((0, 1)) # loop on each field dictionary (ie each time step saved) for field in fieldsList: # find cells with snow - nonZeroIndex = np.where(field['FT'] > 0) + nonZeroIndex = np.where(field["FT"] > 0) xArray = X[nonZeroIndex] yArray = Y[nonZeroIndex] zArray, _ = gT.projectOnGrid(xArray, yArray, demRaster, csz=csz, xllc=xllc, yllc=yllc) - mArray = field['FM'][nonZeroIndex] - particles = {'x': xArray, 'y': yArray, 'z': zArray} + mArray = field["FM"][nonZeroIndex] + particles = {"x": xArray, "y": yArray, "z": zArray} if addVelocityInfo: - uArray = field['FV'][nonZeroIndex] - u2Array = uArray*uArray - kineticEneArray = 0.5*mArray*u2Array - particles['u2'] = u2Array - particles['ekin'] = kineticEneArray + uArray = field["FV"][nonZeroIndex] + u2Array = uArray * uArray + kineticEneArray = 0.5 * mArray * u2Array + particles["u2"] = u2Array + particles["ekin"] = kineticEneArray # mass-averaged path avaProfileMass = appendAverageStd(propList, avaProfileMass, particles, mArray) if addVelocityInfo: - avaProfileMass['totEKin'] = np.append(avaProfileMass['totEKin'], np.nansum(kineticEneArray)) + avaProfileMass["totEKin"] = np.append(avaProfileMass["totEKin"], np.nansum(kineticEneArray)) - avaProfileMass['x'] = avaProfileMass['x'] - xllc - avaProfileMass['y'] = avaProfileMass['y'] - yllc + avaProfileMass["x"] = avaProfileMass["x"] - xllc + avaProfileMass["y"] = avaProfileMass["y"] - yllc # compute s avaProfileMass = gT.computeS(avaProfileMass) @@ -382,9 +398,10 @@ def extendDFAPath(cfg, avaProfile, dem, particlesIni, fieldPFT=None): extended profile at top and bottom (x, y, z). """ # resample the profile - resampleDistance = cfg.getfloat('nCellsResample') * dem['header']['cellsize'] - avaProfile, _ = gT.prepareLine(dem, avaProfile, distance=resampleDistance, Point=None) - avaProfile = extendProfileTop(cfg.getint('extTopOption'), particlesIni, avaProfile) + resampleDistance = cfg.getfloat("nCellsResample") * dem["header"]["cellsize"] + kResample = cfg.getint("kResample") + avaProfile, _ = gT.prepareLine(dem, avaProfile, distance=resampleDistance, Point=None, k=kResample) + avaProfile = extendProfileTop(cfg.getint('extTopOption'), particlesIni, avaProfile, dem, cfg) if cfg.getint('extBottomOption', fallback=0) == 1: if fieldPFT is None: log.warning('extBottomOption is 1 but no peak flow thickness field was provided, ' @@ -397,22 +414,30 @@ def extendDFAPath(cfg, avaProfile, dem, particlesIni, fieldPFT=None): return avaProfile -def extendProfileTop(extTopOption, particlesIni, profile): - """ extend the DFA path at the top (release) +def extendProfileTop(extTopOption, particlesIni, profile, dem=None, cfg=None, considerLLC=False): + """extend the DFA path at the top (release) Either towards the highest point in particlesIni (extTopOption = 0) or the point leading to the longest runout (extTopOption = 1) + or in the upslope direction of the thalweg (extTopOption = 2) Parameters ----------- extTopOption: int decide how to extend towards the top if 0, extrapolate towards the highest point in the release - if 1, extrapolate towards the point leading to the lonest runout + if 1, extrapolate towards the point leading to the longest runout + if 2, extrapolate in upslope direction of the thalweg between a defined distance of the thalweg particlesIni: dict initial particles dict profile: dict profile to extend + dem: dict + contains DEM info (necessary for extTopOption = 2) + cfg: configparser.ConfigParser + settings (necessary for extTopOption = 2) + considerLLC: bool + if True the coordinate of the lower left center is considered (regarding extTopOption = 2) Returns -------- @@ -421,51 +446,80 @@ def extendProfileTop(extTopOption, particlesIni, profile): """ if extTopOption == 0: # get highest particle - indTop = np.argmax(particlesIni['z']) - xExtTop = particlesIni['x'][indTop] - yExtTop = particlesIni['y'][indTop] - zExtTop = particlesIni['z'][indTop] - dx = xExtTop - profile['x'][0] - dy = yExtTop - profile['y'][0] + indTop = np.argmax(particlesIni["z"]) + xExtTop = particlesIni["x"][indTop] + yExtTop = particlesIni["y"][indTop] + zExtTop = particlesIni["z"][indTop] + dx = xExtTop - profile["x"][0] + dy = yExtTop - profile["y"][0] ds = np.sqrt(dx**2 + dy**2) elif extTopOption == 1: # get point with the most important runout gain # get first particle of the path - xFirst = profile['x'][0] - yFirst = profile['y'][0] - zFirst = profile['z'][0] + xFirst = profile["x"][0] + yFirst = profile["y"][0] + zFirst = profile["z"][0] # get last particle of the path - sLast = profile['s'][-1] - zLast = profile['z'][-1] + sLast = profile["s"][-1] + zLast = profile["z"][-1] # compute runout angle for averaged path - tanAngle = (zFirst-zLast)/sLast + tanAngle = (zFirst - zLast) / sLast # compute ds - dx = particlesIni['x'] - xFirst - dy = particlesIni['y'] - yFirst + dx = particlesIni["x"] - xFirst + dy = particlesIni["y"] - yFirst ds = np.sqrt(dx**2 + dy**2) # compute dz - dz = particlesIni['z'] - zFirst + dz = particlesIni["z"] - zFirst # remove the elevation needed to match the runout angle dz1 = dz - tanAngle * ds # get the particle with the highest potential indTop = np.argmax(dz1) - xExtTop = particlesIni['x'][indTop] - yExtTop = particlesIni['y'][indTop] - zExtTop = particlesIni['z'][indTop] + xExtTop = particlesIni["x"][indTop] + yExtTop = particlesIni["y"][indTop] + zExtTop = particlesIni["z"][indTop] ds = ds[indTop] + elif extTopOption == 2: + if len(profile["x"]) <= 1: + # skip computation if thalweg is one point + log.warning("Skip top extension of thalweg since profile contains only one point.") + return profile + + if dem is None: + message = "If extTopOption = 2, the dem needs to be provided" + log.error(message) + raise ValueError(message) + if cfg is None: + message = f"If extTopOption = 2, the cfg needs to be provided" + log.error(message) + raise ValueError(message) + extProfile, _ = extendProfileDirection("top", cfg, dem, profile, considerLLC=considerLLC) + + if "x" in extProfile: + xExtTop = extProfile["x"] + yExtTop = extProfile["y"] + zExtTop = extProfile["z"] + dx = xExtTop - profile["x"][0] + dy = yExtTop - profile["y"][0] + ds = np.sqrt(dx ** 2 + dy ** 2) + + else: + message = f"The extend top option {extTopOption} is not valid, please change!" + log.error(message) + raise ValueError(message) # extend profile - profile['x'] = np.append(xExtTop, profile['x']) - profile['y'] = np.append(yExtTop, profile['y']) - profile['z'] = np.append(zExtTop, profile['z']) - profile['s'] = np.append(0, profile['s'] + ds) + if xExtTop is not None: + profile["x"] = np.append(xExtTop, profile["x"]) + profile["y"] = np.append(yExtTop, profile["y"]) + profile["z"] = np.append(zExtTop, profile["z"]) + profile["s"] = np.append(0, profile["s"] + ds) if debugPlot: debPlot.plotPathExtTop(profile, particlesIni, xFirst, yFirst, zFirst, dz1) return profile -def extendProfileBottom(cfg, dem, profile): - """ extend the DFA path at the bottom (runout area) +def extendProfileBottom(cfg, dem, profile, considerLLC=False): + """extend the DFA path at the bottom (runout area) Find the direction in which to extend considering the last point of the profile and a few previous ones but discarding the ones that are too close @@ -485,69 +539,158 @@ def extendProfileBottom(cfg, dem, profile): dem dict profile: dict profile to extend + considerLLC: bool + If True, the lower left center coordinates are considered when reading z coordinates Returns -------- profile: dict extended profile """ - header = dem['header'] - csz = header['cellsize'] - zRaster = dem['rasterData'] - # get last point - xLast = profile['x'][-1] - yLast = profile['y'][-1] - sLast = profile['s'][-1] - # compute distance from last point: - r = DFAtls.norm(profile['x']-xLast, profile['y']-yLast, 0) + + extProfile, _ = extendProfileDirection("bottom", cfg, dem, profile, considerLLC=considerLLC) + + if "x" in extProfile: + xLast = profile["x"][-1] + yLast = profile["y"][-1] + sLast = profile["s"][-1] + # extend profile + profile["x"] = np.append(profile["x"], extProfile["x"]) + profile["y"] = np.append(profile["y"], extProfile["y"]) + profile["z"] = np.append(profile["z"], extProfile["z"]) + profile["s"] = np.append( + profile["s"], sLast + np.sqrt((xLast - extProfile["x"]) ** 2 + (yLast - extProfile["y"]) ** 2) + ) + + if debugPlot: + _, interestProfile = extendProfileDirection("bottom", cfg, dem, profile, considerLLC=considerLLC) + xInterest = interestProfile["x"] + yInterest = interestProfile["y"] + debPlot.plotPathExtBot(profile, xInterest, yInterest, 0 * yInterest, xLast, yLast) + return profile + + +def extendProfileDirection(direction, cfg, dem, profile, considerLLC, ): + """extend profile to upslope/ downslope direction of thalweg + + Find the direction in which to extend considering the first/last point of the profile + and a few previous ones but discarding the ones that are too close + (nCellsMinExtend* csz < distFromLast <= nCellsMaxExtend * csz). + Extend in this direction for a distance factBottomExt * length of the path (to top or to bottom). + + Parameters + ----------- + cfg: configParser + nCellsMinExtend: int, when extending towards the bottom, take points + at more than nCellsMinExtend*demCellSize from first/last point to get the direction + nCellsMaxExtend: int, when extending towards the bottom, take points at + less than nCellsMaxExtend*demCellSize from first/last point to get the direction + factBottomExt: float, extend the profile from factBottomExt*sMax + dem: dict + DEM header and rasterData + profile: dict + profile to extend + considerLLC: bool + If True, the lower left center coordinates are considered when reading z coordinates + + Returns + -------- + profile: dict + contains the point for extension + """ + + if direction.lower() == "top": + indexInterest = int(0) + elif direction.lower() == "bottom": + indexInterest = int(-1) + else: + message = "direction must be either 'top' or 'bottom'" + log.error(message) + raise ValueError(message) + + header = dem["header"] + csz = header["cellsize"] + if considerLLC: + # compute center coordinates of lower left cell + xllcenter = header["xllcenter"] + yllcenter = header["yllcenter"] + else: + xllcenter = 0 + yllcenter = 0 + zRaster = dem["rasterData"] + # get first/last point + xInt = profile["x"][indexInterest] + yInt = profile["y"][indexInterest] + sLast = profile["s"][-1] + # compute distance from first/last point: + r = DFAtls.norm(profile["x"] - xInt, profile["y"] - yInt, 0) # find the previous points - extendMinDistance = cfg.getfloat('nCellsMinExtend') * csz - extendMaxDistance = cfg.getfloat('nCellsMaxExtend') * csz - pointsOfInterestLast = np.where((r < extendMaxDistance) & (r > extendMinDistance))[0] - xInterest = profile['x'][pointsOfInterestLast] - yInterest = profile['y'][pointsOfInterestLast] + extendMinDistance = cfg.getfloat("nCellsMinExtend") * csz + extendMaxDistance = cfg.getfloat("nCellsMaxExtend") * csz + pointsOfInterest = np.where((r < extendMaxDistance) & (r > extendMinDistance))[0] # check if points are found to compute direction of extension - if len(xInterest) > 0: + if pointsOfInterest.size > 0: + xInterest = profile["x"][pointsOfInterest] + yInterest = profile["y"][pointsOfInterest] # find the direction in which we need to extend the path - vDirX = xLast - xInterest - vDirY = yLast - yInterest - vDirX, vDirY, vDirZ = DFAtls.normalize(np.array([vDirX]), np.array([vDirY]), 0*np.array([vDirY])) - vDirX = np.sum(vDirX) - vDirY = np.sum(vDirY) - vDirZ = np.sum(vDirZ) - vDirX, vDirY, vDirZ = DFAtls.normalize(np.array([vDirX]), np.array([vDirY]), np.array([vDirZ])) + if direction.lower() == "top": + vDirX = xInterest - xInt + vDirY = yInterest - yInt + elif direction.lower() == "bottom": + vDirX = xInt - xInterest + vDirY = yInt - yInterest + + vDirX, vDirY, vDirZ = DFAtls.getAveragedDirection(vDirX, vDirY) # extend in this direction - factExt = cfg.getfloat('factBottomExt') - gamma = factExt * sLast / np.sqrt(vDirX**2 + vDirY**2) - xExtBottom = np.array([xLast + gamma * vDirX]) - yExtBottom = np.array([yLast + gamma * vDirY]) + factExt = cfg.getfloat("factBottomExt") + vNorm = DFAtls.norm(vDirX, vDirY, vDirX * 0) + gamma = factExt * sLast / vNorm + if direction.lower() == "top": + xExt = np.array([xInt - gamma * vDirX]) + yExt = np.array([yInt - gamma * vDirY]) + elif direction.lower() == "bottom": + xExt = np.array([xInt + gamma * vDirX]) + yExt = np.array([yInt + gamma * vDirY]) # project on DEM - zExtBottom, _ = gT.projectOnGrid(xExtBottom, yExtBottom, zRaster, csz=csz) - # Dicothomie method to find the last point on the extention and on the dem - if np.isnan(zExtBottom): - factExt = factExt/2 + zExt, _ = gT.projectOnGrid( + xExt, yExt, zRaster, csz=csz, xllc=xllcenter, yllc=yllcenter + ) + # Dicothomie method to find the first/last point on the extension and on the dem + if np.isnan(zExt): + factExt = factExt / 2 stepSize = factExt isOut = True else: isOut = False stepSize = 0 count = 0 - # remember last point found inside + # remember first/last point found inside factLast = 0 - while count < cfg.getint('maxIterationExtBot') and stepSize * sLast > cfg.getint('nBottomExtPrecision')*csz: + while ( + count < cfg.getint("maxIterationExtBot") + and stepSize * sLast > cfg.getint("nBottomExtPrecision") * csz + ): count = count + 1 - gamma = factExt * sLast / np.sqrt(vDirX**2 + vDirY**2) - xExtBottom = np.array([xLast + gamma * vDirX]) - yExtBottom = np.array([yLast + gamma * vDirY]) + gamma = factExt * sLast / vNorm + + if direction.lower() == "top": + xExt = np.array([xInt - gamma * vDirX]) + yExt = np.array([yInt - gamma * vDirY]) + elif direction.lower() == "bottom": + xExt = np.array([xInt + gamma * vDirX]) + yExt = np.array([yInt + gamma * vDirY]) + # project on DEM - zExtBottom, _ = gT.projectOnGrid(xExtBottom, yExtBottom, zRaster, csz=csz) - stepSize = stepSize/2 - if np.isnan(zExtBottom): + zExt, _ = gT.projectOnGrid( + xExt, yExt, zRaster, csz=csz, xllc=xllcenter, yllc=yllcenter + ) + stepSize = stepSize / 2 + if np.isnan(zExt): factExt = factExt - stepSize isOut = True else: - # remember last point found inside + # remember first/last point found inside factLast = factExt factExt = factExt + stepSize isOut = False @@ -555,26 +698,31 @@ def extendProfileBottom(cfg, dem, profile): if isOut: # the last iteration is not in the domain, fall back to last point in domain factExt = factLast - gamma = factExt * sLast / np.sqrt(vDirX**2 + vDirY**2) - xExtBottom = np.array([xLast + gamma * vDirX]) - yExtBottom = np.array([yLast + gamma * vDirY]) + gamma = factExt * sLast / np.sqrt(vDirX ** 2 + vDirY ** 2) + + if direction.lower() == "top": + xExt = np.array([xInt - gamma * vDirX]) + yExt = np.array([yInt - gamma * vDirY]) + elif direction.lower() == "bottom": + xExt = np.array([xInt + gamma * vDirX]) + yExt = np.array([yInt + gamma * vDirY]) + # project on DEM - zExtBottom, _ = gT.projectOnGrid(xExtBottom, yExtBottom, zRaster, csz=csz) - log.info('found extention after %d iterations, precision is %.2f m' % (count, stepSize * sLast)) + zExt, _ = gT.projectOnGrid( + xExt, yExt, zRaster, csz=csz, xllc=xllcenter, yllc=yllcenter + ) - # 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'], sLast + np.sqrt((xLast-xExtBottom)**2 + (yLast-yExtBottom)**2)) + log.info("found extension after %d iterations, precision is %.2f m" % (count, stepSize * sLast)) + extProfile = {"x": xExt, "y": yExt, "z": zExt} + intProfile = {"x": xInterest, "y": yInterest} + return extProfile, intProfile else: - log.warning('Path not extended at bottom as no point of interest for computing direction \ - of where to extend path is found') - - if debugPlot: - debPlot.plotPathExtBot(profile, xInterest, yInterest, 0*yInterest, xLast, yLast) - return profile + log.warning( + "Path not extended at bottom as no point of interest for computing direction \ + of where to extend path is found" + ) + return {}, {} def extendProfileToFront(cfg, dem, profile, fieldPFT): @@ -799,34 +947,34 @@ def getParabolicFit(cfg, avaProfile, dem): parabolicFit: dict a, b, c coefficients of the parabolic fit (y = a*a*x + b*x + c) """ - s = avaProfile['s'] + s = avaProfile["s"] sE = s[-1] - z = avaProfile['z'] + z = avaProfile["z"] z0 = z[0] zE = z[-1] # same start and end point, minimize distance between curves - if cfg.getfloat('fitOption') == 0: - SumNom = np.sum(s*(s-sE)*((zE-z0)/sE*s+z0-z)) - SumDenom = s*(s-sE) + if cfg.getfloat("fitOption") == 0: + SumNom = np.sum(s * (s - sE) * ((zE - z0) / sE * s + z0 - z)) + SumDenom = s * (s - sE) SumDenom = np.dot(SumDenom, SumDenom) - a = - SumNom/SumDenom - b = (zE-z0)/sE - a*sE - elif cfg.getfloat('fitOption') == 1: + a = -SumNom / SumDenom + b = (zE - z0) / sE - a * sE + elif cfg.getfloat("fitOption") == 1: angleProf, tmpProf, dsProf = gT.prepareAngleProfile(10, avaProfile, raiseWarning=False) - r = avaProfile['s'] - avaProfile['s'][-1] - resampleDistance = cfg.getfloat('nCellsSlope') * dem['header']['cellsize'] - pointsOfInterestLast = np.where(np.abs(r) < resampleDistance) - slope = np.nansum(angleProf[pointsOfInterestLast])/np.size(pointsOfInterestLast) + r = avaProfile["s"] - avaProfile["s"][-1] + resampleDistance = cfg.getfloat("nCellsSlope") * dem["header"]["cellsize"] + pointsOfInterest = np.where(np.abs(r) < resampleDistance) + slope = np.nansum(angleProf[pointsOfInterest]) / np.size(pointsOfInterest) slope = -np.tan(np.radians(slope)) - a = (slope*sE + (z0 - zE))/(sE*sE) - b = -slope - 2*(z0 - zE)/sE + a = (slope * sE + (z0 - zE)) / (sE * sE) + b = -slope - 2 * (z0 - zE) / sE c = z0 - parabolicFit = {'a': a, 'b': b, 'c': c} + parabolicFit = {"a": a, "b": b, "c": c} return parabolicFit def getSplitPoint(cfg, avaProfile, parabolicFit): - """ find the split point corresponding to an avalanche profile, with parabolic fit and the slopeSplitPoint + """find the split point corresponding to an avalanche profile, with parabolic fit and the slopeSplitPoint Parameters ----------- @@ -847,32 +995,48 @@ def getSplitPoint(cfg, avaProfile, parabolicFit): splitPoint: dict (x, y, z, zPra, s) at split point location. """ - indFirst = avaProfile['indStartMassAverage'] - indEnd = avaProfile['indEndMassAverage'] - s0 = avaProfile['s'][indFirst] - sEnd = avaProfile['s'][indEnd] - s = avaProfile['s'] - z = avaProfile['z'] + indFirst = avaProfile["indStartMassAverage"] + indEnd = avaProfile["indEndMassAverage"] + s0 = avaProfile["s"][indFirst] + sEnd = avaProfile["s"][indEnd] + s = avaProfile["s"] + z = avaProfile["z"] sNew = s - s0 - zPara = parabolicFit['a']*sNew*sNew+parabolicFit['b']*sNew+parabolicFit['c'] - parabolicProfile = {'s': sNew, 'z': zPara} + zPara = parabolicFit["a"] * sNew * sNew + parabolicFit["b"] * sNew + parabolicFit["c"] + parabolicProfile = {"s": sNew, "z": zPara} - anglePara, tmpPara, dsPara = gT.prepareAngleProfile(cfg.getfloat('slopeSplitPoint'), parabolicProfile, - raiseWarning=False) + anglePara, tmpPara, dsPara = gT.prepareAngleProfile( + cfg.getfloat("slopeSplitPoint"), parabolicProfile, raiseWarning=False + ) try: - indSplitPoint = gT.findAngleProfile(tmpPara, dsPara, cfg.getfloat('dsMin')) - splitPoint = {'x': avaProfile['x'][indSplitPoint], 'y': avaProfile['y'][indSplitPoint], - 'z': z[indSplitPoint], 'zPara': zPara[indSplitPoint], 's': sNew[indSplitPoint]} + indSplitPoint = gT.findAngleProfile(tmpPara, dsPara, cfg.getfloat("dsMin")) + splitPoint = { + "x": avaProfile["x"][indSplitPoint], + "y": avaProfile["y"][indSplitPoint], + "z": z[indSplitPoint], + "zPara": zPara[indSplitPoint], + "s": sNew[indSplitPoint], + } except IndexError: - noSplitPointFoundMessage = ('Automated split point generation failed as no point where slope is less than %s°' - 'was found, setting split point at the top. Correct split point manually.' - % cfg.getfloat('slopeSplitPoint')) - splitPoint = {'x': avaProfile['x'][0], 'y': avaProfile['y'][0], - 'z': z[0], 'zPara': zPara[0], 's': sNew[0], 'isTopSplitPoint': True} + noSplitPointFoundMessage = ( + "Automated split point generation failed as no point where slope is less than %s°" + "was found, setting split point at the top. Correct split point manually." + % cfg.getfloat("slopeSplitPoint") + ) + splitPoint = { + "x": avaProfile["x"][0], + "y": avaProfile["y"][0], + "z": z[0], + "zPara": zPara[0], + "s": sNew[0], + "isTopSplitPoint": True, + } log.warning(noSplitPointFoundMessage) if debugPlot: - angleProf, tmpProf, dsProf = gT.prepareAngleProfile(cfg.getfloat('slopeSplitPoint'), avaProfile) - debPlot.plotFindAngle(avaProfile, angleProf, parabolicProfile, anglePara, s0, sEnd, splitPoint, indSplitPoint) + angleProf, tmpProf, dsProf = gT.prepareAngleProfile(cfg.getfloat("slopeSplitPoint"), avaProfile) + debPlot.plotFindAngle( + avaProfile, angleProf, parabolicProfile, anglePara, s0, sEnd, splitPoint, indSplitPoint + ) return splitPoint @@ -897,20 +1061,22 @@ def resamplePath(cfg, dem, avaProfile): avaProfile: dict resampled path profile """ - resampleDistance = cfg.getfloat('nCellsResample') * dem['header']['cellsize'] - indFirst = avaProfile['indStartMassAverage'] - indEnd = avaProfile['indEndMassAverage'] - s0 = avaProfile['s'][indFirst] - sEnd = avaProfile['s'][indEnd] - avaProfile, _ = gT.prepareLine(dem, avaProfile, distance=resampleDistance, Point=None) + resampleDistance = cfg.getfloat("nCellsResample") * dem["header"]["cellsize"] + kResample = cfg.getint("kResample") + indFirst = avaProfile["indStartMassAverage"] + indEnd = avaProfile["indEndMassAverage"] + s0 = avaProfile["s"][indFirst] + sEnd = avaProfile["s"][indEnd] + avaProfile, _ = gT.prepareLine(dem, avaProfile, distance=resampleDistance, Point=None, k=kResample) # 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; 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 + indEndCandidates = np.argwhere(avaProfile["s"] >= sEnd + resampleDistance / 3) + indEnd = indEndCandidates[0][0] - 1 if len(indEndCandidates) > 0 else np.size(avaProfile["s"]) - 1 avaProfile['indEndMassAverage'] = indEnd + + avaProfile['indStartMassAverage'] = indFirst return avaProfile @@ -938,33 +1104,38 @@ def saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem): file path to the saved shapefile for the split point """ # put path back in original location - if splitPoint != '': - splitPoint['x'] = splitPoint['x'] + dem['originalHeader']['xllcenter'] - splitPoint['y'] = splitPoint['y'] + dem['originalHeader']['yllcenter'] - avaProfileMass['x'] = avaProfileMass['x'] + dem['originalHeader']['xllcenter'] - avaProfileMass['y'] = avaProfileMass['y'] + dem['originalHeader']['yllcenter'] + if splitPoint != "": + splitPoint["x"] = splitPoint["x"] + dem["originalHeader"]["xllcenter"] + splitPoint["y"] = splitPoint["y"] + dem["originalHeader"]["yllcenter"] + avaProfileMass["x"] = avaProfileMass["x"] + dem["originalHeader"]["xllcenter"] + avaProfileMass["y"] = avaProfileMass["y"] + dem["originalHeader"]["yllcenter"] # get projection from release shp layer - simName = simDFrow['simName'] + simName = simDFrow["simName"] relName = cfgUtils.parseSimName(simName)["releaseName"] - inProjection = pathlib.Path(avalancheDir, 'Inputs', 'REL', relName + '.prj') + inProjection = pathlib.Path(avalancheDir, "Inputs", "REL", relName + ".prj") # save profile in Inputs - pathAB = pathlib.Path(avalancheDir, 'Outputs', 'ana5Utils', 'DFAPath', 'massAvgPath_%s_AB_aimec' % simName) - name = 'massAvaPath' + pathAB = pathlib.Path( + avalancheDir, "Outputs", "ana5Utils", "DFAPath", "massAvgPath_%s_AB_aimec" % simName + ) + name = "massAvaPath" shpConv.writeLine2SHPfile(avaProfileMass, name, pathAB) if inProjection.is_file(): - shutil.copy(inProjection, pathAB.with_suffix('.prj')) + shutil.copy(inProjection, pathAB.with_suffix(".prj")) else: - message = ('No projection layer for shp file %s' % inProjection) + message = "No projection layer for shp file %s" % inProjection log.warning(message) - log.info('Saved path to: %s', pathAB) - if splitPoint != '': - splitAB = pathlib.Path(avalancheDir, 'Outputs', 'ana5Utils', 'DFAPath', 'splitPointParabolicFit_%s_AB_aimec' % simName) - name = 'parabolaSplitPoint' + log.info("Saved path to: %s", pathAB) + if splitPoint != "": + splitAB = pathlib.Path( + avalancheDir, "Outputs", "ana5Utils", "DFAPath", "splitPointParabolicFit_%s_AB_aimec" % simName + ) + name = "parabolaSplitPoint" shpConv.writePoint2SHPfile(splitPoint, name, splitAB) if inProjection.is_file(): - shutil.copy(inProjection, splitAB.with_suffix('.prj')) - log.info('Saved split point to: %s', splitAB) - return pathAB,splitAB + shutil.copy(inProjection, splitAB.with_suffix(".prj")) + log.info("Saved split point to: %s", splitAB) + return pathAB, splitAB + def weightedAvgAndStd(values, weights): """ @@ -974,12 +1145,12 @@ def weightedAvgAndStd(values, weights): """ average = np.average(values, weights=weights) # Fast and numerically precise: - variance = np.average((values-average)**2, weights=weights) + variance = np.average((values - average) ** 2, weights=weights) return (average, math.sqrt(variance)) -def appendAverageStd(propList, avaProfile, particles, weights, naming=''): - """ append averaged to path +def appendAverageStd(propList, avaProfile, particles, weights, naming=""): + """append averaged to path Parameters ----------- @@ -999,9 +1170,9 @@ def appendAverageStd(propList, avaProfile, particles, weights, naming=''): avaProfile: dict averaged profile """ - propListNames = naming if naming != '' else propList + propListNames = naming if naming != "" else propList for prop, propName in zip(propList, propListNames): avg, std = weightedAvgAndStd(particles[propName], weights) avaProfile[prop] = np.append(avaProfile[prop], avg) - avaProfile[prop + 'std'] = np.append(avaProfile[prop + 'std'], std) + avaProfile[prop + "std"] = np.append(avaProfile[prop + "std"], std) return avaProfile diff --git a/avaframe/ana5Utils/DFAPathGenerationCfg.ini b/avaframe/ana5Utils/DFAPathGenerationCfg.ini index c2fb8fbe7..ceace0e27 100644 --- a/avaframe/ana5Utils/DFAPathGenerationCfg.ini +++ b/avaframe/ana5Utils/DFAPathGenerationCfg.ini @@ -15,10 +15,13 @@ addVelocityInfo = False # the path extracted from the DFA simulation is re-sampled # re-sampling step size is defined resampleDistance = nCellsResample x cellSize) nCellsResample = 10 +# Degree of the spline for splprep. +kResample = 3 # extension method at the top # option 0: take the highest particle in the release # option 1: find the point that will lead to the longest runout +# option 2: extend the path in the direction of its first points upslope by factBottomExt x sMax extTopOption = 1 # extension method at the bottom diff --git a/avaframe/ana5Utils/preparePathGeneral.py b/avaframe/ana5Utils/preparePathGeneral.py new file mode 100644 index 000000000..51e259675 --- /dev/null +++ b/avaframe/ana5Utils/preparePathGeneral.py @@ -0,0 +1,140 @@ +""" +generate thalweg from x and y coordinates (including extension to top and bottom and resampling) +""" + +import numpy as np +import logging +import copy +# local imports +import avaframe.in3Utils.geoTrans as gT +from avaframe.ana5Utils import DFAPathGeneration +# create local logger +log = logging.getLogger(__name__) + + +def preparePathGeneralMain(profile, cfgDFAPath, dem): + """ + prepare thalweg from x and y coordinates: + 1. read z coordinates from DEM and compute horizontally projected distance + 2. extend path to bottom and top + 3. resample path points + + Parameters + ------------- + profile: dict + contains x and y coordinates of thalweg location + cfgDFAPath: configparser object + configuration for DFA path generation + dem: dict + dictionary with header and raster data of elevation model + + Returns + ------------- + profileAveraged: dict + s and z coordinates are added (x and y original) to input profile + profileExtended: dict + x, y, s, z of extended and resampled path + """ + # get profile with normalized x and y coordinates and z and s values + profileAveraged = updateSZProfile(profile, dem) + profileExtended = copy.deepcopy(profileAveraged) + + # skip profile that only contains one or two points + if len(profileAveraged["x"]) <= 2: + profileExtended["indStartMassAverage"] = 0 + profileExtended["indEndMassAverage"] = max(len(profileExtended["x"]) - 1, 0) + return profileAveraged, profileExtended + + # if extTopOption == 2, particlesIni are not used!! + profileExtended = pathExtension(profileExtended, dem, cfgDFAPath) + profileExtended = updateSZProfile(profileExtended, dem) + # resample profile/ path and save in an extra dictionary + profileExtended = DFAPathGeneration.resamplePath(cfgDFAPath["PATH"], dem, profileExtended) + + # add input parameters to extended profile if they exist + for inputPara in ["alpha", "exponent", "zDeltaMax"]: + if inputPara in profile.keys(): + profileExtended[inputPara] = profile[inputPara] + + return profileAveraged, profileExtended + + +def pathExtension(profile, demDict, cfgPathGen): + """ + extend thalweg to top and bottom of path + + Parameters + ------------ + profile: dict + thalweg data + demDict: dict + DEM data + cfgPathGen: confiparser object + configuration setup for DFA Path generation + + Returns + ------------- + profile: dict + thalweg data that are extended to top and bottom + """ + + profile["indStartMassAverage"] = 1 + # do not use the last two points because the last points are weird sometimes + profile["indEndMassAverage"] = np.size(profile["x"]) - 2 + + if cfgPathGen["PATH"].getint("extTopOption") != 2: + # TODO: if we provide particlesIni, the other options would also work. + message = "Up to now only top-extension option 2 works!" + log.error(message) + raise ValueError(message) + + profile = DFAPathGeneration.extendProfileTop( + cfgPathGen["PATH"].getint("extTopOption"), + {}, + profile, + dem=demDict, + cfg=cfgPathGen["PATH"], + considerLLC=True, + ) + + # extend the bottom + profile = DFAPathGeneration.extendProfileBottom(cfgPathGen["PATH"], demDict, profile, considerLLC=True) + + return profile + + +def updateSZProfile(profile, dem): + """ + for given coordinates (of the talweg) read z values + from DEM and compute distance between coordinates + + Parameters + ------------ + profile: dict + contains at least x and y coordinates + dem: dict + contains dem data + + Returns + ----------- + profile: dict + profile with added a and z values + """ + x = profile["x"] + y = profile["y"] + + demHeader = dem["header"] + + z, _ = gT.projectOnGrid( + x, + y, + dem["rasterData"], + csz=demHeader["cellsize"], + xllc=demHeader["xllcenter"], + yllc=demHeader["yllcenter"], + ) + s = np.append([0], gT.computeLengthOfLine2D(x, y)) + profile["z"] = z + profile["s"] = s + + return profile diff --git a/avaframe/com1DFA/DFAtools.py b/avaframe/com1DFA/DFAtools.py index e22176c4f..00c40c0fa 100644 --- a/avaframe/com1DFA/DFAtools.py +++ b/avaframe/com1DFA/DFAtools.py @@ -159,3 +159,44 @@ def scalProd(ux, uy, uz, vx, vy, vz): scal = ux*vx + uy*vy + uz*vz return scal + + +def getAveragedDirection(vDirX, vDirY, vDirZ=None): + """ + Compute a single averaged unit direction vector from multiple 2D direction vectors. + + Parameters + ---------- + vDirX: numpy array + x components of the direction vectors to average + vDirY: numpy array + y components of the direction vectors to average + vDirZ: numpy array, optional + z components of the direction vectors to average. If None (default), the + vectors are treated as 2D and z is set to 0. + + Returns + ------- + vDirX: float + x component of the normalized, averaged direction vector + vDirY: float + y component of the normalized, averaged direction vector + vDirZ: float + z component of the normalized, averaged direction vector (0 if input is 2D) + """ + if vDirZ is None: + vDirZ = 0 * np.array([vDirY]) + else: + vDirZ = np.array([vDirZ]) + + vDirX, vDirY, vDirZ = normalize( + np.array([vDirX]), np.array([vDirY]), vDirZ) + + # summed / averaged direction + vDirX = np.sum(vDirX) + vDirY = np.sum(vDirY) + vDirZ = np.sum(vDirZ) + # get unit vector + vDirX, vDirY, vDirZ = normalize(np.array([vDirX]), np.array([vDirY]), np.array([vDirZ])) + + return vDirX, vDirY, vDirZ diff --git a/avaframe/out3Plot/outCom1DFA.py b/avaframe/out3Plot/outCom1DFA.py index c10ca2efa..7b4b9ffd6 100644 --- a/avaframe/out3Plot/outCom1DFA.py +++ b/avaframe/out3Plot/outCom1DFA.py @@ -239,7 +239,7 @@ def addParticles2Plot(particles, ax, dem, whatS="m", whatC="h", colBarResType="" return ax, cb -def addDem2Plot(ax, dem, what="slope", extent="", origHeader=False): +def addDem2Plot(ax, dem, what="slope", extent="", origHeader=False, cmap=None): """Add dem to the background of a plot Parameters @@ -255,6 +255,8 @@ def addDem2Plot(ax, dem, what="slope", extent="", origHeader=False): optional: extent of NonUnifIm plot corresponding coordinates to dem data array at center locations origHeader: bool if True use originalHeader and not header + cmap: matplotlib colormap + colormap for DEM plot """ if origHeader: header = dem["originalHeader"] @@ -267,8 +269,10 @@ def addDem2Plot(ax, dem, what="slope", extent="", origHeader=False): csz = header["cellsize"] xArray = np.linspace(xllc, xllc + (ncols - 1) * csz, ncols) yArray = np.linspace(yllc, yllc + (nrows - 1) * csz, nrows) - cmap = pU.cmapGreys - cmap.set_bad(color="white") + + if cmap is None: + cmap = pU.cmapGreys + cmap.set_bad(color="white") if what == "slope": value = dem["Nz"] / DFAtls.norm(dem["Nx"], dem["Ny"], dem["Nz"]) diff --git a/avaframe/out3Plot/outCom3Plots.py b/avaframe/out3Plot/outCom3Plots.py index d7c18ea9d..37363222b 100644 --- a/avaframe/out3Plot/outCom3Plots.py +++ b/avaframe/out3Plot/outCom3Plots.py @@ -15,82 +15,104 @@ def hybridProfilePlot(avalancheDir, resultsHybrid): """Update profile plot with result of curent iteration""" - fig = plt.figure(figsize=(3*pU.figW, 2*pU.figH)) + fig = plt.figure(figsize=(3 * pU.figW, 2 * pU.figH)) ax = plt.subplot(111) nIter = len(resultsHybrid.keys()) i = 0 cmap, _, ticks, norm = pU.makeColorMap(pU.cmapAvaframeCont, 0, nIter, continuous=pU.contCmap) for key, dict in resultsHybrid.items(): - avaProfileMassExt = dict['path'] + avaProfileMassExt = dict["path"] avaProfileMassExt = geoTrans.computeS(avaProfileMassExt) - alpha = dict['alpha'] - sBetaPoint = dict['sBetaPoint'] + alpha = dict["alpha"] + sBetaPoint = dict["sBetaPoint"] col = cmap(norm(i)) # Plot the whole profile with beta, alpha ... points and lines - ax.plot(avaProfileMassExt['s'], avaProfileMassExt['z'], linestyle='-', color=col, label='profile (iteration %d)' % i) - ax.axvline(x=sBetaPoint, color=col, linestyle=':', linewidth=1, label='Beta point (iteration %d)' % i) - s = avaProfileMassExt['s'][[0, -1]] - z = avaProfileMassExt['z'][0] - s*np.tan(np.deg2rad(alpha)) - ax.plot(s, z, color=col, linestyle='--', label='AlphaLine (iteration %d)' % i) - i = i+1 - - titleText = r'Profiles extracted from the DFA simulations with corresponding $\alpha-\beta$ model results' + ax.plot( + avaProfileMassExt["s"], + avaProfileMassExt["z"], + linestyle="-", + color=col, + label="profile (iteration %d)" % i, + ) + ax.axvline( + x=sBetaPoint, color=col, linestyle=":", linewidth=1, label="Beta point (iteration %d)" % i + ) + s = avaProfileMassExt["s"][[0, -1]] + z = avaProfileMassExt["z"][0] - s * np.tan(np.deg2rad(alpha)) + ax.plot(s, z, color=col, linestyle="--", label="AlphaLine (iteration %d)" % i) + i = i + 1 + + titleText = ( + r"Profiles extracted from the DFA simulations with corresponding $\alpha-\beta$ model results" + ) ax.set_title(titleText) - ax.set_xlabel('projectd length s [m]') - ax.set_ylabel('Height [m]') - ax.set_aspect('equal', adjustable='box') - ax.grid(linestyle=':', color='0.9') + ax.set_xlabel("projectd length s [m]") + ax.set_ylabel("Height [m]") + ax.set_aspect("equal", adjustable="box") + ax.grid(linestyle=":", color="0.9") ax.legend(frameon=False) - title = ('com3HybProfPlot') - l = ax.legend(loc='lower left') + title = "com3HybProfPlot" + l = ax.legend(loc="lower left") l.set_zorder(40) pU.putAvaNameOnPlot(ax, avalancheDir) - path = pathlib.Path(avalancheDir, 'Outputs', 'com3Hybrid') - pU.saveAndOrPlot({'pathResult': path}, title, fig) + path = pathlib.Path(avalancheDir, "Outputs", "com3Hybrid") + pU.saveAndOrPlot({"pathResult": path}, title, fig) def hybridPathPlot(avalancheDir, dem, resultsHybrid, fields, particles, muArray): """Update path plot with result of curent iteration""" - headerOri = dem['originalHeader'] - xllcOri = headerOri['xllcenter'] - yllcOri = headerOri['yllcenter'] - fig = plt.figure(figsize=(3*pU.figW, 2*pU.figH)) + headerOri = dem["originalHeader"] + xllcOri = headerOri["xllcenter"] + yllcOri = headerOri["yllcenter"] + fig = plt.figure(figsize=(3 * pU.figW, 2 * pU.figH)) ax = plt.subplot(111) nIter = len(resultsHybrid.keys()) i = 0 cmap, _, ticks, norm = pU.makeColorMap(pU.cmapAvaframeCont, 0, nIter, continuous=pU.contCmap) for key, dict in resultsHybrid.items(): mu = muArray[i] - avaProfileMassExt = dict['path'] - xAB = dict['xAB'] - yAB = dict['yAB'] + avaProfileMassExt = dict["path"] + xAB = dict["xAB"] + yAB = dict["yAB"] col = cmap(norm(i)) - ax.plot(avaProfileMassExt['x'], avaProfileMassExt['y'], color=col, - label='Center of mass path iteration %d ($\mu$ = %.2f )' % (i, mu), zorder = 20) - ax.plot(xAB - xllcOri, yAB - yllcOri, 'X', color=col, markersize=8, - label=r'com2AB $\alpha$ point iteration %d' % i, zorder = 20) - i = i+1 - titleText = 'Avalanche path for each iteration with peak travel angle field \n and particles flow thickness in final time step' + ax.plot( + avaProfileMassExt["x"], + avaProfileMassExt["y"], + color=col, + label="Center of mass path iteration %d ($\mu$ = %.2f )" % (i, mu), + zorder=20, + ) + ax.plot( + xAB - xllcOri, + yAB - yllcOri, + "X", + color=col, + markersize=8, + label=r"com2AB $\alpha$ point iteration %d" % i, + zorder=20, + ) + i = i + 1 + titleText = "Avalanche path for each iteration with peak travel angle field \n and particles flow thickness in final time step" ax.set_title(titleText) - ax.set_ylabel('x [m]') - ax.set_ylabel('y [m]') - ax, extent, cb, CS = outCom1DFA.addResult2Plot(ax, dem['header'], fields['pta'], 'pta') - cb.ax.set_ylabel(pU.cfgPlotUtils['namepta']) - ax = outCom1DFA.addDem2Plot(ax, dem, what='slope', extent=extent) - ax, cb2 = outCom1DFA.addParticles2Plot(particles, ax, dem, whatS='h') - cb2.ax.set_ylabel('particle ' + pU.cfgPlotUtils['nameFT']) + ax.set_ylabel("x [m]") + ax.set_ylabel("y [m]") + ax, extent, cb, CS = outCom1DFA.addResult2Plot(ax, dem["header"], fields["pta"], "pta") + cb.ax.set_ylabel(pU.cfgPlotUtils["namepta"]) + ax = outCom1DFA.addDem2Plot(ax, dem, what="slope", extent=extent) + ax, cb2 = outCom1DFA.addParticles2Plot(particles, ax, dem, whatS="h") + cb2.ax.set_ylabel("particle " + pU.cfgPlotUtils["nameFT"]) ax.set_ylim(extent[2:]) ax.set_xlim(extent[:2]) - title = ('com3HybRasterPlot') - l = ax.legend(loc='lower left') + title = "com3HybRasterPlot" + l = ax.legend(loc="lower left") l.set_zorder(40) pU.putAvaNameOnPlot(ax, avalancheDir) - path = pathlib.Path(avalancheDir, 'Outputs', 'com3Hybrid') - pU.saveAndOrPlot({'pathResult': path}, title, fig) + path = pathlib.Path(avalancheDir, "Outputs", "com3Hybrid") + pU.saveAndOrPlot({"pathResult": path}, title, fig) def generateCom1DFAPathPlot(avalancheDir, cfgPath, avaProfileMass, dem, parabolicFit, splitPoint, simName): - """ Make energy test analysis and plot results + """Make energy test analysis and plot results Parameters ----------- @@ -108,102 +130,210 @@ def generateCom1DFAPathPlot(avalancheDir, cfgPath, avaProfileMass, dem, paraboli simulation name """ # read field - fieldsList, fieldHeader, timeList = com1DFA.readFields(avalancheDir, ['pta'], simName=simName, - flagAvaDir=True, comModule='com1DFA') + fieldsList, fieldHeader, timeList = com1DFA.readFields( + avalancheDir, ["pta"], simName=simName, flagAvaDir=True, comModule="com1DFA" + ) # compute simulation run out angle - indStart = avaProfileMass['indStartMassAverage'] - indEnd = avaProfileMass['indEndMassAverage'] - runOutAngleRad, runOutAngleDeg = energyLineTest.getRunOutAngle(avaProfileMass, indStart=indStart, indEnd=indEnd) - s0 = avaProfileMass['s'][indStart] - avaProfileMass['s'] = avaProfileMass['s'] - s0 - z0 = avaProfileMass['z'][indStart] + indStart = avaProfileMass["indStartMassAverage"] + indEnd = avaProfileMass["indEndMassAverage"] + runOutAngleRad, runOutAngleDeg = energyLineTest.getRunOutAngle( + avaProfileMass, indStart=indStart, indEnd=indEnd + ) + s0 = avaProfileMass["s"][indStart] + avaProfileMass["s"] = avaProfileMass["s"] - s0 # get parabola - sPara = np.linspace(avaProfileMass['s'][0], avaProfileMass['s'][-1], 500) - zPara = parabolicFit['a']*sPara*sPara+parabolicFit['b']*sPara+parabolicFit['c'] - parabolicProfile = {'s': sPara, 'z': zPara} + sPara = np.linspace(avaProfileMass["s"][0], avaProfileMass["s"][-1], 500) + zPara = parabolicFit["a"] * sPara * sPara + parabolicFit["b"] * sPara + parabolicFit["c"] + parabolicProfile = {"s": sPara, "z": zPara} # get angles of profiles - anglePara, tmpPara, dsPara = geoTrans.prepareAngleProfile(cfgPath.getfloat('slopeSplitPoint'), parabolicProfile, - raiseWarning=False) - angleProf, tmpProf, dsProf = geoTrans.prepareAngleProfile(cfgPath.getfloat('slopeSplitPoint'), avaProfileMass, - raiseWarning=False) + anglePara, tmpPara, dsPara = geoTrans.prepareAngleProfile( + cfgPath.getfloat("slopeSplitPoint"), parabolicProfile, raiseWarning=False + ) + angleProf, tmpProf, dsProf = geoTrans.prepareAngleProfile( + cfgPath.getfloat("slopeSplitPoint"), avaProfileMass, raiseWarning=False + ) # Create figures and plots - fig = plt.figure(figsize=(pU.figW*2, pU.figH*1.5)) + fig = plt.figure(figsize=(pU.figW * 2, pU.figH * 1.5)) # make the top-down view plot ax1 = plt.subplot2grid((2, 2), (1, 0), colspan=1) - rowsMin, rowsMax, colsMin, colsMax = pU.constrainPlotsToData(fieldsList[-1]['pta'], 5, extentOption=True, - constrainedData=False, buffer='') - ax1, extent, cbar0, cs1 = outCom1DFA.addResult2Plot(ax1, dem['header'], fieldsList[-1]['pta'], 'pta') - cbar0.ax.set_ylabel('peak travel angle') - # add DEM hillshade with contour lines - ax1 = outCom1DFA.addDem2Plot(ax1, dem, what='hillshade', extent=extent) - # add path - ax1.plot(avaProfileMass['x'][:indStart+1], avaProfileMass['y'][:indStart+1], '-y.', zorder=20, - label='_top extension', lw=2, path_effects=[pe.Stroke(linewidth=3, foreground='b'), pe.Normal()]) - ax1.plot(avaProfileMass['x'][indEnd:], avaProfileMass['y'][indEnd:], '-y.', zorder=20, - label='_bottom extension', lw=2, path_effects=[pe.Stroke(linewidth=3, foreground='g'), pe.Normal()]) - ax1.plot(avaProfileMass['x'][indStart:indEnd+1], avaProfileMass['y'][indStart:indEnd+1], '-y.', zorder=20, - label='_Center of mass path', lw=2, path_effects=[pe.Stroke(linewidth=3, foreground='k'), pe.Normal()]) - if not splitPoint.get('isTopSplitPoint', False): #if not a top split point - ax1.plot(splitPoint['x'], splitPoint['y'], 'P', color='r', label='_Split point', zorder=20) - ax1.set_xlabel('x [m]') - ax1.set_ylabel('y [m]') - ax1.axis('equal') - ax1.set_ylim([rowsMin, rowsMax]) - ax1.set_xlim([colsMin, colsMax]) - ax1.set_title('Avalanche thalweg') - pU.putAvaNameOnPlot(ax1, avalancheDir) + ax1 = avalancheThalwegPlot(ax1, fieldsList[-1]["pta"], dem, avaProfileMass, splitPoint, avalancheDir) # plot angle of profile and parabola ax2 = plt.subplot2grid((2, 2), (1, 1), colspan=1) # add path - ax2.plot(sPara, anglePara, 'k', lw=1, label='_parabolic fit') - ax2.plot(avaProfileMass['s'][:indStart+1], angleProf[:indStart+1], 'y.', - label='_top extension', lw=2, path_effects=[pe.Stroke(linewidth=3, foreground='b'), pe.Normal()]) - ax2.plot(avaProfileMass['s'][indEnd:], angleProf[indEnd:], 'y.', - label='_bottom extension', lw=2, path_effects=[pe.Stroke(linewidth=3, foreground='g'), pe.Normal()]) - ax2.plot(avaProfileMass['s'][indStart:indEnd+1], angleProf[indStart:indEnd+1], 'y.', - label='_Center of mass path slope', lw=2, path_effects=[pe.Stroke(linewidth=3, foreground='k'), pe.Normal()]) + ax2.plot(sPara, anglePara, "k", lw=1, label="_parabolic fit") + ax2.plot( + avaProfileMass["s"][: indStart + 1], + angleProf[: indStart + 1], + "y.", + label="_top extension", + lw=2, + path_effects=[pe.Stroke(linewidth=3, foreground="b"), pe.Normal()], + ) + ax2.plot( + avaProfileMass["s"][indEnd:], + angleProf[indEnd:], + "y.", + label="_bottom extension", + lw=2, + path_effects=[pe.Stroke(linewidth=3, foreground="g"), pe.Normal()], + ) + ax2.plot( + avaProfileMass["s"][indStart : indEnd + 1], + angleProf[indStart : indEnd + 1], + "y.", + label="_Center of mass path slope", + lw=2, + path_effects=[pe.Stroke(linewidth=3, foreground="k"), pe.Normal()], + ) minY, _ = ax2.get_ylim() minX, _ = ax2.get_xlim() - if not splitPoint.get('isTopSplitPoint', False): #if not a top split point - ax2.axvline(x=splitPoint['s'], color='r', linewidth=1, linestyle='-.', label='_Split point') - ax2.text(splitPoint['s'], minY, "%.2f m" % (splitPoint['s']), color='r', ha="right", va="bottom") - ax2.axhline(y=cfgPath.getfloat('slopeSplitPoint'), color='r', linewidth=1, linestyle='-.', - label='_Split point angle (%.0f°)' % cfgPath.getfloat('slopeSplitPoint')) - ax2.text(minX, cfgPath.getfloat('slopeSplitPoint'), "%.0f°" % (cfgPath.getfloat('slopeSplitPoint')), color='r', ha="left", va="bottom") - ax2.axhline(y=10, color='0.8', linewidth=1, linestyle='-.', label='_Beta angle (10°)') - ax2.text(minX, 10, "%.0f°" % (10), color='0.8', ha="left", va="bottom") - ax2.set_xlabel('$s_{xy}$ [m]') - ax2.set_ylabel('slope angle [°]') - ax2.set_title('Avalanche thalweg profile slope') + if not splitPoint.get("isTopSplitPoint", False): # if not a top split point + ax2.axvline(x=splitPoint["s"], color="r", linewidth=1, linestyle="-.", label="_Split point") + ax2.text(splitPoint["s"], minY, "%.2f m" % (splitPoint["s"]), color="r", ha="right", va="bottom") + ax2.axhline( + y=cfgPath.getfloat("slopeSplitPoint"), + color="r", + linewidth=1, + linestyle="-.", + label="_Split point angle (%.0f°)" % cfgPath.getfloat("slopeSplitPoint"), + ) + ax2.text( + minX, + cfgPath.getfloat("slopeSplitPoint"), + "%.0f°" % (cfgPath.getfloat("slopeSplitPoint")), + color="r", + ha="left", + va="bottom", + ) + ax2.axhline(y=10, color="0.8", linewidth=1, linestyle="-.", label="_Beta angle (10°)") + ax2.text(minX, 10, "%.0f°" % (10), color="0.8", ha="left", va="bottom") + ax2.set_xlabel("$s_{xy}$ [m]") + ax2.set_ylabel("slope angle [°]") + ax2.set_title("Avalanche thalweg profile slope") # make profile plot, zoomed out ax3 = plt.subplot2grid((2, 2), (0, 0), colspan=2) # plot mass averaged center of mass - ax3.plot(avaProfileMass['s'][:indStart+1], avaProfileMass['z'][:indStart+1], '-y.', label='top extension', - lw=1, path_effects=[pe.Stroke(linewidth=3, foreground='b'), pe.Normal()]) - ax3.plot(avaProfileMass['s'][indEnd:], avaProfileMass['z'][indEnd:], '-y.', label='bottom extension', - lw=1, path_effects=[pe.Stroke(linewidth=3, foreground='g'), pe.Normal()]) - ax3.plot(avaProfileMass['s'][indStart:indEnd+1], avaProfileMass['z'][indStart:indEnd+1], '-y.', - label='Center of mass path / profile / angle', - lw=1, path_effects=[pe.Stroke(linewidth=3, foreground='k'), pe.Normal()]) - ax3.plot(sPara, zPara, '-k', label='Parabolic fit') - if not splitPoint.get('isTopSplitPoint', False): #if not a top split point - ax3.axvline(x=splitPoint['s'], color='r', linewidth=1, linestyle='-.', label='Split point') - ax3.axhline(y=splitPoint['z'], color='r', linewidth=1, linestyle='-.', label='_Split point') + ax3.plot( + avaProfileMass["s"][: indStart + 1], + avaProfileMass["z"][: indStart + 1], + "-y.", + label="top extension", + lw=1, + path_effects=[pe.Stroke(linewidth=3, foreground="b"), pe.Normal()], + ) + ax3.plot( + avaProfileMass["s"][indEnd:], + avaProfileMass["z"][indEnd:], + "-y.", + label="bottom extension", + lw=1, + path_effects=[pe.Stroke(linewidth=3, foreground="g"), pe.Normal()], + ) + ax3.plot( + avaProfileMass["s"][indStart : indEnd + 1], + avaProfileMass["z"][indStart : indEnd + 1], + "-y.", + label="Center of mass path / profile / angle", + lw=1, + path_effects=[pe.Stroke(linewidth=3, foreground="k"), pe.Normal()], + ) + ax3.plot(sPara, zPara, "-k", label="Parabolic fit") + if not splitPoint.get("isTopSplitPoint", False): # if not a top split point + ax3.axvline(x=splitPoint["s"], color="r", linewidth=1, linestyle="-.", label="Split point") + ax3.axhline(y=splitPoint["z"], color="r", linewidth=1, linestyle="-.", label="_Split point") minY, _ = ax3.get_ylim() minX, _ = ax3.get_xlim() - ax3.text(splitPoint['s'], minY, "%.2f m" % (splitPoint['s']), color='r', ha="left", va="bottom") - ax3.text(minX, splitPoint['z'], "%.2f m" % (splitPoint['z']), color='r', ha="left", va="bottom") - ax3.set_xlabel('$s_{xy}$ [m]') - ax3.set_ylabel('z [m]') + ax3.text(splitPoint["s"], minY, "%.2f m" % (splitPoint["s"]), color="r", ha="left", va="bottom") + ax3.text(minX, splitPoint["z"], "%.2f m" % (splitPoint["z"]), color="r", ha="left", va="bottom") + ax3.set_xlabel("$s_{xy}$ [m]") + ax3.set_ylabel("z [m]") ax3.legend() - ax3.set_title('Avalanche thalweg profile') + ax3.set_title("Avalanche thalweg profile") - outFileName = '_'.join([simName, 'DFAPath']) - outDir = pathlib.Path(avalancheDir, 'Outputs', 'ana5Utils', 'DFAPath') + outFileName = "_".join([simName, "DFAPath"]) + outDir = pathlib.Path(avalancheDir, "Outputs", "ana5Utils", "DFAPath") plt.tight_layout() - outPath = pU.saveAndOrPlot({'pathResult': outDir}, outFileName, fig) + outPath = pU.saveAndOrPlot({"pathResult": outDir}, outFileName, fig) return outPath + + +def avalancheThalwegPlot(ax, fieldRaster, dem, avaProfileMass, splitPoint=None, avalancheDir="", cmapHS=None): + """ + plots the location of the thalweg (including extensions to top and bottom) + on a raster Field and a hillshade + + Parameters + ---------- + ax: matplotlib.axis + axis on which the thalweg is plotted + fieldRaster: numpy array + raster that is plotted + dem: dict + dictionary containing DEM header and data + avaProfileMass: dict + contains x and y coordinates of thalweg, + and the indices of start and end of the averaged profile + splitPoint: dict or None + contains x and y coordinates of splitpoints that are plotted + avalancheDir: str + name of the avalanche directory to add optionally on the plot + cmapHS: matplotlib.colormap + optionally defined colormap for the plotted hillshade + Returns + ------- + ax: matplotlib.axis + axis containing the plotted thalweg with the raster and hillshade + """ + + indStart = avaProfileMass["indStartMassAverage"] + indEnd = avaProfileMass["indEndMassAverage"] + rowsMin, rowsMax, colsMin, colsMax = pU.constrainPlotsToData( + fieldRaster, dem["header"]["cellsize"], extentOption=True, constrainedData=False, buffer="" + ) + ax, extent, cbar0, cs1 = outCom1DFA.addResult2Plot(ax, dem["header"], fieldRaster, "pta") + cbar0.ax.set_ylabel("peak travel angle") + # add DEM hillshade with contour lines + ax = outCom1DFA.addDem2Plot(ax, dem, what="hillshade", extent=extent, cmap=cmapHS) + # add path + ax.plot( + avaProfileMass["x"][: indStart + 1], + avaProfileMass["y"][: indStart + 1], + "-y.", + zorder=20, + label="top extension", + lw=2, + path_effects=[pe.Stroke(linewidth=3, foreground="b"), pe.Normal()], + ) + ax.plot( + avaProfileMass["x"][indEnd:], + avaProfileMass["y"][indEnd:], + "-y.", + zorder=20, + label="bottom extension", + lw=2, + path_effects=[pe.Stroke(linewidth=3, foreground="g"), pe.Normal()], + ) + ax.plot( + avaProfileMass["x"][indStart : indEnd + 1], + avaProfileMass["y"][indStart : indEnd + 1], + "-y.", + zorder=20, + label="Center of mass path", + lw=2, + path_effects=[pe.Stroke(linewidth=3, foreground="k"), pe.Normal()], + ) + if splitPoint is not None: + if not splitPoint.get("isTopSplitPoint", False): # if not a top split point + ax.plot(splitPoint["x"], splitPoint["y"], "P", color="r", label="_Split point", zorder=20) + ax.set_xlabel("x [m]") + ax.set_ylabel("y [m]") + ax.axis("equal") + ax.set_ylim([rowsMin, rowsMax]) + ax.set_xlim([colsMin, colsMax]) + ax.set_title("Avalanche thalweg") + pU.putAvaNameOnPlot(ax, avalancheDir) + return ax diff --git a/avaframe/out3Plot/plotUtils.py b/avaframe/out3Plot/plotUtils.py index c239b6631..62f222811 100644 --- a/avaframe/out3Plot/plotUtils.py +++ b/avaframe/out3Plot/plotUtils.py @@ -903,6 +903,7 @@ def addHillShadeContours( else: extentPlot = extent + # TODO: restore changes and rebase master! hs = ls.hillshade(data, vert_exag=vertExag, dx=data.shape[1], dy=data.shape[0]) # normalize hillshade and increase the contrast if np.nanmin(hs) != np.nanmax(hs): diff --git a/avaframe/tests/test_DFAPathGeneration.py b/avaframe/tests/test_DFAPathGeneration.py index 4460dd1be..210db26e8 100644 --- a/avaframe/tests/test_DFAPathGeneration.py +++ b/avaframe/tests/test_DFAPathGeneration.py @@ -3,6 +3,7 @@ import math import configparser import pytest +import copy # Local imports import avaframe.ana5Utils.DFAPathGeneration as DFAPathGeneration @@ -66,7 +67,7 @@ def test_extendDFAPath(): cfg['PATH'] = {'nCellsResample': '1', 'extTopOption': '1', 'nCellsMinExtend': '1', 'nCellsMaxExtend': '2', 'factBottomExt': 0.2, 'maxIterationExtBot': 10, 'nBottomExtPrecision': 10, - 'uInterval': '1000'} + 'uInterval': '1000', 'kResample': '3'} # TODO if k=3 for spline needs at least 4 pointsin path avaProfile = {'x': np.array([1, 2, 3, 8]), 'y': np.array([1, 2, 3, 8]), 'z': np.array([40, 30, 20, 0]) @@ -101,7 +102,8 @@ def test_extendDFAPath(): # now use the highest point method cfg = configparser.ConfigParser() cfg['PATH'] = {'nCellsResample': '5', 'extTopOption': '0', 'nCellsMinExtend': '1', - 'nCellsMaxExtend': '2', 'factBottomExt': 0.2, 'maxIterationExtBot': 10, 'nBottomExtPrecision': 10} + 'nCellsMaxExtend': '2', 'factBottomExt': 0.2, 'maxIterationExtBot': 10, 'nBottomExtPrecision': 10, + 'kResample': '3'} avaProfile = {'x': np.array([10, 20, 30]), 'y': np.array([10, 20, 30]), 'z': np.array([40, 30, 20])} particlesIni = {'x': np.array([7., 6.9]), 'y': np.array([10, 20])} @@ -123,7 +125,8 @@ def test_extendDFAPath(): # now If we extend too 1 cfg = configparser.ConfigParser() cfg['PATH'] = {'nCellsResample': '5', 'extTopOption': '0', 'nCellsMinExtend': '2', - 'nCellsMaxExtend': '30', 'factBottomExt': 1, 'maxIterationExtBot': 10, 'nBottomExtPrecision': 1} + 'nCellsMaxExtend': '30', 'factBottomExt': 1, 'maxIterationExtBot': 10, 'nBottomExtPrecision': 1, + 'kResample': '3'} avaProfile = {'x': np.array([10, 20, 30, 70]), 'y': np.array([10, 20, 30, 70]), 'z': np.array([40, 30, 20, 0])} avaProfileExt = DFAPathGeneration.extendDFAPath(cfg['PATH'], avaProfile, dem, particlesIni) @@ -198,7 +201,7 @@ def test_extendProfileToFront(): 'nCellsMinExtend': '1', 'nCellsMaxExtend': '20', 'factBottomExt': 0.2, 'maxIterationExtBot': 10, 'nBottomExtPrecision': 10, 'ftThreshold': 0.01, 'lowFrontFraction': 0.05, - 'upSlopePenalty': 10., 'flowDistPenalty': 5.} + 'upSlopePenalty': 10., 'flowDistPenalty': 5., 'kResample': '3'} 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))} @@ -266,7 +269,7 @@ def test_resamplePath(): """""" # setup required inputs cfg = configparser.ConfigParser() - cfg['PATH'] = {'nCellsResample': '1', 'uInterval': '1000'} + cfg['PATH'] = {'nCellsResample': '1', 'uInterval': '1000', "kResample": "3"} avaProfile = {'x': np.array([5, 15, 20, 25, 30, 35]), 'y': np.array([5, 15, 20, 25, 30, 35]), 'z': np.array([40, 30, 20, 10, 0, 0]), 's': np.array([0, math.sqrt(200), math.sqrt(450), math.sqrt(800), math.sqrt(1250), math.sqrt(1800)]), @@ -465,3 +468,227 @@ def test_getMassAvgPathFromFields_noVelocity(): # Should have one time step assert len(result['x']) == 1 + + +def test_extendProfileTop_option0(): + """ test extending profile at top towards highest point in release (extTopOption = 0) """ + particlesIni = { + "x": np.array([10.0, 20.0, 30.0]), + "y": np.array([5.0, 15.0, 25.0]), + "z": np.array([100.0, 150.0, 120.0]), + } + profile = { + "x": np.array([30.0, 40.0, 50.0]), + "y": np.array([25.0, 35.0, 45.0]), + "z": np.array([120.0, 90.0, 60.0]), + "s": np.array([0.0, 14.14, 28.28]), + } + + result = DFAPathGeneration.extendProfileTop(0, particlesIni, copy.deepcopy(profile), dem=None, cfg=None) + + # highest particle is index 1 (z=150) + xExtTop = 20.0 + yExtTop = 15.0 + zExtTop = 150.0 + dx = xExtTop - 30.0 + dy = yExtTop - 25.0 + dsExpected = np.sqrt(dx ** 2 + dy ** 2) + + assert result["x"][0] == xExtTop + assert result["y"][0] == yExtTop + assert result["z"][0] == zExtTop + assert result["s"][0] == 0.0 + assert result["s"][1] == dsExpected + # rest of profile is preserved, shifted by ds + assert result["x"][1] == 30.0 + assert result["s"][-1] == pytest.approx(profile["s"][-1] + dsExpected) + assert len(result["x"]) == len(profile["x"]) + 1 + + +def test_extendProfileTop_option1(): + """ test extending profile at top towards point giving longest runout (extTopOption = 1) """ + particlesIni = { + "x": np.array([10.0, 20.0, 30.0]), + "y": np.array([5.0, 15.0, 25.0]), + "z": np.array([100.0, 200.0, 120.0]), + } + profile = { + "x": np.array([30.0, 40.0, 50.0]), + "y": np.array([25.0, 35.0, 45.0]), + "z": np.array([120.0, 90.0, 60.0]), + "s": np.array([0.0, 14.14, 28.28]), + } + + result = DFAPathGeneration.extendProfileTop(1, particlesIni, copy.deepcopy(profile), dem=None, cfg=None) + + xFirst, yFirst, zFirst = profile["x"][0], profile["y"][0], profile["z"][0] + sLast, zLast = profile["s"][-1], profile["z"][-1] + tanAngle = (zFirst - zLast) / sLast + + dx = particlesIni["x"] - xFirst + dy = particlesIni["y"] - yFirst + ds = np.sqrt(dx ** 2 + dy ** 2) + dz = particlesIni["z"] - zFirst + dz1 = dz - tanAngle * ds + indTop = np.argmax(dz1) + + assert result["x"][0] == particlesIni["x"][indTop] + assert result["y"][0] == particlesIni["y"][indTop] + assert result["z"][0] == particlesIni["z"][indTop] + assert result["s"][0] == 0.0 + assert result["s"][1] == pytest.approx(ds[indTop] + 0.0) + assert len(result["x"]) == len(profile["x"]) + 1 + + +def test_extendProfileTop_option2(): + """ test extending profile at top in thalweg direction, projected on a flat DEM (extTopOption = 2) """ + cellsize = 1.0 + nrows, ncols = 100, 100 + zRaster = np.zeros((nrows, ncols)) + header = { + "cellsize": cellsize, + "nrows": nrows, + "ncols": ncols, + "xllcenter": 0.0, + "yllcenter": 0.0, + } + dem = {"rasterData": zRaster, "header": header} + + cfg = configparser.ConfigParser() + cfg["GENERAL"] = {} + cfg = cfg["GENERAL"] + cfg["nCellsMinExtend"] = "1" + cfg["nCellsMaxExtend"] = "20" + cfg["factBottomExt"] = "0.1" + cfg["maxIterationExtBot"] = "10" + cfg["nBottomExtPrecision"] = "1" + + # profile running roughly along the x axis, starting well inside the DEM + profile = { + "x": np.array([50.0, 55.0, 60.0, 65.0]), + "y": np.array([50.0, 50.0, 50.0, 50.0]), + "z": np.array([0.0, 0.0, 0.0, 0.0]), + "s": np.array([0.0, 5.0, 10.0, 15.0]), + } + + result = DFAPathGeneration.extendProfileTop(2, None, copy.deepcopy(profile), dem=dem, cfg=cfg, considerLLC=False) + + # a point should have been prepended, extending "backwards" (towards lower x) + assert len(result["x"]) == len(profile["x"]) + 1 + assert result["x"][0] < profile["x"][0] + assert result["s"][0] == 0.0 + # the rest of the s values should be shifted by the added segment length + ds = result["s"][1] + assert result["s"][1] == pytest.approx(ds) + assert result["s"][2] == pytest.approx(profile["s"][1] + ds) + + +def test_extendProfileTop_option2_singlePointProfile(): + """ test that a one-point profile is returned unchanged for extTopOption = 2 """ + profile = { + "x": np.array([50.0]), + "y": np.array([50.0]), + "z": np.array([0.0]), + "s": np.array([0.0]), + } + cfg = configparser.ConfigParser() + cfg["GENERAL"] = {} + dem = {"rasterData": np.zeros((10, 10)), "header": {"cellsize": 1.0, "xllcenter": 0.0, "yllcenter": 0.0}} + + result = DFAPathGeneration.extendProfileTop(2, None, copy.deepcopy(profile), dem=dem, cfg=cfg["GENERAL"]) + + for key in result: + assert result[key] == profile[key] + assert len(result["x"]) == 1 + + +def test_extendProfileTop_option2_missingDemOrCfg(): + """ test that missing dem or cfg raises ValueError for extTopOption = 2 """ + profile = { + "x": np.array([50.0, 55.0]), + "y": np.array([50.0, 50.0]), + "z": np.array([0.0, 0.0]), + "s": np.array([0.0, 5.0]), + } + cfg = configparser.ConfigParser() + cfg["GENERAL"] = {} + + with pytest.raises(ValueError): + DFAPathGeneration.extendProfileTop(2, None, profile, dem=None, cfg=cfg["GENERAL"]) + + dem = {"rasterData": np.zeros((10, 10)), "header": {"cellsize": 1.0, "xllcenter": 0.0, "yllcenter": 0.0}} + with pytest.raises(ValueError): + DFAPathGeneration.extendProfileTop(2, None, profile, dem=dem, cfg=None) + + +def test_extendProfileTop_invalidOption(): + """ test that an invalid extTopOption raises ValueError """ + particlesIni = {"x": np.array([10.0]), "y": np.array([5.0]), "z": np.array([100.0])} + profile = { + "x": np.array([30.0, 40.0]), + "y": np.array([25.0, 35.0]), + "z": np.array([120.0, 90.0]), + "s": np.array([0.0, 14.14]), + } + + with pytest.raises(ValueError): + DFAPathGeneration.extendProfileTop(99, particlesIni, profile) + + +def test_extendProfileTop_option2_exactValues(): + """ test extending profile at top for extTopOption = 2 + + Uses a profile perfectly aligned along the x-axis so the averaged extension + direction is unambiguously (-1, 0), and a flat DEM large enough that the + extension point is found inside the domain on the first try (no dichotomy + iterations needed). + """ + cellsize = 1.0 + nrows, ncols = 100, 100 + zRaster = np.zeros((nrows, ncols)) + header = { + "cellsize": cellsize, + "nrows": nrows, + "ncols": ncols, + "xllcenter": 0.0, + "yllcenter": 0.0, + } + dem = {"rasterData": zRaster, "header": header} + + cfg = configparser.ConfigParser() + cfg["GENERAL"] = {} + cfg = cfg["GENERAL"] + cfg["nCellsMinExtend"] = "1" + cfg["nCellsMaxExtend"] = "20" + cfg["factBottomExt"] = "0.1" + cfg["maxIterationExtBot"] = "10" + cfg["nBottomExtPrecision"] = "1" + + # profile running exactly along the x axis, starting well inside the DEM + profile = { + "x": np.array([50.0, 55.0, 60.0, 65.0]), + "y": np.array([50.0, 50.0, 50.0, 50.0]), + "z": np.array([0.0, 0.0, 0.0, 0.0]), + "s": np.array([0.0, 5.0, 10.0, 15.0]), + } + + result = DFAPathGeneration.extendProfileTop(2, None, copy.deepcopy(profile), dem=dem, cfg=cfg, considerLLC=False) + + # hand-computed expected values: + # xFirst=50, yFirst=50, sTotal=15 + # points of interest (1 < r < 20): x=[55,60,65], y=[50,50,50] + # -> averaged direction is exactly (1, 0) since all vectors point along +x + # factExt = 0.1, vNorm = 1 -> gamma = 0.1 * 15 / 1 = 1.5 + # xExtTop = 50 - 1.5*1 = 48.5, yExtTop = 50 - 1.5*0 = 50 + # zExtTop = 0 (flat raster), so ds = sqrt((48.5-50)^2 + 0^2) = 1.5 + assert result["x"][0] == pytest.approx(50 - 1.5 * 1) + assert result["y"][0] == pytest.approx(50.0) + assert result["z"][0] == pytest.approx(0.0) + assert result["s"][0] == 0.0 + assert result["s"][1] == pytest.approx(1.5) + + # rest of the profile is preserved, each s-value shifted by ds = 1.5 + np.testing.assert_allclose(result["x"][1:], profile["x"]) + np.testing.assert_allclose(result["y"][1:], profile["y"]) + np.testing.assert_allclose(result["z"][1:], profile["z"]) + np.testing.assert_allclose(result["s"][1:], profile["s"] + 1.5) diff --git a/avaframe/tests/test_preparePathGeneral.py b/avaframe/tests/test_preparePathGeneral.py new file mode 100644 index 000000000..62987b1ff --- /dev/null +++ b/avaframe/tests/test_preparePathGeneral.py @@ -0,0 +1,169 @@ +"""Tests for module preparePathGeneral""" + +import numpy as np +import pytest +import configparser + +import avaframe.ana5Utils.preparePathGeneral as prepPathGeneral + + +# --------------------------------------------------------------------------- +# updateSZProfile +# --------------------------------------------------------------------------- + +def test_updateSZProfile_flatDEM(): + """ test that z is read as constant from a flat DEM and s is the cumulative 2D distance """ + cellsize = 1.0 + nrows, ncols = 100, 100 + zRaster = np.zeros((nrows, ncols)) + header = {"cellsize": cellsize, "nrows": nrows, "ncols": ncols, "xllcenter": 0.0, "yllcenter": 0.0} + dem = {"rasterData": zRaster, "header": header} + + profile = { + "x": np.array([10.0, 13.0, 13.0, 20.0]), + "y": np.array([10.0, 10.0, 14.0, 14.0]), + } + expectedS = np.array([0.0, 3.0, 7.0, 14.0]) + + result = prepPathGeneral.updateSZProfile(profile, dem) + + assert "z" in result + assert "s" in result + np.testing.assert_allclose(result["z"], [0.0, 0.0, 0.0, 0.0]) + # segment lengths: 3, 4, 7 (3-4-5 triangle then straight run) -> cumulative distance from 0 + np.testing.assert_allclose(result["s"], expectedS, atol=1e-6) + # x, y are untouched + np.testing.assert_allclose(result["x"], profile["x"]) + np.testing.assert_allclose(result["y"], profile["y"]) + + +def test_updateSZProfile_slopedDEM(): + """ test that z values are correctly read (bilinearly interpolated) from a DEM with a linear x-gradient """ + cellsize = 1.0 + nrows, ncols = 100, 100 + # elevation ramp: z = x-coordinate everywhere (constant along y) + zRaster = np.tile(np.arange(ncols) * cellsize, (nrows, 1)).astype(float) + header = {"cellsize": cellsize, "nrows": nrows, "ncols": ncols, "xllcenter": 0.0, "yllcenter": 0.0} + dem = {"rasterData": zRaster, "header": header} + + profile = { + "x": np.array([10.0, 20.0, 35.5]), + "y": np.array([50.0, 50.0, 50.0]), + } + + result = prepPathGeneral.updateSZProfile(profile, dem) + + # since the ramp is linear, bilinear interpolation should exactly recover z == x + np.testing.assert_allclose(result["z"], profile["x"], atol=1e-6) + assert result["s"][0] == 0.0 + np.testing.assert_allclose(result["s"], [0.0, 10.0, 25.5], atol=1e-6) + + +# --------------------------------------------------------------------------- +# pathExtension +# --------------------------------------------------------------------------- + +def _makeFlatDemAndCfg(extTopOption=2): + cellsize = 1.0 + nrows, ncols = 200, 200 + zRaster = np.zeros((nrows, ncols)) + header = {"cellsize": cellsize, "nrows": nrows, "ncols": ncols, "xllcenter": 0.0, "yllcenter": 0.0} + dem = {"rasterData": zRaster, "header": header} + + cfg = configparser.ConfigParser() + cfg["PATH"] = { + "extTopOption": str(extTopOption), + "nCellsMinExtend": "1", + "nCellsMaxExtend": "20", + "factBottomExt": "0.1", + "maxIterationExtBot": "10", + "nBottomExtPrecision": "1", + } + return dem, cfg + + +def test_pathExtension_option2_extendsAndSetsIndices(): + """ test that pathExtension sets the mass-average indices and extends the profile at both ends """ + dem, cfg = _makeFlatDemAndCfg(extTopOption=2) + + profile = { + "x": np.array([100.0, 105.0, 110.0, 115.0, 120.0]), + "y": np.array([100.0, 100.0, 100.0, 100.0, 100.0]), + "z": np.array([0.0, 0.0, 0.0, 0.0, 0.0]), + "s": np.array([0.0, 5.0, 10.0, 15.0, 20.0]), + } + origLen = len(profile["x"]) + + result = prepPathGeneral.pathExtension(profile, dem, cfg) + + # indices are set based on the ORIGINAL (pre-extension) length + assert result["indStartMassAverage"] == 1 + assert result["indEndMassAverage"] == origLen - 2 + + # profile grew by one point at the top and one at the bottom + assert len(result["x"]) == origLen + 2 + # top extension moves in -x direction, bottom extension moves in +x direction + assert result["x"][0] < profile["x"][1] # extended top point is further "back" + assert result["x"][-1] > profile["x"][-2] # extended bottom point is further "forward" + assert result["y"][0] == profile["y"][0] + + +def test_pathExtension_invalidTopOptionRaises(): + """ test that pathExtension raises ValueError for any extTopOption other than 2 """ + dem, cfg = _makeFlatDemAndCfg(extTopOption=0) + + profile = { + "x": np.array([100.0, 105.0, 110.0]), + "y": np.array([100.0, 100.0, 100.0]), + "z": np.array([0.0, 0.0, 0.0]), + "s": np.array([0.0, 5.0, 10.0]), + } + + with pytest.raises(ValueError): + prepPathGeneral.pathExtension(profile, dem, cfg) + + +# --------------------------------------------------------------------------- +# preparePathGeneralMain +# --------------------------------------------------------------------------- + +def test_preparePathGeneralMain_shortProfile(): + """ test the short-circuit branch for a profile with 2 or fewer points """ + dem, cfg = _makeFlatDemAndCfg(extTopOption=2) + + profile = { + "x": np.array([10.0, 20.0]), + "y": np.array([10.0, 10.0]), + } + + profileAveraged, profileExtended = prepPathGeneral.preparePathGeneralMain(profile, cfg, dem) + + # z and s were added by updateSZProfile + assert "z" in profileAveraged + assert "s" in profileAveraged + np.testing.assert_allclose(profileAveraged["s"], [0.0, 10.0]) + + # profileExtended is a copy with mass-average indices set to the short-profile defaults + assert profileExtended["indStartMassAverage"] == 0 + assert profileExtended["indEndMassAverage"] == 1 + np.testing.assert_allclose(profileExtended["x"], profileAveraged["x"]) + + # pathExtension/resamplePath must NOT have been invoked -> no extension happened, + # i.e. length is unchanged from the input + assert len(profileExtended["x"]) == 2 + + +def test_preparePathGeneralMain_shortProfile_singlePoint(): + """ test the short-circuit branch for a profile with only 1 point """ + dem, cfg = _makeFlatDemAndCfg(extTopOption=2) + + profile = { + "x": np.array([10.0]), + "y": np.array([10.0]), + } + + profileAveraged, profileExtended = prepPathGeneral.preparePathGeneralMain(profile, cfg, dem) + + assert profileExtended["indStartMassAverage"] == 0 + # max(len(x), 1) == 1 for a single-point profile + assert profileExtended["indEndMassAverage"] == 0 diff --git a/docs/moduleAna5Utils.rst b/docs/moduleAna5Utils.rst index 697ccfb08..5ee541a48 100644 --- a/docs/moduleAna5Utils.rst +++ b/docs/moduleAna5Utils.rst @@ -118,6 +118,15 @@ There are two options available to extend the mass-averaged path profile in the line. :math:`\Delta z` and :math:`\Delta s` represent the vertical and horizontal distance between a point in the release and the first point of the mass-averaged path profile. +2. Extend the path in the direction of the thalweg, upwards. The extension direction is + found from the points of the profile lying between ``nCellsMinExtend`` * cellSize and + ``nCellsMaxExtend`` * cellSize of the first point (same logic as used for the bottom extension, + ``extBottomOption = 0``). The path is then extended upwards, opposite to this direction, by a + length of ``factBottomExt`` * :math:`s_{total}`, where :math:`s_{total}` is the total length of the + profile, and the new point is projected onto the DEM. If this point falls outside the DEM, a + bisection search (up to ``maxIterationExtBot`` iterations, with a precision of + ``nBottomExtPrecision`` * cellSize) is used to find the furthest point along the extension + direction that still lies within the DEM. We also extend the path at the bottom, to have some buffer in the runout area. Two options exist (``extBottomOption``): @@ -167,6 +176,40 @@ profile. This parabolic fit determines the split point location. It is the first point for which the slope is lower than the ``slopeSplitPoint`` angle. This point is then projected on the avalanche path profile. +Path preparation from given x, y coordinates +============================================== +In addition to generating a path from a DFA simulation, it is also possible to prepare a path +starting from a simple set of x, y coordinates (for example a thalweg digitized by hand or +imported from another source). This is handled by :py:func:`ana5Utils.preparePathGeneral.preparePathGeneralMain`, +which brings a raw x, y profile through the same steps used for the mass-averaged path: + +1. **Compute z and s:** for every x, y point, the elevation is read from the DEM and the horizontal + distance travelled along the profile is computed, giving a complete (x, y, s, z) profile + (:py:func:`ana5Utils.preparePathGeneral.updateSZProfile`). + +2. **Extend to top and bottom:** the profile is lengthened using the same + :py:func:`ana5Utils.DFAPathGeneration.extendProfileTop` and + :py:func:`ana5Utils.DFAPathGeneration.extendProfileBottom` functions used for the mass-averaged + path (see :ref:`moduleAna5Utils:Path extension`), via :py:func:`ana5Utils.preparePathGeneral.pathExtension`. + At present, only ``extTopOption = 2`` is supported for this top extension; other values raise an + error. After extension, z and s values are recomputed for the new points. + +3. **Resample:** the extended profile is resampled at an approximate spacing of + ``nCellsResample`` * cellSize using :py:func:`ana5Utils.DFAPathGeneration.resamplePath` + (see :ref:`moduleAna5Utils:Resampling`). + +The function returns two profiles: + +* ``profileAveraged``: the original x, y coordinates with z and s added, left otherwise unmodified. +* ``profileExtended``: the extended and resampled version of the path, ready to be used as input + for modules such as :ref:`moduleCom2AB:com2AB: Alpha Beta Model` or :ref:`moduleAna3AIMEC:ana3AIMEC: Aimec`. + + +.. Note:: + This entry point is useful when a path profile already exists (e.g. from manual digitization) + and only the extension and resampling steps of the automated path generation are needed, + without running a DFA simulation to derive a mass-averaged path. + Distance-Time Analysis ----------------------