internship/data/processing/projection.py

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

D = np.diff(np.stack((artha, buoy)), axis=0)
x = np.arange(
    -150,
    np.sqrt((D**2).sum()) + 150,
    config.getfloat("bathy", "step", fallback=1),
)
theta = np.angle(D.dot((1, 1j)))

coords = artha + (x * np.stack((np.cos(theta), np.sin(theta)))).T

log.info("Interpolating bathymetry in 1D")
z = interpolate.griddata(bathy[:, :2], bathy[:, 2], coords)
log.debug(z)