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 b1fa40a8d..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,13 +145,13 @@ 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 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,18 +160,18 @@ 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) 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,11 +230,11 @@ 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 - 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 @@ -530,8 +612,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) @@ -584,7 +668,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 +677,10 @@ def saveSplitAndPath(avalancheDir, simDFrow, splitPoint, avaProfileMass, dem): com1DFA simulation dictionary Returns -------- - avaProfile: dict - resampled path profile + 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 != '': @@ -607,7 +693,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(): @@ -617,13 +703,13 @@ 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(): shutil.copy(inProjection, splitAB.with_suffix('.prj')) log.info('Saved split point to: %s', splitAB) - + return pathAB,splitAB def weightedAvgAndStd(values, weights): """ diff --git a/avaframe/ana5Utils/DFAPathGenerationCfg.ini b/avaframe/ana5Utils/DFAPathGenerationCfg.ini index 523849cd3..35ef554b4 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 = True + [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,6 +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 = True +[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 + +# 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 + +# 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/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index 3f8624802..f86d2ad33 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.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) + 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 diff --git a/avaframe/out3Plot/outCom3Plots.py b/avaframe/out3Plot/outCom3Plots.py index aac6622d0..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,84 +145,65 @@ 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): #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]) - l1 = 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 splitPoint != '': - 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 splitPoint != '': - 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', 'DFAPath') + outDir = pathlib.Path(avalancheDir, 'Outputs', 'ana5Utils', 'DFAPath') plt.tight_layout() outPath = pU.saveAndOrPlot({'pathResult': outDir}, outFileName, fig) return outPath diff --git a/avaframe/out3Plot/plotUtils.py b/avaframe/out3Plot/plotUtils.py index 98d37ed02..03cbc7b65 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/runAna5DFAPathGeneration.py b/avaframe/runAna5DFAPathGeneration.py new file mode 100644 index 000000000..39f3a11c3 --- /dev/null +++ b/avaframe/runAna5DFAPathGeneration.py @@ -0,0 +1,94 @@ +""" + 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 time +import argparse + +# local imports +from avaframe.in3Utils import cfgUtils +from avaframe.in3Utils import logUtils +from avaframe.in3Utils import initializeProject as initProj +from avaframe.ana5Utils import DFAPathGeneration + +def runAna5DFAPathGeneration(avalancheDir='', runDFAModule=''): + """ + 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 (optional) + whether to run the DFA module (True) or load existing results (False); overrides setting in ini + + 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 + + """ + + # Time the whole routine + startTime = time.time() + + # log file name; leave empty to use default runLog.log + logName = 'runAna5DFAPathGeneration' + + # 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(DFAPathGeneration) + + # 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.cleanModuleFiles(avalancheDir, DFAPathGeneration, deleteOutput=True) + + # Run or load DFA results depending on runDFAModule bool in command call + if args.runDFA: + runDFAModule = args.runDFA.lower() == 'true' + else: + 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 %.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', type=str, choices=['true', 'false'], default='', + help='add to override ini setting to run DFA module') + + args = parser.parse_args() + runAna5DFAPathGeneration(str(args.avadir), str(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) 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