diff --git a/data/config.ini b/data/config.ini new file mode 100644 index 0000000..e0e702f --- /dev/null +++ b/data/config.ini @@ -0,0 +1,12 @@ +[bathy] +inp=Database_20220224.xyz +sub=bathy_sub.npy +out=bathy.npy + +[artha] +lat=43.398450 +lon=-1.673097 + +[buoy] +lat=43.408333 +lon=-1.681667 diff --git a/data/processing/lambert.py b/data/processing/lambert.py new file mode 100644 index 0000000..87d1bc4 --- /dev/null +++ b/data/processing/lambert.py @@ -0,0 +1,100 @@ +import numpy as np + + +class Lambert: + def __init__( + self, + X0=1.7e6, + Y0=2.2e6, + lambda0=3, + phi0=43, + phi1=42.25, + phi2=43.75, + a=6378137, + e=0.081819190842622, + ): + self._X0 = X0 + self._Y0 = Y0 + + self._lambda0 = np.radians(lambda0) + self._phi0 = np.radians(phi0) + self._phi1 = np.radians(phi1) + self._phi2 = np.radians(phi2) + + self._a = a + self._e = e + + self._set_n() + self._set_rho0() + + def cartesian(self, lam, phi): + lam = np.radians(lam) + phi = np.radians(phi) + + theta = self._n * (lam - self._lambda0) + rho = self._rho(phi) + rho0 = self._rho(self._phi0) + + X = self._X0 + rho * np.sin(theta) + Y = self._Y0 + rho0 - rho * np.cos(theta) + + return X, Y + + def _set_n(self): + self._n = ( + np.log(np.cos(self._phi2) / np.cos(self._phi1)) + + ( + np.log( + (1 - (self._e * np.sin(self._phi1) ** 2)) + / (1 - (self._e * np.sin(self._phi2) ** 2)) + ) + / 2 + ) + ) / ( + np.log( + ( + np.tan(self._phi1 / 2 + np.pi / 4) + * ( + (1 - self._e * np.sin(self._phi2)) + * (1 + self._e * np.sin(self._phi1)) + ) + ** (self._e / 2) + ) + / ( + np.tan(self._phi2 / 2 + np.pi / 4) + * ( + (1 + self._e * np.sin(self._phi2)) + * (1 - self._e * np.sin(self._phi1)) + ) + ** (self._e / 2) + ) + ) + ) + + def _set_rho0(self): + self._rho0 = ( + self._a + * np.cos(self._phi1) + / (self._n * np.sqrt(1 - self._e**2 * np.sin(self._phi1) ** 2)) + * ( + np.tan(self._phi1 / 2 + np.pi / 4) + * ( + (1 - self._e * np.sin(self._phi1)) + / (1 + self._e * np.sin(self._phi1)) + ) + ** (self._e / 2) + ) + ** self._n + ) + + def _rho(self, phi): + return ( + self._rho0 + * ( + 1 + / np.tan(phi / 2 + np.pi / 4) + * ((1 + self._e * np.sin(phi)) / (1 - self._e * np.sin(phi))) + ** (self._e / 2) + ) + ** self._n + ) diff --git a/data/processing/projection.py b/data/processing/projection.py new file mode 100644 index 0000000..4f63591 --- /dev/null +++ b/data/processing/projection.py @@ -0,0 +1,82 @@ +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) diff --git a/data/processing/subdomain.py b/data/processing/subdomain.py new file mode 100644 index 0000000..88c3d96 --- /dev/null +++ b/data/processing/subdomain.py @@ -0,0 +1,50 @@ +import argparse +import configparser +import logging +import pathlib + +import numpy as np + +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") + +artha = np.array( + (config.getfloat("artha", "lon"), config.getfloat("artha", "lat")) +) +buoy = np.array( + (config.getfloat("buoy", "lon"), config.getfloat("buoy", "lat")) +) +log.debug(f"artha: {artha}") +log.debug(f"buoy: {buoy}") + +domain = np.stack((artha, buoy)) +domain.sort(axis=0) +log.debug(f"domain: {domain}") + +domain[0] -= 0.002 +domain[1] += 0.002 +log.debug(f"domain: {domain}") + +bathy_inp = pathlib.Path(config.get("bathy", "inp")) +bathy_out = pathlib.Path(config.get("bathy", "sub")) +log.info(f"Reading bathymetry from '{bathy_inp}'") +raw_bathy = np.genfromtxt(bathy_inp) +log.debug(f"Initial size: {raw_bathy.shape}") + +bathy = raw_bathy[ + ((raw_bathy[:, :2] > domain[0]) & (raw_bathy[:, :2] < domain[1])).all( + axis=1 + ) +] +log.debug(f"Final size: {bathy.shape}") + +log.info(f"Saving subdomain to 'bathy'") +np.save(bathy_out, bathy)