Add bathy without breakwater

swash/config.ini

@@ -9,6 +9,7 @@ hstru=Hstru.dat
 poro=Poro.dat
 psize=Psize.dat
 out=out_data
+out_nb=out_data_nb
 
 [swash]
 input=sws/INPUT.sws

swash/processing/bathy_nb.py

@@ -0,0 +1,94 @@
+import argparse
+import configparser
+import logging
+import pathlib
+import sys
+
+import numpy as np
+import pandas as pd
+
+try:
+    import matplotlib.pyplot as plt
+except ImportError:
+    plt = None
+
+
+parser = argparse.ArgumentParser(description="Pre-process bathymetry")
+parser.add_argument("-v", "--verbose", action="count", default=0)
+args = parser.parse_args()
+
+logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
+log = logging.getLogger("bathy")
+
+log.info("Starting bathymetry pre-processing")
+config = configparser.ConfigParser()
+config.read("config.ini")
+
+root = pathlib.Path(config.get("data", "root"))
+
+log.info(f"Reading input data from '{root}'")
+bathy_hires = np.loadtxt(root.joinpath(config.get("data", "hires")))
+bathy_lores = np.loadtxt(root.joinpath(config.get("data", "bathy")))
+hstru = np.loadtxt(root.joinpath(config.get("data", "hstru")))
+poro = np.loadtxt(root.joinpath(config.get("data", "poro")))
+psize = np.loadtxt(root.joinpath(config.get("data", "psize")))
+
+log.info("Generating grid")
+x_hires = -np.arange(0, 0.5 * bathy_hires.size, 0.5)[::-1]
+x_lores = -np.arange(0, 1 * bathy_lores.size, 1)[::-1]
+x_hstru = -np.arange(0, 0.5 * hstru.size, 0.5)[::-1]
+
+log.info("Generating output data")
+bathy_hires_pd = pd.Series(bathy_hires.copy(), index=x_hires)
+bathy_lores_pd = pd.Series(bathy_lores.copy(), index=x_lores)
+
+bathy = pd.DataFrame(
+    index=bathy_lores_pd.index.union(bathy_hires_pd.index),
+    columns=("z", "hstru", "poro", "psize"),
+    data=0,
+)
+bathy.z[bathy_lores_pd.index] = bathy_lores_pd
+bathy.z[bathy_hires_pd.index] = bathy_hires_pd
+bathy.z = np.minimum(bathy.z, -15)
+
+#bathy.loc[x_hstru, ("hstru", "poro", "psize")] = np.array(
+#    (hstru, poro, psize)
+#).T
+
+bathy = bathy.reindex(bathy_lores_pd.index)
+log.debug(f"Bathymetry:\n{bathy}")
+log.info(
+    f"xmin: {bathy.index.min()}, "
+    f"xmax: {bathy.index.max()}, "
+    f"n: {bathy.index.size}"
+)
+
+if config.has_option("data", "out_nb"):
+    out = pathlib.Path(config.get("data", "out_nb"))
+    log.info(f"Writing output data to '{out}'")
+    out.mkdir(exist_ok=True)
+    np.savetxt(out.joinpath("bathy.dat"), bathy.z, newline=" ")
+    np.savetxt(out.joinpath("hstru.dat"), bathy.hstru, newline=" ")
+    np.savetxt(out.joinpath("poro.dat"), bathy.poro, newline=" ")
+    np.savetxt(out.joinpath("psize.dat"), bathy.psize, newline=" ")
+
+    bathy.to_hdf(out.joinpath("bathy.h5"), "bathy", mode="w")
+
+if config.getboolean("proc", "plot", fallback=False):
+    if plt is None:
+        log.error("Could not import PyPlot")
+        sys.exit(1)
+
+    log.info("Plotting data")
+
+    fig, ax = plt.subplots()
+    ax.plot(x_hires, bathy_hires, label="High-res")
+    ax.plot(x_lores, bathy_lores, label="Low-res")
+    ax.plot(bathy.index, bathy.z, ls="-.", c="k", label="Combined")
+    ax.plot(bathy.index, bathy.z + bathy.hstru, label="Hstru")
+
+    ax.grid()
+    ax.legend()
+    plt.show(block=True)
+
+log.info("Processing finished")