import argparse import configparser import logging import pathlib import matplotlib.pyplot as plt import numpy as np from scipy import interpolate from .lambert import Lambert 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") bathy_inp = pathlib.Path(config.get("bathy", "sub")) bathy_out = pathlib.Path(config.get("bathy", "out")) 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)) def display(): x = np.linspace(bathy[:, 0].min(), bathy[:, 0].max()) y = np.linspace(bathy[:, 1].min(), bathy[:, 1].max()) X, Y = np.meshgrid(x, y) Z = interpolate.griddata( bathy[:, :2], bathy[:, 2], (X, Y), method="nearest" ) fix, ax = plt.subplots() ax.pcolormesh(X, Y, Z) ax.scatter(*artha, c="k") ax.scatter(*buoy, c="k") ax.axis("equal") return ax D = np.diff(np.stack((artha, buoy)), axis=0) x = np.arange(-150, np.sqrt((D**2).sum()) + 150) theta = np.angle(D.dot((1, 1j))) coords = artha + (x * np.stack((np.cos(theta), np.sin(theta)))).T z = interpolate.griddata(bathy[:,:2], bathy[:,2], coords) ax = display() ax.scatter(*coords.T, c="k", marker=".") fig_1d, ax_1d = plt.subplots() ax_1d.plot(x, z) plt.show(block=True)