Route A: Wood-Anderson (response-removed) amplitude pipeline for magnitudes

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 <picks.csv> --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 <picks.csv> --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.

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()

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()