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")