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

142 lines
4.1 KiB
Python
Raw Normal View History

import argparse
import configparser
import logging
import pathlib
import numpy as np
from scipy import interpolate
2022-03-17 09:50:21 +01:00
try:
import matplotlib.pyplot as plt
except ImportError:
plt = None
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")
2022-03-11 14:55:16 +01:00
inp_root = pathlib.Path(config.get("inp", "root"))
out_root = pathlib.Path(config.get("out", "root"))
bathy_inp = out_root.joinpath(config.get("out", "sub"))
hires_inp = inp_root.joinpath(config.get("inp", "hires"))
2022-03-17 09:50:21 +01:00
hstru_inp = inp_root.joinpath(config.get("inp", "hstru"))
poro_inp = inp_root.joinpath(config.get("inp", "poro"))
psize_inp = inp_root.joinpath(config.get("inp", "psize"))
2022-03-11 14:55:16 +01:00
bathy_out = inp_root.joinpath(config.get("out", "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(
2022-03-14 13:34:21 +01:00
config.getfloat("out", "left", fallback=0),
np.sqrt((D**2).sum()) + config.getfloat("out", "right", fallback=0),
2022-03-11 14:55:16 +01:00
config.getfloat("out", "step", fallback=1),
2022-03-11 14:11:48 +01:00
)
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(f"z: {z}")
2022-03-11 16:02:36 +01:00
_hires = np.loadtxt(hires_inp)[::-1]
bathy_hires = np.stack(
(
np.linspace(
0,
2022-03-11 14:55:16 +01:00
(_hires.size - 1) * config.getfloat("inp", "hires_step"),
_hires.size,
),
_hires,
),
axis=1,
)
del _hires
log.debug(f"Bathy hires: {bathy_hires}")
z_cr = 5
hires_crossing = np.diff(np.signbit(bathy_hires[:, 1] - z_cr)).nonzero()[0][-1]
log.debug(f"Hires crossing: {hires_crossing}")
z_crossing = np.diff(np.signbit(z - z_cr)).nonzero()[0][-1]
log.debug(f"Z crossing: {z_crossing}")
x_min_hires = x[z_crossing] + (
bathy_hires[:, 0].min() - bathy_hires[hires_crossing, 0]
)
x_max_hires = x[z_crossing] + (
bathy_hires[:, 0].max() - bathy_hires[hires_crossing, 0]
)
log.debug(f"Replacing range: [{x_min_hires},{x_max_hires}]")
flt_x = (x > x_min_hires) & (x < x_max_hires)
z[flt_x] = interpolate.griddata(
(bathy_hires[:, 0],),
bathy_hires[:, 1],
(x[flt_x] - x[z_crossing] + bathy_hires[hires_crossing, 0]),
)
2022-03-11 14:55:16 +01:00
2022-03-17 09:50:21 +01:00
hstru_in = np.loadtxt(hstru_inp)[::-1]
hstru = np.zeros(z.shape)
hstru[flt_x] = interpolate.griddata(
(bathy_hires[:,0],),
hstru_in,
(x[flt_x] - x[z_crossing] + bathy_hires[hires_crossing, 0]),
)
poro_in = np.loadtxt(poro_inp)[::-1]
poro = np.zeros(z.shape)
poro[flt_x] = interpolate.griddata(
(bathy_hires[:,0],),
poro_in,
(x[flt_x] - x[z_crossing] + bathy_hires[hires_crossing, 0]),
)
psize_in = np.loadtxt(psize_inp)[::-1]
psize = np.zeros(z.shape)
psize[flt_x] = interpolate.griddata(
(bathy_hires[:,0],),
psize_in,
(x[flt_x] - x[z_crossing] + bathy_hires[hires_crossing, 0]),
)
2022-03-14 14:16:25 +01:00
np.savetxt(out_root.joinpath("bathy.dat"), z[::-1], newline=" ")
2022-03-17 09:50:21 +01:00
np.savetxt(out_root.joinpath("hstru.dat"), hstru[::-1], newline=" ")
np.savetxt(out_root.joinpath("poro.dat"), poro[::-1], newline=" ")
np.savetxt(out_root.joinpath("psize.dat"), psize[::-1], newline=" ")
if plt is not None and config.getboolean("out", "plot", fallback=False):
fig, ax = plt.subplots()
ax.plot(-x, z, color="k")
ax.fill_between(-x, z+hstru, z, color="k", alpha=.2)
plt.show(block=True)