[in1Utils/out1Plot/runScripts] Function for generating generic debris-flow topography

debrisframe/in1Utils/generateTopo.py

@@ -0,0 +1,172 @@
+"""
+  Create generic/idealised topographies
+"""
+
+# load modules
+import logging
+import numpy as np
+from scipy.stats import norm
+import pathlib
+
+# local imports
+import avaframe.in3Utils.generateTopo as genTop 
+
+# create local logger
+# change log level in calling module to DEBUG to see log messages
+log = logging.getLogger("avaframe.debrisframe.in1Utils")
+
+def debrisFlowTopoAverage(cfg):
+    """
+    Compute coordinates of an average parabolic-shaped slope as a generic topography for debris-flow simulations
+    defined by a 2nd-degree polynomial: ax**2 + bx + c
+    The parameters of the polynomial function are derived from real watershed profiles (Kessler 2019)
+
+    Parameters
+    ----------------------
+    cfg: configparser Object
+        configuration setup for topo generation
+    """
+    # input parameters
+    C = float(cfg["TOPO"]["C"]) # C = 709 m
+    cff = float(cfg["CHANNELS"]["cff"])
+    cRadius = float(cfg["CHANNELS"]["cRadius"])
+    cInit = float(cfg["CHANNELS"]["cInit"])
+    cMustart = float(cfg["CHANNELS"]["cMustart"])
+    cMuend = float(cfg["CHANNELS"]["cMuend"])
+
+    # Get grid definitons
+    dx, xEnd, yEnd = genTop.getGridDefs(cfg)
+
+    # Compute coordinate grid
+    xv, yv, zv, x, y, nRows, nCols = genTop.computeCoordGrid(dx, xEnd, yEnd)
+
+    # If a channel shall be introduced
+    # Get parabola Parameters
+    [A, B, fLen] = genTop.getParabolaParams(cfg) # fLen = x0 = 1679 m
+
+    # Set surface elevation
+    mask = np.zeros(np.shape(xv))
+    mask[np.where(xv < fLen)] = 1 
+    zv = (A * xv ** 2 + B * xv + C) * mask
+
+    # If a channel shall be introduced
+    if cfg["TOPO"].getboolean("channel"):
+        # Compute cumulative distribution function - c1 for upper part (start)
+        # of channel and c2 for lower part (end) of channel
+        c1 = norm.cdf(xv, cMustart * fLen, cff)
+        c2 = 1.0 - norm.cdf(xv, cMuend * fLen, cff) # cMuend = 0.62, cff = 120
+
+        # combine both into one function separated at the the middle of
+        #  the channel longprofile location
+        mask_c1 = np.zeros(np.shape(xv))
+        mask_c1[np.where(xv < (fLen * (0.5 * (cMustart + cMuend))))] = 1
+        c0 = c1 * mask_c1
+
+        mask_c2 = np.zeros(np.shape(xv))
+        mask_c2[np.where(xv >= (fLen * (0.5 * (cMustart + cMuend))))] = 1
+        c0 = c0 + c2 * mask_c2
+
+        # Is the channel of constant width or narrowing
+        if cfg["TOPO"].getboolean("narrowing"):
+            # upper part of channel: constant width
+            mask_c1 = np.zeros(np.shape(xv))
+            mask_c1[np.where(xv < (fLen * (0.5 * (cMustart + cMuend))))] = 1
+            cExtent_c1 = np.zeros(np.shape(xv)) + cRadius
+
+            # lower part of channel: narrowing
+            mask_c2 = np.zeros(np.shape(xv))
+            mask_c2[np.where(xv >= (fLen * (0.5 * (cMustart + cMuend))))] = 1
+            cExtent_c2 = cInit * (1 - c0[:]) + (c0[:] * cRadius)
+
+            cExtent = cExtent_c1 * mask_c1 + cExtent_c2 * mask_c2
+        else:
+            cExtent = np.zeros(np.shape(xv)) + cRadius
+
+        # Set surface elevation
+        mask = np.zeros(np.shape(y))
+        mask[np.where(abs(y) < cExtent)] = 1
+        # Add surface elevation modification introduced by channel
+        if cfg["TOPO"].getboolean("topoAdd"):
+            zv = zv + cRadius * c0 * (1.0 - np.sqrt(np.abs(1.0 - (np.square(y) / (cExtent ** 2))))) * mask # changed from cExtent to cRadius
+            # outside of the channel, add layer of channel thickness
+            mask = np.zeros(np.shape(y))
+            mask[np.where(abs(y) >= cExtent)] = 1
+            mask_c2 = np.ones(np.shape(xv)) #added, smooth transition from upper to lower part of channel
+            c0 = c2 * mask_c2 #added to extend lower distribution to upper edge of topography
+            zv = zv + cRadius * mask * c0 # changed from cExtent to cRadius
+        else:
+            zv = zv - cExtent * c0 * np.sqrt(np.abs(1.0 - (np.square(y) / (cExtent ** 2)))) * mask
+
+    # Log info here
+    log.info("Generic debris-flow topography is computed")
+
+    return x, y, zv
+
+def generateTopo(cfg, avalancheDir):
+    """  
+    Compute coordinates of desired topography with given inputs  
+
+    Parameters  
+    ----------------------  
+    cfg: configparser Object  
+        configuration setup for topo generation  
+    debrisDir: string  
+        directory to data folder  
+    """ 
+
+    # Which DEM type
+    demType = cfg["TOPO"]["demType"]
+
+    log.info("DEM type is set to: %s" % demType)
+
+    # Set Output directory
+    outDir = pathlib.Path(avalancheDir, "Inputs")
+    if outDir.is_dir():
+        log.info("The new DEM is saved to %s" % (outDir))
+    else:
+        log.error(
+            "Required folder structure: NameOfAvalanche/Inputs missing! \
+                    Run runInitializeProject first!"
+        )
+
+    # Call topography type
+    if demType == "FP":
+        [x, y, z] = genTop.flatplane(cfg)
+
+    elif demType == "IP":
+        [x, y, z] = genTop.inclinedplane(cfg)
+
+    elif demType == "PF":
+        [x, y, z] = genTop.parabola(cfg)
+
+    elif demType == "TPF":
+        [x, y, z] = genTop.parabolaRotation(cfg)
+
+    elif demType == "HS":
+        [x, y, z] = genTop.hockey(cfg)
+
+    elif demType == "BL":
+        [x, y, z] = genTop.bowl(cfg)
+
+    elif demType == "HX":
+        [x, y, z] = genTop.helix(cfg)
+
+    elif demType == "PY":
+        [x, y, z] = genTop.pyramid(cfg)
+
+    elif demType == "DFTA":
+        [x, y, z] = debrisFlowTopoAverage(cfg)
+
+    # If a drop shall be introduced
+    if cfg["TOPO"].getboolean("drop"):
+        z = genTop.addDrop(cfg, x, y, z)
+
+    # moves topo in z direction
+    if cfg["DEMDATA"]["zEdit"] != "":
+        z = z + cfg["DEMDATA"].getfloat("zEdit")
+        log.info("Changed topo elevation by %.2f" % cfg["DEMDATA"].getfloat("zEdit"))
+
+    # Write DEM to file
+    genTop.writeDEM(cfg, z, outDir)
+
+    return (z, demType, outDir)

