diff --git a/debrisframe/c1TIF/c1TIFCfg.ini b/debrisframe/c1TIF/c1TIFCfg.ini index 829fd31..0212dde 100644 --- a/debrisframe/c1TIF/c1TIFCfg.ini +++ b/debrisframe/c1TIF/c1TIFCfg.ini @@ -3,8 +3,11 @@ [GENERAL] -# so far no general parameters to be set as debris flow option is fully implemented in com1DFA -# hence checkout the com1DFA_override section + +# if inputHydrograph is True, the in2TopoHyd module is computed before the c1TIF computation +# Hydrograph input data need to be provided and the flag timeDependentRelease needs to be set to TRue +inputHydrograph = False + [com1DFA_com1DFA_override] @@ -23,7 +26,7 @@ simTypeList = available modelType = dfa #+++++++++++++ Output++++++++++++ -# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, particles) - separated by | +# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, dmDet, sfcChange, demAdapted, timeInfo, entrDepth, particles) - separated by | resType = pft|pfv|ppr #+++++++++Debris flow properties @@ -38,6 +41,11 @@ rhoEnt = 2000 # mass per particle (if MPPDIR is used) [kg] massPerPart = 280000. +#+++++++++++++Mesh +# remesh the input rasters or look for remeshed rasters +# expected mesh size [m], default in AvaFrame: 5 m +meshCellSize = 2 + #+++++++++++++general start conditions: time dependent release # if timeDependentRelease is True, provide the the timesteps, thickness and velocity # for a releases in a csv-table in the REL folder @@ -141,7 +149,33 @@ entEroEnergy = 5000 adaptSfcStopped = 0 adaptSfcDetrainment = 0 adaptSfcEntrainment = 0 +# deposition is only entrainable if it is an entrainment simulation and if adaptSfcStopped = 1 and adaptSfcEntrainment = 1 +# then mass that is stopped can be entrained (also if the entrainment area is not in this stopped area) +entrainableDeposition = False # only adapt topography if changing height in at least one cell is > thresholdAdaptSfc [m] thresholdAdaptSfc = 0.1 # use the adapted topography as background in the report plots adaptDemPlot = False + +# specify a particular release area scenario, provide name of shapefile with or without extension .shp (optional) +releaseScenario = + + +[in2TopoHyd_in2TopoHyd_override] +# use default in2TopoHyd config as base configuration (True) and override following parameters +# if False and local_in2TopoHydCfg is available use local +defaultConfig = True + +#++++++++++++++++ Hydraulic Boundary Conditions +# Elevation steps in vertical direction to calclules the flow area +dElev = 0.1 +# method how the flow velocity as initial condition is calculated +# currently only the method after Rickenmann (1999) is available +velType = rickenmann +# mean slope for the cross section in flow direction [m/m] +# optional; if left blank, the slope is computed automatically +slope = + +# export cell centers along the release line as points in a csv-file; for plausability check +# True or False +exportCrossSectionCells = True diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/HYDR/hydrograph.csv b/debrisframe/data/debrisParabolaChanneled/Inputs/HYDR/hydrograph.csv new file mode 100644 index 0000000..29d702d --- /dev/null +++ b/debrisframe/data/debrisParabolaChanneled/Inputs/HYDR/hydrograph.csv @@ -0,0 +1,18 @@ +timestep,discharge +0,14.17 +5,28.33 +10,42.5 +15,56.67 +20,70.83 +25,85.0 +30,77.92 +35,70.83 +40,63.75 +45,56.67 +50,49.58 +55,42.5 +60,35.42 +65,28.33 +70,21.25 +75,14.17 +80,7.08 diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.cpg b/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.cpg new file mode 100644 index 0000000..3ad133c --- /dev/null +++ b/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.dbf b/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.dbf new file mode 100644 index 0000000..c2aa4c1 Binary files /dev/null and b/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.dbf differ diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.prj b/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.prj new file mode 100644 index 0000000..10a9c72 --- /dev/null +++ b/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.prj @@ -0,0 +1 @@ +PROJCS["ETRS_1989_Austria_Lambert",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",400000.0],PARAMETER["Central_Meridian",13.3333333333333],PARAMETER["Standard_Parallel_1",49.0],PARAMETER["Standard_Parallel_2",46.0],PARAMETER["Latitude_Of_Origin",47.5],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.shp b/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.shp new file mode 100644 index 0000000..977cc0b Binary files /dev/null and b/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.shp differ diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.shx b/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.shx new file mode 100644 index 0000000..2ff0354 Binary files /dev/null and b/debrisframe/data/debrisParabolaChanneled/Inputs/LEVEE/crossSectionLevee.shx differ diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/REL/releaseParabolaChanneledTimeDep.csv b/debrisframe/data/debrisParabolaChanneled/Inputs/REL/releaseParabolaChanneledTimeDep.csv deleted file mode 100644 index 36d3842..0000000 --- a/debrisframe/data/debrisParabolaChanneled/Inputs/REL/releaseParabolaChanneledTimeDep.csv +++ /dev/null @@ -1,13 +0,0 @@ -timestep,thickness,velocity -0,0.965,3.7 -20,7.33,7.4 -40,10.328,8.3 -60,9.295,8 -80,8.64,7.7 -100,7.29,7.3 -120,5.979,6.9 -140,5.148,6.5 -160,3.805,6 -180,2.718,5.3 -200,1.819,4.5 -220,0.544,3.1 \ No newline at end of file diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.cpg b/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.cpg new file mode 100644 index 0000000..3ad133c --- /dev/null +++ b/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.dbf b/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.dbf new file mode 100644 index 0000000..44b9091 Binary files /dev/null and b/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.dbf differ diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.prj b/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.prj new file mode 100644 index 0000000..10a9c72 --- /dev/null +++ b/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.prj @@ -0,0 +1 @@ +PROJCS["ETRS_1989_Austria_Lambert",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",400000.0],PARAMETER["Central_Meridian",13.3333333333333],PARAMETER["Standard_Parallel_1",49.0],PARAMETER["Standard_Parallel_2",46.0],PARAMETER["Latitude_Of_Origin",47.5],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.shp b/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.shp new file mode 100644 index 0000000..418d6bf Binary files /dev/null and b/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.shp differ diff --git a/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.shx b/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.shx new file mode 100644 index 0000000..ec321e6 Binary files /dev/null and b/debrisframe/data/debrisParabolaChanneled/Inputs/XSECT/crossSection.shx differ diff --git a/debrisframe/debrisframeCfg.ini b/debrisframe/debrisframeCfg.ini index cc45b8f..fd1187a 100644 --- a/debrisframe/debrisframeCfg.ini +++ b/debrisframe/debrisframeCfg.ini @@ -5,7 +5,7 @@ [MAIN] # Path to avalanche directory #TODO: change avalancheDir to debrisDir -avalancheDir = data/debrisGenericTopoTimeDep +avalancheDir = data/debrisParabolaChanneled # number of CPU cores to use for the computation of com1DFA # possible values are: diff --git a/debrisframe/in1Utils/fileUtils.py b/debrisframe/in1Utils/fileUtils.py new file mode 100644 index 0000000..f5e6460 --- /dev/null +++ b/debrisframe/in1Utils/fileUtils.py @@ -0,0 +1,35 @@ +""" +Directory and file handling helper functions +""" + +import pathlib +import shutil +import logging + +log = logging.getLogger("avaframe.debrisframe.in1Utils.fileUtils") + + +def copyHydrToInput(debrisDir): + """ + copy the output file from in2TopoHyd into Inputs/REL + to serve as input data for c1TIF + + Parameters + ----------- + debrisDir: str or pathlib path + path to debris-flow directory + """ + + debrisDir = pathlib.Path(debrisDir) + + outputFile = debrisDir / "Outputs" / "in2TopoHyd" / "initCondHyd.csv" + inputsDir = debrisDir / "Inputs" / "REL" + inputsDir.mkdir(parents=True, exist_ok=True) + + inputsHydFile = inputsDir / outputFile.name + + # TODO: How should we handle if the file already exists? + if inputsHydFile.exists(): + raise FileExistsError(f"File already exists: {inputsHydFile}") + + shutil.copy2(outputFile, inputsHydFile) diff --git a/debrisframe/in2TopoHyd/__init__.py b/debrisframe/in2TopoHyd/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/debrisframe/in2TopoHyd/in2TopoHyd.py b/debrisframe/in2TopoHyd/in2TopoHyd.py new file mode 100644 index 0000000..f0bb819 --- /dev/null +++ b/debrisframe/in2TopoHyd/in2TopoHyd.py @@ -0,0 +1,873 @@ +""" +Get initial conditions for hydrograph +""" + +# Load modules +import pathlib +import math +import logging +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + +# local imports +from avaframe.in1Data import getInput as gI +from avaframe.in3Utils import geoTrans +import avaframe.in2Trans.shpConversion as shpConv +import avaframe.com1DFA.DFAtools as DFAtls +import avaframe.in3Utils.fileHandlerUtils as fU + +# create local logger under avaframe namespace to use its logging configuration +log = logging.getLogger("avaframe.debrisframe.in2TopoHyd") + + +def assignCrossSectionCoords(leveePoints, crossSection): + """ + Assign the nearest cross-section coordinates to the levee points + + Parameters + ----------- + leveePoints: dict + dictionary containing x,y-coordinates of levee points + crossSection: dict + dictionary containing x,y-coordinates of cross-section cells + + Returns + -------- + crossSection: dict + see input parameter + additional levee information about x,y-coordinates, + elevation, distance and indices projected on nearest + cross-section cell + """ + # assign levee points to cross section cells + xcoordLevee = leveePoints["x"] + ycoordLevee = leveePoints["y"] + + idPoint = [] + for x, y in zip(xcoordLevee, ycoordLevee): + diffX = crossSection["x"] - x + diffY = crossSection["y"] - y + dist = np.sqrt(diffX * diffX + diffY * diffY) + id = min(dist) + id = np.where(dist == id) + idPoint.append(id[0][0]) + + xLevee = crossSection["x"][idPoint] + yLevee = crossSection["y"][idPoint] + elevLevee = crossSection["elevation"][idPoint] + sLevee = crossSection["s"][idPoint] + + if np.min(sLevee) == 0.0 or np.max(sLevee) == crossSection["s"][-1]: + message = """At least one levee point lies on the starting or ending point of your release line! + Ensure that you keep a distance of at least one raster cell interior to the edges""" + log.error(message) + raise IndexError(message) + + crossSection["idLevee"] = idPoint + crossSection["xLevee"] = xLevee + crossSection["yLevee"] = yLevee + crossSection["elevLevee"] = elevLevee + crossSection["sLevee"] = sLevee + + return crossSection + + +def assignRasterCoords(cellSize, releaseLine): + """ + Assign nearest raster coordinates to the release line + + Parameters + ----------- + cellSize: int + dem cell size + releaseLine: dict + dictionary including starting and ending point of release line + + Returns + -------- + releaseLine: dict + dictionary including new coordinates for start and ending point + """ + + # assign raster cell coord to starting and ending point of release line + xcoord = [] + ycoord = [] + for x, y in zip(releaseLine["x"], releaseLine["y"]): + dx = x % cellSize + dy = y % cellSize + if dx < cellSize / 2: + xcoord.append(x - x % cellSize) + else: + xcoord.append(x + (cellSize - x % cellSize)) + if dy < cellSize / 2: + ycoord.append(y - y % cellSize) + else: + ycoord.append(y + (cellSize - y % cellSize)) + + releaseLine["xRaster"] = np.array(xcoord) + releaseLine["yRaster"] = np.array(ycoord) + + return releaseLine + + +def computeSubArea(elevation, distance, surfElev, wetCellIdx, idx): + """ + Computes flow subareas at the boundaries of a wetted cross section + Handles the first (0) and last element (-1) of an index array + + Parameters + ----------- + elevation: 1D-array + the elevation of the cross-section cells + distance: 1D-array + the path along (distance) along the cross-section cells + surfElev: float + surface elevation of debris flow + wetCellIdx: 1D-array + indices of wetted cells + idx: int + 0 or -1 + Starting or ending index + + Returns + -------- + subarea: float + flow subarea at the boundary of a wetted cross section + dx: float + horizontal distance that includes the flow subarea + """ + + if idx not in (0, -1): + raise ValueError("idx must be 0 or -1") + + cellIdx = wetCellIdx[idx] + + if idx == 0: + neighborIdx = cellIdx - 1 + elif idx == -1: + neighborIdx = cellIdx + 1 + + ds = abs(distance[neighborIdx] - distance[cellIdx]) + diffElev = elevation[neighborIdx] - elevation[cellIdx] + slope = diffElev / ds + th = surfElev - elevation[cellIdx] + if slope == 0.0: + dx = ds + else: + dx = th / slope + + subArea = 0.5 * th * dx + + return subArea, dx + + +def computeRatingCurve(crossSection, topoHydCfg): + """ + compute relation between given discharge and flow thickness + in a given topographic cross section + + Parameters + ----------- + crossSection: dict + dictionary containing x,y-coordinates of cross-section cells, + the elevation of the cross-section cells, + the path (distance) along the cross-section cells + topoHydCfg: configparser object + configuration settings for the in2TopoHyd-module + + Returns + -------- + ratingCurve: dict + dictionary containing + thickness values/surface elevations and corresponding flow area + minimum elevation of cross section (channel) + min and max values of horizontal distance of flow area + """ + # get the location of the lowest point of the cross section (channel) + # channel is the area between the levee points + idLevee = crossSection["idLevee"] + elevLevee = crossSection["elevLevee"] + elevation = crossSection["elevation"] + minElev = np.min(elevation[min(idLevee) : max(idLevee) + 1]) + # horizontal distance between neighbouring cells + distance = crossSection["s"] + # get elevation steps + dElev = topoHydCfg["GENERAL"].getfloat("dElev") + + # initialize lists + thickness = [] + flowArea = [] + surfaceLevel = [] + xmin = [] + xmax = [] + + # number of iteration steps + numIt = (np.min(elevLevee) - minElev) / dElev + numIt = math.floor(numIt) + + # routine loop: calculate flow area for given flow thickness + for step in range(1, numIt + 1): + # get surface level + thick = step * dElev + surfElev = minElev + thick + # find indices of cells that are equal to or lie below surface level + # search only in area between the two levee points + idx = np.where(elevation[min(idLevee) : max(idLevee) + 1] <= surfElev)[0] + idx = idx + min(idLevee) + + if len(idx) == 0: + message = "Debris-flow surface level lies below the lowest channel elevation point!" + log.error(message) + raise ValueError(message) + + # np.trapz() only calculates the area between the vertices (cell centers) + # if the thickness at the very left and the very right cell is still > 0, + # there are remaining subareas on both sides that have to be considered + areaLeft, dxLeft = computeSubArea( + elevation=elevation, distance=distance, surfElev=surfElev, wetCellIdx=idx, idx=0 + ) + areaRight, dxRight = computeSubArea( + elevation=elevation, distance=distance, surfElev=surfElev, wetCellIdx=idx, idx=-1 + ) + + # get cell elevations and distances between cells + elev = elevation[idx] + dist = distance[idx] + # compute flow area + thickCells = surfElev - elev + # TODO: function np.trapz was changed to np.trapezoid in later numpy versions + flowArea.append(np.trapz(np.maximum(thickCells, 0), dist) + np.sum([areaLeft, areaRight])) + + # save results + thickness.append(thick) + surfaceLevel.append(surfElev) + xmin.append(np.min(dist) - dxLeft) + xmax.append(np.max(dist) + dxRight) + + ratingCurve = { + "thickness": np.array(thickness), + "flowArea": np.array(flowArea), + "minElevation": minElev, + "surfElev": np.array(surfaceLevel), + "xmin": np.array(xmin), + "xmax": np.array(xmax), + } + + return ratingCurve + + +def computeSlope(dem, crossSection): + """ + Computation of the slope of every DEM cell along the release line + between levee points + + Parameters + ----------- + dem: dict + dictionary with dem header and rasterData (numpy nd array of z values) + crossSection: dict + dictionary containing + x,y-coordinates of cross-section cells + indices of cross-section cells along dem-raster + + Returns + -------- + slope: float + mean channel slope [m/m] between levee points + """ + + # the slope has to be in flow direction + # get starting and ending points of cross section + xcoordStart = crossSection["x"][0] + ycoordStart = crossSection["y"][0] + xcoordEnd = crossSection["x"][-1] + ycoordEnd = crossSection["y"][-1] + # get cell size of dem + csz = dem["header"]["cellsize"] + + + # get direction of cross section + dx = xcoordEnd - xcoordStart + dy = ycoordEnd - ycoordStart + + # calculate unit normal vectors on left- and right-handside + # of direction vector + nx1, ny1, _ = DFAtls.normalize(-dy, dx, 0) + nx2, ny2, _ = DFAtls.normalize(dy, -dx, 0) + + # compute two additional cross sections on both handsides of the release line, respectively, + # in distance d of original cross section + d = 2 * np.sqrt(csz**2 + csz**2) + # move starting and ending points of release line in direction of n1 and n2 + xcoordStart1 = xcoordStart + d * nx1 + ycoordStart1 = ycoordStart + d * ny1 + xcoordEnd1 = xcoordEnd + d * nx1 + ycoordEnd1 = ycoordEnd + d * ny1 + xcoordStart2 = xcoordStart + d * nx2 + ycoordStart2 = ycoordStart + d * ny2 + xcoordEnd2 = xcoordEnd + d * nx2 + ycoordEnd2 = ycoordEnd + d * ny2 + # save results in a dictionary + coord1 = {"x": np.array([xcoordStart1, xcoordEnd1]), "y": np.array([ycoordStart1, ycoordEnd1])} + coord2 = {"x": np.array([xcoordStart2, xcoordEnd2]), "y": np.array([ycoordStart2, ycoordEnd2])} + # get elevation for new cross sections + crossSect1 = getCrossSectionCells(dem, coord1) + crossSect2 = getCrossSectionCells(dem, coord2) + # between levee points + start = min(crossSection["idLevee"]) + end = max(crossSection["idLevee"]) + # get mean elevation on each handside + meanElev1 = np.mean(crossSect1["elevation"][start : end + 1]) + meanElev2 = np.mean(crossSect2["elevation"][start : end + 1]) + # get slope as central difference + slope = abs((meanElev1 - meanElev2) / (2 * d)) + + return slope + + +def computeRelFlowThVel(discharge, topoHydCfg, ratingCurve, dem, crossSection): + """ + For a given discharge this function computes flow thickness (starting from the lowest point) + and the mean velocity at a certain cross section. + These flow thicknesses and velocities are used as starting condition. + + Parameters + ----------- + discharge: 1D-array + discharge value + topoHydCfg: configparser object + configuration settings for the in2TopoHyd-module + ratingCurve: dict + dictionary containing + thickness values and corresponding flow area + dem: dict + dictionary with dem header and rasterData (numpy nd array of z values) + crossSection: dict + dictionary containing indices of cross-section cells along dem-raster + + Returns + -------- + thickness: 1D-array + flow thickness + flowVel: 1D-array + flow velocity + """ + # TODO: same for other functions + slope = topoHydCfg["GENERAL"].get("slope", fallback="") + velType = topoHydCfg["GENERAL"].get("velType", fallback="rickenmann") + + # definition of the average slope of the cross section in flow direction + if slope == "": + slope = computeSlope(dem, crossSection) + else: + slope = float(slope) + + log.info(f"chosen slope: {slope:.2f}") + + flowVel = [] + # compute flow velocity + if velType == "rickenmann": + # flow velocity after Rickenmann (1999) + for q in discharge: + v = 2.1 * math.pow(q, 0.33) * math.pow(slope, 0.33) + flowVel.append(v) + + # TODO: add additional methods for calculating the flow velocity + + else: + message = "No velType defined!" + log.error(message) + raise ValueError(message) + + flowVel = np.round(np.array(flowVel), decimals=1) + + # calculate corresponding flow area + flowArea = discharge / flowVel + + # fetch rating curve + thicknessRC = ratingCurve["thickness"] + flowAreaRC = ratingCurve["flowArea"] + # check if the discharge is overtopping the channel + idx = np.where(flowArea > max(flowAreaRC))[0] + if len(idx) != 0: + qOver = min(discharge[idx]) + message = "Discharge of at least %.02f is overtopping the debris-flow channel!" % qOver + log.error(message) + raise ValueError(message) + # interpolate start flow thickness + thickness = np.interp(flowArea, flowAreaRC, thicknessRC) + thickness = np.round(thickness, decimals=2) + + return thickness, flowVel + + +def getCrossSectionCells(dem, releaseLine): + """ + Get cells along the release line to generate a terrain cross section + + Parameters + ----------- + dem: dict + elevation raster data + releaseLine: dict + dictionary including starting and ending point of release line + + Returns + -------- + crossSection: dict + x,y-coordinates and elevation of cells along release line + indices of cross-section cells along dem-raster + """ + + # get origin and cell size of DEM + xllcenter = dem["header"]["xllcenter"] + yllcenter = dem["header"]["yllcenter"] + csz = dem["header"]["cellsize"] + # get elevation data + elevation = dem["rasterData"] + + # assign raster coords to starting and ending point of release line + releaseLine = assignRasterCoords(csz, releaseLine) + xcoordStart = releaseLine["xRaster"][0] + ycoordStart = releaseLine["yRaster"][0] + xcoordEnd = releaseLine["xRaster"][1] + ycoordEnd = releaseLine["yRaster"][1] + + # get direction of release line + dx = xcoordEnd - xcoordStart + dy = ycoordEnd - ycoordStart + + # exception if the release line is a parallel to the + # horizontal or the vertical of the DEM grid + if dx == 0.0: + num = abs(int(dy / csz)) + 1 + crossSectionY = np.linspace(ycoordStart, ycoordEnd, num) + crossSectionX = np.array([xcoordStart] * len(crossSectionY)) + elif dy == 0.0: + num = abs(int(dx / csz)) + 1 + crossSectionX = np.linspace(xcoordStart, xcoordEnd, num) + crossSectionY = np.array([ycoordStart] * len(crossSectionX)) + else: + # get x-coords of raster cells + num = abs(int(dx / csz)) + 1 + crossSectionX = np.linspace(xcoordStart, xcoordEnd, num) + # get distance along x-direction + x = np.cumsum(np.diff(crossSectionX)) + x = np.insert(x, 0, 0) + # get y-coords of raster cells + crossSectionY = dy / dx * x + ycoordStart + # get modulus of division by cell size + modulus = crossSectionY % csz + # round to values that are divisible by cell size + crossSectionY = np.where(modulus < csz / 2, crossSectionY - modulus, crossSectionY + (csz - modulus)) + + # get indices of cross-section cells + col = np.int32(np.round((crossSectionX - xllcenter) / csz)) + row = np.int32(np.round((crossSectionY - yllcenter) / csz)) + # get elevation of cross-section cells + elevCrossSection = elevation[row, col] + + crossSectIdx = np.array([row, col]) + + crossSection = { + "x": crossSectionX, + "y": crossSectionY, + "elevation": elevCrossSection, + "crossSectIdx": crossSectIdx, + } + + return crossSection + + +def getFlowDirection(crossSection, dem): + """ + This function computes the flow direction. + Normal to the cross section. + + Parameters + ----------- + crossSection: dict + dictionary containing x,y-coordinates of cross-section cells, + the elevation of the cross-section cells, + the path (distance) along the cross-section cells + dem: dict + elevation raster data + + + Returns + -------- + flwDir: 1D-array + flow direction including x,y,z-components as unit normal vector of cross section + + """ + + # get starting and ending points of cross section + xcoordStart = crossSection["x"][0] + ycoordStart = crossSection["y"][0] + xcoordEnd = crossSection["x"][-1] + ycoordEnd = crossSection["y"][-1] + # get cell size of dem + csz = dem["header"]["cellsize"] + + # get direction of cross section + dx = xcoordEnd - xcoordStart + dy = ycoordEnd - ycoordStart + + # calculate unit normal vectors on left- and right-handside + # of direction vector; x-y-plane + nx1, ny1, _ = DFAtls.normalize(-dy, dx, 0) + nx2, ny2, _ = DFAtls.normalize(dy, -dx, 0) + + # compute two additional cross sections on both handsides of the release line, respectively, + # in distance d of original cross section + d = 2 * np.sqrt(csz**2 + csz**2) + # move starting and ending points of release line in direction of n1 and n2 + xcoordStart1 = xcoordStart + d * nx1 + ycoordStart1 = ycoordStart + d * ny1 + xcoordEnd1 = xcoordEnd + d * nx1 + ycoordEnd1 = ycoordEnd + d * ny1 + xcoordStart2 = xcoordStart + d * nx2 + ycoordStart2 = ycoordStart + d * ny2 + xcoordEnd2 = xcoordEnd + d * nx2 + ycoordEnd2 = ycoordEnd + d * ny2 + # save results in a dictionary + coord1 = {"x": np.array([xcoordStart1, xcoordEnd1]), "y": np.array([ycoordStart1, ycoordEnd1])} + coord2 = {"x": np.array([xcoordStart2, xcoordEnd2]), "y": np.array([ycoordStart2, ycoordEnd2])} + # get elevation for new cross sections + crossSect1 = getCrossSectionCells(dem, coord1) + crossSect2 = getCrossSectionCells(dem, coord2) + # get mean elevation on each handside + meanElev1 = np.mean(crossSect1["elevation"]) + meanElev2 = np.mean(crossSect2["elevation"]) + # get flow direction in x-y-plane + if meanElev1 > meanElev2: + flwDir = np.array([nx2, ny2]) + else: + flwDir = np.array([nx1, ny1]) + + # get z-component + # get normal vector of the grid mesh + demDict = {"header": dem["header"].copy(), "rasterData": dem["rasterData"]} + demDict = geoTrans.getNormalMesh(demDict) + NxNormed, NyNormed, NzNormed = DFAtls.normalize(demDict["Nx"], demDict["Ny"], demDict["Nz"]) + # compute z-component: 2D-vector must be rotated in 3D-space + # normal vector of grid and normal vector of x-y-plane have to be perpendicular -> n * g = 0 + nz = -(flwDir[0] * NxNormed + flwDir[1] * NyNormed) / NzNormed + + # get normal vector along cross section + row = crossSection["crossSectIdx"][0] + col = crossSection["crossSectIdx"][1] + nz = nz[row, col] + # between levee points + start = min(crossSection["idLevee"]) + end = max(crossSection["idLevee"]) + nz = nz[start : end + 1] + nz = np.mean(nz) + # combine to 3D-unit-normal vector + flwDir = DFAtls.normalize(flwDir[0], flwDir[1], nz) + flwDir = np.array(flwDir) + + return flwDir + + +def assignRelFlowTh(crossSection, ratingCurve, releaseThickness): + """ + This function assigns a release flow thickness to any wet cell + Starting from the lowest point of the channel (releaseThickness), + the corresponding thicknesses are assigned to the other wet cells. + + Parameters + ----------- + crossSection: dict + dictionary containing x,y-coordinates of cross-section cells, + the elevation of the cross-section cells, + the path along (distance) along the cross-section cells + ratingCurve: dict + dictionary containing the minimum elevation of cross section (channel) + releaseThickness: 1D-array + array containing the release thickness values for each discharge value + + Returns + -------- + wetCells: dict + dictionary containing x,y-coordinates and thicknesses for any wet cell + """ + + xCoords = crossSection["x"] + yCoords = crossSection["y"] + elevation = crossSection["elevation"] + idLevee = crossSection["idLevee"] + minElev = ratingCurve["minElevation"] + + # only consider cross section between levee points + elevation = elevation[min(idLevee) : max(idLevee) + 1] + + # identify wet cells + # any cell that lies below debris-flow surface table is considered as wet + thicknessCells = [] + wetXcoords = [] + wetYcoords = [] + + for th in releaseThickness: + # get surface level + surfElev = minElev + th + # find indices of cells that are equal to or lie below surface level + idx = np.where(elevation <= surfElev)[0] + # get cell flow thicknesses + elev = elevation[idx] + thickCells = np.round(surfElev - elev, decimals=2) + thicknessCells.append(thickCells) + # get x,y-coordinates of wet cells + idx = min(idLevee) + idx + wetXcoords.append(xCoords[idx]) + wetYcoords.append(yCoords[idx]) + + wetCells = {"thicknessCells": thicknessCells, "wetXcoords": wetXcoords, "wetYcoords": wetYcoords} + + return wetCells + + +def plotCrossSection(crossSection, outputDir): + """ + This function plots the terrain cross section + + Parameters + ----------- + crossSection: dict + dictionary containing + the elevation of the cross-section cells and levee points, + the path along (distance) along the cross-section cells + outputDir: str or Path + path to output directory of in2TopoHyd module + + Returns + -------- + Plot: .png + saves plot in debris-flow directory + """ + + # plot cross section + fig, ax = plt.subplots(figsize=(10, 8)) + + ax.plot(crossSection["s"], crossSection["elevation"]) + ax.scatter(crossSection["sLevee"], crossSection["elevLevee"], color="red", label="Levee points") + ax.set_xlabel("distance [m]"), ax.set_ylabel("elevation [m]") + ax.grid(color="gray", linestyle="--", linewidth=0.5, alpha=0.6) + ax.set_title("Cross Section") + + plt.legend() + plt.tight_layout() + + # check if directory already exists + path = outputDir / "Plots" + fU.makeADir(path) + + fig.savefig(path / "crossSection.png") + + +def plotRatingCurve(ratingCurve, crossSection, outputDir): + """ + This function plots the terrain cross section + and the corresponding rating curve + + Parameters + ----------- + ratingCurve: dict + dictionary containing + thickness values/surface elevations and corresponding flow area + min and max values of horizontal distance of flow area + crossSection: dict + dictionary containing + the elevation of the cross-section cells, + the path along (distance) along the cross-section cells + outDir: str or Path + path to output directory of in2TopoHyd module + + Returns + -------- + Plot: .png + saves plot in debris-flow directory + """ + + thickness = ratingCurve["thickness"] + flowArea = ratingCurve["flowArea"] + surfElev = ratingCurve["surfElev"] + xmin = ratingCurve["xmin"] + xmax = ratingCurve["xmax"] + + # plot cross section + fig, ax = plt.subplots(ncols=1, nrows=2, figsize=(10, 8)) + + ax[0].plot(crossSection["s"], crossSection["elevation"]) + ax[0].scatter(crossSection["sLevee"], crossSection["elevLevee"], color="red", label="Levee points") + ax[0].hlines( + surfElev, xmin=xmin, xmax=xmax, linestyles="--", colors="grey", lw=0.5, label="elevation increments" + ) + ax[0].set_xlabel("distance [m]"), ax[0].set_ylabel("elevation [m]") + ax[0].grid(color="gray", linestyle="--", linewidth=0.5, alpha=0.6) + ax[0].set_title("Cross Section") + ax[0].legend() + + ax[1].plot(thickness, flowArea) + ax[1].set_xlabel("flow thickness [m]"), ax[1].set_ylabel("flow area [m²]") + ax[1].grid(color="gray", linestyle="--", linewidth=0.5, alpha=0.6) + ax[1].set_title("Rating Curve") + + plt.tight_layout() + + path = outputDir / "Plots" / "ratingCurve.png" + fig.savefig(path) + + +def in2TopoHydMain(debrisDir, topoHydCfg, debrisCfg): + """ + Main script to get the initial conditions for a release line as a csv-file + + Parameters + ----------- + debrisDir: str or pathlib path + path to debris-flow directory + topoHydCfg: configparser.ConfigParser + configuration file for the in2TopoHyd-module + debrisCfg: configparser.ConfigParser + configuration file for the c1TIF-module + + Returns + -------- + csv-file with initial conditions at the release line + Plots of cross section and rating curve + cross-section cell centers as csv-file - optional + + """ + + # create output directory + outputDir = pathlib.Path(debrisDir, "Outputs", "in2TopoHyd") + fU.makeADir(outputDir) + + # +++ 1. read input data + log.info("Read input data") + + # get dem + dem = gI.initializeDEM(debrisDir) + + # get file name of release line + # first, check if name is provided in the c1TIF-config file + inputDir = pathlib.Path(debrisDir, "Inputs") + fname, *_ = gI.getAndCheckInputFiles(inputDir=inputDir, inputType="cross section", folder="XSECT") + + # get release line + releaseLine = shpConv.readLine(fname, "release1", dem) + + log.info(f"release line used: {fname}") + + # fetch number of points + nPoints = len(releaseLine["x"]) + # check if line includes only two points: starting and ending point + if nPoints != 2: + message = "Release line consists of more/less than 2 points! Only starting and ending point allowed!" + log.error(message) + raise ValueError(message) + + # get file name of levee points + inputDir = pathlib.Path(debrisDir, "Inputs") + fname, *_ = gI.getAndCheckInputFiles(inputDir=inputDir, inputType="Levee", folder="LEVEE") + + # get levee points + leveePoints = shpConv.readLine(fname, "release1", dem) + + log.info(f"levee points used: {fname}") + + # +++ 2. get all cells along the release line + log.info("Get all cells along the release line") + + # get cross section cells + crossSection = getCrossSectionCells(dem, releaseLine) + # get distance between cells + crossSection = geoTrans.computeS(crossSection) + + # assign levee points to neares cross-section coordinates + crossSection = assignCrossSectionCoords(leveePoints, crossSection) + + # export cross section cell centers as points for plausibility check + if topoHydCfg["EXPORTS"].getboolean("exportCrossSectionCells"): + file = pd.DataFrame( + {"x": crossSection["x"], "y": crossSection["y"], "elev": crossSection["elevation"]} + ) + path = outputDir / "crossSectionCells.csv" + file.to_csv(path, sep=",", decimal=".", header=True, index=False) + + plotCrossSection(crossSection=crossSection, outputDir=outputDir) + + # +++ 4. calculate hydraulic boundary conditions + log.info("Calculate hydraulic boundary conditions") + + # compute rating curve + ratingCurve = computeRatingCurve(crossSection, topoHydCfg) + + # plot cross section and rating curve for plausibility check + plotRatingCurve(ratingCurve=ratingCurve, crossSection=crossSection, outputDir=outputDir) + + # compute release flow thicknesses and velocities + # get file name of hydrograph + fname, *_ = gI.getAndCheckInputFiles( + inputDir=inputDir, folder="HYDR", inputType="Hydrograph", fileExt="csv" + ) + # read hydrograph + hydrograph = pd.read_csv(fname, sep=",", decimal=".", header=0) + discharge = np.array(hydrograph["discharge"]) + timestep = np.array(hydrograph["timestep"]) + + log.info(f"hydrograph used: {fname}") + + # get release flow thicknesses and velocities + relTh, vel = computeRelFlowThVel( + discharge, topoHydCfg, ratingCurve=ratingCurve, dem=dem, crossSection=crossSection + ) + + # distribute mean flow thickness over wetted cells + wetCells = assignRelFlowTh(crossSection, ratingCurve, relTh) + + # +++ 5. get flow direction + log.info("Get flow direction") + + # normal to release line + # Vx, Vy, Vz -> vel magnitude multiplied by unit vector of direction + flwDirection = getFlowDirection(crossSection=crossSection, dem=dem) + vx = np.round(vel * flwDirection[0], decimals=2) + vy = np.round(vel * flwDirection[1], decimals=2) + vz = np.round(vel * flwDirection[2], decimals=2) + + # +++ 6. export x,y-coords, timesteps, thickness and x,y,z(?)-velocities to csv-file + log.info("Export csv-file") + + time = [] + th = [] + velx = [] + vely = [] + velz = [] + x = [] + y = [] + for t, i in enumerate(timestep): + for j in range(len(wetCells["wetXcoords"][t])): + time.append(i) + th.append(wetCells["thicknessCells"][t][j]) + velx.append(vx[t]) + vely.append(vy[t]) + velz.append(vz[t]) + x.append(wetCells["wetXcoords"][t][j]) + y.append(wetCells["wetYcoords"][t][j]) + + output = pd.DataFrame( + { + "timestep": time, + "thickness": th, + "velocityX": velx, + "velocityY": vely, + "velocityZ": velz, + "x": x, + "y": y, + } + ) + + output.to_csv(outputDir / "initCondHyd.csv", sep=",", decimal=".", header=True, index=False) diff --git a/debrisframe/in2TopoHyd/in2TopoHydCfg.ini b/debrisframe/in2TopoHyd/in2TopoHydCfg.ini new file mode 100644 index 0000000..5743a6d --- /dev/null +++ b/debrisframe/in2TopoHyd/in2TopoHydCfg.ini @@ -0,0 +1,21 @@ +### Config File - This file contains the main settings for the in2TopoHyd module run +## Copy to local_in2TopoHydCfg.ini and set your parameters +# This file is part of DebrisFrame + +[GENERAL] + +#++++++++++++++++ Hydraulic Boundary Conditions +# Elevation steps in vertical direction to calclules the flow area +dElev = 0.1 +# method how the flow velocity as initial condition is calculated +# currently only the method after Rickenmann (1999) is available +velType = rickenmann +# mean slope for the cross section in flow direction [m/m] +# optional; if left blank, the slope is computed automatically +slope = + +[EXPORTS] +# export cell centers along the release line as points in a csv-file; for plausability check +# True or False +exportCrossSectionCells = True + diff --git a/debrisframe/runC1TIF.py b/debrisframe/runC1TIF.py index 93bd998..9e78fd8 100644 --- a/debrisframe/runC1TIF.py +++ b/debrisframe/runC1TIF.py @@ -10,21 +10,27 @@ # import config and init tools from avaframe.in3Utils import cfgUtils from avaframe.in3Utils import logUtils +from avaframe.in3Utils import cfgHandling import avaframe.in3Utils.initializeProject as initProj from avaframe.in3Utils import fileHandlerUtils as fU # import computation modules import debrisframe as debf from debrisframe.c1TIF import c1TIF +from debrisframe.in2TopoHyd import in2TopoHyd +from debrisframe.in1Utils import fileUtils -def runC1TIF(debrisDir=""): +def runC1TIF(debrisDir="", inHydr=False): """Run com1DFA with debris flow parameters with only an avalanche/ debris flow directory as input Parameters ---------- debrisDir: str path to debris flow directory (setup e.g. with init scripts) + inHydr: bool + if inHydr is True, the initial conditions for c1TIF + are computed from a hydrograph first by executing in2TopoHyd Returns ------- @@ -61,6 +67,37 @@ def runC1TIF(debrisDir=""): # load debris flow config DebrisCfg = cfgUtils.getModuleConfig(c1TIF) + # --------------------- + # check if in2TopoHyd computes input data for c1TIF + # TODO: should we also override the flag, in Expert mode? + if inHydr: + DebrisCfg["GENERAL"]["inputHydrograph"] = "True" + else: + inHydr = DebrisCfg["GENERAL"].getboolean("inputHydrograph") + + if inHydr: + # TODO: put this in an separate function? -> where? + if DebrisCfg["com1DFA_com1DFA_override"].getboolean("timeDependentRelease") is False: + message = "If input data are computed from hydrograph, timeDependentRelease needs to be set to True." + log.error(message) + raise ValueError(message) + + in2TopoHydCfg = cfgUtils.getModuleConfig( + in2TopoHyd, + fileOverride="", + modInfo=False, + toPrint=False, + onlyDefault=DebrisCfg["in2TopoHyd_in2TopoHyd_override"].getboolean("defaultConfig"), + ) + in2TopoHydCfg, debrisCfg = cfgHandling.applyCfgOverride( + in2TopoHydCfg, DebrisCfg, in2TopoHyd, addModValues=False + ) + + in2TopoHyd.in2TopoHydMain(debrisDir, in2TopoHydCfg, DebrisCfg) + + # copy in2TopoHyd output into Inputs folder for c1TIF + fileUtils.copyHydrToInput(debrisDir) + # perform com1DFA simulation with debris flow settings _, plotDict, reportDictList, _ = c1TIF.c1TIFMain(cfgMain, DebrisCfg) @@ -86,6 +123,12 @@ def runC1TIF(debrisDir=""): default="", help="the avalanche/ debris directory", ) - print(parser) + parser.add_argument( + "-inHydr", + "--inputHydrograph", + action="store_true", + help="If set, input data is computed from a hydrograph. " + + "If omitted, the default/ini configuration is used." + ) args = parser.parse_args() - runC1TIF(str(args.debrisdir)) + runC1TIF(str(args.debrisdir), args.inputHydrograph) diff --git a/debrisframe/runIn2TopoHyd.py b/debrisframe/runIn2TopoHyd.py new file mode 100644 index 0000000..c5df959 --- /dev/null +++ b/debrisframe/runIn2TopoHyd.py @@ -0,0 +1,86 @@ +""" +Run the In2TopoHyd module to get the initial conditions at the release line +""" + +# load modules +import time +import pathlib +import argparse + +# local imports +from avaframe.in3Utils import cfgUtils +from avaframe.in3Utils import logUtils + +# import computation modules +import debrisframe as debf +from debrisframe.c1TIF import c1TIF +from debrisframe.in2TopoHyd import in2TopoHyd + + +def runIn2TopoHyd(debrisDir=""): + """ + Run in2TopoHyd with only a debris flow directory as input + + Parameters + ---------- + debrisDir: str + path to debris flow directory (setup e.g. with init scripts) + + Returns + ------- + csv-file with initial conditions at the release line + Plots of cross section and rating curve + cross-section cell centers as csv-file - optional + """ + + # Time the whole routine + startTime = time.time() + + # log file name; leave empty to use default runLog.log + logName = "runIn2TopoHyd" + + # Load debris flow directory from general configuration file + modPath = pathlib.Path(debf.__file__).resolve().parent + cfgNameFile = modPath / "debrisframeCfg.ini" + cfgMain = cfgUtils.getGeneralConfig(nameFile=cfgNameFile) + if debrisDir != "": + cfgMain["MAIN"]["avalancheDir"] = debrisDir + # TODO: change avalancheDir to debrisDir + else: + debrisDir = cfgMain["MAIN"]["avalancheDir"] + + # Start logging + log = logUtils.initiateLogger(debrisDir, logName) + log.info("MAIN SCRIPT") + log.info("Current debris flow: %s", debrisDir) + + # load module config + # topoHydCfg + topoHydCfg = cfgUtils.getModuleConfig(in2TopoHyd, debrisDir, toPrint=False) + # c1TIF + debrisCfg = cfgUtils.getModuleConfig(c1TIF, debrisDir, toPrint=False) + + # ---------------- + # Run in2TopoHyd + in2TopoHyd.in2TopoHydMain(debrisDir, topoHydCfg, debrisCfg) + + # Print time needed + endTime = time.time() + log.info("Took %6.1f seconds to calculate." % (endTime - startTime)) + + return + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="Run debris flow workflow") + parser.add_argument( + "debrisdir", + metavar="debrisdir", + type=str, + nargs="?", + default="", + help="the avalanche/ debris directory", + ) + print(parser) + args = parser.parse_args() + runIn2TopoHyd(str(args.debrisdir)) diff --git a/docs/requirements.txt b/docs/requirements.txt index 53fc1f3..5d17f79 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,2 +1,3 @@ sphinx==7.1.2 sphinx-rtd-theme==1.3.0rc1 +sphinxcontrib-bibtex diff --git a/docs/source/_static/css/custom.css b/docs/source/_static/css/custom.css new file mode 100644 index 0000000..014eb0b --- /dev/null +++ b/docs/source/_static/css/custom.css @@ -0,0 +1,18 @@ +.math { + text-align: left; +} +.eqno { + float: right; +} + +/* .wy-nav-content { + max-width: 100%; +} */ + +/* disable horizontal scroll bar for tables */ +.wy-table-responsive table td { + white-space: normal !important; +} +.wy-table-responsive { + overflow: visible !important; +} diff --git a/docs/source/_static/ratingCurve.png b/docs/source/_static/ratingCurve.png new file mode 100644 index 0000000..3533b0b Binary files /dev/null and b/docs/source/_static/ratingCurve.png differ diff --git a/docs/source/conf.py b/docs/source/conf.py index 3a5e885..e011b54 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -18,8 +18,11 @@ "sphinx.ext.autosummary", "sphinx.ext.intersphinx", "sphinx.ext.autosectionlabel", + "sphinxcontrib.bibtex", ] +bibtex_bibfiles = ["references_all.bib"] + intersphinx_mapping = { "python": ("https://docs.python.org/3/", None), "sphinx": ("https://www.sphinx-doc.org/en/master/", None), @@ -45,9 +48,23 @@ # of the sidebar. html_logo = "_static/logo.png" +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ["_static"] + # -- Options for LaTeX output --------------------------------------------- latex_logo = "_static/logo.png" # -- Options for EPUB output epub_show_urls = "footnote" + +# -- Options for referencing ------------------------------------------- +numfig = True +math_numfig = True +math_eqref_format = "Eq.{number}" + + +def setup(app): + app.add_css_file("css/custom.css") \ No newline at end of file diff --git a/docs/source/index.rst b/docs/source/index.rst index a28f681..cf4e2e0 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -44,6 +44,21 @@ Computational modules moduleC1TIF.rst +Input modules +===================== + +.. only:: html + + * :doc:`moduleIn2TopoHyd` + +.. toctree:: + :maxdepth: 3 + :hidden: + :caption: Input modules + + moduleIn2TopoHyd.rst + + Indices and tables ================== @@ -51,4 +66,16 @@ Indices and tables * :ref:`modindex` * :ref:`search` +References +========== + + * :doc:`zreferences` + +.. toctree:: + :maxdepth: 0 + :hidden: + :caption: References + + zreferences.rst + diff --git a/docs/source/moduleC1TIF.rst b/docs/source/moduleC1TIF.rst index ff5b2d9..7b45bea 100644 --- a/docs/source/moduleC1TIF.rst +++ b/docs/source/moduleC1TIF.rst @@ -85,6 +85,14 @@ To run pixi run python runC1TIF.py +Run from hydrograph +-------- + +If in ``local_c1TIFCfg.ini`` ``inputHydrograph`` is set to ``True``, the :py:mod:`c1TIF` input data is computed from +a hydrograph by executing the :ref:`moduleIn2TopoHyd:in2TopoHyd: Initial Hydraulic Conditions Module` module. +From in2TopoHyd :ref:`moduleIn2TopoHyd:Output`, the resulting time dependent release (csv) file is copied into the +``Inputs/REL`` directory for the :py:mod:`c1TIF` simulation. + Theory -------- diff --git a/docs/source/moduleIn2TopoHyd.rst b/docs/source/moduleIn2TopoHyd.rst new file mode 100644 index 0000000..cb3c4f9 --- /dev/null +++ b/docs/source/moduleIn2TopoHyd.rst @@ -0,0 +1,274 @@ +in2TopoHyd: Initial Hydraulic Conditions Module +================================================ + +:py:mod:`in2TopoHyd` is a module for deriving initial hydraulic conditions for the :py:mod:`c1TIF` computation module +for a debris-flow hydrograph at a prescribed release line. The module combines the topography of the release area with +a time-dependent discharge hydrograph (Topo-Hydrograph) to calculate the corresponding flow thickness and flow velocity at the release line. + +The hydraulic conditions are derived from the terrain cross section along the release line. +A rating curve is calculated from the cross-sectional geometry, relating flow thickness to flow area. +For each discharge value of the input hydrograph, the module determines the corresponding mean flow velocity and flow thickness. +The flow thickness is then distributed over the wetted cells of the release cross section. + +The resulting initial conditions contain the flow thickness as well as the velocity components in x, y and z direction for every wet cell and every hydrograph timestep. + +.. Note:: + The module is still under development and the functions are not fully tested yet! + +Theory +------- + +The hydraulic boundary conditions are derived in several steps. + +1. Terrain cross section +^^^^^^^^^^^^^^^^^^^^^ + +The release line constists of a starting and an ending point and is first mapped to the DEM raster grid. +The two release-line points are assigned to the nearest raster-cell coordinates. All DEM cells along the resulting line are then extracted. +For each cross-section cell, the module stores its x and y coordinates and its terrain elevation. The distance along the cross section is calculated and used for the subsequent hydraulic calculations. + +2. Levee points +^^^^^^^^^^^^^^^^ + +Two so called levee points define the left and right top of the channel embankment which constrain the flow to the main channel. +Only this part of the cross section is considered for the calculation of the flow area. + +3. Rating curve +^^^^^^^^^^^^ + +A rating curve is calculated from the terrain cross section between the levee points. +The lowest elevation within this channel section is taken as the minimum channel elevation. Starting from this elevation, the module increases the flow surface elevation in increments of ``dElev``. +For each surface elevation, all cross-section cells below the flow surface are considered wetted. The corresponding flow area is calculated by integrating the flow thickness along the cross section. +It is subsequently used to convert the flow area required for a given discharge into a corresponding flow thickness. + +.. figure:: _static/ratingCurve.png + :width: 90% + + Channel cross section and corresponding rating curve + +4. Flow velocity +^^^^^^^^^^^^ + +.. Note:: + The current implementation provides the approach after :cite:`Rickenmann1999` as the only available method for calculating the mean flow velocity. + +For each discharge value, the mean flow velocity is calculated. +The currently implemented approach :cite:`Rickenmann1999` calculates the veloctiy [m/s] as: + +.. math:: + v = 2.1 \cdot Q^{0.33} \cdot S^{0.33} + :label: velRickenmann + +where :math:`Q` is the discharge [m³/s] and :math:`S` is the mean channel slope [m/m] in flow direction. +The slope can either be supplied through the configuration file via :math:`slope` or calculated automatically. +For the automatic calculation, two additional cross sections are generated in a predefined normal distance on each side of the release line. +The mean elevations between the levee points of this virtual cross sections are used to estimate the slope using a central difference: + +.. math:: + S = \frac{Elev_j - Elev_i}{2d} + :label: slope + +where `d` is the normal distance between the release line and the virtual cross section. + +5. Flow area and flow thickness +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Once the mean velocity is known, the required flow area :math:`A` [m²] is obtained from the mass balance equation: + +.. math:: + A = \frac{Q}{v} + :label: massBalance + +The required flow area is then compared with the previously calculated rating curve and +the flow thickness corresponding to the required flow area is obtained by interpolation of the rating curve. +If the required flow area exceeds the maximum flow area represented by the rating curve, +the module stops and reports that the discharge is overtopping the debris-flow channel. + +6. Distribution over the release cross section +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The calculated flow thickness represents the thickness measured from the lowest point of the channel. +For each timestep, the corresponding surface elevation :math:`Elev_{surf}` [m] is calculated as: + +.. math:: + Elev_{surf} = Elev_{min} + \text{flow thickness}. + +All cross-section cells with an elevation below this surface are considered wet. For each wet cell, the local flow thickness is calculated as +the difference between the flow surface elevation and the terrain elevation. +Consequently, the resulting initial condition is spatially distributed across the cross section +rather than assigning one constant thickness to every cell. + +7. Flow direction +^^^^^^^^^^^^^^ + +The flow direction is by defintion perpendicular to the release line. + +.. Note:: + It is in the responsibility of the user to create a release line normal to the channel flow direction! + Otherwise correct results cannot be guaranteed. + +In a predefined normal distance the program generates two additional cross sections on each side of the release line (compare to calculation of slope). +The mean elevations of this virtual cross sections are then compared to determine which of the two possible normal directions corresponds to the downslope direction. + +Input +--------- + +The :py:mod:`in2TopoHyd` module requires a + +* **digital elevation model as raster file**, + +* **a time-dependent discharge hydrograph**, + +* **a release line**, + +* **and levee points**. + +:py:mod:`in2TopoHyd` calculations are performed within a process directory, organized with the +folder structure described below. + +.. Note:: + + :: + + NameOfDebrisFlow/ + Inputs/ + DEM raster file + CFGs/ - expert configuration files (optional) + REL/ - release line + LEVEE/ - levee points + Outputs/ + in2TopoHyd/ + initCondHyd.csv + crossSectionCells.csv (optional) + Plots/ + crossSection.png + ratingCurve.png + Work/ + + +Digital elevation model +^^^^^^^^^^^^^^^^^^^^^^^ + +The module uses the provided DEM in ``Inputs`` to obtain the terrain elevation along the release line. +The DEM cell size and raster origin are also used to map the release line to raster cells. + + +Release line +^^^^^^^^^^^^ + +The release line is read from the ``Inputs/REL`` directory. It must contain **exactly two points**, representing the starting and ending point of the release line. +If a ``releaseScenario`` is specified in ``(local_)c1TIFCfg.ini``, this file is used (**with** extension .shp). Otherwise, the module searches the ``REL`` directory for a release file. +If no unique release file can be identified, the module stops with an error. +The release line defines the terrain cross section used for the hydraulic calculations. + +Levee points +^^^^^^^^^^^^ + +The module requires exactly one point shapefile in ``Inputs/LEVEE`` whose filename ends with ``*levee.shp``. +The levee points define the lateral limits of the debris-flow channel. Each levee point is assigned to the nearest DEM cell of the release cross section. + +Discharge hydrograph +^^^^^^^^^^^^^^^^^^^^ + +The time-dependent discharge is read from the hydrograph file (.csv) in ``Inputs/REL``. + +The hydrograph must contain the columns: + +* `timestep` — timestep of the discharge value [s] +* `discharge` — discharge value [m³/s] + +The discharge values are used to derive the hydraulic conditions for every timestep of the hydrograph. + +Model configuration +-------------------- + +The model configuration is read from ``in2TopoHydCfg.ini``. A local copy can be created and modified for individual process directories. + +The available parameters are: + +* ``dElev`` + + Elevation increment used to calculate the cross-sectional flow area. For low discharge values it might be necessary + to decrease ``dElev`` + +* ``velType`` + + Method used to calculate the mean flow velocity. + +* ``slope`` + + Mean slope in flow direction [m/m]. If this parameter is left empty, the slope is calculated automatically from the DEM. + +* ``exportCrossSectionCells`` + + If set to ``True``, the coordinates and elevations of the cross-section cell centers are exported as a CSV file for plausibility checking. + +The default parameter values are stored in ``in2TopoHydCfg.ini``. + +Output +------- + +Initial hydraulic conditions +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The main output is ``initCondHyd.csv``. + +For every wet cross-section cell and hydrograph timestep, the file contains: + +* `timestep` — hydrograph timestep [s] +* `thickness` — local flow thickness [m] +* `velocityX` — x-component of flow velocity [m/s] +* `velocityY` — y-component of flow velocity [m/s] +* `velocityZ` — z-component of flow velocity [m/s] +* `x` — coordinate of the wet cell +* `y` — coordinate of the wet cell + +The file therefore provides spatially distributed hydraulic initial conditions along the release line. + +Cross-section cells +^^^^^^^^^^^^^^^^^^^ + +If ``exportCrossSectionCells = True``, the module additionally writes ``crossSectionCells.csv``. + +The file contains: + +* `x` — coordinate +* `y` — coordinate +* `elev` — DEM elevation [m] + +This output can be used to check whether the release line has been correctly mapped to the DEM and whether the resulting terrain cross section is plausible. + +Plots +^^^^^ + +Two plots are generated for plausibility checks. + +* ``crossSection.png`` + + Shows the terrain elevation along the release cross section and the location of the levee points. + +* ``ratingCurve.png`` + + Contains both the terrain cross section with the calculated flow-surface elevation steps and the resulting relationship between flow thickness and flow area. + +These plots are intended to help verify the geometric and hydraulic calculations. + +To run +------- + +* first go to ``DebrisFrame/debrisframe`` + +* copy ``debrisframeCfg.ini`` to ``local_debrisframeCfg.ini`` and set your desired process directory name + +* create a process directory with the required DEM, release line, levee points and hydrograph + +* copy ``in2TopoHyd/in2TopoHydCfg.ini`` to ``in2TopoHyd/local_in2TopoHydCfg.ini`` and, if desired, change the configuration settings + +* run: + :: + + pixi run python runIn2TopoHyd.py + + +The complete workflow reads the DEM and release information, generates the release cross section, calculates the rating curve, +derives the hydraulic conditions from the discharge hydrograph and writes the resulting initial conditions to ``Outputs/in2TopoHyd/initCondHyd.csv``. diff --git a/docs/source/references_all.bib b/docs/source/references_all.bib index 07d63da..0045765 100644 --- a/docs/source/references_all.bib +++ b/docs/source/references_all.bib @@ -17,4 +17,15 @@ @article{MeMcMuMuBa2023 pages = {2405–2420} } +@article{Rickenmann1999, + title = {Empirical {Relationships} for {Debris} {Flows}}, + volume = {19}, + doi = {10.1023/A:1008064220727}, + number = {1}, + journal = {Natural Hazards}, + author = {Rickenmann, Dieter}, + month = jan, + year = {1999}, + pages = {47--77}, +}