seismo-relay/sfm/shape_metrics.py

"""Waveform-shape metrics for false-trigger detection.

A false trigger is an isolated impulse (quiet → spike → quiet); a real event
rings for many cycles. Two numbers separate them: crest factor (how far the
peak stands above the typical sample) and how many samples sit near the peak.
"""
from __future__ import annotations
import numpy as np

_GEO_CHANNELS = ("Tran", "Vert", "Long")
NEAR_PEAK_FRACTION = 0.5   # a sample "near the peak" is >= this * peak amplitude


def channel_shape(x) -> dict | None:
    x = np.asarray(x, dtype=float)
    if x.size < 2:
        return None
    peak = float(np.max(np.abs(x)))
    if peak <= 0:
        return None
    rms = float(np.sqrt(np.mean(x ** 2)))
    if rms <= 0:
        return None
    near = int(np.sum(np.abs(x) >= NEAR_PEAK_FRACTION * peak))
    return {"crest_factor": peak / rms, "near_peak_count": near,
            "sample_count": int(x.size)}


def shape_from_samples(chans: dict) -> dict | None:
    best_axis, best_peak, best_x = None, -1.0, None
    for ax in _GEO_CHANNELS:
        x = chans.get(ax)
        if x is None:
            continue
        x = np.asarray(x, dtype=float)
        if x.size < 2:
            continue
        p = float(np.max(np.abs(x)))
        if p > best_peak:
            best_axis, best_peak, best_x = ax, p, x
    if best_axis is None:
        return None
    s = channel_shape(best_x)
    if s is None:
        return None
    s["axis"] = best_axis
    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:
        with h5py.File(path, "r") as f:
            chans = {ax: f[f"samples/{ax}"][:] for ax in _GEO_CHANNELS
                     if f"samples/{ax}" in f}
    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)