Bathymetry processing (subdomain + 1d interpolation)

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

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

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)

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)