From ef55613d93d091ced1dd9df8bd28461591b6dc2f Mon Sep 17 00:00:00 2001 From: leon-wagner Date: Tue, 22 Oct 2024 13:41:36 +0200 Subject: [PATCH 1/9] fixes #1038, adds extra checks for what type is expected (point, float or other) --- avaframe/in2Trans/shpConversion.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/avaframe/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index 3f8624802..9ae91440a 100644 --- a/avaframe/in2Trans/shpConversion.py +++ b/avaframe/in2Trans/shpConversion.py @@ -475,12 +475,20 @@ def writePoint2SHPfile(pointDict, pointName, fileName): """ fileName = str(fileName) w = shapefile.Writer(fileName) - w.field("name", "C") - if len(pointDict["x"]) > 1 or len(pointDict["y"]) > 1: - message = "Length of pointDict is not allowed to exceed one" + w.field('name', 'C') + if isinstance(pointDict['x'], (float, np.float64)) and isinstance(pointDict['y'], (float, np.float64)): + w.point(pointDict['x'], pointDict['y']) + elif isinstance(pointDict['x'], (np.array)) and isinstance(pointDict['y'], (np.array)): + if len(pointDict['x']) > 1 or len(pointDict['y']) > 1: + message = 'Length of pointDict is not allowed to exceed one' + log.error(message) + raise ValueError(message) + else: + w.point(pointDict['x'][0], pointDict['y'][0]) + else: + message = 'Format of point is neither float nor array of length 1' log.error(message) - raise ValueError(message) - w.point(pointDict["x"][0], pointDict["y"][0]) + raise TypeError(message) w.record(pointName) w.close() return fileName From 3e29f5030652b0f0ca83098ff43e21cd02c6c3c8 Mon Sep 17 00:00:00 2001 From: leon-wagner Date: Wed, 23 Oct 2024 10:59:49 +0200 Subject: [PATCH 2/9] rework of runScript, and corresponding adjustments in DFAPathGeneration.py - com1 overrides moved from runScript to ini - new --runDFA flag logic - path is generated from fields by default (resType changed accordingly) - moved runScript to main avaframe folder (for QGIS) - clean DFAPath output folder in avadir before generating path (necessary for QGIS, so if script is run multiple times, paths that are previously generated are not returned to QGIS) --- avaframe/ana5Utils/DFAPathGeneration.py | 15 +- avaframe/ana5Utils/DFAPathGenerationCfg.ini | 20 ++- avaframe/runComputeDFAPath.py | 149 ++++++++++++++++++++ avaframe/runScripts/runComputeDFAPath.py | 108 -------------- 4 files changed, 176 insertions(+), 116 deletions(-) create mode 100644 avaframe/runComputeDFAPath.py delete mode 100644 avaframe/runScripts/runComputeDFAPath.py diff --git a/avaframe/ana5Utils/DFAPathGeneration.py b/avaframe/ana5Utils/DFAPathGeneration.py index b1fa40a8d..420e4a422 100644 --- a/avaframe/ana5Utils/DFAPathGeneration.py +++ b/avaframe/ana5Utils/DFAPathGeneration.py @@ -69,7 +69,7 @@ def generateAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInf # read field fieldName = ['FT', 'FM'] if addVelocityInfo: - fieldName.append['FV'] + fieldName.append('FV') fieldsList, fieldHeader, timeList = com1DFA.readFields(avalancheDir, fieldName, simName=simName, flagAvaDir=True, comModule='com1DFA') # get fields header @@ -78,7 +78,7 @@ def generateAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInf 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]['FD'] > 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, _ = gT.projectOnRaster(dem, particlesIni) @@ -152,7 +152,7 @@ def getDFAPathFromField(fieldsList, fieldHeader, dem): """ compute mass averaged path from fields Also returns the averaged velocity and kinetic energy associated - The dem and fieldsList (FT, FV and FM) need to have identical dimentions and cell size. + The dem and fieldsList (FT, FM and FV) need to have identical dimensions and cell size. If FV is not provided, information about velocity and kinetic energy is not computed Parameters @@ -584,7 +584,7 @@ def saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem): avalancheDir: pathlib avalanche directory pathlib path simDFrow: pandas series - row of the siulation dataframe coresponding to the current simulation analyzed + row of the simulation dataframe corresponding to the current simulation analyzed splitPoint: dict point dictionary avaProfileMass: dict @@ -593,8 +593,10 @@ def saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem): com1DFA simulation dictionary Returns -------- - avaProfile: dict - resampled path profile + pathShpFile: str + File path to the saved shapefile for the mass-averaged path. + pointShpFile: str + File path to the saved shapefile for the split point. """ # put path back in original location if splitPoint != '': @@ -623,6 +625,7 @@ def saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem): if inProjection.is_file(): shutil.copy(inProjection, splitAB.with_suffix('.prj')) log.info('Saved split point to: %s', splitAB) + return str(pathAB), str(splitAB) def weightedAvgAndStd(values, weights): diff --git a/avaframe/ana5Utils/DFAPathGenerationCfg.ini b/avaframe/ana5Utils/DFAPathGenerationCfg.ini index 523849cd3..3d019cace 100644 --- a/avaframe/ana5Utils/DFAPathGenerationCfg.ini +++ b/avaframe/ana5Utils/DFAPathGenerationCfg.ini @@ -42,5 +42,21 @@ slopeSplitPoint = 20 dsMin = 20 # True to get the path from particles, False to get the path from fields (FT and FM) -pathFromPart = True - +pathFromPart = False + +[com1DFA_com1DFA_override] +# use default com1DFA config as base configuration (True) and override following parameters +# if False and local is available use local +# Especially one needs to save multiple intermediate time steps in order to +# compute the path. One also needs to save pta, and the correct particle or fields variables. +defaultConfig = True +tSteps = 0:5 +resType = pta|FT|FM +simTypeList = null +sphKernelRadiusTimeStepping = True +cMax = 0.01 +sphOption = 2 +massPerParticleDeterminationMethod = MPPKR +explicitFriction = 1 +frictModel = Coulomb +mucoulomb = 0.42 diff --git a/avaframe/runComputeDFAPath.py b/avaframe/runComputeDFAPath.py new file mode 100644 index 000000000..79e0c8e39 --- /dev/null +++ b/avaframe/runComputeDFAPath.py @@ -0,0 +1,149 @@ +""" + Create an avalanche path from DFA simulation results and create a split point + Configuration should be specified in DFAPathGenerationCfg.ini + or the local version of it. + It is possible to generate a path from particles or from fields. + From particles, you need to save the particles dictionaries at + multiple time steps (first, some in the middle and last) + From fields, you need the FT, FM at multiple time steps +""" +# import general python modules +import pathlib +import numpy as np +import time +import argparse + +# local imports +from avaframe.in3Utils import cfgUtils +from avaframe.in3Utils import logUtils +import avaframe.in3Utils.initializeProject as initProj +import avaframe.ana5Utils.DFAPathGeneration as DFAPath +from avaframe.in1Data import getInput as gI +from avaframe.in3Utils import fileHandlerUtils as fU +from avaframe.com1DFA import com1DFA +import avaframe.in3Utils.geoTrans as gT +from avaframe.in3Utils import cfgHandling +from avaframe.out3Plot import outCom3Plots + + +def runComputeDFAPath(avalancheDir='', runDFAModule=False): + """ + Run DFA path generation in the default configuration with only an avalanche directory as input. + + Parameters + ---------- + avalancheDir: str + path to avalanche directory (setup e.g. with init scripts) + runDFAModule: bool + if com1DFA should be run to create required DFA results; overrides ini setting + + Returns + ------- + massAvgPath: str + file path to the mass-averaged path result saved as a shapefile. + splitPoint: str + file path to the split point result saved as a shapefile. + """ + + # Time the whole routine + startTime = time.time() + + # log file name; leave empty to use default runLog.log + logName = 'runComputeDFAPath' + + # Load avalanche directory from general configuration file + # More information about the configuration can be found here + # on the Configuration page in the documentation + cfgMain = cfgUtils.getGeneralConfig() + if avalancheDir != '': + cfgMain['MAIN']['avalancheDir'] = avalancheDir + else: + avalancheDir = cfgMain['MAIN']['avalancheDir'] + + # Start logging + log = logUtils.initiateLogger(avalancheDir, logName) + log.info("MAIN SCRIPT") + log.info('Current avalanche: %s', avalancheDir) + + # load config for path generation (from DFAPathGenerationCfg.ini or its local) + cfgDFAPath = cfgUtils.getModuleConfig(DFAPath) + + # Clean input directory(ies) of old work and output files + # If you just created the ``avalancheDir`` this one should be clean but if you + # already did some calculations you might want to clean it:: + initProj.cleanSingleAvaDir(avalancheDir, deleteOutput=False) + + if runDFAModule: + # Clean input 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, fileOverride='', modInfo=False, 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', 'DFAPath') + fU.makeADir(outDir) + # write configuration to file + 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=com1DFACfgFile) + else: + # read simulation dem + demOri = gI.readDEM(avalancheDir) + dem = com1DFA.setDEMoriginToZero(demOri) + 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 particles or FV, FT and FM results + # create dataFrame of results (read DFA data) + simDF, _ = cfgUtils.readAllConfigurationInfo(avalancheDir) + + #Clean DFAPath output in avalanche directory + initProj.cleanModuleFiles(avalancheDir, DFAPath, deleteOutput=True) + + 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'] + # get the mass average path + avaProfileMass, particlesIni = DFAPath.generateAveragePath(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']) + # make the parabolic fit + parabolicFit = DFAPath.getParabolicFit(cfgDFAPath['PATH'], avaProfileMass, dem) + # here the avaProfileMass given in input is overwritten and returns only an x, y, z extended profile + avaProfileMass = DFAPath.extendDFAPath(cfgDFAPath['PATH'], avaProfileMass, dem, particlesIni) + # resample path and keep track of start and end of mass averaged part + avaProfileMass = DFAPath.resamplePath(cfgDFAPath['PATH'], dem, avaProfileMass) + # get split point + splitPoint = DFAPath.getSplitPoint(cfgDFAPath['PATH'], avaProfileMass, parabolicFit) + # make analysis and generate plots + _ = outCom3Plots.generateCom1DFAPathPlot(avalancheDir, cfgDFAPath['PATH'], avaProfileMass, dem, + parabolicFit, splitPoint, simName) + # now save the path and split point as shapefiles + massAvgPath,splitPoint = DFAPath.saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem) + + log.info("DFA path generation completed") + + # Print time needed + endTime = time.time() + log.info('Took %6.1f seconds to calculate.' % (endTime - startTime)) + + return massAvgPath,splitPoint + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description='Run computeDFAPath workflow') + parser.add_argument('avadir', metavar='avalancheDir', type=str, nargs='?', default='', + help='the avalanche directory') + parser.add_argument('--runDFA', action="store_true", + help='run Com1DFA with adjusted parameters to generate usable example results') + + args = parser.parse_args() + runComputeDFAPath(str(args.avadir), bool(args.runDFA)) diff --git a/avaframe/runScripts/runComputeDFAPath.py b/avaframe/runScripts/runComputeDFAPath.py deleted file mode 100644 index a9258de55..000000000 --- a/avaframe/runScripts/runComputeDFAPath.py +++ /dev/null @@ -1,108 +0,0 @@ -""" - Create an avalanche path from DFA simulation results and create a split point - Configuration should be specified in DFAPathGenerationCfg.ini - or the local version of it. - It is possible to generate a path from particles or from fields. - From particles, you need to save the particles dictionaries at - multiple time steps (first, some in the middle and last) - From fields, you need the FT, FM at multiple time steps -""" -import pathlib -import numpy as np -# Local imports -import avaframe.in3Utils.initializeProject as initProj -from avaframe.in3Utils import cfgUtils -from avaframe.in3Utils import logUtils -import avaframe.ana5Utils.DFAPathGeneration as DFAPath -from avaframe.in1Data import getInput as gI -from avaframe.in3Utils import fileHandlerUtils as fU -from avaframe.com1DFA import com1DFA -import avaframe.in3Utils.geoTrans as gT -# import plotting tools -from avaframe.out3Plot import outCom3Plots - -# set to true if you want to run com1DFA first to create the DFA input needed for the path generation -runDFAModule = True - -# +++++++++SETUP CONFIGURATION++++++++++++++++++++++++ -# log file name; leave empty to use default runLog.log -logName = 'runComputeDFAPath' - -# Load avalanche directory from general configuration file -cfgMain = cfgUtils.getGeneralConfig() -avaList = [cfgMain['MAIN']['avalancheDir']] -# name of avalanche directory as list, multiple possible -# avaList = ['avaAlr', 'avaHit', 'avaKot', 'avaMal', 'avaWog', 'avaGar', 'avaHelix', -# 'avaBowl', 'avaParabola', 'avaHockeySmall', 'avaHockeyChannel'] -# avaList = ['data/' + name for name in avaList] - -for avaName in avaList: - # set avaDir - avalancheDir = pathlib.Path(avaName) - # Start logging - log = logUtils.initiateLogger(avalancheDir, logName) - log.info('MAIN SCRIPT') - log.info('Current avalanche: %s', avalancheDir) - - if runDFAModule: - # Clean input directory of old work and output files from module - initProj.cleanModuleFiles(avalancheDir, com1DFA, deleteOutput=True) - # create com1DFA configuration - # read the default one (and only this one, no local is read) - com1DFACfg = cfgUtils.getModuleConfig(com1DFA, fileOverride='', toPrint=False, onlyDefault=True) - # and adapt settings. Especially one needs to save multiple intermediate time steps in order to compute the - # path. One also needs to save the correct particle or fields variables in order to compute the path - com1DFACfg['GENERAL']['tSteps'] = '0:5' - com1DFACfg['GENERAL']['resType'] = 'pta|particles' - com1DFACfg['GENERAL']['simTypeList'] = 'null' - com1DFACfg['GENERAL']['sphKernelRadiusTimeStepping'] = 'True' - com1DFACfg['GENERAL']['cMax'] = '0.01' - com1DFACfg['GENERAL']['sphOption'] = '2' - com1DFACfg['GENERAL']['massPerParticleDeterminationMethod'] = 'MPPKR' - com1DFACfg['GENERAL']['explicitFriction'] = '1' - com1DFACfg['GENERAL']['frictModel'] = 'Coulomb' - com1DFACfg['GENERAL']['mucoulomb'] = '0.42' - outDir = pathlib.Path(avalancheDir, 'Outputs', 'DFAPath') - fU.makeADir(outDir) - # write configuration to file - 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=com1DFACfgFile) - else: - # read simulation dem - demOri = gI.readDEM(avalancheDir) - dem = com1DFA.setDEMoriginToZero(demOri) - 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 particles or FV, FD and FM results - # create dataFrame of results (read DFA data) - simDF, _ = cfgUtils.readAllConfigurationInfo(avalancheDir) - - # load config for path generation (fron DFAPathGenerationCfg.ini or its local) - DFAPathCfg = cfgUtils.getModuleConfig(DFAPath) - # loop over all simulations - for simName, simDFrow in simDF.iterrows(): - log.info('Computing avalanche path from simulation: %s', simName) - pathFromPart = DFAPathCfg['PATH'].getboolean('pathFromPart') - resampleDistance = DFAPathCfg['PATH'].getfloat('nCellsResample') * dem['header']['cellsize'] - # get the mass average path - avaProfileMass, particlesIni = DFAPath.generateAveragePath(avalancheDir, pathFromPart, simName, dem) - avaProfileMass, _ = gT.prepareLine(dem, avaProfileMass, distance=resampleDistance, Point=None) - avaProfileMass['indStartMassAverage'] = 1 - avaProfileMass['indEndMassAverage'] = np.size(avaProfileMass['x']) - # make the parabolic fit - parabolicFit = DFAPath.getParabolicFit(DFAPathCfg['PATH'], avaProfileMass, dem) - # here the avaProfileMass given in input is overwriten and returns only a x, y, z extended profile - avaProfileMass = DFAPath.extendDFAPath(DFAPathCfg['PATH'], avaProfileMass, dem, particlesIni) - # resample path and keep track of start and end of mass averaged part - avaProfileMass = DFAPath.resamplePath(DFAPathCfg['PATH'], dem, avaProfileMass) - # get split point - splitPoint = DFAPath.getSplitPoint(DFAPathCfg['PATH'], avaProfileMass, parabolicFit) - # make analysis and generate plots - _ = outCom3Plots.generateCom1DFAPathPlot(avalancheDir, DFAPathCfg['PATH'], avaProfileMass, dem, - parabolicFit, splitPoint, simName) - # now save the path and split point as shape file - DFAPath.saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem) From 3c5234d1e4ae3b6f3c675812e5069e1367d15243 Mon Sep 17 00:00:00 2001 From: leon-wagner Date: Thu, 24 Oct 2024 15:48:26 +0200 Subject: [PATCH 3/9] adjust config - removed coulomb friction and other com1 overrides, path is now generated from standard run - added option to add velocity info when FV is saved --- avaframe/ana5Utils/DFAPathGenerationCfg.ini | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/avaframe/ana5Utils/DFAPathGenerationCfg.ini b/avaframe/ana5Utils/DFAPathGenerationCfg.ini index 3d019cace..bb5d06cf8 100644 --- a/avaframe/ana5Utils/DFAPathGenerationCfg.ini +++ b/avaframe/ana5Utils/DFAPathGenerationCfg.ini @@ -44,19 +44,17 @@ dsMin = 20 # True to get the path from particles, False to get the path from fields (FT and FM) pathFromPart = False +# If path from fields, True to add velocity info (needs FV saved) +addVelocityInfo = False + [com1DFA_com1DFA_override] -# use default com1DFA config as base configuration (True) and override following parameters +# use default com1DFA config as base configuration (True) # if False and local is available use local +defaultConfig = True +# and override following parameters # Especially one needs to save multiple intermediate time steps in order to # compute the path. One also needs to save pta, and the correct particle or fields variables. -defaultConfig = True tSteps = 0:5 resType = pta|FT|FM -simTypeList = null -sphKernelRadiusTimeStepping = True -cMax = 0.01 -sphOption = 2 -massPerParticleDeterminationMethod = MPPKR -explicitFriction = 1 -frictModel = Coulomb -mucoulomb = 0.42 +#resType = pta|particles +simTypeList = null \ No newline at end of file From d0fa1b89ce0d783f1a28508a8e7be63252133904 Mon Sep 17 00:00:00 2001 From: leon-wagner Date: Wed, 6 Nov 2024 13:22:52 +0100 Subject: [PATCH 4/9] set splitpoint to be generated at the top of path if it is not found, instead of not at all (avoids issues with QGIS) - adjusted corresponding error message - accordingly adjusted plot generation to exclude splitPoints set at the top --- avaframe/ana5Utils/DFAPathGeneration.py | 6 ++++-- avaframe/out3Plot/outCom3Plots.py | 6 +++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/avaframe/ana5Utils/DFAPathGeneration.py b/avaframe/ana5Utils/DFAPathGeneration.py index 420e4a422..3f4e035f1 100644 --- a/avaframe/ana5Utils/DFAPathGeneration.py +++ b/avaframe/ana5Utils/DFAPathGeneration.py @@ -530,8 +530,10 @@ def getSplitPoint(cfg, avaProfile, parabolicFit): '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, provide the split point manually.' % cfg.getfloat('slopeSplitPoint')) - splitPoint = '' + '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) diff --git a/avaframe/out3Plot/outCom3Plots.py b/avaframe/out3Plot/outCom3Plots.py index aac6622d0..dd18cd996 100644 --- a/avaframe/out3Plot/outCom3Plots.py +++ b/avaframe/out3Plot/outCom3Plots.py @@ -146,7 +146,7 @@ def generateCom1DFAPathPlot(avalancheDir, cfgPath, avaProfileMass, dem, paraboli 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 splitPoint != '': + if not splitPoint.get('isTopSplitPoint', False): ax1.plot(splitPoint['x'], splitPoint['y'], 'P', color='r', label='_Split point', zorder=20) ax1.set_xlabel('x [m]') ax1.set_ylabel('y [m]') @@ -172,7 +172,7 @@ def generateCom1DFAPathPlot(avalancheDir, cfgPath, avaProfileMass, dem, paraboli label='_Center of mass path slope', lw=2, path_effects=[pe.Stroke(linewidth=3, foreground='k'), pe.Normal()]) minY, _ = ax3.get_ylim() minX, _ = ax3.get_xlim() - if splitPoint != '': + if not splitPoint.get('isTopSplitPoint', False): ax3.axvline(x=splitPoint['s'], color='r', linewidth=1, linestyle='-.', label='_Split point') ax3.text(splitPoint['s'], minY, "%.2f m" % (splitPoint['s']), color='r', ha="right", va="bottom") ax3.axhline(y=cfgPath.getfloat('slopeSplitPoint'), color='r', linewidth=1, linestyle='-.', @@ -199,7 +199,7 @@ def generateCom1DFAPathPlot(avalancheDir, cfgPath, avaProfileMass, dem, paraboli label='Center of mass path / profile / angle', lw=1, path_effects=[pe.Stroke(linewidth=3, foreground='k'), pe.Normal()]) ax2.plot(sPara, zPara, '-k', label='Parabolic fit') - if splitPoint != '': + if not splitPoint.get('isTopSplitPoint', False): ax2.axvline(x=splitPoint['s'], color='r', linewidth=1, linestyle='-.', label='Split point') ax2.axhline(y=splitPoint['z'], color='r', linewidth=1, linestyle='-.', label='_Split point') minY, _ = ax2.get_ylim() From dc799fd91d5dad061933e6e8e251f6e79837a0a0 Mon Sep 17 00:00:00 2001 From: leon-wagner Date: Thu, 7 Nov 2024 15:59:53 +0100 Subject: [PATCH 5/9] general changes to DFAPath output plots MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - add missing discrete color levels for Travel Angle and changed it from a mix of lapaz and batlow to nuuk color map - grey out faulty legend plotting line thereby removing user warning: “no artists found with label”… in log and removes empty legend in top right corner of birds eye view panel - more descriptive runDFA flag comment in runScript --- avaframe/out3Plot/outCom3Plots.py | 10 +++++----- avaframe/out3Plot/plotUtils.py | 14 +++++++------- avaframe/runComputeDFAPath.py | 3 ++- 3 files changed, 14 insertions(+), 13 deletions(-) diff --git a/avaframe/out3Plot/outCom3Plots.py b/avaframe/out3Plot/outCom3Plots.py index dd18cd996..3779329bc 100644 --- a/avaframe/out3Plot/outCom3Plots.py +++ b/avaframe/out3Plot/outCom3Plots.py @@ -153,8 +153,8 @@ def generateCom1DFAPathPlot(avalancheDir, cfgPath, avaProfileMass, dem, paraboli ax1.axis('equal') ax1.set_ylim([rowsMin, rowsMax]) ax1.set_xlim([colsMin, colsMax]) - l1 = ax1.legend() - l1.set_zorder(40) + #ax1.legend() + #l1.set_zorder(40) ax1.set_title('Avalanche thalweg') pU.putAvaNameOnPlot(ax1, avalancheDir) @@ -181,11 +181,11 @@ def generateCom1DFAPathPlot(avalancheDir, cfgPath, avaProfileMass, dem, paraboli ax3.text(minX, cfgPath.getfloat('slopeSplitPoint'), "%.0f°" % (cfgPath.getfloat('slopeSplitPoint')), color='r', ha="left", va="bottom") ax3.axhline(y=10, color='0.8', linewidth=1, linestyle='-.', label='_Beta angle (10°)') ax3.text(minX, 10, "%.0f°" % (10), color='0.8', ha="left", va="bottom") - # ax3.axvline(x=avaProfileMass['s'][indStart], color='b', linewidth=1, linestyle='-.', label='Start of profile') - # ax3.axvline(x=avaProfileMass['s'][indEnd], color='g', linewidth=1, linestyle='-.', label='End of profile') + #ax3.axvline(x=avaProfileMass['s'][indStart], color='b', linewidth=1, linestyle='-.', label='Start of profile') + #ax3.axvline(x=avaProfileMass['s'][indEnd], color='g', linewidth=1, linestyle='-.', label='End of profile') ax3.set_xlabel('$s_{xy}$ [m]') ax3.set_ylabel('slope angle [°]') - ax3.legend() + #ax3.legend() ax3.set_title('Avalanche thalweg profile slope') # make profile plot, zoomed out diff --git a/avaframe/out3Plot/plotUtils.py b/avaframe/out3Plot/plotUtils.py index 24a90c4fd..78f370a4a 100644 --- a/avaframe/out3Plot/plotUtils.py +++ b/avaframe/out3Plot/plotUtils.py @@ -155,8 +155,8 @@ # for the choice of the colormaps, check https://www.fabiocrameri.ch/colourmaps/ # and http://hclwizard.org:3000/hclwizard/ -# multi sequential colormap for pressure +# multi sequential colormap for pressure levP = list(fU.splitIniValueToArraySteps(cfgPlotUtils["pressureColorLevels"])) # Hawaii color map colorsP = ["#B0F4FA", "#75C165", "#A96C00", "#8B0069"] @@ -168,17 +168,17 @@ colorsT = ["#FCFFC9", "#EBCE7B", "#DE9529", "#BE5A32", "#7F2B3F", "#1D0B14"] cmapT = copy.copy(cmapCrameri.lajolla) -# multi sequential colormap for speed +# multi sequential colormap for velocity levS = list(fU.splitIniValueToArraySteps(cfgPlotUtils["speedColorLevels"])) -# Batflow color map +# Batlow color map colorsS = ["#FFCEF4", "#FFA7A8", "#C19A1B", "#578B21", "#007054", "#004960", "#201158"] cmapS = copy.copy(cmapCrameri.batlow.reversed()) # multi sequential colormap for Travel Angle levTA = list(fU.splitIniValueToArraySteps(cfgPlotUtils["travelAngleColorLevels"])) -# Batflow color map -colorsTA = ["#FFCEF4", "#FFA7A8", "#C19A1B", "#578B21", "#007054", "#004960", "#201158"] -cmapTA = copy.copy(cmapCrameri.lapaz) +# Nuuk color map +colorsTA = ['#fefeb2', '#d2d184', '#bab98d', '#a0a598', '#6f878d', '#386982', '#05598c'] +cmapTA = copy.copy(cmapCrameri.nuuk.reversed()) # colormap used if no resType provided cmapNN = copy.copy(cmapCrameri.imola.reversed()) @@ -216,7 +216,7 @@ cmapSpeed = {"cmap": cmapS, "colors": colorsS, "levels": levS} -cmapTravelAngle = {"cmap": cmapTA, "colors": colorsTA} +cmapTravelAngle = {"cmap": cmapTA, "colors": colorsTA, "levels": levTA} cmapProb = {"cmap": cmapProbmap, "colors": colorsProb, "levels": levProb} diff --git a/avaframe/runComputeDFAPath.py b/avaframe/runComputeDFAPath.py index 79e0c8e39..0fde76bcf 100644 --- a/avaframe/runComputeDFAPath.py +++ b/avaframe/runComputeDFAPath.py @@ -143,7 +143,8 @@ def runComputeDFAPath(avalancheDir='', runDFAModule=False): parser.add_argument('avadir', metavar='avalancheDir', type=str, nargs='?', default='', help='the avalanche directory') parser.add_argument('--runDFA', action="store_true", - help='run Com1DFA with adjusted parameters to generate usable example results') + help='set flag to run com1DFA with cfg override to generate usable example results before' + ' generating path') args = parser.parse_args() runComputeDFAPath(str(args.avadir), bool(args.runDFA)) From c24a916e431aa8a8470fa051c66b8e12cf35e696 Mon Sep 17 00:00:00 2001 From: leon-wagner Date: Mon, 11 Nov 2024 10:06:13 +0100 Subject: [PATCH 6/9] return path and splitpoint as strings --- avaframe/runComputeDFAPath.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/avaframe/runComputeDFAPath.py b/avaframe/runComputeDFAPath.py index 0fde76bcf..8f0b29348 100644 --- a/avaframe/runComputeDFAPath.py +++ b/avaframe/runComputeDFAPath.py @@ -35,7 +35,7 @@ def runComputeDFAPath(avalancheDir='', runDFAModule=False): avalancheDir: str path to avalanche directory (setup e.g. with init scripts) runDFAModule: bool - if com1DFA should be run to create required DFA results; overrides ini setting + True if com1DFA should be run to create required DFA results before path generation; overrides ini setting Returns ------- @@ -135,7 +135,7 @@ def runComputeDFAPath(avalancheDir='', runDFAModule=False): endTime = time.time() log.info('Took %6.1f seconds to calculate.' % (endTime - startTime)) - return massAvgPath,splitPoint + return str(massAvgPath), str(splitPoint) if __name__ == "__main__": From fe7f05e1df81e142f4a2586a6ab6b74e11885c0d Mon Sep 17 00:00:00 2001 From: leon-wagner Date: Thu, 21 Nov 2024 09:20:54 +0100 Subject: [PATCH 7/9] PR feedback - set dt to 1 --- avaframe/ana5Utils/DFAPathGeneration.py | 4 ++-- avaframe/ana5Utils/DFAPathGenerationCfg.ini | 3 ++- avaframe/out3Plot/outCom3Plots.py | 2 +- .../{runComputeDFAPath.py => runAna5DFAPathGeneration.py} | 8 ++++---- 4 files changed, 9 insertions(+), 8 deletions(-) rename avaframe/{runComputeDFAPath.py => runAna5DFAPathGeneration.py} (96%) diff --git a/avaframe/ana5Utils/DFAPathGeneration.py b/avaframe/ana5Utils/DFAPathGeneration.py index 3f4e035f1..8d6c09566 100644 --- a/avaframe/ana5Utils/DFAPathGeneration.py +++ b/avaframe/ana5Utils/DFAPathGeneration.py @@ -611,7 +611,7 @@ def saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem): relName = simName.split('_')[0] inProjection = pathlib.Path(avalancheDir, 'Inputs', 'REL', relName + '.prj') # save profile in Inputs - pathAB = pathlib.Path(avalancheDir, 'Outputs', 'DFAPath', 'massAvgPath_%s_AB_aimec' % simName) + pathAB = pathlib.Path(avalancheDir, 'Outputs', 'ana5Utils', 'DFAPath', 'massAvgPath_%s_AB_aimec' % simName) name = 'massAvaPath' shpConv.writeLine2SHPfile(avaProfileMass, name, pathAB) if inProjection.is_file(): @@ -621,7 +621,7 @@ def saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem): log.warning(message) log.info('Saved path to: %s', pathAB) if splitPoint != '': - splitAB = pathlib.Path(avalancheDir, 'Outputs', 'DFAPath', 'splitPointParabolicFit_%s_AB_aimec' % simName) + splitAB = pathlib.Path(avalancheDir, 'Outputs', 'ana5Utils', 'DFAPath', 'splitPointParabolicFit_%s_AB_aimec' % simName) name = 'parabolaSplitPoint' shpConv.writePoint2SHPfile(splitPoint, name, splitAB) if inProjection.is_file(): diff --git a/avaframe/ana5Utils/DFAPathGenerationCfg.ini b/avaframe/ana5Utils/DFAPathGenerationCfg.ini index bb5d06cf8..b7793356d 100644 --- a/avaframe/ana5Utils/DFAPathGenerationCfg.ini +++ b/avaframe/ana5Utils/DFAPathGenerationCfg.ini @@ -57,4 +57,5 @@ defaultConfig = True tSteps = 0:5 resType = pta|FT|FM #resType = pta|particles -simTypeList = null \ No newline at end of file +simTypeList = null +dt = 1 \ No newline at end of file diff --git a/avaframe/out3Plot/outCom3Plots.py b/avaframe/out3Plot/outCom3Plots.py index 3779329bc..cf346e7fe 100644 --- a/avaframe/out3Plot/outCom3Plots.py +++ b/avaframe/out3Plot/outCom3Plots.py @@ -223,7 +223,7 @@ def generateCom1DFAPathPlot(avalancheDir, cfgPath, avaProfileMass, dem, paraboli ax2.set_title('Avalanche thalweg profile') outFileName = '_'.join([simName, 'DFAPath']) - outDir = pathlib.Path(avalancheDir, 'Outputs', 'DFAPath') + outDir = pathlib.Path(avalancheDir, 'Outputs', 'ana5Utils', 'DFAPath') plt.tight_layout() outPath = pU.saveAndOrPlot({'pathResult': outDir}, outFileName, fig) return outPath diff --git a/avaframe/runComputeDFAPath.py b/avaframe/runAna5DFAPathGeneration.py similarity index 96% rename from avaframe/runComputeDFAPath.py rename to avaframe/runAna5DFAPathGeneration.py index 8f0b29348..47654d9bb 100644 --- a/avaframe/runComputeDFAPath.py +++ b/avaframe/runAna5DFAPathGeneration.py @@ -26,7 +26,7 @@ from avaframe.out3Plot import outCom3Plots -def runComputeDFAPath(avalancheDir='', runDFAModule=False): +def runAna5DFAPathGeneration(avalancheDir='', runDFAModule=False): """ Run DFA path generation in the default configuration with only an avalanche directory as input. @@ -49,7 +49,7 @@ def runComputeDFAPath(avalancheDir='', runDFAModule=False): startTime = time.time() # log file name; leave empty to use default runLog.log - logName = 'runComputeDFAPath' + logName = 'runAna5DFAPathGeneration' # Load avalanche directory from general configuration file # More information about the configuration can be found here @@ -83,7 +83,7 @@ def runComputeDFAPath(avalancheDir='', runDFAModule=False): # and override with settings from DFAPath config com1DFACfg, cfgDFAPath = cfgHandling.applyCfgOverride(com1DFACfg, cfgDFAPath, com1DFA, addModValues=False) - outDir = pathlib.Path(avalancheDir, 'Outputs', 'DFAPath') + outDir = pathlib.Path(avalancheDir, 'Outputs', 'ana5Utils', 'DFAPath') fU.makeADir(outDir) # write configuration to file com1DFACfgFile = outDir / 'com1DFAPathGenerationCfg.ini' @@ -147,4 +147,4 @@ def runComputeDFAPath(avalancheDir='', runDFAModule=False): ' generating path') args = parser.parse_args() - runComputeDFAPath(str(args.avadir), bool(args.runDFA)) + runAna5DFAPathGeneration(str(args.avadir), bool(args.runDFA)) From fdf3c6b7b6c73018605abf1f8c0b2b90fc31ba89 Mon Sep 17 00:00:00 2001 From: leon-wagner Date: Tue, 22 Oct 2024 13:41:36 +0200 Subject: [PATCH 8/9] fixes #1038, adds extra checks for what type is expected (point, float or other) - adjustment to writePoint2SHPfile to account for pytest. tested and it still fixes the addressed problem. --- avaframe/in2Trans/shpConversion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/avaframe/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index 9ae91440a..f86d2ad33 100644 --- a/avaframe/in2Trans/shpConversion.py +++ b/avaframe/in2Trans/shpConversion.py @@ -478,7 +478,7 @@ def writePoint2SHPfile(pointDict, pointName, fileName): w.field('name', 'C') if isinstance(pointDict['x'], (float, np.float64)) and isinstance(pointDict['y'], (float, np.float64)): w.point(pointDict['x'], pointDict['y']) - elif isinstance(pointDict['x'], (np.array)) and isinstance(pointDict['y'], (np.array)): + elif isinstance(pointDict['x'], (np.ndarray)) and isinstance(pointDict['y'], (np.ndarray)): if len(pointDict['x']) > 1 or len(pointDict['y']) > 1: message = 'Length of pointDict is not allowed to exceed one' log.error(message) From 23a123d67add3c341b5d2248bec83fe86147c804 Mon Sep 17 00:00:00 2001 From: leon-wagner Date: Mon, 2 Dec 2024 13:20:32 +0100 Subject: [PATCH 9/9] Incorporate feedback from PR review - adjusted runDFAModule flag, added to ini - renamed functions in DFAPathGeneration.py to reflect better what they do --- avaframe/ana1Tests/energyLineTest.py | 4 +- avaframe/ana5Utils/DFAPathGeneration.py | 109 +++++++++++++++++--- avaframe/ana5Utils/DFAPathGenerationCfg.ini | 32 ++++-- avaframe/com3Hybrid/com3Hybrid.py | 2 +- avaframe/out3Plot/outCom3Plots.py | 94 +++++++---------- avaframe/runAna5DFAPathGeneration.py | 105 +++++-------------- avaframe/tests/test_DFAPathGeneration.py | 4 +- 7 files changed, 184 insertions(+), 166 deletions(-) diff --git a/avaframe/ana1Tests/energyLineTest.py b/avaframe/ana1Tests/energyLineTest.py index cc6a28db6..4a3b7da5f 100644 --- a/avaframe/ana1Tests/energyLineTest.py +++ b/avaframe/ana1Tests/energyLineTest.py @@ -33,8 +33,8 @@ def mainEnergyLineTest(avalancheDir, energyLineTestCfg, com1DFACfg, simName, dem """ log.info('Energy line test for simulation: %s' % simName) pathFromPart = energyLineTestCfg['energyLineTest'].getboolean('pathFromPart') - avaProfileMass, particlesIni = DFAPath.generateAveragePath(avalancheDir, pathFromPart, simName, dem, - addVelocityInfo=True) + avaProfileMass, particlesIni = DFAPath.generateMassAveragePath(avalancheDir, pathFromPart, simName, dem, + addVelocityInfo=True) # read pfv for plot fieldsList, fieldHeader, timeList = com1DFA.readFields(avalancheDir, ['pfv'], simName=simName, flagAvaDir=True, diff --git a/avaframe/ana5Utils/DFAPathGeneration.py b/avaframe/ana5Utils/DFAPathGeneration.py index 8d6c09566..4b38b6c47 100644 --- a/avaframe/ana5Utils/DFAPathGeneration.py +++ b/avaframe/ana5Utils/DFAPathGeneration.py @@ -10,23 +10,105 @@ import shutil # Local imports +from avaframe.in1Data import getInput as gI import avaframe.in2Trans.shpConversion as shpConv from avaframe.in3Utils import cfgUtils +from avaframe.in3Utils import fileHandlerUtils as fU +from avaframe.in3Utils import cfgHandling +import avaframe.in3Utils.initializeProject as initProj import avaframe.in3Utils.geoTrans as gT import avaframe.com1DFA.particleTools as particleTools import avaframe.com1DFA.DFAtools as DFAtls from avaframe.com1DFA import com1DFA import avaframe.out3Plot.outDebugPlots as debPlot +from avaframe.out3Plot import outCom3Plots + # create local logger # change log level in calling module to DEBUG to see log messages log = logging.getLogger(__name__) cfgAVA = cfgUtils.getGeneralConfig() debugPlot = cfgAVA['FLAGS'].getboolean('debugPlot') - -def generateAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInfo=False, flagAvaDir=True, - comModule='com1DFA'): - """ extract path from fileds or particles +def generatePathAndSplitpoint(avalancheDir, cfgDFAPath, cfgMain, runDFAModule): + """ + Parameters: + avalancheDir (str): + path to the avalanche directory + cfgDFAPath (configParser object): + DFAPath configuration + cfgMain (configParser object): + main avaframe configuration + runDFAModule (bool): + whether to run the DFA module (True) or load existing results (False) + Returns + ------- + massAvgPath: pathlib + file path to the mass-averaged path result saved as a shapefile + 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, fileOverride='', modInfo=False, 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 + 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=com1DFACfgFile) + else: # read existing simulation results + # read simulation dem + demOri = gI.readDEM(avalancheDir) + dem = com1DFA.setDEMoriginToZero(demOri) + 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 + # create dataFrame of results (read DFA data) + simDF, _ = cfgUtils.readAllConfigurationInfo(avalancheDir) + if simDF is None: + message = "Did not find any com1DFA simulations in %s/Outputs/com1DFA/" % avalancheDir + log.error(message) + raise FileExistsError(message) + + 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'] + # 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']) + # 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) + # resample path and keep track of start and end of mass averaged part + avaProfileMass = resamplePath(cfgDFAPath['PATH'], dem, avaProfileMass) + # get split point + splitPoint = getSplitPoint(cfgDFAPath['PATH'], avaProfileMass, parabolicFit) + # make analysis and generate plots + _ = 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) + + return avaPath, splitPoint + +def generateMassAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInfo=False, flagAvaDir=True, + comModule='com1DFA'): + """ extract path from fields or particles Parameters ----------- @@ -63,7 +145,7 @@ def generateAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInf particlesIni = particlesList[0] log.info('Using particles to generate avalanche path profile') # postprocess to extract path and energy line - avaProfileMass = getDFAPathFromPart(particlesList, addVelocityInfo=addVelocityInfo) + avaProfileMass = getMassAvgPathFromPart(particlesList, addVelocityInfo=addVelocityInfo) else: particlesList = '' # read field @@ -84,12 +166,12 @@ def generateAveragePath(avalancheDir, pathFromPart, simName, dem, addVelocityInf particlesIni, _ = gT.projectOnRaster(dem, particlesIni) log.info('Using fields to generate avalanche path profile') # postprocess to extract path and energy line - avaProfileMass = getDFAPathFromField(fieldsList, fieldHeader, dem) + avaProfileMass = getMassAvgPathFromFields(fieldsList, fieldHeader, dem) return avaProfileMass, particlesIni -def getDFAPathFromPart(particlesList, addVelocityInfo=False): +def getMassAvgPathFromPart(particlesList, addVelocityInfo=False): """ compute mass averaged path from particles Also returns the averaged velocity and kinetic energy associated @@ -148,7 +230,7 @@ def getDFAPathFromPart(particlesList, addVelocityInfo=False): return avaProfileMass -def getDFAPathFromField(fieldsList, fieldHeader, dem): +def getMassAvgPathFromFields(fieldsList, fieldHeader, dem): """ compute mass averaged path from fields Also returns the averaged velocity and kinetic energy associated @@ -595,10 +677,10 @@ def saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem): com1DFA simulation dictionary Returns -------- - pathShpFile: str - File path to the saved shapefile for the mass-averaged path. - pointShpFile: str - File path to the saved shapefile for the split point. + pathAB: pathlib + file path to the saved shapefile for the mass-averaged path + splitAB: pathlib + file path to the saved shapefile for the split point """ # put path back in original location if splitPoint != '': @@ -627,8 +709,7 @@ def saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem): if inProjection.is_file(): shutil.copy(inProjection, splitAB.with_suffix('.prj')) log.info('Saved split point to: %s', splitAB) - return str(pathAB), str(splitAB) - + return pathAB,splitAB def weightedAvgAndStd(values, weights): """ diff --git a/avaframe/ana5Utils/DFAPathGenerationCfg.ini b/avaframe/ana5Utils/DFAPathGenerationCfg.ini index b7793356d..0efe3ff24 100644 --- a/avaframe/ana5Utils/DFAPathGenerationCfg.ini +++ b/avaframe/ana5Utils/DFAPathGenerationCfg.ini @@ -1,4 +1,17 @@ +[GENERAL] +# flag to generate usable DFA results prior to path generation. Warning: generating results like this will delete any +# previously generated results from the com1DFA Output directory +# set to 'True' to always run DFA simulations +# set to 'False' to load existing simulation results from the "Outputs" folder +runDFAModule = False + [PATH] +# True to get the path from particles, False to get the path from fields (FT and FM) +pathFromPart = False + +# If path from fields, True to add velocity info (needs FV saved in addition to FT and FM) +addVelocityInfo = False + # the path extracted from the DFA simulation is re-sampled # re-sampling step size is defined resampleDistance = nCellsResample x cellSize) nCellsResample = 10 @@ -41,21 +54,22 @@ slopeSplitPoint = 20 # dsMin meters after the split point also have an angle bellow slopeSplitPoint dsMin = 20 -# True to get the path from particles, False to get the path from fields (FT and FM) -pathFromPart = False - -# If path from fields, True to add velocity info (needs FV saved) -addVelocityInfo = False - [com1DFA_com1DFA_override] # use default com1DFA config as base configuration (True) # if False and local is available use local defaultConfig = True # and override following parameters -# Especially one needs to save multiple intermediate time steps in order to -# compute the path. One also needs to save pta, and the correct particle or fields variables. + +# saving time step, i.e. time in seconds (first and last time step are always saved). Multiple intermediate time steps +# are required in order to compute the path 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 -#resType = pta|particles + +# list of simulations that shall be performed (null, ent, res, entres, available (use all available input data)) simTypeList = null + +# simulation time step [s]. a larger time step is applied here to reduce processing time dt = 1 \ No newline at end of file diff --git a/avaframe/com3Hybrid/com3Hybrid.py b/avaframe/com3Hybrid/com3Hybrid.py index f425a9d7c..0c942050a 100644 --- a/avaframe/com3Hybrid/com3Hybrid.py +++ b/avaframe/com3Hybrid/com3Hybrid.py @@ -99,7 +99,7 @@ def maincom3Hybrid(cfgMain, cfgHybrid): # ++++++++++ GENERATE PATH +++++++++++ # postprocess to extract path and energy line - avaProfileMass = DFAPath.getDFAPathFromPart(particlesList, addVelocityInfo=True) + avaProfileMass = DFAPath.getMassAvgPathFromPart(particlesList, addVelocityInfo=True) # make a copy because extendDFAPathKernel might modify avaProfileMass avaProfileMassExt = DFAPath.extendDFAPath(DFAPathGenerationCfg['PATH'], avaProfileMass.copy(), dem, particlesList[0]) diff --git a/avaframe/out3Plot/outCom3Plots.py b/avaframe/out3Plot/outCom3Plots.py index cf346e7fe..d7c18ea9d 100644 --- a/avaframe/out3Plot/outCom3Plots.py +++ b/avaframe/out3Plot/outCom3Plots.py @@ -127,16 +127,15 @@ def generateCom1DFAPathPlot(avalancheDir, cfgPath, avaProfileMass, dem, paraboli angleProf, tmpProf, dsProf = geoTrans.prepareAngleProfile(cfgPath.getfloat('slopeSplitPoint'), avaProfileMass, raiseWarning=False) - - # create figures and plots + # Create figures and plots fig = plt.figure(figsize=(pU.figW*2, pU.figH*1.5)) - # make the bird view plot + + # 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 @@ -146,81 +145,62 @@ def generateCom1DFAPathPlot(avalancheDir, cfgPath, avaProfileMass, dem, paraboli 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 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.legend() - #l1.set_zorder(40) ax1.set_title('Avalanche thalweg') pU.putAvaNameOnPlot(ax1, avalancheDir) - # plot angle of profile and parabola - ax3 = plt.subplot2grid((2, 2), (1, 1), colspan=1) + ax2 = plt.subplot2grid((2, 2), (1, 1), colspan=1) # add path - ax3.plot(sPara, anglePara, 'k', lw=1, label='_parabolic fit') - - ax3.plot(avaProfileMass['s'][:indStart+1], angleProf[:indStart+1], 'y.', + 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()]) - ax3.plot(avaProfileMass['s'][indEnd:], angleProf[indEnd:], 'y.', + ax2.plot(avaProfileMass['s'][indEnd:], angleProf[indEnd:], 'y.', label='_bottom extension', lw=2, path_effects=[pe.Stroke(linewidth=3, foreground='g'), pe.Normal()]) - ax3.plot(avaProfileMass['s'][indStart:indEnd+1], angleProf[indStart:indEnd+1], 'y.', + 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, _ = ax3.get_ylim() - minX, _ = ax3.get_xlim() - if not splitPoint.get('isTopSplitPoint', False): - ax3.axvline(x=splitPoint['s'], color='r', linewidth=1, linestyle='-.', label='_Split point') - ax3.text(splitPoint['s'], minY, "%.2f m" % (splitPoint['s']), color='r', ha="right", va="bottom") - ax3.axhline(y=cfgPath.getfloat('slopeSplitPoint'), color='r', linewidth=1, linestyle='-.', + 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')) - - ax3.text(minX, cfgPath.getfloat('slopeSplitPoint'), "%.0f°" % (cfgPath.getfloat('slopeSplitPoint')), color='r', ha="left", va="bottom") - ax3.axhline(y=10, color='0.8', linewidth=1, linestyle='-.', label='_Beta angle (10°)') - ax3.text(minX, 10, "%.0f°" % (10), color='0.8', ha="left", va="bottom") - #ax3.axvline(x=avaProfileMass['s'][indStart], color='b', linewidth=1, linestyle='-.', label='Start of profile') - #ax3.axvline(x=avaProfileMass['s'][indEnd], color='g', linewidth=1, linestyle='-.', label='End of profile') - ax3.set_xlabel('$s_{xy}$ [m]') - ax3.set_ylabel('slope angle [°]') - #ax3.legend() - ax3.set_title('Avalanche thalweg profile slope') + 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 - ax2 = plt.subplot2grid((2, 2), (0, 0), colspan=2) + ax3 = plt.subplot2grid((2, 2), (0, 0), colspan=2) # plot mass averaged center of mass - ax2.plot(avaProfileMass['s'][:indStart+1], avaProfileMass['z'][:indStart+1], '-y.', label='top extension', + 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()]) - ax2.plot(avaProfileMass['s'][indEnd:], avaProfileMass['z'][indEnd:], '-y.', label='bottom extension', + ax3.plot(avaProfileMass['s'][indEnd:], avaProfileMass['z'][indEnd:], '-y.', label='bottom extension', lw=1, path_effects=[pe.Stroke(linewidth=3, foreground='g'), pe.Normal()]) - ax2.plot(avaProfileMass['s'][indStart:indEnd+1], avaProfileMass['z'][indStart:indEnd+1], '-y.', + 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()]) - ax2.plot(sPara, zPara, '-k', label='Parabolic fit') - if not splitPoint.get('isTopSplitPoint', False): - ax2.axvline(x=splitPoint['s'], color='r', linewidth=1, linestyle='-.', label='Split point') - ax2.axhline(y=splitPoint['z'], color='r', linewidth=1, linestyle='-.', label='_Split point') - minY, _ = ax2.get_ylim() - minX, _ = ax2.get_xlim() - ax2.text(splitPoint['s'], minY, "%.2f m" % (splitPoint['s']), color='r', ha="left", va="bottom") - ax2.text(minX, splitPoint['z'], "%.2f m" % (splitPoint['z']), color='r', ha="left", va="bottom") - - # add runout line - s = avaProfileMass['s'][[indStart, indEnd]] - # ax2.plot(s, z0-np.tan(runOutAngleRad)*s, '-b', label='com1dfa center of mass runout line (%.2f°)' % runOutAngleDeg) - # ax2.axvline(x=s[0], color='b', linewidth=1, linestyle='-.', label='Start of profile') - # ax2.axvline(x=s[-1], color='g', linewidth=1, linestyle='-.', label='End of profile') - zLim = ax2.get_ylim() - sLim = ax2.get_xlim() - - ax2.set_xlabel('$s_{xy}$ [m]') - ax2.set_ylabel('z [m]') - ax2.set_xlim(sLim) - ax2.set_ylim(zLim) - ax2.legend() - ax2.set_title('Avalanche thalweg profile') + 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.legend() + ax3.set_title('Avalanche thalweg profile') outFileName = '_'.join([simName, 'DFAPath']) outDir = pathlib.Path(avalancheDir, 'Outputs', 'ana5Utils', 'DFAPath') diff --git a/avaframe/runAna5DFAPathGeneration.py b/avaframe/runAna5DFAPathGeneration.py index 47654d9bb..25676e2af 100644 --- a/avaframe/runAna5DFAPathGeneration.py +++ b/avaframe/runAna5DFAPathGeneration.py @@ -9,24 +9,16 @@ """ # import general python modules import pathlib -import numpy as np import time import argparse # local imports from avaframe.in3Utils import cfgUtils from avaframe.in3Utils import logUtils -import avaframe.in3Utils.initializeProject as initProj -import avaframe.ana5Utils.DFAPathGeneration as DFAPath -from avaframe.in1Data import getInput as gI -from avaframe.in3Utils import fileHandlerUtils as fU -from avaframe.com1DFA import com1DFA -import avaframe.in3Utils.geoTrans as gT -from avaframe.in3Utils import cfgHandling -from avaframe.out3Plot import outCom3Plots +from avaframe.in3Utils import initializeProject as initProj +from avaframe.ana5Utils import DFAPathGeneration - -def runAna5DFAPathGeneration(avalancheDir='', runDFAModule=False): +def runAna5DFAPathGeneration(avalancheDir='', runDFAModule=''): """ Run DFA path generation in the default configuration with only an avalanche directory as input. @@ -34,15 +26,15 @@ def runAna5DFAPathGeneration(avalancheDir='', runDFAModule=False): ---------- avalancheDir: str path to avalanche directory (setup e.g. with init scripts) - runDFAModule: bool - True if com1DFA should be run to create required DFA results before path generation; overrides ini setting + runDFAModule: bool (optional) + whether to run the DFA module (True) or load existing results (False); overrides setting in ini Returns ------- - massAvgPath: str - file path to the mass-averaged path result saved as a shapefile. - splitPoint: str - file path to the split point result saved as a shapefile. + massAvgPath: pathlib + file path to the mass-averaged path result saved as a shapefile + splitPoint: pathlib + file path to the split point result saved as a shapefile """ # Time the whole routine @@ -66,85 +58,36 @@ def runAna5DFAPathGeneration(avalancheDir='', runDFAModule=False): log.info('Current avalanche: %s', avalancheDir) # load config for path generation (from DFAPathGenerationCfg.ini or its local) - cfgDFAPath = cfgUtils.getModuleConfig(DFAPath) + cfgDFAPath = cfgUtils.getModuleConfig(DFAPathGeneration) - # Clean input directory(ies) of old work and output files - # If you just created the ``avalancheDir`` this one should be clean but if you + # Clean DFAPath output in avalanche directory + # If you just created the directory this one should be clean but if you # already did some calculations you might want to clean it:: - initProj.cleanSingleAvaDir(avalancheDir, deleteOutput=False) - - if runDFAModule: - # Clean input 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, fileOverride='', modInfo=False, 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 - 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=com1DFACfgFile) + initProj.cleanModuleFiles(avalancheDir, DFAPathGeneration, deleteOutput=True) + + # Run or load DFA results depending on runDFAModule=True or False + if runDFAModule != '': + runDFAModule = bool(runDFAModule) else: - # read simulation dem - demOri = gI.readDEM(avalancheDir) - dem = com1DFA.setDEMoriginToZero(demOri) - 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 particles or FV, FT and FM results - # create dataFrame of results (read DFA data) - simDF, _ = cfgUtils.readAllConfigurationInfo(avalancheDir) - - #Clean DFAPath output in avalanche directory - initProj.cleanModuleFiles(avalancheDir, DFAPath, deleteOutput=True) - - 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'] - # get the mass average path - avaProfileMass, particlesIni = DFAPath.generateAveragePath(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']) - # make the parabolic fit - parabolicFit = DFAPath.getParabolicFit(cfgDFAPath['PATH'], avaProfileMass, dem) - # here the avaProfileMass given in input is overwritten and returns only an x, y, z extended profile - avaProfileMass = DFAPath.extendDFAPath(cfgDFAPath['PATH'], avaProfileMass, dem, particlesIni) - # resample path and keep track of start and end of mass averaged part - avaProfileMass = DFAPath.resamplePath(cfgDFAPath['PATH'], dem, avaProfileMass) - # get split point - splitPoint = DFAPath.getSplitPoint(cfgDFAPath['PATH'], avaProfileMass, parabolicFit) - # make analysis and generate plots - _ = outCom3Plots.generateCom1DFAPathPlot(avalancheDir, cfgDFAPath['PATH'], avaProfileMass, dem, - parabolicFit, splitPoint, simName) - # now save the path and split point as shapefiles - massAvgPath,splitPoint = DFAPath.saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem) + runDFAModule = cfgDFAPath['GENERAL'].getboolean('runDFAModule') + + # Call main function to generate path and splitPoint + avaPath, splitPoint = DFAPathGeneration.generatePathAndSplitpoint(avalancheDir, cfgDFAPath, cfgMain, runDFAModule) log.info("DFA path generation completed") # Print time needed endTime = time.time() - log.info('Took %6.1f seconds to calculate.' % (endTime - startTime)) - - return str(massAvgPath), str(splitPoint) + log.info('Took %.1f seconds to calculate.' % (endTime - startTime)) + return avaPath,splitPoint if __name__ == "__main__": parser = argparse.ArgumentParser(description='Run computeDFAPath workflow') parser.add_argument('avadir', metavar='avalancheDir', type=str, nargs='?', default='', help='the avalanche directory') parser.add_argument('--runDFA', action="store_true", - help='set flag to run com1DFA with cfg override to generate usable example results before' - ' generating path') + help='add to override ini setting and set runDFAModule as True') args = parser.parse_args() runAna5DFAPathGeneration(str(args.avadir), bool(args.runDFA)) diff --git a/avaframe/tests/test_DFAPathGeneration.py b/avaframe/tests/test_DFAPathGeneration.py index 665163db2..720eba791 100644 --- a/avaframe/tests/test_DFAPathGeneration.py +++ b/avaframe/tests/test_DFAPathGeneration.py @@ -41,13 +41,13 @@ def test_getDFAPathFromPart(): particlesList = [{'nPart': 5, 'm': weights, 'x': values, 'y': values, 'z': values, 'trajectoryLengthXY': values, 'trajectoryLengthXYCor': values, 'ux': values, 'uy': values, 'uz': values}] - avaProfile = DFAPathGeneration.getDFAPathFromPart(particlesList, addVelocityInfo=False) + avaProfile = DFAPathGeneration.getMassAvgPathFromPart(particlesList, addVelocityInfo=False) # print(avaProfile) for prop in ['x', 'y', 'z', 's', 'sCor']: assert avaProfile[prop] == 3 assert avaProfile[prop + 'std'] == 1.5 - avaProfile = DFAPathGeneration.getDFAPathFromPart(particlesList, addVelocityInfo=True) + avaProfile = DFAPathGeneration.getMassAvgPathFromPart(particlesList, addVelocityInfo=True) # print(avaProfile) for prop in ['x', 'y', 'z', 's', 'sCor']: assert avaProfile[prop] == 3