diff --git a/debrisframe/in1Utils/__init__.py b/debrisframe/in1Utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/debrisframe/in1Utils/generateTopo.py b/debrisframe/in1Utils/generateTopo.py new file mode 100644 index 0000000..bf85a63 --- /dev/null +++ b/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 = cfg["TOPO"].getfloat("C") + cff = cfg["CHANNELS"].getfloat("cff") + cRadius = cfg["CHANNELS"].getfloat("cRadius") + cInit = cfg["CHANNELS"].getfloat("cInit") + cMustart = cfg["CHANNELS"].getfloat("cMustart") + cMuend = cfg["CHANNELS"].getfloat("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) + + # 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) + + # 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, debrisDir): + """ + 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(debrisDir, "Inputs") + if outDir.is_dir(): + log.info("The new DEM is saved to %s" % (outDir)) + else: + log.error( + "Required folder structure: NameOfDebrisFlow/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) diff --git a/debrisframe/in1Utils/generateTopoCfg.ini b/debrisframe/in1Utils/generateTopoCfg.ini new file mode 100644 index 0000000..3922097 --- /dev/null +++ b/debrisframe/in1Utils/generateTopoCfg.ini @@ -0,0 +1,54 @@ +### Config File - This file contains the main settings for the topography generation +## Set your parameters +# This file is part of DebrisFrame. + +# General Topography parameters ---------------------- +[in3Utils_generateTopo_override] + +# use default generateTopo config as base configuration (True) and override following parameters +# if False and local_generateTopoCfg is available use local +defaultConfig = True + +# total horizontal extent of the domain [m] +xEnd = 2000 + +# 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 + +# 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 + +#------------------------------------------------------ + +# Channel parameters ----------------------------------- + +# standard channel radius - width of channel +cRadius = 20 + + # 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 + +# standard deviation sigma +cff = 120 + + + diff --git a/debrisframe/runScripts/runGenerateTopo.py b/debrisframe/runScripts/runGenerateTopo.py new file mode 100644 index 0000000..2dc27c5 --- /dev/null +++ b/debrisframe/runScripts/runGenerateTopo.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python +""" + Run script for generateTopo in module in3Utils +""" +import pathlib + +# Local imports +from debrisframe.in1Utils import generateTopo as gT +from avaframe.out3Plot import outTopo as oT +from avaframe.in3Utils import cfgUtils, logUtils, generateTopo, cfgHandling + +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 default input parameters from configuration file and update with override parameters + debGenTopoCfg = cfgUtils.getModuleConfig(gT) + defaultCfg = cfgUtils.getModuleConfig( + generateTopo, onlyDefault=debGenTopoCfg["in3Utils_generateTopo_override"].getboolean("defaultConfig") + ) + defaultCfg, debrisCfg = cfgHandling.applyCfgOverride(defaultCfg, debGenTopoCfg, generateTopo, addModValues=False) + + # Call main function to generate DEMs + [z, name_ext, outDir] = gT.generateTopo(defaultCfg, debrisDir) + + # Plot new topogrpahy + oT.plotGeneratedDEM(z, name_ext, defaultCfg, outDir, cfgMain)