diff --git a/swash/config.ini b/swash/config.ini index f815c62..560cfa6 100644 --- a/swash/config.ini +++ b/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 diff --git a/swash/processing/bathy_nb.py b/swash/processing/bathy_nb.py new file mode 100644 index 0000000..dace50a --- /dev/null +++ b/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")