#!/usr/bin/env python3 """ Wave-height analysis of a Cleggan -> Inishbofin ferry crossing. Inputs (as recorded by strixie_mobile_service): eb7e97c7.csv.gz IMU/environment log (sensor,timestamp_ns,v0..v4) timestamp_ns is SensorEvent.timestamp = elapsedRealtimeNanos 9466feee.nmea.gz GNSS NMEA stream Method 1. Resample accelerometer + gyroscope onto a uniform 20 Hz grid. 2. Track the world "up" direction in the device frame with a Mahony complementary filter (gyro integration, slow gravity correction, tau ~ 50 s so that wave-band tilt is not absorbed into the gravity estimate) and project acceleration onto it -> vertical acceleration. 3. Double-integrate in the frequency domain with a band-limited transfer function -> heave displacement; take spectral moments -> Hm0, Tp, Tz. 4. Align the IMU record (boot-relative clock) to the GNSS record (UTC) by cross-correlating magnetometer/accelerometer heading against COG. Outputs: wave_height.png, wave_summary.json """ import gzip, json import numpy as np from scipy.signal import butter, filtfilt, welch, spectrogram HERE = __file__.rsplit('/', 1)[0] IMU_FILE = f'{HERE}/eb7e97c7.csv.gz' NMEA_FILE = f'{HERE}/9466feee.nmea.gz' FS = 20.0 # analysis sample rate (Hz) F_LO, F_HI = 0.07, 0.50 # wave band (2-14 s). Below F_LO the 1/f^4 # integration is noise-dominated; see README. TAU_ATT = 50.0 # gravity-tracking time constant (s) WIN = 240.0 # sliding-window length for Hm0 (s) STEP = 30.0 # ---------------------------------------------------------------- loading def load_imu(path): cols = {'accelerometer': 3, 'gyroscope': 3, 'magnetometer': 3, 'barometer': 1} rows = {k: [] for k in cols} with gzip.open(path, 'rt') as f: f.readline() for line in f: p = line.rstrip('\n').split(',') n = cols.get(p[0]) if n is None: continue rows[p[0]].append([float(p[1])] + [float(x) for x in p[2:2 + n]]) out = {} for k, v in rows.items(): a = np.asarray(v) a[:, 0] /= 1e9 # ns -> s (boot-relative) out[k] = a return out def load_nmea(path): def dm2deg(s, h): if not s: return np.nan i = s.index('.') - 2 v = float(s[:i]) + float(s[i:]) / 60.0 return -v if h in ('S', 'W') else v rmc = [] with gzip.open(path, 'rt', errors='replace') as f: for line in f: if not line.startswith('$GNRMC'): continue p = line.strip().split('*')[0].split(',') if len(p) < 10 or p[2] != 'A': continue t = p[1] sec = int(t[0:2]) * 3600 + int(t[2:4]) * 60 + float(t[4:]) rmc.append([sec, dm2deg(p[3], p[4]), dm2deg(p[5], p[6]), float(p[7]) * 0.514444 if p[7] else np.nan, float(p[8]) if p[8] else np.nan]) return np.asarray(rmc) # utc_s, lat, lon, sog_ms, cog_deg # ------------------------------------------------------- signal utilities def antialias(a, fc=8.0): fs = 1.0 / np.median(np.diff(a[:, 0])) b, c = butter(4, fc / (fs / 2), 'low') o = a.copy() o[:, 1:] = filtfilt(b, c, a[:, 1:], axis=0) return o def resample(a, tg): return np.column_stack([np.interp(tg, a[:, 0], a[:, i]) for i in range(1, a.shape[1])]) def track_up(A, W, fs, tau=TAU_ATT): """Mahony filter: world 'up' unit vector expressed in the device frame. A world-fixed vector seen from the rotating device obeys du/dt = -w x u, so the gyro carries the fast (wave-band) attitude and the accelerometer only trims the slow drift -- at tau=50 s the trim sits an octave below the longest wave of interest, so wave-band roll/pitch is not absorbed into the gravity estimate and does not leak horizontal acceleration into vertical. """ Kp, Ki, dt = 1.0 / tau, 1.0 / (3 * tau * tau), 1.0 / fs u = A[0] / np.linalg.norm(A[0]) bias = np.zeros(3) U = np.empty_like(A) for k in range(len(A)): am = A[k] / np.linalg.norm(A[k]) e = np.cross(u, am) # rotation vector carrying u -> am bias += Ki * e * dt # gyro bias estimate w = W[k] - bias - Kp * e u -= np.cross(w, u) * dt u /= np.linalg.norm(u) U[k] = u return U def accel_to_eta(a, fs, flo=F_LO, fhi=F_HI, order=4): """Band-limited double integration, acceleration -> displacement.""" N = len(a) F = np.fft.rfft(a) f = np.fft.rfftfreq(N, 1 / fs) H = np.zeros_like(f) nz = f > 0 hp = 1 / (1 + (flo / f[nz]) ** (2 * order)) lp = 1 / (1 + (f[nz] / fhi) ** (2 * order)) H[nz] = -hp * lp / (2 * np.pi * f[nz]) ** 2 return np.fft.irfft(F * H, N) def eta_spectrum(a, fs, nperseg): f, Sa = welch(a, fs=fs, nperseg=int(nperseg * fs), noverlap=int(0.75 * nperseg * fs), detrend='linear') S = np.zeros_like(Sa) nz = f > 0 S[nz] = Sa[nz] / (2 * np.pi * f[nz]) ** 4 return f, S def moments(f, S, flo=F_LO, fhi=F_HI): s = (f >= flo) & (f <= fhi) m = {n: np.trapezoid(S[s] * f[s] ** n, f[s]) for n in (0, 1, 2)} fp = f[s][np.argmax(S[s])] return dict(Hm0=4 * np.sqrt(m[0]), m0=m[0], Tp=1 / fp, Tm01=m[0] / m[1], Tm02=np.sqrt(m[0] / m[2])) def zero_upcrossing(eta, t): neg = np.signbit(eta) idx = np.where(~neg[1:] & neg[:-1])[0] H, T, tc = [], [], [] for a, b in zip(idx[:-1], idx[1:]): seg = eta[a:b] H.append(seg.max() - seg.min()) T.append(t[b] - t[a]) tc.append(t[a]) H, T, tc = np.array(H), np.array(T), np.array(tc) ok = (T > 1 / F_HI) & (T < 1 / F_LO) return H[ok], T[ok], tc[ok] def heading_from_imu(A, M, U): """Tilt-compensated heading of the phone's +y axis (top of the device).""" mh = M - np.sum(M * U, axis=1, keepdims=True) * U north = mh / np.linalg.norm(mh, axis=1, keepdims=True) east = np.cross(north, U) # ENU is right-handed: E = N x U east /= np.linalg.norm(east, axis=1, keepdims=True) yd = np.tile([0.0, 1.0, 0.0], (len(U), 1)) yh = yd - np.sum(yd * U, axis=1, keepdims=True) * U yh /= np.linalg.norm(yh, axis=1, keepdims=True) return np.arctan2(np.sum(yh * east, axis=1), np.sum(yh * north, axis=1)) def align_to_gnss(t_imu, imu_heading_rad, rmc): """Offset (s) of the IMU record within the GNSS record, by COG correlation. The phone's top pointed aft, so its heading tracks COG + 180 deg. """ rel = t_imu - t_imu[0] course = np.unwrap(imu_heading_rad + np.pi) ut = rmc[:, 0] - rmc[0, 0] ok = ~np.isnan(rmc[:, 4]) cogu = np.unwrap(np.radians(rmc[ok, 4])) offs = np.arange(0, ut[-1] - rel[-1], 0.5) score = np.array([np.std(np.angle(np.exp( 1j * (course - np.interp(rel + o, ut[ok], cogu))))) for o in offs]) return offs[np.argmin(score)], offs, np.degrees(score) # ------------------------------------------------------------------ main def analyse(): imu = load_imu(IMU_FILE) rmc = load_nmea(NMEA_FILE) acc, gyr, mag = imu['accelerometer'], imu['gyroscope'], imu['magnetometer'] tg = np.arange(max(acc[0, 0], gyr[0, 0]), min(acc[-1, 0], gyr[-1, 0]), 1 / FS) A = resample(antialias(acc), tg) W = resample(antialias(gyr), tg) M = resample(antialias(mag), tg) U = track_up(A, W, FS) g_est = np.sum(A * U, axis=1) az = g_est - g_est.mean() eta = accel_to_eta(az, FS) f, S = eta_spectrum(az, FS, nperseg=256) mom = moments(f, S) H, T, tc = zero_upcrossing(eta, tg) H13 = np.mean(np.sort(H)[-len(H) // 3:]) H110 = np.mean(np.sort(H)[-len(H) // 10:]) offset, offs, score = align_to_gnss(tg, heading_from_imu(A, M, U), rmc) utc0 = rmc[0, 0] + offset # UTC seconds at tg[0] # sliding-window sea state slide = [] for s0 in np.arange(0, tg[-1] - tg[0] - WIN, STEP): seg = az[int(s0 * FS):int((s0 + WIN) * FS)] ff, SS = eta_spectrum(seg, FS, nperseg=120) mm = moments(ff, SS) sel = (tc - tg[0] >= s0) & (tc - tg[0] < s0 + WIN) slide.append([s0 + WIN / 2, mm['Hm0'], mm['Tp'], H[sel].max() if sel.any() else np.nan]) slide = np.asarray(slide) # band energy split bands = [] for lo, hi, name in [(0.07, 0.14, 'swell 7-14 s'), (0.14, 0.25, 'wind sea 4-7 s'), (0.25, 0.50, 'short < 4 s')]: s = (f >= lo) & (f < hi) m0 = np.trapezoid(S[s], f[s]) bands.append((name, lo, hi, m0, 4 * np.sqrt(m0), 100 * m0 / mom['m0'])) return dict(tg=tg, az=az, eta=eta, f=f, S=S, mom=mom, H=H, T=T, tc=tc, H13=H13, H110=H110, slide=slide, bands=bands, rmc=rmc, offset=offset, offs=offs, score=score, utc0=utc0, g_est=g_est, U=U, A=A, W=W) # --------------------------------------------------------------- plotting import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from matplotlib.gridspec import GridSpec from matplotlib.colors import ListedColormap, Normalize from matplotlib.collections import LineCollection SURFACE = '#fcfcfb' INK = '#0b0b0b' INK2 = '#52514e' MUTED = '#8a8880' GRID = '#e6e5e1' S1 = '#2a78d6' # blue S2 = '#eb6834' # orange SPECTRAL_RAMP = ListedColormap(['#2c7bb6', '#218abf', '#00a6ca', '#00b8c3', '#00ccbc', '#48d7ad', '#90eb9d', '#c7f38a', '#ffff8c', '#fbd45b', '#f29e2e', '#e76818', '#d7191c']) BLUE_RAMP = ListedColormap(['#cde2fb', '#b7d3f6', '#9ec5f4', '#86b6ef', '#6da7ec', '#5598e7', '#3987e5', '#2a78d6', '#256abf', '#1c5cab', '#184f95', '#104281', '#0d366b']) def hhmm(sec): sec = sec % 86400 return '%02d:%02d' % (sec // 3600, (sec % 3600) // 60) def style(ax, xlabel=None, ylabel=None, title=None): ax.set_facecolor(SURFACE) for s in ('top', 'right'): ax.spines[s].set_visible(False) for s in ('left', 'bottom'): ax.spines[s].set_color(GRID) ax.tick_params(colors=INK2, labelsize=9, length=3, color=GRID) ax.grid(True, color=GRID, lw=0.7, zorder=0) ax.set_axisbelow(True) if xlabel: ax.set_xlabel(xlabel, color=INK2, fontsize=9.5) if ylabel: ax.set_ylabel(ylabel, color=INK2, fontsize=9.5) if title: ax.set_title(title, color=INK, fontsize=12, fontweight='600', loc='left', pad=10) def make_figure(R, path): tg, eta, slide = R['tg'], R['eta'], R['slide'] trel = tg - tg[0] utc = R['utc0'] + trel mom = R['mom'] fig = plt.figure(figsize=(15.5, 17.5), facecolor=SURFACE) gs = GridSpec(5, 2, height_ratios=[0.52, 2.6, 1.5, 1.5, 1.35], hspace=0.46, wspace=0.20, left=0.062, right=0.965, top=0.915, bottom=0.068) fig.text(0.062, 0.982, 'Sea state on the Cleggan – Inishbofin crossing', fontsize=22, fontweight='700', color=INK, va='top') fig.text(0.062, 0.9655, '18 September 2026 · vertical motion from a hand-held phone IMU, ' '13:14:46–13:34:20 UTC (19.6 min) · 90 Hz accelerometer + 52 Hz gyroscope', fontsize=11, color=INK2, va='top') # ---- hero stat tiles ------------------------------------------------- axh = fig.add_subplot(gs[0, :]); axh.axis('off') tiles = [('Significant wave height', '%.1f m' % mom['Hm0'], 'Hₘ₀ = 4√m₀'), ('Highest single wave', '%.1f m' % R['H'].max(), 'trough to crest'), ('Mean zero-crossing period', '%.1f s' % np.mean(R['T']), 'as encountered'), ('Peak energy period', '%.1f s' % mom['Tp'], 'as encountered'), ('Vertical acceleration', '%.2f g' % (R['az'].std() / 9.81), 'RMS')] for i, (lab, val, sub) in enumerate(tiles): x = i / len(tiles) axh.text(x, 0.72, lab.upper(), fontsize=8.6, color=MUTED, transform=axh.transAxes, fontweight='600') axh.text(x, 0.16, val, fontsize=29, color=INK, transform=axh.transAxes, fontweight='700') axh.text(x, -0.14, sub, fontsize=8.8, color=MUTED, transform=axh.transAxes) # ---- map ------------------------------------------------------------- ax = fig.add_subplot(gs[1, 0]) rmc = R['rmc'] ut = rmc[:, 0] - rmc[0, 0] lat, lon = rmc[:, 1], rmc[:, 2] ax.plot(lon, lat, color='#d8d7d2', lw=5, solid_capstyle='round', zorder=1) # colour the instrumented stretch by windowed Hm0 seg_t = R['offset'] + slide[:, 0] slat = np.interp(seg_t, ut, lat); slon = np.interp(seg_t, ut, lon) pts = np.column_stack([slon, slat]).reshape(-1, 1, 2) segs = np.concatenate([pts[:-1], pts[1:]], axis=1) norm = Normalize(vmin=1.2, vmax=3.0) lc = LineCollection(segs, cmap=SPECTRAL_RAMP, norm=norm, lw=6.5, capstyle='round', zorder=3) lc.set_array(slide[:-1, 1]) ax.add_collection(lc) cb = fig.colorbar(lc, ax=ax, orientation='horizontal', fraction=0.048, pad=0.11, aspect=34) cb.set_label('Hₘ₀ in 4-minute windows (m)', color=INK2, fontsize=9.5) cb.ax.tick_params(colors=INK2, labelsize=9); cb.outline.set_visible(False) ax.plot(lon[0], lat[0], 'o', ms=9, mfc=SURFACE, mec=INK2, mew=2, zorder=4) ax.plot(lon[-1], lat[-1], 'o', ms=9, mfc=INK, mec=SURFACE, mew=2, zorder=4) ax.annotate('Cleggan ' + hhmm(rmc[0, 0]), (lon[0], lat[0]), textcoords='offset points', xytext=(-12, 6), fontsize=10, color=INK2, ha='right') ax.annotate('Inishbofin ' + hhmm(rmc[-1, 0]), (lon[-1], lat[-1]), textcoords='offset points', xytext=(14, -2), fontsize=10, color=INK2, ha='left') ax.annotate('IMU recording starts', (slon[0], slat[0]), textcoords='offset points', xytext=(10, 10), fontsize=9, color=MUTED) ax.set_aspect(1 / np.cos(np.radians(lat.mean()))) style(ax, 'longitude (°E)', 'latitude (°N)', 'Where the motion was measured') ax.margins(0.10) ax.text(0.015, 0.02, 'grey: full GNSS track \u00b7 colour: IMU coverage', transform=ax.transAxes, fontsize=9, color=MUTED, va='bottom') # ---- spectrum -------------------------------------------------------- ax = fig.add_subplot(gs[1, 1]) f, S = R['f'], R['S'] m = f > 0.012 ax.axvspan(0.012, F_LO, color='#f0efec', zorder=0) ax.text(0.019, 0.62, 'noise-dominated\n(1/f⁴ integration)', fontsize=8.6, color=MUTED, transform=ax.get_xaxis_transform(), va='top') bandcol = {'swell 7-14 s': '#eaf2fd', 'wind sea 4-7 s': '#f4f9ff'} for name, lo, hi, _m0, hb, pct in R['bands']: if name in bandcol: ax.axvspan(lo, hi, color=bandcol[name], zorder=0) ax.text(np.sqrt(lo * hi), 2.6e-3, '%s\n%.1f m \u00b7 %.0f%%' % (name.split()[0] if 'short' not in name else 'short', hb, pct), fontsize=9, color=INK2, ha='center', va='bottom', linespacing=1.5) ax.plot(f[m], S[m], color=S1, lw=2, zorder=3) ax.axvline(1 / mom['Tp'], color=S2, lw=2, ls=(0, (4, 3)), zorder=4) ax.annotate('peak %.1f s' % mom['Tp'], (1 / mom['Tp'], 60), textcoords='offset points', xytext=(9, 0), fontsize=9.5, color=S2, fontweight='600') ax.set_xscale('log'); ax.set_yscale('log') ax.set_xlim(0.012, 1.0); ax.set_ylim(1e-3, 3e3) style(ax, 'encounter frequency (Hz)', 'heave energy density (m²/Hz)', 'Wave energy spectrum') sec = ax.secondary_xaxis('top', functions=(lambda x: 1 / np.maximum(x, 1e-9), lambda x: 1 / np.maximum(x, 1e-9))) sec.set_xlabel('period (s)', color=INK2, fontsize=9.5) sec.set_xticks([60, 30, 14, 7, 4, 2, 1]) sec.xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(lambda v, _: '%g' % v)) sec.xaxis.set_minor_locator(matplotlib.ticker.NullLocator()) sec.tick_params(colors=INK2, labelsize=9); sec.spines['top'].set_color(GRID) ax.text(0.985, 0.955, 'shaded: band used for Hₘ₀ (0.07–0.50 Hz)', transform=ax.transAxes, fontsize=9, color=MUTED, ha='right') # ---- heave time series ---------------------------------------------- ax = fig.add_subplot(gs[2, :]) ax.fill_between(trel, eta, 0, color=S1, alpha=0.18, lw=0, zorder=2) ax.plot(trel, eta, color=S1, lw=1.0, zorder=3) ax.axhline(0, color=MUTED, lw=0.8, zorder=2) imax = np.argmax(R['H']); tmax = R['tc'][imax] - tg[0] ax.annotate('largest wave %.1f m' % R['H'][imax], (tmax, eta[int(tmax * FS)]), textcoords='offset points', xytext=(0, 62), ha='center', fontsize=10, color=S2, fontweight='600', arrowprops=dict(arrowstyle='-', color=S2, lw=1.5, shrinkB=6)) for y, lab in ((mom['Hm0'] / 2, None), (-mom['Hm0'] / 2, None)): ax.axhline(y, color=S2, lw=1.2, ls=(0, (5, 4)), zorder=4) ax.text(0.998, mom['Hm0'] / 2 + 0.1, '±Hₘ₀/2', color=S2, fontsize=9.5, va='bottom', ha='right', fontweight='600', transform=ax.get_yaxis_transform()) ax.set_xlim(0, trel[-1]); ax.set_ylim(-2.6, 2.6) style(ax, 'UTC', 'heave (m)', 'Vertical motion, band-passed to 2–14 s') ticks = np.arange(0, trel[-1], 120) ax.set_xticks(ticks); ax.set_xticklabels([hhmm(R['utc0'] + t) for t in ticks]) # ---- Hm0 over time --------------------------------------------------- ax = fig.add_subplot(gs[3, :]) ax.scatter(R['tc'] - tg[0], R['H'], s=13, color=MUTED, alpha=0.55, lw=0, zorder=2, label='individual waves (zero up-crossing)') ax.plot(slide[:, 0], slide[:, 3], color=S2, lw=2, zorder=3, label='largest wave per 4-min window') ax.plot(slide[:, 0], slide[:, 1], color=S1, lw=2.6, zorder=4, label='Hₘ₀, 4-min window') ax.set_xlim(0, trel[-1]); ax.set_ylim(0, 4.6) style(ax, 'UTC', 'wave height (m)', 'Sea state through the crossing') ax.set_xticks(ticks); ax.set_xticklabels([hhmm(R['utc0'] + t) for t in ticks]) lg = ax.legend(loc='lower right', bbox_to_anchor=(1, 1.004), frameon=False, fontsize=9.5, ncol=3, handlelength=1.6) for t in lg.get_texts(): t.set_color(INK2) # ---- spectrogram ----------------------------------------------------- ax = fig.add_subplot(gs[4, :]) ff, tt, Sxx = spectrogram(R['az'], fs=FS, nperseg=int(150 * FS), noverlap=int(135 * FS), detrend='linear') nz = ff > 0 Se = Sxx[nz] / (2 * np.pi * ff[nz])[:, None] ** 4 fmask = (ff[nz] >= 0.05) & (ff[nz] <= 0.5) pm = ax.pcolormesh(tt, ff[nz][fmask], np.log10(Se[fmask] + 1e-4), cmap=BLUE_RAMP, shading='gouraud', vmin=-1.6, vmax=0.9, zorder=2) cb = fig.colorbar(pm, ax=ax, fraction=0.028, pad=0.012) cb.set_label('log₁₀ energy density (m²/Hz)', color=INK2, fontsize=9.5) cb.ax.tick_params(colors=INK2, labelsize=9); cb.outline.set_visible(False) ax.set_xlim(tt[0], tt[-1]) style(ax, 'UTC', 'encounter frequency (Hz)', 'Where the energy sat, minute by minute') ax.set_xticks(ticks); ax.set_xticklabels([hhmm(R['utc0'] + t) for t in ticks]) ax.grid(False) fig.text(0.062, 0.010, 'Heave = vertical acceleration (gravity removed via a gyro-led attitude filter) ' 'double-integrated in the frequency domain over 0.07–0.50 Hz. ' 'Periods are as encountered by a vessel making 11 kn, not true wave periods.\n' 'The IMU clock is boot-relative; it was aligned to GNSS UTC by correlating ' 'magnetometer heading against course over ground (\u00b115 s \u2248 \u00b185 m along track).', fontsize=8.8, color=MUTED, linespacing=1.7, va='bottom') fig.savefig(path, dpi=130, facecolor=SURFACE) print('wrote', path) if __name__ == '__main__': R = analyse() m = R['mom'] print('IMU offset within GNSS record: %.1f s (UTC start %s)' % (R['offset'], hhmm(R['utc0']))) print('gravity estimate %.4f m/s2' % R['g_est'].mean()) print('Hm0 %.2f m H1/3 %.2f m H1/10 %.2f m Hmax %.2f m' % (m['Hm0'], R['H13'], R['H110'], R['H'].max())) print('Tp %.1f s Tm01 %.1f s Tm02 %.1f s Tz %.1f s (%d waves)' % (m['Tp'], m['Tm01'], m['Tm02'], np.mean(R['T']), len(R['H']))) for n, lo, hi, m0, h, p in R['bands']: print(' %-16s Hm0 %.2f m %2.0f%%' % (n, h, p)) make_figure(R, f'{HERE}/wave_height.png') json.dump(dict( utc_start=hhmm(R['utc0']), duration_s=float(R['tg'][-1] - R['tg'][0]), imu_offset_in_gnss_s=float(R['offset']), Hm0_m=float(m['Hm0']), H13_m=float(R['H13']), H110_m=float(R['H110']), Hmax_m=float(R['H'].max()), Tp_s=float(m['Tp']), Tm01_s=float(m['Tm01']), Tm02_s=float(m['Tm02']), Tz_s=float(np.mean(R['T'])), n_waves=int(len(R['H'])), rms_vertical_accel_g=float(R['az'].std() / 9.81), Hm0_window_min_m=float(R['slide'][:, 1].min()), Hm0_window_max_m=float(R['slide'][:, 1].max()), bands=[dict(name=n, f_lo=lo, f_hi=hi, Hm0_m=float(h), pct_energy=float(p)) for n, lo, hi, _, h, p in R['bands']], band_hz=[F_LO, F_HI], note='heave of a hand-held phone; periods are encounter periods', ), open(f'{HERE}/wave_summary.json', 'w'), indent=2)