Data processing + photo analysis

data/photos.md

@@ -0,0 +1,65 @@
+# Times
+## Block Displacement
+18:28:12
+## Previous
+18:27:39
+## Next
+18:28:28
+## Other Waves
+18:25:01
+18:23:13
+18:22:55
+18:22:00
+18:09:08
+18:08:43
+18:08:27
+18:07:18
+18:04:37
+17:59:08
+17:58:31
+17:54:15
+17:53:55
+17:53:39
+17:47:53
+17:39:40
+17:38:45
+17:38:28
+17:32:06
+17:31:46
+17:26:06
+17:25:12
+17:24:47
+17:23:55
+17:23:18
+
+# Periods
+## Block Displacement
+### Before
+33
+### After
+16
+## Other Waves
+
+18
+
+
+25
+16
+
+
+
+37
+
+20
+16
+
+
+
+17
+
+20
+
+
+25
+
+37

data/processing/spec.py

@@ -45,7 +45,7 @@ if cycle is None:
     f = inp["f"]
     S = inp["S"] * Sm
 else:
-    f = np.arange(inp["f"].min(), inp["f"].max() + 1/cycle, 1/cycle)
+    f = np.arange(inp["f"].min(), inp["f"].max() + 1 / cycle, 1 / cycle)
     S = griddata(inp["f"], inp["S"] * Sm, f)
 
 with out_spec.open("w") as out:

data/processing/ts.py

@@ -25,12 +25,14 @@ out_ts = out_root.joinpath("ts.dat")
 
 raw_ts = []
 for tsi in config.get("inp", "raw_ts").split(","):
-    raw_ts.append(np.loadtxt(
+    raw_ts.append(
-        inp_root.joinpath(tsi),
+        np.loadtxt(
-        dtype=[("state", int), ("z", float), ("y", float), ("x", float)],
+            inp_root.joinpath(tsi),
-        delimiter=",",
+            dtype=[("state", int), ("z", float), ("y", float), ("x", float)],
-        max_rows=2304,
+            delimiter=",",
-    ))
+            max_rows=2304,
+        )
+    )
 n = len(raw_ts)
 raw_ts = np.concatenate(raw_ts)
 log.debug(f"{raw_ts=}")
@@ -39,11 +41,11 @@ 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]
+t = np.linspace(0, 30 * 60 * n, 2304 * n + 1)[:-1]
 log.debug(f"{t=}")
 
 log.info(f"Saving timeseries to '{out_ts}'")
-np.savetxt(out_ts, np.stack((t, raw_ts["z"]/100), axis=1))
+np.savetxt(out_ts, np.stack((t, raw_ts["z"] / 100), axis=1))
 
 fig, ax = plt.subplots()
 ax.plot(t, raw_ts["z"])

data/processing/wavelet.py

