[ana5]: add general thalweg extension & resampling

avaframe/ana5Utils/DFAPathGenerationCfg.ini

@@ -15,10 +15,13 @@ addVelocityInfo = False
 # the path extracted from the DFA simulation is re-sampled
 # re-sampling step size is defined resampleDistance = nCellsResample x cellSize)
 nCellsResample = 10
+# Degree of the spline for splprep.
+kResample = 3
 
 # extension method at the top
 # option 0: take the highest particle in the release
 # option 1: find the point that will lead to the longest runout
+# option 2: extend the path in the direction of its first points upslope by factBottomExt x sMax
 extTopOption = 1
 
 # extension method at the bottom

avaframe/ana5Utils/preparePathGeneral.py

@@ -0,0 +1,140 @@
+"""
+generate thalweg from x and y coordinates (including extension to top and bottom and resampling)
+"""
+
+import numpy as np
+import logging
+import copy
+# local imports
+import avaframe.in3Utils.geoTrans as gT
+from avaframe.ana5Utils import DFAPathGeneration
+# create local logger
+log = logging.getLogger(__name__)
+
+
+def preparePathGeneralMain(profile, cfgDFAPath, dem):
+    """
+     prepare thalweg from x and y coordinates:
+     1. read z coordinates from DEM and compute horizontally projected distance
+     2. extend path to bottom and top
+     3. resample path points
+
+     Parameters
+     -------------
+     profile: dict
+        contains x and y coordinates of thalweg location
+    cfgDFAPath: configparser object
+        configuration for DFA path generation
+    dem: dict
+        dictionary with header and raster data of elevation model
+
+    Returns
+    -------------
+    profileAveraged: dict
+        s and z coordinates are added (x and y original) to input profile
+    profileExtended: dict
+        x, y, s, z of extended and resampled path
+    """
+    # get profile with normalized x and y coordinates and z and s values
+    profileAveraged = updateSZProfile(profile, dem)
+    profileExtended = copy.deepcopy(profileAveraged)
+
+    # skip profile that only contains one or two points
+    if len(profileAveraged["x"]) <= 2:
+        profileExtended["indStartMassAverage"] = 0
+        profileExtended["indEndMassAverage"] = max(len(profileExtended["x"]) - 1, 0)
+        return profileAveraged, profileExtended
+
+    # if extTopOption == 2, particlesIni are not used!!
+    profileExtended = pathExtension(profileExtended, dem, cfgDFAPath)
+    profileExtended = updateSZProfile(profileExtended, dem)
+    # resample profile/ path and save in an extra dictionary
+    profileExtended = DFAPathGeneration.resamplePath(cfgDFAPath["PATH"], dem, profileExtended)
+
+    # add input parameters to extended profile if they exist
+    for inputPara in ["alpha", "exponent", "zDeltaMax"]:
+        if inputPara in profile.keys():
+            profileExtended[inputPara] = profile[inputPara]
+
+    return profileAveraged, profileExtended
+
+
+def pathExtension(profile, demDict, cfgPathGen):
+    """
+    extend thalweg to top and bottom of path
+
+    Parameters
+    ------------
+    profile: dict
+        thalweg data
+    demDict: dict
+        DEM data
+    cfgPathGen: confiparser object
+        configuration setup for DFA Path generation
+
+    Returns
+    -------------
+    profile: dict
+        thalweg data that are extended to top and bottom
+    """
+
+    profile["indStartMassAverage"] = 1
+    # do not use the last two points because the last points are weird sometimes
+    profile["indEndMassAverage"] = np.size(profile["x"]) - 2
+
+    if cfgPathGen["PATH"].getint("extTopOption") != 2:
+        # TODO: if we provide particlesIni, the other options would also work.
+        message = "Up to now only top-extension option 2 works!"
+        log.error(message)
+        raise ValueError(message)
+
+    profile = DFAPathGeneration.extendProfileTop(
+        cfgPathGen["PATH"].getint("extTopOption"),
+        {},
+        profile,
+        dem=demDict,
+        cfg=cfgPathGen["PATH"],
+        considerLLC=True,
+    )
+
+    # extend the bottom
+    profile = DFAPathGeneration.extendProfileBottom(cfgPathGen["PATH"], demDict, profile, considerLLC=True)
+
+    return profile
+
+
+def updateSZProfile(profile, dem):
+    """
+    for given coordinates (of the talweg) read z values
+    from DEM and compute distance between coordinates
+
+    Parameters
+    ------------
+    profile: dict
+        contains at least x and y coordinates
+    dem: dict
+        contains dem data
+
+    Returns
+    -----------
+    profile: dict
+        profile with added a and z values
+    """
+    x = profile["x"]
+    y = profile["y"]
+
+    demHeader = dem["header"]
+
+    z, _ = gT.projectOnGrid(
+        x,
+        y,
+        dem["rasterData"],
+        csz=demHeader["cellsize"],
+        xllc=demHeader["xllcenter"],
+        yllc=demHeader["yllcenter"],
+    )
+    s = np.append([0], gT.computeLengthOfLine2D(x, y))
+    profile["z"] = z
+    profile["s"] = s
+
+    return profile

