From c443947e3c86f42765bd188379c42953abf8361b Mon Sep 17 00:00:00 2001 From: "Edgar P. Burkhart" Date: Wed, 6 Jul 2022 07:53:59 +0200 Subject: [PATCH] Outputs for article --- data/processing/plot.py | 76 +++++++++++++++++++++++++++++++++++ data/processing/projection.py | 2 +- data/processing/zero_cross.py | 16 ++++---- 3 files changed, 86 insertions(+), 8 deletions(-) create mode 100644 data/processing/plot.py diff --git a/data/processing/plot.py b/data/processing/plot.py new file mode 100644 index 0000000..852682a --- /dev/null +++ b/data/processing/plot.py @@ -0,0 +1,76 @@ +import argparse +import configparser +import logging +import pathlib + +import numpy as np +from scipy import interpolate + +import matplotlib.pyplot as plt + +from .lambert import Lambert + +parser = argparse.ArgumentParser(description="Pre-process bathymetry") +parser.add_argument("-v", "--verbose", action="count", default=0) +parser.add_argument("-c", "--config", default="config.ini") +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(args.config) + +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")) +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")) +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], + ), + 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)) + +bathy[:, :2] = bathy[:, :2] - artha + +fig, ax = plt.subplots(figsize=(6 / 2.54, 5 / 2.54), constrained_layout=True, dpi=200) + +c = ax.tricontourf( + bathy[:, 0], + bathy[:, 1], + bathy[:, 2], + cmap="plasma", + levels=np.arange(-30, 10, 5), + extend="both", +) +ax.plot(*(np.stack((artha, buoy)) - artha).T, lw=1, ls="-.", c="k", marker="x") +ax.set(xlim=(bathy[np.argmax(bathy[:, 1]), 0], bathy[np.argmin(bathy[:, 1]), 0])) +ax.set(ylim=(bathy[np.argmin(bathy[:, 0]), 1], bathy[np.argmax(bathy[:, 0]), 1])) +ax.set(xlabel="x (m)", ylabel="y (m)") +fig.colorbar(c, label="z (m)") +ax.set_aspect("equal") +ax.set_rasterization_zorder(1.5) + +fig.savefig("bathy2d.pdf") +plt.show() diff --git a/data/processing/projection.py b/data/processing/projection.py index 2c22998..650a395 100644 --- a/data/processing/projection.py +++ b/data/processing/projection.py @@ -135,7 +135,7 @@ fig, ax = plt.subplots(figsize=(16 / 2.54, 2 / 3 * 10 / 2.54), constrained_layou ax.plot(-x, z, color="k") ax.fill_between(-x, z + hstru, z, color="k", alpha=0.2) #ax.set_title(f"N={z.size-1}, x=[{-x.max()};{-x.min()}]") -ax.set(ylim=(-30, 15)) +ax.set(ylim=(-40, 15)) ax.set(xlabel="x (m)", ylabel="z (m)") ax.autoscale(True, "x", True) ax.grid() diff --git a/data/processing/zero_cross.py b/data/processing/zero_cross.py index ce76b19..e010b34 100644 --- a/data/processing/zero_cross.py +++ b/data/processing/zero_cross.py @@ -121,10 +121,10 @@ if bigw.size > 32: log.warning(f"Number of large waves: {bigw.size}") sys.exit() -for w in bigw: - fig, (ax2, ax) = plt.subplots(2, figsize=(15/2.54, 2/3*10/2.54), constrained_layout=True) - i0 = cr0[w] - int(1200 / dt) - i1 = cr0[w + 2] + int(1200 / dt) +fig, ax_ = plt.subplots(2 * (bigw.size // 2), 2, figsize=(15/2.54, 4/3*10/2.54), constrained_layout=True) +for w, ax2, ax in zip(bigw, ax_[::2].flatten(), ax_[1::2].flatten()): + i0 = cr0[w] - int(400 / dt) + i1 = cr0[w + 2] + int(400 / dt) # a = [t0[i0], t0[i1], *Mlims] # c = ax2.imshow(M[:, i0:i1], extent=a, aspect="auto", cmap="Spectral", vmin=-v, vmax=+v) ws = np.ptp(raw_ts["z"][cr0[w]:cr0[w+2]]) * 1e-2 @@ -143,13 +143,15 @@ for w in bigw: #ax.plot(t[i0:i1], z[i0:i1], c="k", lw=1, alpha=0.2, ls="-.") # ax.vlines(t[raw_ts["state"] != 0], -20, 20, colors=c[np.where(st != 777, st, 0)]) ax.set(xlim=(t[i0], t[i1 - 1]), ylim=(-ym, ym)) + ax2.set(ylim=(2, 200)) ax2.set(ylabel="T (s)") ax2.grid(c="k", alpha=0.2) ax2.semilogy() ax.grid(c="k", alpha=.2) #ax.axhline(0, c="k", alpha=0.2, lw=1, ls="-.") #ax.set(zorder=1, frame_on=False) - ax.set(xlabel="t (s)", ylabel="z (m)") + ax.set_xlabel("t (s)", loc="left") + ax.set_ylabel("z (m)") ax.axvspan(t[cr0[w]], t[cr0[w+2]], color="k", alpha=.1) locator = mdates.AutoDateLocator(minticks=3, maxticks=7) @@ -162,8 +164,8 @@ for w in bigw: ax2.set_rasterization_zorder(1.5) - fig.savefig(out_root.joinpath(f"wavelet{w}.pdf"), dpi=300) - fig.savefig(out_root.joinpath(f"wavelet{w}.png"), dpi=200) +fig.savefig(out_root.joinpath(f"wavelet.pdf"), dpi=300) +fig.savefig(out_root.joinpath(f"wavelet.png"), dpi=200) #fig, ax = plt.subplots(constrained_layout=True) ## ax.plot(fft.rfftfreq(raw_ts["z"].size, dt), np.abs(fft.rfft(raw_ts["z"])), c="k", alpha=.2, lw=1)