@@ -0,0 +1,61 @@
+import argparse
+import configparser
+import logging
+import pathlib
+
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.signal as sgl
+
+parser = argparse.ArgumentParser(description="Pre-process time-series")
+parser.add_argument("-v", "--verbose", action="count", default=0)
+parser.add_argument("-c", "--config", default="config.ini")
+args = parser.parse_args()
+
+logging.basicConfig()
+log = logging.getLogger("bathy")
+log.setLevel(max((10, 20 - 10 * args.verbose)))
+
+log.info("Starting time-series pre-processing")
+config = configparser.ConfigParser()
+config.read(args.config)
+
+inp_root = pathlib.Path(config.get("inp", "root"), "cerema/raw")
+out_root = pathlib.Path(config.get("out", "root"))
+
+raw_ts = []
+for tsi in sorted(inp_root.glob("2017022817*.raw")):
+#for tsi in sorted(inp_root.glob("*.raw")):
+    raw_ts.append(
+        np.loadtxt(
+            tsi,
+            dtype=[("state", int), ("z", float), ("y", float), ("x", float)],
+            delimiter=",",
+            max_rows=2304,
+        )
+    )
+    log.debug(f"Loading <{tsi}>")
+n = len(raw_ts)
+raw_ts = np.concatenate(raw_ts)
+log.debug(f"{raw_ts=}")
+
+# t = np.linspace(0, 30 * 60 * n * 1e3, 2304 * n + 1)[:-1].astype("timedelta64[ms]") + np.datetime64("2017-02-28T00:00")
+t = np.linspace(0, 30 * 60 * n, 2304 * n, endpoint=False)
+
+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)))}")
+    # log.debug(f"{t[raw_ts['state'] != 0]}")
+
+z = raw_ts["z"]
+# z = np.cos(2 * np.pi * 7 * t) + sgl.gausspulse(t - 0.4, fc=2)
+
+M = sgl.cwt(z, sgl.morlet, np.arange(1, 30 / (30 * 60 / 2304)))
+print(M)
+
+fig, ax = plt.subplots()
+c = ax.imshow(M, aspect="auto", cmap="spring", vmin=0)
+ax2 = ax.twinx()
+ax2.plot(z, c="k")
+fig.colorbar(c)
+plt.show()

data/processing/zero_cross.py

