diff --git a/data/config.ini b/data/config.ini index 955650f..3daa05e 100644 --- a/data/config.ini +++ b/data/config.ini @@ -5,7 +5,7 @@ hires=bathyhires.dat hstru=Hstru.dat poro=Poro.dat psize=Psize.dat -raw_ts=201702281700.raw,201702281730.raw +raw_ts=cerema/raw/201702281700.raw,cerema/raw/201702281730.raw hires_step=0.5 [out] diff --git a/data/data/.gitignore b/data/data/.gitignore index 47293d7..702fc79 100644 --- a/data/data/.gitignore +++ b/data/data/.gitignore @@ -1,2 +1,3 @@ *.xyz *.raw +cerema \ No newline at end of file diff --git a/data/processing/orbitals.py b/data/processing/orbitals.py new file mode 100644 index 0000000..63eafb5 --- /dev/null +++ b/data/processing/orbitals.py @@ -0,0 +1,97 @@ +import argparse +import configparser +import logging +import pathlib + +import matplotlib.pyplot as plt +import numpy as np + +parser = argparse.ArgumentParser(description="Plot orbitals") +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 time-series 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")) + +out_ts = out_root.joinpath("ts.dat") + +raw_ts = [] +for tsi in config.get("inp", "raw_ts").split(","): + raw_ts.append( + np.loadtxt( + inp_root.joinpath(tsi), + dtype=[("state", int), ("z", float), ("y", float), ("x", float)], + delimiter=",", + max_rows=2304, + ) + ) +n = len(raw_ts) +raw_ts = np.concatenate(raw_ts) +log.debug(f"{raw_ts=}") + +if (errs := np.count_nonzero(raw_ts["state"])) != 0: + log.warning(f"{errs} transmission errors!") + log.debug(f"{dict(zip(*np.unique(raw_ts['state'], return_counts=True)))}") + +t = np.linspace(0, 30 * 60 * n, 2304 * n + 1)[:-1] +log.debug(f"{t=}") + +flt = (t > 1370) & (t < 1405) + +figt, axt = plt.subplots(3) +axt[0].plot(t, raw_ts["x"]) +axt[1].plot(t, raw_ts["y"]) +axt[2].plot(t, raw_ts["z"]) +for ax in axt: + ax.axvline(t[flt].min(), c="k") + ax.axvline(t[flt].max(), c="k") + ax.grid() + ax.set(xlim=(t.min(), t.max())) + +ts_flt = raw_ts[flt] +z0 = ts_flt["z"] +figtz, axtz = plt.subplots(3) +axtz[0].plot(t[flt], ts_flt["x"]) +axtz[1].plot(t[flt], ts_flt["y"]) +axtz[2].plot(t[flt], z0) +for ax in axtz: + ax.grid() + ax.set(xlim=(t[flt].min(), t[flt].max())) + +fig3d = plt.figure() +ax3d = fig3d.add_subplot(projection="3d") +ax3d.plot(ts_flt["x"], ts_flt["y"], z0, c="#0066ff") +ax3d.quiver3D( + ts_flt["x"][:-1], + ts_flt["y"][:-1], + z0[:-1], + np.diff(ts_flt["x"])[:], + np.diff(ts_flt["y"])[:], + np.diff(z0)[:], + color="#0066ff", +) +ax3d.set(xlabel="x (cm)", ylabel="y (cm)", zlabel="z (cm)") + +theta = np.angle(raw_ts["x"] + 1j * raw_ts["y"]).mean() +fig2dv, ax2dv = plt.subplots() +x0 = ts_flt["x"] * np.cos(theta) + ts_flt["y"] * np.sin(theta) +ax2dv.plot(x0, z0, c="#0066ff", lw=1) +ax2dv.quiver( + x0[:-1], + z0[:-1], + np.diff(x0)[:], + np.diff(z0)[:], + color="#0066ff", +) +ax2dv.grid() +ax2dv.set(aspect="equal") + +plt.show() diff --git a/swash/processing/r_test.py b/swash/processing/r_test.py new file mode 100644 index 0000000..5a67f1c --- /dev/null +++ b/swash/processing/r_test.py @@ -0,0 +1,49 @@ +import matplotlib.pyplot as plt +import numpy as np +from numpy import random +import scipy.signal as sgl + +yi = random.normal(size=2**20) +yr = np.roll(yi, -(2**10)) + +figy, axy = plt.subplots() +axy.plot(np.arange(2**10, 2**11), yi[2**10 : 2**11]) +axy.plot(np.arange(2**10, 2**11), yr[2**10 : 2**11]) + +figf, axf = plt.subplots() +axf.plot(*sgl.welch(yi)) +axf.plot(*sgl.welch(yr)) + +eta = lambda r: yi + r * yr +u = lambda r: -yi + r * yr + + +def puv(eta, u): + f, phi_eta = sgl.welch(eta) + phi_u = sgl.welch(u)[1] + phi_eta_u = np.abs(sgl.csd(eta, u)[1].real) + + return f, np.sqrt( + (phi_eta + phi_u - 2 * phi_eta_u) / (phi_eta + phi_u + 2 * phi_eta_u) + ) + + +figr, axr = plt.subplots() + +for r in np.arange(0, 1.1, 0.1): + axr.plot(*puv(eta(r), u(r)), c="k") + Rn = puv( + eta(r) + 0.4 * random.normal(size=2**20), + u(r) + 0.4 * random.normal(size=2**20), + ) + axr.plot( + *Rn, + c="#ff6600", + ) + axr.annotate(f"{r=:.1f}", (Rn[0][0], Rn[1][0]), bbox={"boxstyle": "square", "facecolor": "w"}) +axr.grid() +axr.autoscale(True, "x", tight=True) +axr.set(ylim=(0, 1), ylabel="R", xlabel="f") +axr.legend(("No noise", "40% noise"), loc="lower left") + +plt.show() diff --git a/tasks.md b/tasks.md index bfb637a..b2f9c99 100644 --- a/tasks.md +++ b/tasks.md @@ -1,10 +1,10 @@ - Vérification méthodes calcul réflection avec données forcées en Python +* Vérification méthodes calcul réflection avec données forcées en Python Vérification vague incidente possible avec digue et sans digue sur temps long (4h) avec spectre Jonswap -Tracer trajectoires bouées (3d, 2d) autour de vague 15m +* Tracer trajectoires bouées (3d, 2d) autour de vague 15m Swash output into binary matlab files -> input with scipy.io.matlab