diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index e928523d3..f63b90b6f 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -2209,6 +2209,8 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # setup a result fields info data frame to save max values of fields and avalanche front resultsDF = setupresultsDF(resTypes, cfg["VISUALISATION"].getboolean("createRangeTimeDiagram")) + if cfg["VISUALISATION"].getboolean("outputHydrograph"): + outHydDF, outHydIdx, NzNormed = debF.setupOutHydrographDF(cfgGen, dem) # Add different time stepping options here log.debug("Use standard time stepping") @@ -2310,6 +2312,8 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si else: rangeValue = "" resultsDF = addMaxValuesToDF(resultsDF, fields, t, resTypes, rangeValue=rangeValue) + if cfg["VISUALISATION"].getboolean("outputHydrograph"): + outHydDF = debF.collectOutputHydrograph(outHydDF, outHydIdx, NzNormed, fields, t) tCPU["nSave"] = nSave particles["t"] = t @@ -2490,7 +2494,14 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # save resultsDF to file resultsDFPath = pathlib.Path(cfgGen["avalancheDir"], "Outputs", "com1DFA", "resultsDF_%s.csv" % simHash) resultsDF.to_csv(resultsDFPath) - + if cfg["VISUALISATION"].getboolean("outputHydrograph"): + # save outputHydrographDF to file + outHydPath = pathlib.Path(cfgGen["avalancheDir"], "Outputs", "out1HelpDFHyd") + fU.makeADir(outHydPath) + outHydDF.to_csv(outHydPath / f"outputHydrographDF_{simHash}.csv") + # plot output hydrograph + debF.outputHydrographPlot(outHydDF, cfgGen, simHash=simHash) + if cfg["EXPORTS"].getboolean("exportData"): exportFields(cfg, t, fields, dem, outDir, cuSimName, TSave="final") diff --git a/avaframe/com1DFA/debrisFunctions.py b/avaframe/com1DFA/debrisFunctions.py index f91f1c1e3..1367f5890 100644 --- a/avaframe/com1DFA/debrisFunctions.py +++ b/avaframe/com1DFA/debrisFunctions.py @@ -6,14 +6,19 @@ import logging import numpy as np import copy +import pandas as pd +import pathlib +import matplotlib.pyplot as plt import avaframe.com1DFA.DFAfunctionsCython as DFAfunC import avaframe.com1DFA.particleTools as particleTools import avaframe.in3Utils.geoTrans as geoTrans import avaframe.com1DFA.com1DFA as com1DFA -import avaframe.in2Trans.shpConversion as shpConv from avaframe.in1Data import getInput as gI import avaframe.in3Utils.cfgUtils as cfgUtils +import avaframe.in3Utils.fileHandlerUtils as fU +import avaframe.out3Plot.plotUtils as pU +import avaframe.com1DFA.DFAtools as DFAtls # create local logger # change log level in calling module to DEBUG to see log messages @@ -78,9 +83,7 @@ def initializeTimeDepRelease(cfg, inputSimLines, particles, fields, dem, zPartAr return particles, fields, zPartArray0 -def addReleaseParticles( - cfg, particles, inputSimLines, timeDepRelValues, dem, zPartArray0, timeDepRelIndex -): +def addReleaseParticles(cfg, particles, inputSimLines, timeDepRelValues, dem, zPartArray0, timeDepRelIndex): """ add new particles initialized by a time dependent release to particles that are in the flow already @@ -131,39 +134,38 @@ def addReleaseParticles( # check if already existing particles are within the release polygon # it's possible that there are still a few particles in the polygon with low velocities # TODO: could think of a threshold of number of particles that are still allowed in the polygons? - mask = geoTrans.getParticlesInPolygon(particles, relLine, cfg["GENERAL"].getfloat("thresholdPointInRel")) + mask = geoTrans.getParticlesInPolygon( + particles, relLine, cfg["GENERAL"].getfloat("thresholdPointInRel") + ) if np.sum(mask) > 0: message = ( - "Already existing particles are within the release polygon, which can cause numerical instabilities (at timestep: %02f s)" - % (particles["t"] + particles["dt"]) + "Already existing particles are within the release polygon, which can cause numerical instabilities (at timestep: %02f s)" + % (particles["t"] + particles["dt"]) ) log.error(message) raise ValueError(message) particlesRelease = com1DFA.initializeParticles( - cfg["GENERAL"], - relLine, - dem, - relThField=relThField, - timestep=particles["t"] + particles["dt"] + cfg["GENERAL"], relLine, dem, relThField=relThField, timestep=particles["t"] + particles["dt"] ) if "velocity" in timeDepRelValues and "x" not in timeDepRelValues: particlesRelease = DFAfunC.updateInitialVelocity(cfg["GENERAL"], particlesRelease, dem, velocityMag) - - elif "velocityX" in timeDepRelValues and "x" in timeDepRelValues: for uComp, timedepParameter in zip(["ux", "uy", "uz"], ["velocityX", "velocityY", "velocityZ"]): - raster = gI.timeDepRelCoordsToRaster(relLine["timeDepRelValues"], timeDepRelIndex, - dem["originalHeader"], parameter=timedepParameter - ) + raster = gI.timeDepRelCoordsToRaster( + relLine["timeDepRelValues"], + timeDepRelIndex, + dem["originalHeader"], + parameter=timedepParameter, + ) rasterDict = {"header": dem["header"], "rasterData": raster} - particlesRelease, _ = geoTrans.projectOnRaster(rasterDict, - particlesRelease, outData=uComp) + particlesRelease, _ = geoTrans.projectOnRaster(rasterDict, particlesRelease, outData=uComp) particlesRelease["uMag"] = np.sqrt( - particlesRelease["ux"] ** 2 + particlesRelease["uy"] ** 2 + particlesRelease["uz"] ** 2) + particlesRelease["ux"] ** 2 + particlesRelease["uy"] ** 2 + particlesRelease["uz"] ** 2 + ) particles = particleTools.mergeParticleDict(particles, particlesRelease) # save initial z position for travel angle computation @@ -191,14 +193,11 @@ def defineReleaseLineFromCoordinates(relFile, timeDepRelValues, demHeader): relThFieldData: numpy array release thickness raster data """ - releaseLine = { - "file": relFile, - "Name": [relFile.stem], - "initializedFrom": "csvfile" - } + releaseLine = {"file": relFile, "Name": [relFile.stem], "initializedFrom": "csvfile"} initialIndex = np.where(timeDepRelValues["timeStep"] == 0)[0] - releaseLine["rasterData"] = gI.timeDepRelCoordsToRaster(timeDepRelValues, initialIndex, demHeader, - parameter="thickness") + releaseLine["rasterData"] = gI.timeDepRelCoordsToRaster( + timeDepRelValues, initialIndex, demHeader, parameter="thickness" + ) relThFieldData = releaseLine["rasterData"] # TODO: define thickness for output report, now mean of thickness values releaseLine["thickness"] = np.nanmean(np.where(relThFieldData == 0, np.nan, relThFieldData)) @@ -270,3 +269,166 @@ def checkTravelledDistance(cfgGen, timeDepRelValues, timeDepRelCsv): # TODO: error or warning? log.error(message) raise ValueError(message) + + +def setupOutHydrographDF(cfgGen, demDictOrig): + """ + setup result fields data frame for initial time step for generating an output hydrograph and + reading raster indices and normal vectors of hydrograph lines and distances of cross-section cell centers + computes the z-component of the normal vector of the grid mesh for computing the flow depth + + Parameters + ---------- + cfgGen: configparser + general configuration for DFA simulation + demDictOrig: dict + dictionary with dem header and rasterData (numpy nd array of z values) + + Returns + -------- + outHydrographDF: data frame + data frame with one line for initial time step for discharges + outHydIdx: data frame + data frame with raster indices, normal vectors of output hydrograph cross section lines + distances of cross-section cell centers + NzNormed: 2D-array + field with z-components of the normal vector of the grid mesh + """ + + # read csv file with ouput hydrograph cross section indices + inputDir = pathlib.Path(cfgGen.get("avalancheDir"), "Inputs") + fileType = "Output hydrograph" + folder = "OUTHYDR" + fname, *_ = gI.getAndCheckInputFiles(inputDir, folder, fileType, fileExt="csv") + if fname is None: + message = f"No {fileType} csv-file in Inputs/{folder}!" + log.error(message) + raise FileNotFoundError(message) + outHydIdx = pd.read_csv(fname, sep=",", decimal=".", header=0) + + # create dataframe with result fields + sId = np.unique(np.array(outHydIdx["station"])) + initialTime = {"timestep": [0.0]} + station = {f"XS {i}": [0.0] for i in sId} + outHydrographDF = pd.DataFrame({**initialTime, **station}).set_index("timestep") + + # get z-component of normal vector of the grid mesh + demDict = {"header": demDictOrig["header"].copy(), "rasterData": demDictOrig["rasterData"]} + demDict = geoTrans.getNormalMesh(demDict) + _, _, NzNormed = DFAtls.normalize(demDict["Nx"], demDict["Ny"], demDict["Nz"]) + + return outHydrographDF, outHydIdx, NzNormed + + +def collectOutputHydrograph(outputHydrographDF, outHydIdx, NzNormed, fields, t): + """ + Read discharge in the cross section cells for the current time step. + + Parameters + ---------- + outputHydrographDF: data frame + data frame with one line for initial time step for discharges + outHydIdx: dataframe + data frame with station names, raster indices, normal vectors of output hydrograph cross section lines + and projected cell sizes of cross-section cell centers + NzNormed: 2D-array + field with z-components of the normal vector of the grid mesh + fields: dict + com1DFA fields (FM) + t: float + current computation time + + Returns + -------- + outputHydrographDF: data frame + updated data frame + """ + + station = np.array(outHydIdx["station"]) + sId = np.unique(station) + row = np.array(outHydIdx["row"]) + col = np.array(outHydIdx["column"]) + Nx = np.array(outHydIdx["Nx"]) + Ny = np.array(outHydIdx["Ny"]) + prjCsz= np.array(outHydIdx["prjCsz"]) + + newLineDischarge = [] + + for i in sId: + idx = np.where(station == i) + r = row[idx] + c = col[idx] + # read field data + velX = fields["Vx"][r, c] + velY = fields["Vy"][r, c] + # compute flow depth + depthField = fields["FT"][r, c] / NzNormed[r, c] + # determination of the flow direction with respect to the cross-section line + # there are two normal vectors of the cross section line + # if the scalar product of the mean velocity vector and normal vector is > 0, both show in the same direction + #TODO: define flow direction in the beginning of the simulation and keep it constant -> backflow + nx = Nx[idx] + ny = Ny[idx] + scalarProduct = velX * nx + velY * ny + if not np.mean(scalarProduct) > 0.0: + nx = -nx + ny = -ny + + # get cell size + csz = prjCsz[idx] + # compute discharge per cell and sum up + discharge = (velX * nx + velY * ny) * csz * depthField + newLineDischarge.append(np.sum(discharge)) + # save to data frame + outputHydrographDF.loc[t] = newLineDischarge + + return outputHydrographDF + + +def outputHydrographPlot(outputHydrographDF, cfgGen, simHash=""): + """ + Plot time series of discharge and cumulative volume for every output cross section (one figure per cross section) + + Parameters + ----------- + outputHydrographDF: data frame + data frame with index column 'timestep' and additional columns for every cross-section labeled 'XS_'. + The entries of the cross-section columns include the discharge values. + cfgGen: configparser + general configuration for DFA simulation + simHash: str + simulation hash for the file name + + Returns + -------- + Plots in outDir-directory + """ + + colnames = list(outputHydrographDF.columns) + t = list(outputHydrographDF.index) + dt = np.diff(t) + + outDir = pathlib.Path(cfgGen["avalancheDir"], "Outputs", "out1HelpDFHyd", "Plots") + fU.makeADir(outDir) + + # create figures + for station in colnames: + + discharge = np.array(outputHydrographDF[station]) + dV = 0.5 * (discharge[:-1] + discharge[1:]) * dt + volume = np.cumsum(dV) + + fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(pU.figW, pU.figH)) + ax[0].plot(t, discharge, color="brown") + ax[1].plot(t[1:], volume, color="brown") + + ax[0].set_xlabel("time [s]") + ax[0].set_ylabel("discharge [m³/s]") + ax[1].set_xlabel("time [s]") + ax[1].set_ylabel("cumulative volume [m³]") + ax[1].set_title(f"total cumulative volume = {np.max(volume):.02f} m³") + plt.suptitle(f"Output hydrograph - {station}") + + # save plot + plotName = "outputHydrograph_%s_%s.png" % (simHash, station) + pU.saveAndOrPlot({"pathResult": outDir}, plotName, fig) diff --git a/avaframe/in1Data/getInput.py b/avaframe/in1Data/getInput.py index 5911c70ee..8f7ad1e7a 100644 --- a/avaframe/in1Data/getInput.py +++ b/avaframe/in1Data/getInput.py @@ -350,7 +350,7 @@ def getAndCheckInputFiles(inputDir, folder, inputType, fileExt="shp", fileSuffix If fileExt is empty, only suffix will be checked. - If a file is found, it is checked to be of an allowed extension; currently .shp .asc and .tif are supported. + If a file is found, it is checked to be of an allowed extension; currently .shp .asc, .tif and .csv are supported. Parameters ---------- @@ -415,7 +415,7 @@ def getAndCheckInputFiles(inputDir, folder, inputType, fileExt="shp", fileSuffix if OutputFile.suffix not in supportedFileFormats: message = ( - "Unsupported file format found for OutputFile %s; shp, asc, tif are allowed" % OutputFile + "Unsupported file format found for OutputFile %s; shp, asc, tif, csv are allowed" % OutputFile ) log.error(message) raise AssertionError(message) diff --git a/avaframe/tests/test_debrisFunctions.py b/avaframe/tests/test_debrisFunctions.py index 696ddc1c6..b0a53c098 100644 --- a/avaframe/tests/test_debrisFunctions.py +++ b/avaframe/tests/test_debrisFunctions.py @@ -3,6 +3,8 @@ import numpy as np import pytest import configparser +import pandas as pd +from pathlib import Path import avaframe.com1DFA.debrisFunctions as debF @@ -147,7 +149,9 @@ def test_addReleaseParticles(): particles["y"] = np.array([5, 3, 30]) with pytest.raises(ValueError): - debF.addReleaseParticles(cfg, particles, inputSimLines, timeDepRelValues, dem, zPartArray0, timeDepRelIndex=0) + debF.addReleaseParticles( + cfg, particles, inputSimLines, timeDepRelValues, dem, zPartArray0, timeDepRelIndex=0 + ) def test_addReleaseParticlesFromCoordinateCsv(): @@ -250,6 +254,128 @@ def test_addReleaseParticlesFromCoordinateCsv(): assert np.array_equal(zPartArray0, np.ones(5)) +@pytest.fixture +def hydrographCfg(tmp_path): + """Provide a temporary avalanche directory with a hydrograph input folder.""" + (tmp_path / "Inputs" / "OUTHYDR").mkdir(parents=True) + cfg = configparser.ConfigParser() + cfg["GENERAL"] = {"avalancheDir": str(tmp_path)} + return cfg + + +def test_setupOutHydrographDF(hydrographCfg): + """Read cross sections and initialize one zero-discharge column per station.""" + inputDir = hydrographCfg["GENERAL"]["avalancheDir"] + + csvFile = Path(inputDir) / "Inputs" / "OUTHYDR" / "crossSections.csv" + crossSections = pd.DataFrame( + { + "station": [1, 1, 1, 2, 2, 2], + "row": [397, 396, 395, 397, 396, 395], + "column": [75, 75, 75, 163, 164, 165], + "Nx": [1.0, 1.0, 1.0, 0.76975, 0.76975, 0.76975], + "Ny": [0.0, 0.0, 0.0, 0.63833, 0.63833, 0.63833], + "prjCsz": [0.86705, 1.86705, 1.99998, 1.31306, 2.72134, 2.17712], + } + ) + crossSections.to_csv(csvFile, index=False) + + # A horizontal DEM must have a normalized z component of 1 everywhere. + dem = { + "header": { + "nrows": 751, + "ncols": 1001, + "cellsize": 2.0, + "xllcenter": 999.000000000000, + "yllcenter": -5001.000000000000, + "nodata_value": -9999.0, + }, + "rasterData": np.ones((751, 1001)), + } + + outputDf, outIdx, nzNormed = debF.setupOutHydrographDF(hydrographCfg["GENERAL"], dem) + + pd.testing.assert_frame_equal(outIdx, crossSections) + assert outputDf.index.name == "timestep" + np.testing.assert_allclose(outputDf.index.to_numpy(), [0.0]) + assert list(outputDf.columns) == ["XS 1", "XS 2"] + np.testing.assert_allclose(outputDf.iloc[0].to_numpy(), [0.0, 0.0]) + np.testing.assert_allclose(nzNormed, np.ones((751, 1001))) + assert "Nx" not in dem # Original DEM remains unchanged. + + +def test_setupOutHydrographDFMissingFile(hydrographCfg): + """A missing input CSV must be reported explicitly.""" + dem = {"header": {}, "rasterData": np.zeros((751, 1001))} + + with pytest.raises(FileNotFoundError, match="No Output hydrograph csv-file"): + debF.setupOutHydrographDF(hydrographCfg["GENERAL"], dem) + + +def test_collectOutputHydrograph(): + """Integrate flux over multiple cells and stations, incl. reversed normals.""" + outIdx = pd.DataFrame( + { + "station": [1, 1, 1, 2, 2, 2], + "row": [397, 396, 395, 397, 396, 395], + "column": [75, 75, 75, 163, 164, 165], + "Nx": [1.0, 1.0, 1.0, 0.76975, 0.76975, 0.76975], + "Ny": [0.0, 0.0, 0.0, 0.63833, 0.63833, 0.63833], + "prjCsz": [0.86705, 1.86705, 1.99998, 1.31306, 2.72134, 2.17712], + } + ) + + shape = (751, 1001) + fields = { + "Vx": np.ones(shape), + "Vy": np.zeros(shape), + "FT": np.ones(shape), + } + + nzNormed = np.ones(shape) + outputDf = pd.DataFrame({"XS 1": [0.0], "XS 2": [0.0]}, index=[0.0]) + outputDf.index.name = "timestep" + # q = (Vx * Nx + Vy * Ny) * FT / nzNormed * distance + qXS1 = 1 * 1 * 0.86705 + 1 * 1 * 1.86705 + 1 * 1 * 1.99998 + qXS2 = ( + (1 * 0.76975 + 0 * 0.63833) * 1 * 1.31306 + + (1 * 0.76975 + 0 * 0.63833) * 1 * 2.72134 + + (1 * 0.76975 + 0 * 0.63833) * 1 * 2.17712 + ) + + t = 1.0 + result = debF.collectOutputHydrograph(outputDf, outIdx, nzNormed, fields, t) + + assert result is outputDf + np.testing.assert_allclose(result.index.to_numpy(), [0.0, 1.0]) + assert result.loc[0.0, "XS 1"] == 0.0 + assert result.loc[0.0, "XS 2"] == 0.0 + assert result.loc[1.0, "XS 1"] == pytest.approx(qXS1) + assert result.loc[1.0, "XS 2"] == pytest.approx(qXS2) + + +def test_collectOutputHydrographZeroVelocity(): + """A stationary flow must produce zero discharge even with nonzero thickness.""" + outIdx = pd.DataFrame( + {"station": [1], "row": [0], "column": [0], "Nx": [1.0], "Ny": [0.0], "prjCsz": [2.0]} + ) + + outputDf = pd.DataFrame({"XS 1": [0.0]}, index=[0.0]) + nzNormed = np.ones((1, 1)) + fields = {"Vx": np.zeros((1, 1)), "Vy": np.zeros((1, 1)), "FT": np.ones((1, 1))} + t = 1.0 + + result = debF.collectOutputHydrograph( + outputDf, + outIdx, + nzNormed, + fields, + t, + ) + + assert result.loc[1.0, "XS 1"] == pytest.approx(0.0) + + """ Test does not work because: When calling pytest, executing DFAfunctionsCython.upfateFieldsC() raises an error ("Fatal Python error: Aborted") (see issue #1002?)