diff --git a/data/photos.md b/data/photos.md new file mode 100644 index 0000000..5d33b18 --- /dev/null +++ b/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 diff --git a/data/processing/wavelet.py b/data/processing/wavelet.py new file mode 100644 index 0000000..c49b50d --- /dev/null +++ b/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() diff --git a/data/processing/zero_cross.py b/data/processing/zero_cross.py new file mode 100644 index 0000000..e4eb9f5 --- /dev/null +++ b/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()