update DFA path runscript [ana5; runscripts]

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,

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):
     """

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,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

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])

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