Merge branch 'master' of ssh://git.edgarpierre.fr:39529/m2cce/internship

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]

data/data/.gitignore

@@ -1,2 +1,3 @@
 *.xyz
 *.raw
+cerema

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

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

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