diff --git a/scripts/backfill_event_shape.py b/scripts/backfill_event_shape.py index 975c716..4b7c308 100644 --- a/scripts/backfill_event_shape.py +++ b/scripts/backfill_event_shape.py @@ -1,12 +1,12 @@ #!/usr/bin/env python3 -"""Backfill events.shape_* from each event's .h5 waveform samples. Idempotent.""" +"""Backfill events.shape_* and shape_offset_* from each event's .h5 samples. Idempotent.""" from __future__ import annotations import argparse, logging, sys from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) from sfm.database import SeismoDb from sfm.waveform_store import WaveformStore -from sfm.shape_metrics import shape_from_h5 +from sfm.shape_metrics import shape_from_h5, offset_from_h5 log = logging.getLogger("backfill_event_shape") @@ -21,6 +21,7 @@ def backfill_shape(db: SeismoDb, store: WaveformStore, *, dry_run: bool = False) if not h5_path.exists(): counts["skipped_no_h5"] += 1; continue shape = shape_from_h5(h5_path) + offset = offset_from_h5(h5_path) if shape is None: # The .h5 can no longer yield a shape (fewer than 2 samples, or a # flat trace). Clear any previously stored value rather than @@ -28,22 +29,32 @@ def backfill_shape(db: SeismoDb, store: WaveformStore, *, dry_run: bool = False) # from and silently feeds the false-trigger detector. Seen after # a decoder fix shrinks an event: 493 rows in the prod snapshot # were carrying metrics from a superseded decode (2026-08-25). - if row.get("shape_crest_factor") is not None: + if (row.get("shape_crest_factor") is not None + or row.get("shape_offset") is not None): if not dry_run: with db._connect() as conn: conn.execute( "UPDATE events SET shape_crest_factor=NULL, " "shape_near_peak_count=NULL, shape_sample_count=NULL, " - "shape_axis=NULL WHERE id=?", (row["id"],)) + "shape_axis=NULL, shape_offset=NULL, shape_offset_axis=NULL, " + "shape_offset_pre=NULL, shape_offset_spread=NULL WHERE id=?", + (row["id"],)) counts["cleared_stale"] += 1 counts["skipped_no_samples"] += 1; continue if not dry_run: with db._connect() as conn: conn.execute( "UPDATE events SET shape_crest_factor=?, shape_near_peak_count=?, " - "shape_sample_count=?, shape_axis=? WHERE id=?", + "shape_sample_count=?, shape_axis=?, shape_offset=?, " + "shape_offset_axis=?, shape_offset_pre=?, shape_offset_spread=? " + "WHERE id=?", (shape["crest_factor"], shape["near_peak_count"], - shape["sample_count"], shape["axis"], row["id"])) + shape["sample_count"], shape["axis"], + (1 if offset["offset"] else 0) if offset else None, + offset["axis"] if offset else None, + offset["pre"] if offset else None, + offset["spread"] if offset else None, + row["id"])) counts["updated"] += 1 log.info("backfill_shape: %s", counts) return counts diff --git a/sfm/database.py b/sfm/database.py index 04c5d24..1096bbe 100644 --- a/sfm/database.py +++ b/sfm/database.py @@ -99,6 +99,10 @@ CREATE TABLE IF NOT EXISTS events ( shape_near_peak_count INTEGER, -- samples >= 0.5 * peak (FT: few; real: many) shape_sample_count INTEGER, -- total samples (to normalize near_peak_count) shape_axis TEXT, -- geophone channel measured ("Tran"/"Vert"/"Long") + shape_offset INTEGER, -- 1 = DC-offset false trigger (pre-trigger baseline off zero + flat). Meaningful for waveforms only. + shape_offset_axis TEXT, -- geo channel the offset was measured on + shape_offset_pre REAL, -- pre-trigger baseline median (in/s) + shape_offset_spread REAL, -- max(pre,mid,end) - min(...) in in/s; small = constant/DC created_at TEXT NOT NULL DEFAULT (strftime('%Y-%m-%dT%H:%M:%SZ', 'now')), UNIQUE(serial, timestamp) ); @@ -225,6 +229,10 @@ class SeismoDb: ("shape_near_peak_count", "INTEGER"), ("shape_sample_count", "INTEGER"), ("shape_axis", "TEXT"), + ("shape_offset", "INTEGER"), + ("shape_offset_axis", "TEXT"), + ("shape_offset_pre", "REAL"), + ("shape_offset_spread", "REAL"), ("reviewed_real", "INTEGER NOT NULL DEFAULT 0"), ): if col not in existing_cols: @@ -430,9 +438,11 @@ class SeismoDb: tran_zc_above_range, vert_zc_above_range, long_zc_above_range, mic_zc_above_range, shape_crest_factor, shape_near_peak_count, - shape_sample_count, shape_axis) + shape_sample_count, shape_axis, + shape_offset, shape_offset_axis, + shape_offset_pre, shape_offset_spread) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, - ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) + ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) """, ( self._new_id(), serial, key, session_id, ts, @@ -464,6 +474,10 @@ class SeismoDb: rec.get("shape_near_peak_count"), rec.get("shape_sample_count"), rec.get("shape_axis"), + rec.get("shape_offset"), + rec.get("shape_offset_axis"), + rec.get("shape_offset_pre"), + rec.get("shape_offset_spread"), ), ) inserted += 1 @@ -517,7 +531,11 @@ class SeismoDb: shape_crest_factor = COALESCE(?, shape_crest_factor), shape_near_peak_count = COALESCE(?, shape_near_peak_count), shape_sample_count = COALESCE(?, shape_sample_count), - shape_axis = COALESCE(?, shape_axis) + shape_axis = COALESCE(?, shape_axis), + shape_offset = COALESCE(?, shape_offset), + shape_offset_axis = COALESCE(?, shape_offset_axis), + shape_offset_pre = COALESCE(?, shape_offset_pre), + shape_offset_spread = COALESCE(?, shape_offset_spread) WHERE serial = ? AND timestamp = ? """, ( @@ -549,6 +567,10 @@ class SeismoDb: rec.get("shape_near_peak_count") if rec else None, rec.get("shape_sample_count") if rec else None, rec.get("shape_axis") if rec else None, + rec.get("shape_offset") if rec else None, + rec.get("shape_offset_axis") if rec else None, + rec.get("shape_offset_pre") if rec else None, + rec.get("shape_offset_spread") if rec else None, serial, ts, ), diff --git a/sfm/shape_metrics.py b/sfm/shape_metrics.py index 05c9dd5..cd02b80 100644 --- a/sfm/shape_metrics.py +++ b/sfm/shape_metrics.py @@ -47,6 +47,60 @@ def shape_from_samples(chans: dict) -> dict | None: return s +# ── Offset (DC-baseline) detection ──────────────────────────────────────────── +# A DC offset is a false trigger where the geophone baseline sits at a constant +# non-zero floor (sensor bumped / settled / drifted) instead of oscillating +# around zero. Brian's method (validated in scratch/offset_scan3.py): the +# pre-trigger window is definitionally quiet, so a true offset shows |pre| off +# zero AND stays flat across the record (pre ≈ mid ≈ end). A transient moves one +# third relative to the others and is rejected by the spread test. +# Thresholds are in in/s (the .h5 samples are already range-scaled); validated at +# Normal range (10 in/s) — the only range in the fleet. +OFFSET_FLOOR = 0.025 # |pre| at/above this reads as an off-zero baseline (5 A/D counts) +OFFSET_MAX_SPREAD = 0.02 # max(pre,mid,end) - min(...) at/below this reads as flat/constant + + +def _channel_offset(x, pretrig_n): + """Return (pre, spread, is_offset) for one channel, or None if unusable.""" + x = np.asarray(x, dtype=float) + n = x.size + if n < 3: + return None + t = n // 3 + pre = x[:pretrig_n] if (pretrig_n and 0 < pretrig_n < n) else x[:t] + mid, end = x[t:2 * t], x[2 * t:] + if pre.size == 0 or mid.size == 0 or end.size == 0: + return None + vals = [float(np.median(seg)) for seg in (pre, mid, end)] + spread = max(vals) - min(vals) + is_offset = abs(vals[0]) >= OFFSET_FLOOR and spread <= OFFSET_MAX_SPREAD + return vals[0], spread, is_offset + + +def offset_from_samples(chans: dict, pretrig_n) -> dict | None: + """Detect a DC-offset false trigger across the geophone channels. + + An event is offset if ANY geo channel's pre-trigger baseline is off zero and + flat across the record. Reports the tripping axis (or, if none trips, the + most-offset-like axis) with its ``pre``/``spread`` for transparency + tuning. + Returns None when no geo channel is usable. + """ + results = [] + for ax in _GEO_CHANNELS: + x = chans.get(ax) + if x is None: + continue + r = _channel_offset(x, pretrig_n) + if r is not None: + results.append((ax, r[0], r[1], r[2])) + if not results: + return None + offenders = [r for r in results if r[3]] + ax, pre, spread, _ = max(offenders or results, key=lambda r: abs(r[1])) + return {"offset": bool(offenders), "axis": ax, + "pre": round(pre, 6), "spread": round(spread, 6)} + + def shape_from_h5(path) -> dict | None: import h5py try: @@ -56,3 +110,18 @@ def shape_from_h5(path) -> dict | None: except Exception: return None return shape_from_samples(chans) + + +def offset_from_h5(path) -> dict | None: + """offset_from_samples fed from an event's .h5 (float32 in/s geo samples + + the pretrig_samples attribute).""" + import h5py + try: + with h5py.File(path, "r") as f: + chans = {ax: f[f"samples/{ax}"][:] for ax in _GEO_CHANNELS + if f"samples/{ax}" in f} + pretrig_n = f.attrs.get("pretrig_samples") + except Exception: + return None + pretrig_n = int(pretrig_n) if pretrig_n is not None else 0 + return offset_from_samples(chans, pretrig_n) diff --git a/sfm/waveform_store.py b/sfm/waveform_store.py index 6b25e14..44d7020 100644 --- a/sfm/waveform_store.py +++ b/sfm/waveform_store.py @@ -41,7 +41,7 @@ from minimateplus.blastware_file import blastware_filename, write_blastware_file from minimateplus.framing import S3Frame from minimateplus.models import Event from sfm import event_hdf5 -from sfm.shape_metrics import shape_from_h5 +from sfm.shape_metrics import shape_from_h5, offset_from_h5 log = logging.getLogger("sfm.waveform_store") @@ -270,6 +270,13 @@ class WaveformStore: "shape_sample_count": _shape["sample_count"], "shape_axis": _shape["axis"], } if _shape else {} + _offset = offset_from_h5(hdf5_path) if hdf5_filename else None + _offset_rec = { + "shape_offset": 1 if _offset["offset"] else 0, + "shape_offset_axis": _offset["axis"], + "shape_offset_pre": _offset["pre"], + "shape_offset_spread": _offset["spread"], + } if _offset else {} return { "filename": filename, "filesize": filesize, @@ -278,6 +285,7 @@ class WaveformStore: "hdf5_filename": hdf5_filename, "sidecar_filename": sidecar_path.name, **_shape_rec, + **_offset_rec, } def save_imported_bw( @@ -461,6 +469,13 @@ class WaveformStore: "shape_sample_count": _shape["sample_count"], "shape_axis": _shape["axis"], } if _shape else {} + _offset = offset_from_h5(hdf5_path) if hdf5_filename else None + _offset_rec = { + "shape_offset": 1 if _offset["offset"] else 0, + "shape_offset_axis": _offset["axis"], + "shape_offset_pre": _offset["pre"], + "shape_offset_spread": _offset["spread"], + } if _offset else {} return ev, { "filename": filename, "filesize": filesize, @@ -470,6 +485,7 @@ class WaveformStore: "sidecar_filename": sidecar_path.name, "serial": serial, **_shape_rec, + **_offset_rec, } def save_imported_idf( @@ -751,6 +767,13 @@ class WaveformStore: "shape_sample_count": _shape["sample_count"], "shape_axis": _shape["axis"], } if _shape else {} + _offset = offset_from_h5(hdf5_path) if hdf5_filename else None + _offset_rec = { + "shape_offset": 1 if _offset["offset"] else 0, + "shape_offset_axis": _offset["axis"], + "shape_offset_pre": _offset["pre"], + "shape_offset_spread": _offset["spread"], + } if _offset else {} return ev, { "filename": filename, "filesize": filesize, @@ -760,6 +783,7 @@ class WaveformStore: "sidecar_filename": sidecar_path.name, "serial": serial, **_shape_rec, + **_offset_rec, } def load_a5(self, serial: str, filename: str) -> Optional[list[S3Frame]]: diff --git a/tests/test_offset_metrics.py b/tests/test_offset_metrics.py new file mode 100644 index 0000000..e349a82 --- /dev/null +++ b/tests/test_offset_metrics.py @@ -0,0 +1,94 @@ +import numpy as np +import h5py +from sfm.shape_metrics import offset_from_samples, offset_from_h5 + + +def test_flags_constant_dc_floor(): + # A geophone channel sitting at a constant +0.05 in/s across the whole record + # is a DC offset: baseline off zero AND flat across pre/mid/end thirds. + n = 300 + chans = {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)} + r = offset_from_samples(chans, pretrig_n=50) + assert r["offset"] is True + assert r["axis"] == "Tran" + assert abs(r["pre"] - 0.05) < 1e-6 + assert r["spread"] < 0.02 + + +def test_transient_rejected_by_spread(): + # Off-zero pre-trigger but the baseline SETTLES back over the record — a + # transient, not a constant offset. The spread test must reject it. + x = np.concatenate([np.full(100, 0.05), np.full(100, 0.025), np.zeros(100)]) + chans = {"Tran": x, "Vert": np.zeros(300), "Long": np.zeros(300)} + r = offset_from_samples(chans, pretrig_n=100) + assert r["offset"] is False + + +def test_clean_oscillation_not_offset(): + t = np.arange(300) + x = 0.4 * np.sin(2 * np.pi * t / 20) # oscillates around zero — baseline IS zero + chans = {"Tran": x, "Vert": np.zeros(300), "Long": np.zeros(300)} + r = offset_from_samples(chans, pretrig_n=50) + assert r["offset"] is False + + +def test_below_floor_not_offset_but_reports_pre(): + # A flat baseline below the floor is not an offset; still report the axis/pre + # for tuning transparency. + n = 300 + chans = {"Tran": np.full(n, 0.01), "Vert": np.zeros(n), "Long": np.zeros(n)} + r = offset_from_samples(chans, pretrig_n=50) + assert r["offset"] is False + assert r["axis"] == "Tran" + assert abs(r["pre"] - 0.01) < 1e-6 + + +def test_none_when_no_geo_channels(): + assert offset_from_samples({"MicL": np.full(300, 0.05)}, pretrig_n=50) is None + + +def test_pretrig_fallback_when_invalid(): + # pretrig_n of 0 (missing/unusable) falls back to the first third. + n = 300 + chans = {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)} + r = offset_from_samples(chans, pretrig_n=0) + assert r["offset"] is True + + +def test_flags_offset_on_any_axis(): + # Offset on Vert alone still flags the event, and Vert is reported. + n = 300 + chans = {"Tran": np.zeros(n), "Vert": np.full(n, -0.06), "Long": np.zeros(n)} + r = offset_from_samples(chans, pretrig_n=50) + assert r["offset"] is True + assert r["axis"] == "Vert" + + +def _write_h5(path, chans, pretrig_n): + with h5py.File(path, "w") as f: + g = f.create_group("samples") + for k, v in chans.items(): + g.create_dataset(k, data=np.asarray(v, dtype="float32")) + if pretrig_n is not None: + f.attrs["pretrig_samples"] = pretrig_n + + +def test_offset_from_h5_reads_pretrig_attr(tmp_path): + p = tmp_path / "ev.h5" + n = 300 + _write_h5(p, {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)}, + pretrig_n=50) + r = offset_from_h5(str(p)) + assert r["offset"] is True and r["axis"] == "Tran" + + +def test_offset_from_h5_missing_pretrig_attr_falls_back(tmp_path): + p = tmp_path / "noattr.h5" + n = 300 + _write_h5(p, {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)}, + pretrig_n=None) + assert offset_from_h5(str(p))["offset"] is True # falls back to first-third + + +def test_offset_from_h5_missing_file_is_none(tmp_path): + assert offset_from_h5(str(tmp_path / "nope.h5")) is None diff --git a/tests/test_offset_storage.py b/tests/test_offset_storage.py new file mode 100644 index 0000000..73cf069 --- /dev/null +++ b/tests/test_offset_storage.py @@ -0,0 +1,64 @@ +from __future__ import annotations +from pathlib import Path + +import numpy as np, h5py + +from sfm.database import SeismoDb +from sfm.waveform_store import WaveformStore +from scripts.backfill_event_shape import backfill_shape +from minimateplus.models import Event, Timestamp, PeakValues + +_FIX = Path(__file__).parent / "fixtures/histogram-extension-re/events-5-21-26/K558LL8B.7I0W" + + +def _event(waveform_key="0111abcd"): + ev = Event(index=0) + ev._waveform_key = bytes.fromhex(waveform_key) + ev.timestamp = Timestamp(raw=b"", flag=0x10, year=2026, unknown_byte=0, + month=6, day=25, hour=8, minute=50, second=0) + ev.record_type = "Waveform" + ev.peak_values = PeakValues(tran=0.075, vert=0.220, long=0.045, + peak_vector_sum=0.231, micl=0.01) + return ev + + +def test_insert_stores_offset_from_record(tmp_path: Path): + db = SeismoDb(tmp_path / "s.db") + ev = _event() + rec = {ev._waveform_key.hex(): { + "filename": "F.CE0W", "filesize": 10, + "shape_offset": 1, "shape_offset_axis": "Tran", + "shape_offset_pre": 0.05, "shape_offset_spread": 0.001}} + db.insert_events([ev], serial="BE1", waveform_records=rec) + row = db.query_events(serial="BE1")[0] + assert row["shape_offset"] == 1 + assert row["shape_offset_axis"] == "Tran" + assert abs(row["shape_offset_pre"] - 0.05) < 1e-6 + assert abs(row["shape_offset_spread"] - 0.001) < 1e-6 + + +def test_save_imported_bw_attaches_offset(tmp_path: Path): + store = WaveformStore(tmp_path / "waveforms") + ev, rec = store.save_imported_bw(_FIX.read_bytes(), source_path=_FIX, serial_hint="BE9558") + assert rec["shape_offset"] in (0, 1) + assert rec["shape_offset_axis"] in ("Tran", "Vert", "Long") + assert "shape_offset_pre" in rec and "shape_offset_spread" in rec + + +def test_backfill_updates_offset(tmp_path: Path): + db = SeismoDb(tmp_path / "s.db") + store = WaveformStore(tmp_path / "waveforms") + ev = Event(index=0); ev._waveform_key = bytes.fromhex("0111abcd") + db.insert_events([ev], serial="BE1", + waveform_records={ev._waveform_key.hex(): {"filename": "F.CE0W", "filesize": 10}}) + p = store.hdf5_path_for("BE1", "F.CE0W") + with h5py.File(p, "w") as f: + g = f.create_group("samples") + g.create_dataset("Tran", data=np.full(300, 0.05, "float32")) + g.create_dataset("Vert", data=np.zeros(300, "float32")) + g.create_dataset("Long", data=np.zeros(300, "float32")) + f.attrs["pretrig_samples"] = 50 + backfill_shape(db, store) + row = db.query_events(serial="BE1")[0] + assert row["shape_offset"] == 1 + assert row["shape_offset_axis"] == "Tran"