import argparse import configparser import logging import pathlib import numpy as np from scipy import interpolate import matplotlib.pyplot as plt from .lambert import Lambert parser = argparse.ArgumentParser(description="Pre-process bathymetry") parser.add_argument("-v", "--verbose", action="count", default=0) parser.add_argument("-c", "--config", default="config.ini") 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(args.config) inp_root = pathlib.Path(config.get("inp", "root")) out_root = pathlib.Path(config.get("out", "root")) bathy_inp = out_root.joinpath(config.get("out", "sub")) log.info(f"Loading bathymetry from {bathy_inp}") bathy_curvi = np.load(bathy_inp) projection = Lambert() bathy = np.stack( ( *projection.cartesian(bathy_curvi[:, 0], bathy_curvi[:, 1]), bathy_curvi[:, 2], ), axis=1, ) log.debug(f"Cartesian bathy: {bathy}") artha_curvi = np.array( (config.getfloat("artha", "lon"), config.getfloat("artha", "lat")) ) buoy_curvi = np.array((config.getfloat("buoy", "lon"), config.getfloat("buoy", "lat"))) artha = np.asarray(projection.cartesian(*artha_curvi)) buoy = np.asarray(projection.cartesian(*buoy_curvi)) bathy[:, :2] = bathy[:, :2] - artha fig, ax = plt.subplots(figsize=(6 / 2.54, 5 / 2.54), constrained_layout=True, dpi=200) c = ax.tricontourf( bathy[:, 0], bathy[:, 1], bathy[:, 2], cmap="plasma", levels=np.arange(-30, 10, 5), extend="both", ) ax.plot(*(np.stack((artha, buoy)) - artha).T, lw=1, ls="-.", c="k", marker="x") ax.set(xlim=(bathy[np.argmax(bathy[:, 1]), 0], bathy[np.argmin(bathy[:, 1]), 0])) ax.set(ylim=(bathy[np.argmin(bathy[:, 0]), 1], bathy[np.argmax(bathy[:, 0]), 1])) ax.set(xlabel="x (m)", ylabel="y (m)") fig.colorbar(c, label="z (m)") ax.set_aspect("equal") ax.set_rasterization_zorder(1.5) fig.savefig(out_root.joinpath("bathy2d.pdf")) plt.show()