debrisframe/in1Utils/generateTopoCfg.ini

@@ -0,0 +1,139 @@
+### Config File - This file contains the main settings for the topography generation
+## Set your parameters
+# This file is part of DebrisFrame.
+
+
+# General Topography parameters ----------------------
+[TOPO]
+# DEM spatial resolution [m]
+dx = 5.
+
+# total horizontal extent of the domain [m]
+xEnd = 2000
+
+# total horizontal extent of the domain [m]
+yEnd = 1500
+
+# topography type
+# demType - topography type options:
+# FP (Flat plane), IP (Inclined plane) [dx, xEnd, yEnd, zElev]
+# PF (Parabolic slope with flat foreland) [dx, xEnd, yEnd, fLens or meanAlpha, C, optional:channel, dam]
+# TPF (Triple parabolic slope with flat foreland) [dx, xEnd, yEnd, fLens, fFlat, C]
+# HS (Hockeystick with linear slope and flat foreland and smooth transition) [dx, xEnd, yEnd, meanAlpha, z0, rCirc, optional:channel]
+# BL (Bowl-shaped topography) [dx, xEnd, yEnd, rBowl]
+# HX (Helix-shaped topography) [dx, xEnd, yEnd, flens or meanAlpha, C, rHelix, optional:channel]
+# PY (pyramid-shaped topography, optional with flat foreland) [dx, xEnd, yEnd, meanAlpha, z0, optional:flatx, flaty, phi]
+# DFTA (generic debris-flow topography) [dx, xEnd, yEnd, fLens, meanAlpha = 0, channel, narrowing]
+demType = DFTA
+
+# distance to point where slope transitions into flat plane [m] - required for PF, HX, DFTA if meanAlpha is not provided
+fLens = 1679
+
+# legnth [m] of the flat foreground for the triple parabola for TPF
+fFlat = 500
+
+# slope angle from max. elevation to start flat plane [°] - or slope of inclined plane [°]
+# this parameter required for IP, HS, PY, (PF, HX - if not fLens is used)
+meanAlpha = 0
+
+# total fall height [m] - required for PF, HX, DFTA
+C = 709
+
+# bowl radius [m] - required for BL
+rBowl = 500
+
+# radius for helix [m] - required for HX
+rHelix = 1250
+
+# max elevation of inclined plane [m] - required for  IP, HS, PY
+z0 = 2200
+
+# elevation of flat plane [m] - required for FP
+zElev = 0
+
+# radius of smoothing circle [m] - required for HS
+rCirc = 200
+
+# flatland in front of pyramid - required for PY - optional
+flatx = 2500
+flaty = 2500
+
+# flags to set channel and narrowing of channel and if channel is cut into
+# topo or a layer is added use topoAdd, True = additional topo is superimposed, False = additional topo is cut out
+channel = True
+narrowing = True
+topoAdd = True
+# dam option for PF topo
+dam = False
+# Add a step/drop in the topography (add a step in the x direction)
+drop = False
+
+# flag if pyramid topography is rotated along z-axis
+flagRot = True
+# rotation angle along z-axis [degree]
+phi = 25.
+
+# flag if not full bowl but only curved slope in x-direction
+curvedSlope = False
+#------------------------------------------------------
+
+
+# Channel parameters -----------------------------------
+[CHANNELS]
+# standard channel radius - Gerinnebreite
+cRadius = 20
+
+# start and end half width of channel that is narrowing in the middle part
+cInit = 250
+ # mean mu - represents upper part of the channel (e.g. 10% of sloping topography part)
+cMustart = 0.1
+
+# mean mu - represents lower part of the channel (e.g. 62% of sloping topography part)
+cMuend = 0.62
+
+# mean mu - represents lower part of the channel (86% of sloping topography part) - flat plane
+cMuendFP = 0.86
+
+# standard deviation sigma
+cff = 120
+#--------------------------------------------------------------
+
+# DAM parameters -----------------------------------
+[DAMS]
+# relative position with respect to where flat outrun starts
+damPos = 0.6
+
+damWidth = 75
+
+damHeight = 75
+
+# Add a drop in the topography (in x direction)
+[DROP]
+# angle of slope of the step in [°]
+alphaDrop = 80
+
+# x position of the start of the drop
+xStartDrop = 1000
+# length of the drop in the x direction
+dxDrop = 20
+
+#DEM outputfile parameters---------------------------------
+[DEMDATA]
+# x coordinate of lower left center
+xl = 1000.0
+
+# y-coordinate of lower left center
+yl = -5000.0
+
+# Prefix of DEM file name
+demName = DEM
+
+# no data value
+nodata_value = -9999
+#---------------------------------------------------------
+
+# moves topography in z direction e.g. +100m or -100m 
+zEdit = 
+
+
+

