1
Fork 0
internship/data/processing/projection.py

64 lines
1.6 KiB
Python
Raw Normal View History

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],
),
2022-03-11 14:11:48 +01:00
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)
2022-03-11 14:11:48 +01:00
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
2022-03-11 14:11:48 +01:00
log.info("Interpolating bathymetry in 1D")
z = interpolate.griddata(bathy[:, :2], bathy[:, 2], coords)
log.debug(z)