avaframe/com1DFA/DFAtools.py

@@ -159,3 +159,44 @@ def scalProd(ux, uy, uz, vx, vy, vz):
     scal = ux*vx + uy*vy + uz*vz
 
     return scal
+
+
+def getAveragedDirection(vDirX, vDirY, vDirZ=None):
+    """
+    Compute a single averaged unit direction vector from multiple 2D direction vectors.
+
+    Parameters
+    ----------
+    vDirX: numpy array
+        x components of the direction vectors to average
+    vDirY: numpy array
+        y components of the direction vectors to average
+    vDirZ: numpy array, optional
+        z components of the direction vectors to average. If None (default), the
+        vectors are treated as 2D and z is set to 0.
+
+    Returns
+    -------
+    vDirX: float
+        x component of the normalized, averaged direction vector
+    vDirY: float
+        y component of the normalized, averaged direction vector
+    vDirZ: float
+        z component of the normalized, averaged direction vector (0 if input is 2D)
+    """
+    if vDirZ is None:
+        vDirZ = 0 * np.array([vDirY])
+    else:
+        vDirZ = np.array([vDirZ])
+
+    vDirX, vDirY, vDirZ = normalize(
+        np.array([vDirX]), np.array([vDirY]), vDirZ)
+
+    # summed / averaged direction
+    vDirX = np.sum(vDirX)
+    vDirY = np.sum(vDirY)
+    vDirZ = np.sum(vDirZ)
+    # get unit vector
+    vDirX, vDirY, vDirZ = normalize(np.array([vDirX]), np.array([vDirY]), np.array([vDirZ]))
+
+    return vDirX, vDirY, vDirZ

avaframe/out3Plot/outCom1DFA.py

@@ -239,7 +239,7 @@ def addParticles2Plot(particles, ax, dem, whatS="m", whatC="h", colBarResType=""
     return ax, cb
 
 
-def addDem2Plot(ax, dem, what="slope", extent="", origHeader=False):
+def addDem2Plot(ax, dem, what="slope", extent="", origHeader=False, cmap=None):
     """Add dem to the background of a plot
 
     Parameters
@@ -255,6 +255,8 @@ def addDem2Plot(ax, dem, what="slope", extent="", origHeader=False):
         optional: extent of NonUnifIm plot corresponding coordinates to dem data array at center locations
     origHeader: bool
         if True use originalHeader and not header
+    cmap: matplotlib colormap
+        colormap for DEM plot
     """
     if origHeader:
         header = dem["originalHeader"]
@@ -267,8 +269,10 @@ def addDem2Plot(ax, dem, what="slope", extent="", origHeader=False):
     csz = header["cellsize"]
     xArray = np.linspace(xllc, xllc + (ncols - 1) * csz, ncols)
     yArray = np.linspace(yllc, yllc + (nrows - 1) * csz, nrows)
-    cmap = pU.cmapGreys
-    cmap.set_bad(color="white")
+
+    if cmap is None:
+        cmap = pU.cmapGreys
+        cmap.set_bad(color="white")
 
     if what == "slope":
         value = dem["Nz"] / DFAtls.norm(dem["Nx"], dem["Ny"], dem["Nz"])