debrisframe/out1Plot/outTopo.py

@@ -0,0 +1,41 @@
+"""
+    Simple plotting for DEMs
+"""
+
+import logging
+
+# local imports
+from avaframe.in3Utils import geoTrans
+import avaframe.out3Plot.plotUtils as pU
+import avaframe.out3Plot.outTopo as oT
+
+# create local logger
+log = logging.getLogger("avaframe.debrisframe.out1Plot")
+
+def plotGeneratedDEM(z, nameExt, cfg, outDir, cfgMain):
+    """ Plot DEM with given information on the origin of the DEM """
+
+    cfgTopo = cfg['TOPO']
+    cfgDEM = cfg['DEMDATA']
+
+    # input parameters
+    dx = float(cfgTopo['dx'])
+    xl = float(cfgDEM['xl'])
+    yl = float(cfgDEM['yl'])
+    demName = cfgDEM['demName']
+
+    # Set coordinate grid with given origin
+    nrows, ncols = z.shape
+    X, Y = geoTrans.makeCoordinateGrid(xl, yl, dx, ncols, nrows)
+
+    topoNames = {'IP': 'inclined Plane', 'FP': 'flat plane', 'PF': 'parabola flat', 'TPF': 'triple parabola flat',
+                 'HS': 'Hockeystick smoothed', 'BL': 'bowl', 'HX': 'Helix', 'PY': 'Pyramid', 'DFTA': 'generic debris-flow topography'}
+
+    ax, fig = oT._generateDEMPlot(X, Y, z, topoNames[nameExt])
+
+    # Save figure to file
+    outName = '%s_%s_plot' % (demName, nameExt)
+    # save and or show figure
+    plotPath = pU.saveAndOrPlot({'pathResult': outDir}, outName, fig)
+    log.info('Saving plot to: %s', plotPath)
+

debrisframe/runScripts/runGenerateTopo.py

@@ -0,0 +1,31 @@
+#!/usr/bin/env python
+"""
+    Run script for generateTopo in module in3Utils
+"""
+import pathlib
+
+# Local imports
+from debrisframe.in1Utils import generateTopo as gT
+from debrisframe.out1Plot import outTopo as oT
+from avaframe.in3Utils import cfgUtils, logUtils
+
+if __name__ == '__main__':
+    # log file name; leave empty to use default runLog.log
+    logName = 'generateTopo'
+
+    # Load avalanche directory from general configuration file
+    cfgMain = cfgUtils.getGeneralConfig(nameFile=pathlib.Path("debrisframeCfg.ini"))
+    debrisDir = cfgMain["MAIN"]["avalancheDir"]
+
+    # Start logging
+    log = logUtils.initiateLogger(debrisDir, logName)
+    log.info('MAIN SCRIPT')
+
+    # Load input parameters from configuration file
+    cfg = cfgUtils.getModuleConfig(gT, debrisDir)
+
+    # Call main function to generate DEMs
+    [z, name_ext, outDir] = gT.generateTopo(cfg, debrisDir)
+
+    # Plot new topogrpahy
+    oT.plotGeneratedDEM(z, name_ext, cfg, outDir, cfgMain)