@@ -0,0 +1,164 @@
+import argparse
+import configparser
+import logging
+import pathlib
+
+import matplotlib.pyplot as plt
+import matplotlib.dates as mdates
+import numpy as np
+import scipy.signal as sgl
+from scipy import fft
+
+parser = argparse.ArgumentParser(description="Pre-process time-series")
+parser.add_argument("-v", "--verbose", action="count", default=0)
+parser.add_argument("-c", "--config", default="config.ini")
+args = parser.parse_args()
+
+logging.basicConfig()
+log = logging.getLogger("bathy")
+log.setLevel(max((10, 20 - 10 * args.verbose)))
+
+log.info("Starting time-series pre-processing")
+config = configparser.ConfigParser()
+config.read(args.config)
+
+inp_root = pathlib.Path(config.get("inp", "root"), "cerema/raw")
+out_root = pathlib.Path(config.get("out", "root"))
+out_root.mkdir(exist_ok=True)
+
+raw_ts = []
+#for tsi in sorted(inp_root.glob("2017022817*.raw")):
+for tsi in sorted(inp_root.glob("*.raw")):
+    raw_ts.append(
+        np.loadtxt(
+            tsi,
+            dtype=[("state", int), ("z", float), ("y", float), ("x", float)],
+            delimiter=",",
+            max_rows=2304,
+        )
+    )
+    log.debug(f"Loading <{tsi}>")
+n = len(raw_ts)
+raw_ts = np.concatenate(raw_ts)
+log.debug(f"{raw_ts=}")
+
+t0 = np.linspace(0, 30 * 60 * n, 2304 * n, endpoint=False)
+t = (t0 * 1e3).astype("timedelta64[ms]") + np.datetime64("2017-02-28T00:00")
+
+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)))}")
+    # log.debug(f"{t[raw_ts['state'] != 0]}")
+
+sos = sgl.butter(1, 0.2, btype="lowpass", fs=2305 / (30 * 60), output="sos")
+z = sgl.sosfiltfilt(sos, raw_ts["z"]*1e-2)
+cr0 = np.where(np.diff(np.sign(z)))[0]
+
+wave = np.fromiter(
+    (
+        np.max(np.abs(z[cr0[i - 1] : cr0[i]])) + np.max(np.abs(z[cr0[i] : cr0[i + 1]]))
+        for i in range(1, len(cr0) - 1)
+    ),
+    dtype=np.single,
+)
+
+log.debug(f"{wave=}")
+log.debug(f"{t=}")
+
+# plt.plot(t[cr0[1:-1]], wave)
+
+nw = len(wave) / 2
+nlw = np.sum(wave > 12)
+H13 = np.quantile(wave, 2 / 3)
+log.info(f"Number of waves: {nw}")
+log.info(f"Number of waves >m: {nlw}")
+log.info(f"Proportion: {nlw/nw:e}")
+log.info(f"H1/3: {H13}m")
+
+dt = 30 * 60 / 2304
+# Mlims = (int(5 / dt), int(30 / dt))
+N = t.size // 24
+s0 = 2 * dt
+dj = 0.5
+J = 1 / dj * np.log2(N * dt / s0)
+j = np.arange(0, J)
+sj = s0 * 2 ** (j * dj)
+# sj = s0 * np.arange(1, 2 ** (J * dj))
+Mw = sj / dt
+Mlims = sj[[0, -1]]
+M = (np.abs(sgl.cwt(raw_ts["z"]*1e-2, sgl.morlet2, Mw))/np.var(raw_ts["z"]*1e-2))**2
+# M = np.abs(sgl.cwt(z, sgl.morlet, Mw))
+v = np.max(np.abs(M))
+
+fig, ax = plt.subplots()
+# ax2 = ax.twinx()
+# ax.plot(t0, raw_ts["z"], lw=.5, c="k", alpha=.2)
+# ax.plot(t0, z, ls="-.", lw=.25, alpha=.2, c="k")
+st = raw_ts["state"][raw_ts["state"] != 0]
+c = np.asarray(["g", "b", "r"])
+# ax.vlines(t0[raw_ts["state"] != 0], -20, 20, colors=c[np.where(st != 777, st, 0)])
+# ax.set(xlabel="t (s)", ylabel="z (cm)")
+# ax.set(xlim=(17 * 3600 + 20 * 60, 17 * 3600 + 30 * 60))
+ax.grid(c="k", alpha=0.2)
+ax.set(zorder=1, frame_on=False)
+ax.semilogy()
+a = [t0[0], t0[-1], *Mlims]
+# c = ax.imshow(M, extent=a, aspect="auto", cmap="plasma", vmin=0)
+c = ax.contourf(t, sj, M, cmap="Greys", vmin=0, vmax=v)
+fig.colorbar(c)
+
+bigw = np.where(wave > 12)[0]
+ym = 1.1 * np.max(np.abs(z))
+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)
+    # a = [t0[i0], t0[i1], *Mlims]
+    # c = ax2.imshow(M[:, i0:i1], extent=a, aspect="auto", cmap="Spectral", vmin=-v, vmax=+v)
+    c = ax2.contourf(
+        t[i0:i1],
+        sj,
+        M[:, i0:i1],
+        cmap="Greys",
+        vmin=0,
+        levels=[1, 2.5, 5, 10, 20, 40],
+        extend="both",
+    )
+    fig.colorbar(c, ax=ax2, label="NWPS")
+    ax.plot(t[i0:i1], (raw_ts["z"]*1e-2)[i0:i1], c="k", lw=1)
+    #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(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.axvspan(t[cr0[w]], t[cr0[w+2]], color="k", alpha=.1)
+
+    locator = mdates.AutoDateLocator(minticks=3, maxticks=7)
+    formatter = mdates.ConciseDateFormatter(locator)
+    ax.xaxis.set_major_locator(locator)
+    ax.xaxis.set_major_formatter(formatter)
+    ax2.xaxis.set_major_locator(locator)
+    ax2.xaxis.set_major_formatter(formatter)
+    ax2.axes.set_xticklabels([])
+
+    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, 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)
+#ax.plot(*sgl.welch(raw_ts["z"], 1 / dt), c="k", alpha=0.2, label="PSD")
+#ax.plot(1 / sj, N * np.mean(M, axis=1), c="k", label="CWT")
+## ax.grid(color="k", alpha=.2)
+#ax.set(xlabel="T (s)", ylabel="PSD")
+## ax2.set(ylabel="Average Wavelet Transform")
+#ax.set(xlim=1 / Mlims)
+#ax.legend()
+
+plt.show()