From 7b44e2c48d8f97b067b7d4658ce8b1444a53f035 Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Sun, 12 Jul 2026 09:54:17 -0700 Subject: [PATCH] Route A: Wood-Anderson (response-removed) amplitude pipeline + runbook Scripts to re-measure amplitudes as Wood-Anderson displacement on a pnwstore host, addressing the Method-B seismological caveats (counts scalar response, OBS redeployment station-term stability, fixed short window, no SNR gate, missing NC/BK): - route_a_build_station_inventory.py: coords + response + operational epochs for all networks (NC/BK via NCEDC); StationXML + epoch CSV. - route_a_wa_amplitudes.py: per P and S pick -> remove response, simulate IASPEI Wood-Anderson, peak WA displacement in a distance-scaled window (all components, vertical fallback), pre-signal SNR, response-epoch tag. Resumable/shardable. - route_a_build_dataset.py: SNR gate + epoch-keyed station ids -> dataset feeding the existing phase3->phase2->phase4/5 pipeline unchanged (--suffix _routeA). - ROUTE_A_RUNBOOK.md: environment, steps, validation vs Method B, parameter notes. Design (per PI): Wood-Anderson; P and S; distance-scaled window; all components + Z fallback. Untested against pnwstore -- validate on a small --limit first. Co-Authored-By: Claude Opus 4.8 (1M context) --- 4_relocation/magnitude/ROUTE_A_RUNBOOK.md | 93 +++++++++ .../magnitude/route_a_build_dataset.py | 67 +++++++ .../route_a_build_station_inventory.py | 91 +++++++++ .../magnitude/route_a_wa_amplitudes.py | 177 ++++++++++++++++++ 4 files changed, 428 insertions(+) create mode 100644 4_relocation/magnitude/ROUTE_A_RUNBOOK.md create mode 100644 4_relocation/magnitude/route_a_build_dataset.py create mode 100644 4_relocation/magnitude/route_a_build_station_inventory.py create mode 100644 4_relocation/magnitude/route_a_wa_amplitudes.py diff --git a/4_relocation/magnitude/ROUTE_A_RUNBOOK.md b/4_relocation/magnitude/ROUTE_A_RUNBOOK.md new file mode 100644 index 00000000..9be3a847 --- /dev/null +++ b/4_relocation/magnitude/ROUTE_A_RUNBOOK.md @@ -0,0 +1,93 @@ +# Route A runbook — Wood-Anderson magnitudes (response-removed) + +Route A re-measures amplitudes as **Wood-Anderson displacement** (instrument response +removed) to address the main seismological caveats of the counts-based Method B: + +| Method-B caveat | Route-A fix | +|---|---| +| counts + scalar station term can't represent a frequency-dependent response | remove response → WA displacement; station term becomes ~pure site | +| one station term per reused OBS code across yearly redeployments | per-**epoch** station terms (inventory epochs → station id `NET.STA@epoch`) | +| fixed 2.5 s window misses the delayed S/Lg peak | **distance-scaled** window per phase | +| no SNR gate (noisy OBS) | pre-signal **SNR** measured and gated | +| Wood-Anderson product missing NC/BK (old IRIS-only run) | NC/BK fetched via **NCEDC** (`utils/data_client.py`) | + +Decisions in force: **Wood-Anderson** amplitude; measured for **both P and S**; +**distance-scaled** window; **all components, vertical as fallback**. + +## Where to run +A UW-internal host with pnwstore + FDSN/NCEDC access: +``` +pixi install --environment internal # adds pnwstore +git pull # get this branch +``` +Step 1 is one network request per pick (~10^6 picks) — slow. It appends on the fly, +supports `--start-index` (resume / shard), and should be validated on a small +`--limit` first. Inputs are the pick-assignment CSV used by Method B +(`Cascadia_updated_catalog_picks_assignment_ver_3.csv`, which already carries event +and station coordinates). + +## Steps + +**0. Station inventory** (coords + response + epochs; any machine with internet): +``` +python route_a_build_station_inventory.py \ + --picks /wd1/.../Cascadia_updated_catalog_picks_assignment_ver_3.csv \ + --out-xml station_inventory.xml --out-csv station_epochs.csv +``` +Check the reported count of "stations with multiple epochs" — those are the OBS +redeployments the epoch-keyed station terms will separate. + +**1. Wood-Anderson amplitudes** (pnwstore host; validate small first): +``` +# smoke test on 500 picks — inspect wa_amp_mm (expect ~1e-4..10 mm) and snr +python route_a_wa_amplitudes.py --picks --inventory station_inventory.xml \ + --out raw_wa_amplitudes.csv --source pnwstore --limit 500 +# full run (optionally shard by --start-index across processes/hosts) +python route_a_wa_amplitudes.py --picks --inventory station_inventory.xml \ + --out raw_wa_amplitudes.csv --source pnwstore +``` +Output `raw_wa_amplitudes.csv`: one row per pick with `wa_amp_mm`, `snr`, `n_comp`, +`epoch`, `dist_hypo_km`, event/station coordinates, and a `reason` field. + +**2. Build the analysis dataset** (SNR gate + epoch-keyed station ids; any machine): +``` +python route_a_build_dataset.py --raw raw_wa_amplitudes.csv \ + --out ../../data/magnitude/amp_distance_dataset_routeA.csv --min-snr 3 --epoch-station +``` + +**3. Inversion → absolute ML → QC** (reuse the existing pipeline, `--suffix _routeA`): +``` +python phase3_route_b_relative_magnitude.py \ + --dataset ../../data/magnitude/amp_distance_dataset_routeA.csv \ + --outdir ../../data/magnitude --fix-n 1.0 --suffix _routeA +python phase2_anchor_comcat_ml.py \ + --events ../../data/magnitude/route_b_event_relative_mag_routeA.csv \ + --catalog ../../data/Cascadia_relocated_catalog_ver_3.csv \ + --outdir ../../data/magnitude --suffix _routeA +python phase4_qc_and_gr.py --catalog ../../data/magnitude/cascadia_catalog_ML_routeA.csv --tag routeA +python phase5_pygmt_map.py --catalog ../../data/magnitude/cascadia_catalog_ML_routeA.csv \ + --out ../../data/magnitude/cascadia_ML_map_routeA.png +``` +Because the epoch is folded into the station id, phase3's per-`(station, phase)` +terms are automatically per-deployment — no change to the inversion code. + +## Validation (Route A vs Method B) +Join `cascadia_catalog_ML_routeA.csv` and `cascadia_catalog_ML_kpos.csv` on +`event_id` and compare ML. They should agree within the ~0.3-mag station scatter; +systematic offsets flag response/gain problems (e.g. a station whose counts-based +term was wrong because of a redeployment). Expect Route A to (a) fill NC/BK, (b) +tighten the station scatter, and (c) reduce the low-magnitude / offshore calibration +bias, since the amplitudes are now physical. + +## Method / parameter notes +- **Wood-Anderson**: IASPEI constants (T0=0.8 s, h=0.7, static gain 2080), applied to + ground velocity after `remove_response(output="VEL")`. The absolute gain is absorbed + by the ComCat-ML calibration, so it does not affect final ML. +- **Window**: P `[t_P−0.3, t_P + min(1+0.03r, 15, 0.8·(t_S−t_P))]`; + S `[t_S−0.3, t_S + min(2+0.06r, 60)]` s (r = hypocentral distance, km). +- **Components**: peak over all available components; vertical used when only a + vertical channel exists. +- **SNR**: peak signal / RMS of a pre-signal noise window (default 10 s); gate at 3. +- **Epochs**: the response epoch containing the pick time tags each measurement, so a + reused OBS station code deployed with different instruments gets separate station + terms. diff --git a/4_relocation/magnitude/route_a_build_dataset.py b/4_relocation/magnitude/route_a_build_dataset.py new file mode 100644 index 00000000..28c4537c --- /dev/null +++ b/4_relocation/magnitude/route_a_build_dataset.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 +""" +Route A, step 2: assemble the analysis-ready amplitude-distance dataset from the raw +Wood-Anderson measurements. Applies the SNR gate and (optionally) folds the response +epoch into the station identifier so the inversion solves an epoch-specific station +term for redeployed OBS. The output matches the schema consumed by +phase3_route_b_relative_magnitude.py, so the rest of the pipeline +(phase3 -> phase2 -> phase4/5/6) runs unchanged with --suffix _routeA. + +Usage: + python route_a_build_dataset.py --raw raw_wa_amplitudes.csv \ + --out ../../data/magnitude/amp_distance_dataset_routeA.csv \ + --min-snr 3 --epoch-station +""" +from __future__ import annotations + +import argparse +import os + +import numpy as np +import pandas as pd + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--raw", required=True, help="raw_wa_amplitudes.csv") + ap.add_argument("--out", default="../../data/magnitude/amp_distance_dataset_routeA.csv") + ap.add_argument("--min-snr", type=float, default=3.0) + ap.add_argument("--epoch-station", action="store_true", + help="append @epoch to the station id (per-deployment station terms)") + args = ap.parse_args(argv) + + df = pd.read_csv(os.path.expanduser(args.raw)) + n0 = len(df) + if "reason" in df.columns: + df = df[df["reason"] == "ok"] + amp = pd.to_numeric(df["wa_amp_mm"], errors="coerce") + snr = pd.to_numeric(df.get("snr"), errors="coerce") + keep = amp.notna() & (amp > 0) & ((snr >= args.min_snr) | snr.isna()) + n_lowsnr = int(((snr < args.min_snr)).sum()) + df = df[keep].copy() + df["amp"] = amp[keep].to_numpy() + df["log10A"] = np.log10(df["amp"]) + + if args.epoch_station: # per-deployment station terms + ep = df["epoch"].fillna("").astype(str) + df["station"] = np.where(ep.ne(""), df["station"].astype(str) + "@" + ep, + df["station"].astype(str)) + + cols = ["arid", "event_id", "station", "network", "phase", "evla", "evlo", "evdp", + "stla", "stlo", "stel_m", "dist_hypo_km", "amp", "log10A", "snr", "epoch"] + cols = [c for c in cols if c in df.columns] + out = os.path.expanduser(args.out) + os.makedirs(os.path.dirname(out), exist_ok=True) + df[cols].to_csv(out, index=False) + + print(f"raw rows : {n0:,}") + print(f"dropped low-SNR (<{args.min_snr}) : {n_lowsnr:,}") + print(f"clean observations : {len(df):,} -> {out}") + print(f"events : {df['event_id'].nunique():,}") + print(f"station terms (ids) : {df['station'].nunique():,}" + + (" (epoch-keyed)" if args.epoch_station else "")) + print("phase counts :", df["phase"].value_counts().to_dict()) + + +if __name__ == "__main__": + main() diff --git a/4_relocation/magnitude/route_a_build_station_inventory.py b/4_relocation/magnitude/route_a_build_station_inventory.py new file mode 100644 index 00000000..b1bc7def --- /dev/null +++ b/4_relocation/magnitude/route_a_build_station_inventory.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python3 +""" +Route A, step 0: build a station inventory (coordinates + instrument response + +operational epochs) for every network in the pick catalog. + +The inventory serves two purposes downstream: + * response removal for Wood-Anderson simulation (route_a_wa_amplitudes.py); + * exposing per-station response *epochs* so ocean-bottom seismometers that were + redeployed with different instruments under a reused station code get an + epoch-specific station term in the inversion (the main Method-B soundness fix). + +NC/BK are served by NCEDC (see utils/data_client.py); all others by the chosen FDSN +datacenter. Run on a host with FDSN + NCEDC reachable (any machine with internet). + +Outputs: + station_inventory.xml StationXML (input to route_a_wa_amplitudes.py) + station_epochs.csv one row per (network, station, channel, epoch) + +Usage: + python route_a_build_station_inventory.py \ + --picks /path/Cascadia_updated_catalog_picks_assignment_ver_3.csv \ + --out-xml station_inventory.xml --out-csv station_epochs.csv +""" +from __future__ import annotations + +import argparse + +import pandas as pd +from obspy import Inventory, UTCDateTime +from obspy.clients.fdsn import Client + +NCEDC_NETWORKS = frozenset(["NC", "BK"]) + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--picks", required=True, help="picks CSV with a NET.STA 'station' column") + ap.add_argument("--out-xml", default="station_inventory.xml") + ap.add_argument("--out-csv", default="station_epochs.csv") + ap.add_argument("--t0", default="2010-01-01") + ap.add_argument("--t1", default="2016-01-01") + ap.add_argument("--fdsn", default="IRIS") + ap.add_argument("--channels", default="?H?,?N?") + args = ap.parse_args(argv) + + picks = pd.read_csv(args.picks) + picks.columns = [c.strip() for c in picks.columns] + networks = sorted({str(s).split(".")[0].strip() + for s in picks["station"].dropna()}) + t0, t1 = UTCDateTime(args.t0), UTCDateTime(args.t1) + print("networks:", networks) + + inv = Inventory(networks=[], source="route_a_build_station_inventory") + rows = [] + for net in networks: + client = Client("NCEDC") if net in NCEDC_NETWORKS else Client(args.fdsn) + try: + sub = client.get_stations(network=net, station="*", location="*", + channel=args.channels, starttime=t0, endtime=t1, + level="response") + except Exception as e: + print(f"{net}: get_stations failed: {e}") + continue + inv += sub + n_ch = 0 + for n in sub: + for s in n: + for c in s.channels: + n_ch += 1 + rows.append(dict( + network=n.code, station=s.code, location=c.location_code, + channel=c.code, latitude=c.latitude, longitude=c.longitude, + elevation=c.elevation, depth=c.depth, + start_date=str(c.start_date), end_date=str(c.end_date), + sample_rate=c.sample_rate, + has_response=c.response is not None)) + print(f"{net}: {n_ch} channel-epochs") + + inv.write(args.out_xml, format="STATIONXML") + df = pd.DataFrame(rows) + df.to_csv(args.out_csv, index=False) + # flag stations with >1 epoch per channel band (redeployments -> epoch-keyed terms) + if len(df): + multi = (df.groupby(["network", "station"])["start_date"].nunique() > 1).sum() + print(f"wrote {args.out_xml} and {args.out_csv}: " + f"{len(df)} channel-epochs, {df.groupby(['network','station']).ngroups} stations, " + f"{multi} stations with multiple epochs (redeployments)") + + +if __name__ == "__main__": + main() diff --git a/4_relocation/magnitude/route_a_wa_amplitudes.py b/4_relocation/magnitude/route_a_wa_amplitudes.py new file mode 100644 index 00000000..5158173c --- /dev/null +++ b/4_relocation/magnitude/route_a_wa_amplitudes.py @@ -0,0 +1,177 @@ +#!/usr/bin/env python3 +""" +Route A, step 1: Wood-Anderson amplitudes for each P and S pick. + +For every pick this script fetches the waveform (NC/BK -> NCEDC via +utils/data_client.py), removes the instrument response, simulates a Wood-Anderson +seismograph, and measures the peak WA displacement in a *distance-scaled* window, +taken as the maximum over all available components (vertical used as a fallback). +It also records a pre-signal signal-to-noise ratio and the response *epoch* the +pick falls in (so redeployed OBS get an epoch-specific station term downstream). + +This addresses the main Method-B seismological caveats: + * counts -> physical Wood-Anderson displacement (instrument response removed); + * distance-scaled measurement window (captures the delayed S/Lg peak); + * SNR for quality gating; per-epoch tagging for OBS redeployments; + * NC/BK filled via NCEDC (missing from the old IRIS-only Wood-Anderson run). + +MUST run on a host with pnwstore + FDSN/NCEDC access +(`pixi install --environment internal`). It is slow (one request per pick); it +appends to --out on the fly and supports --start-index for resume / sharding. +VALIDATE on a small --limit first and sanity-check wa_amp_mm before the full run. + +Usage: + python route_a_wa_amplitudes.py \ + --picks /path/Cascadia_updated_catalog_picks_assignment_ver_3.csv \ + --inventory station_inventory.xml \ + --out raw_wa_amplitudes.csv --source pnwstore +""" +from __future__ import annotations + +import argparse +import csv +import os +import sys +import time + +import numpy as np +import pandas as pd +from obspy import UTCDateTime, read_inventory + +sys.path.append(os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "..", "utils")) +from data_client import get_waveforms # noqa: E402 (routes NC/BK -> NCEDC) + +R_EARTH_KM = 6371.0 +VP, VS = 6.0, 3.5 # km/s, for S-P timing and window scaling + +# IASPEI-standard Wood-Anderson (T0 = 0.8 s, damping h = 0.7, static gain 2080), +# applied to GROUND VELOCITY (one zero at the origin). The absolute gain is later +# absorbed by the ComCat-ML calibration, so 2080 vs 2800 does not affect final ML. +PAZ_WA = {"poles": [-5.49779 - 5.60886j, -5.49779 + 5.60886j], + "zeros": [0j], "gain": 1.0, "sensitivity": 2080.0} + +OUT_COLS = ["arid", "event_id", "network", "station", "phase", "evla", "evlo", "evdp", + "stla", "stlo", "stel_m", "dist_hypo_km", "wa_amp_mm", "snr", "n_comp", + "epoch", "reason"] + + +def haversine_km(la1, lo1, la2, lo2): + la1, lo1, la2, lo2 = map(np.radians, (la1, lo1, la2, lo2)) + a = (np.sin((la2 - la1) / 2) ** 2 + + np.cos(la1) * np.cos(la2) * np.sin((lo2 - lo1) / 2) ** 2) + return 2 * R_EARTH_KM * np.arcsin(np.sqrt(np.clip(a, 0, 1))) + + +def phase_window(phase, r_km): + """(pre, post) seconds around the pick. Post-window grows with distance to + capture the delayed peak; the P window is capped short of the S arrival.""" + tsp = r_km * (1.0 / VS - 1.0 / VP) # S-P time (s) + if phase == "P": + post = min(1.0 + 0.03 * r_km, 15.0, 0.8 * max(tsp, 1.25)) + return 0.3, max(post, 1.0) + post = min(2.0 + 0.06 * r_km, 60.0) # S/Lg coda grows with distance + return 0.3, post + + +def epoch_id(inv, net, sta, t): + """Index of the response epoch containing time t for net.sta (redeployment tag).""" + try: + chans = inv.select(network=net, station=sta, time=t).get_contents()["channels"] + starts = sorted({inv.get_channel_metadata(cid, t).get("starttime") + for cid in chans} - {None}) + return f"{starts[0].date}" if starts else "" + except Exception: + return "" + + +def _phase(v): + s = str(v).strip().upper() + return "P" if s in ("P", "0") else "S" if s in ("S", "1") else s + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--picks", required=True) + ap.add_argument("--inventory", required=True) + ap.add_argument("--out", default="raw_wa_amplitudes.csv") + ap.add_argument("--source", default="pnwstore") + ap.add_argument("--sample-rate", type=int, default=100) + ap.add_argument("--highpass", type=float, default=1.0, help="post-WA high-pass (Hz)") + ap.add_argument("--noise-win", type=float, default=10.0, help="pre-signal noise window (s)") + ap.add_argument("--start-index", type=int, default=0, help="resume / shard start row") + ap.add_argument("--limit", type=int, default=None, help="process at most N picks (testing)") + args = ap.parse_args(argv) + + inv = read_inventory(args.inventory) + picks = pd.read_csv(args.picks) + picks.columns = [c.strip() for c in picks.columns] + ev_col = "idx" if "idx" in picks.columns else "event_id" + + sl = slice(args.start_index, args.start_index + args.limit if args.limit else None) + write_header = (args.start_index == 0) or (not os.path.exists(args.out)) + fh = open(args.out, "w" if write_header else "a", newline="") + w = csv.DictWriter(fh, fieldnames=OUT_COLS) + if write_header: + w.writeheader() + + for _, row in picks.iloc[sl].iterrows(): + net, sta = str(row["station"]).split(".")[0].strip(), str(row["station"]).split(".")[1].strip() + phase = _phase(row["phase"]) + rec = {k: "" for k in OUT_COLS} + rec.update(arid=row.get("arid"), event_id=row.get(ev_col), network=net, + station=f"{net}.{sta}", phase=phase, n_comp=0) + try: + tp = UTCDateTime(str(row["time_pick"])) + evla, evlo, evdp = float(row["latitude"]), float(row["longitude"]), float(row["depth"]) + stla, stlo, stel = float(row["slatitude"]), float(row["slongitude"]), float(row["selevation"]) + except Exception as e: + rec["reason"] = f"row:{e}"; w.writerow(rec); continue + r = float(np.hypot(haversine_km(evla, evlo, stla, stlo), evdp + stel / 1000.0)) + pre, post = phase_window(phase, r) + rec.update(evla=evla, evlo=evlo, evdp=evdp, stla=stla, stlo=stlo, stel_m=stel, + dist_hypo_km=round(r, 3)) + + try: + st = get_waveforms(net, sta, "*H*", tp - (args.noise_win + 5), tp + post + 5, + source=args.source) + except Exception as e: + rec["reason"] = f"fetch:{str(e)[:80]}"; w.writerow(rec); continue + if len(st) == 0: + rec["reason"] = "no_data"; w.writerow(rec); continue + + try: + st.merge(method=1, fill_value="interpolate") + st.resample(args.sample_rate) + st.detrend("demean"); st.taper(0.05) + st.remove_response(inventory=inv, output="VEL", water_level=60) + st.simulate(paz_simulate=PAZ_WA) # ground velocity -> WA displacement (m) + st.filter("highpass", freq=args.highpass) + except Exception as e: + rec["reason"] = f"resp:{str(e)[:80]}"; w.writerow(rec); continue + + rec["epoch"] = epoch_id(inv, net, sta, tp) + sig_lo, sig_hi = tp - pre, tp + post + noi_lo, noi_hi = tp - args.noise_win - 1, tp - 1 + amp, noise, ncomp = 0.0, 0.0, 0 + for tr in st: + d = tr.slice(sig_lo, sig_hi).data + n = tr.slice(noi_lo, noi_hi).data + if len(d): + amp = max(amp, float(np.max(np.abs(d)))); ncomp += 1 + if len(n): + noise = max(noise, float(np.sqrt(np.mean(n.astype(float) ** 2)))) + if ncomp == 0: + rec["reason"] = "no_window_data"; w.writerow(rec); continue + + rec.update(wa_amp_mm=amp * 1000.0, # m -> mm + snr=(amp / noise) if noise > 0 else np.nan, + n_comp=ncomp, reason="ok") + w.writerow(rec); fh.flush() + time.sleep(0.05) + + fh.close() + print(f"done -> {args.out}") + + +if __name__ == "__main__": + main()