Phase pickers on three sequences outside the United States¶
The global campaign will run one PhaseNet weight set over every network
EarthScope holds, most of it outside the United States and Japan. The
fine-tuned quakescope2026 was meant to be the more general of the two
candidates; the western-states benchmark says its jma_wc parent still
recovers more arrivals. Neither of those tests left the United States. This
one does.
Three sequences, each with the network operator's own analyst picks published through an FDSN event service, chosen for large magnitude, a dense reviewed aftershock catalogue in the hours after the mainshock, and open waveforms.
| Sequence | Date | M | Operator | Setting |
|---|---|---|---|---|
| Kaikōura, New Zealand | 2016-11-13 | 7.8 | GeoNet (GNS Science) | oblique thrust and strike-slip through more than a dozen crustal faults; sparse permanent network, nearest station 12 km, most reviewed picks 80–120 km out |
| Norcia, central Italy | 2016-10-30 | 6.5 | INGV | normal faulting in the Apennines, the largest shock of the Amatrice–Visso–Norcia sequence; permanent network plus the temporary stations deployed after August |
| Thessaly, central Greece | 2021-03-03 | 6.3 | NOA | normal-faulting doublet (M6.3, then M6.0 the next day); dense HL/HT/HP coverage with a station 5 km from the epicentre |
Kaiser et al. (2017, SRL 88, doi:10.1785/0220170018) is the seismological report for Kaikōura and Lanza et al. (2019, GRL 46, doi:10.1029/2019GL082780) relocated its aftershocks. Chiaraluce et al. (2017, SRL 88, doi:10.1785/0220160221) describe the central Italy sequence; Michele et al. (2020, Sci. Rep. 10, doi:10.1038/s41598-019-43393-2) hand-picked its early aftershocks. Karakostas et al. (2021, Bull. Geol. Soc. Greece 58) and Kassaras et al. (2022, J. Geodyn. 150, doi:10.1016/j.jog.2022.101898) analyse the Thessaly doublet. None of those papers' own pick files is used here: the reference is the operator's routine bulletin, harvested live, which is what a global campaign will be judged against.
Weight sets¶
The same four as the US comparison, so the two reports read together.
quakescope2026 |
the v7 fine-tune, the production candidate |
jma_wc |
the SeisBench Japanese model v7 was fine-tuned from, the baseline that matters |
original |
Zhu & Beroza (2019), the published reference |
instance |
trained on the Italian INSTANCE dataset: in-domain for Norcia, which makes it the control there |
How it is scored¶
The method is the one in
phasenet_sequence_comparison.ipynb:
an aftershock window starting ten minutes after the mainshock, six stations,
every model on identical waveforms, recall against manual analyst picks
within 0.5 s. Recall is the metric; precision is not computed, because an
aftershock catalogue is not exhaustive and an unmatched model pick may be an
arrival nobody had time to mark. Section 7 holds the pick budget fixed rather
than the threshold, which the US comparison showed is the only fair way to
rank weight sets whose probabilities sit on different scales.
Two things differ from the US notebook. Stations are chosen by where the analysts' picks are, not by distance: nearest-first at Kaikōura would put the single close station in and leave the reference thin. And there is no curated event list, because none of these operators serves a revisable catalogue for a closed sequence in the way ComCat does; the event set is pinned instead by caching the harvest beside the notebook, so a re-run scores against the same picks.
Where the picks come from, and why it is three different queries¶
Every operator publishes analyst arrivals through an FDSN event service, and
every one does it differently. GeoNet rejects includearrivals outright but
carries picks and arrivals in the per-event QuakeML, so each event is fetched
by id. INGV accepts includearrivals only with an eventid, so it is the
same loop. NOA answers region queries with arrivals in one go but refuses
anything larger than about an hour of a busy sequence, so the window is asked
for in hourly pieces. Section 4a reports what each service's QuakeML actually
says about its picks.
import io, json, time
from collections import Counter, defaultdict
from pathlib import Path
import boto3
import matplotlib.pyplot as plt
import numpy as np
import obspy
import pandas as pd
import seisbench.models as sbm
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.clients.fdsn.header import FDSNException
from obspy.geodetics import locations2degrees
%matplotlib inline
1. Configuration¶
# Three sequences. `radius` bounds both the event harvest and the station search;
# `min_mag` is the operator's completeness in the hours after the mainshock,
# roughly - a lower floor only adds events with few picks.
SEQUENCES = {
"Kaikoura 2016": dict(
time=UTCDateTime("2016-11-13T11:02:56"), lat=-42.69, lon=173.02, mag=7.8,
agency="GEONET", radius=2.0, min_mag=2.5, window_min=180,
note="sparse permanent network; GeoNet reviews every located event",
),
"Norcia 2016": dict(
time=UTCDateTime("2016-10-30T06:40:18"), lat=42.83, lon=13.11, mag=6.5,
agency="INGV", radius=0.8, min_mag=2.0, window_min=120,
note="largest of the Amatrice-Visso-Norcia sequence; permanent plus temporary stations",
),
"Thessaly 2021": dict(
time=UTCDateTime("2021-03-03T10:16:08"), lat=39.75, lon=22.20, mag=6.3,
agency="NOA", radius=1.0, min_mag=2.0, window_min=180,
note="doublet with an M6.0 the next day; dense Greek network with a station 5 km away",
),
}
WEIGHTS = ["quakescope2026", "jma_wc", "original", "instance"]
WINDOW_START = 600 # s after the origin: past the mainshock coda
ORIGIN_LEAD = 180 # s: an origin this far before the window can still put arrivals in it
N_STATIONS = 6
CHANNEL_PREFERENCE = ["HH", "EH", "BH"]
DETECT_FLOOR = 0.02 # run once low, threshold offline
REPORT_THRESHOLD = 0.3
THRESHOLD_SWEEP = [0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7]
MATCH_TOL = 0.5 # s, an analyst pick counts as recovered within this
CACHE = Path("global_sequences_cache") # gitignored: harvested picks, so the reference is pinned
CACHE.mkdir(exist_ok=True)
COLORS = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100"]
C_P, C_S = "#2a78d6", "#eb6834"
plt.rcParams.update({"figure.dpi": 110, "axes.grid": True, "grid.alpha": 0.25,
"grid.linewidth": 0.5, "axes.axisbelow": True})
2. Data access¶
Waveforms and picks both come from the operator's FDSN services, routed by agency. Nothing here touches the campaign's S3 readers; section 8 checks whether the campaign has reached these networks yet.
_clients = {}
def client_for(agency):
if agency not in _clients:
_clients[agency] = Client(agency, timeout=300)
return _clients[agency]
def _event_id(ev):
rid = str(ev.resource_id.id)
return rid.split("eventId=")[-1].split("eventid=")[-1].split("/")[-1]
def _pick_rows(ev, label):
"""One row per arrival of the preferred origin, with the pick's own metadata."""
origin = ev.preferred_origin() or (ev.origins[0] if ev.origins else None)
if origin is None:
return []
mag = ev.preferred_magnitude() or (ev.magnitudes[0] if ev.magnitudes else None)
by_id = {p.resource_id.id: p for p in ev.picks}
rows = []
for arr in origin.arrivals:
pick = by_id.get(arr.pick_id.id)
if pick is None or not arr.phase:
continue
phase = arr.phase[0].upper()
if phase not in ("P", "S"):
continue
w = pick.waveform_id
rows.append(dict(
sequence=label, event=_event_id(ev), origin=origin.time.datetime,
mag=mag.mag if mag else np.nan,
station=f"{w.network_code}.{w.station_code}", channel=w.channel_code or "",
phase=phase, time=pick.time.datetime,
mode=str(pick.evaluation_mode), status=str(pick.evaluation_status),
method=str(pick.method_id.id).split("/")[-1] if pick.method_id else "",
agency=pick.creation_info.agency_id if pick.creation_info else "",
time_weight=arr.time_weight if arr.time_weight is not None else np.nan,
onset=str(pick.onset), uncertainty=pick.time_errors.uncertainty
if pick.time_errors and pick.time_errors.uncertainty is not None else np.nan,
))
return rows
def harvest(label, seq, t0, t1):
"""Every arrival the operator publishes for events in the window, all stations.
Cached beside the notebook. The event services serve a revisable catalogue,
so the cache is what pins the reference between runs - delete it to re-harvest.
"""
cache = CACHE / f"{label.replace(' ', '_')}_picks.parquet"
if cache.exists():
return pd.read_parquet(cache)
client = client_for(seq["agency"])
region = dict(latitude=seq["lat"], longitude=seq["lon"], maxradius=seq["radius"],
minmagnitude=seq["min_mag"])
rows, seen = [], set()
t_start = time.time()
if seq["agency"] == "NOA":
# region queries carry arrivals, but only for about an hour of a busy
# sequence at a time; split further if the service still refuses
edges = np.arange(float(t0 - ORIGIN_LEAD), float(t1) + 1, 3600.0)
pieces = [(UTCDateTime(a), UTCDateTime(min(a + 3600.0, float(t1)))) for a in edges if a < float(t1)]
while pieces:
a, b = pieces.pop(0)
try:
cat = client.get_events(starttime=a, endtime=b, includearrivals=True, **region)
except FDSNException as exc:
if "413" in str(exc) or "too much" in str(exc).lower():
m = a + (b - a) / 2
pieces = [(a, m), (m, b)] + pieces
continue
if "204" in str(exc) or "No data" in str(exc):
continue
raise
for ev in cat:
if _event_id(ev) in seen:
continue
seen.add(_event_id(ev)); rows += _pick_rows(ev, label)
else:
cat = client.get_events(starttime=t0 - ORIGIN_LEAD, endtime=t1, **region)
print(f" {seq['agency']}: {len(cat)} events, fetching each one's arrivals", flush=True)
for ev in cat:
eid = _event_id(ev)
try:
if seq["agency"] == "INGV":
full = client.get_events(eventid=eid, includearrivals=True)
else: # GeoNet: the per-event QuakeML carries picks
full = client.get_events(eventid=eid)
except FDSNException:
continue
for e in full:
rows += _pick_rows(e, label)
df = pd.DataFrame(rows)
if len(df):
df["origin"] = pd.to_datetime(df.origin, utc=True)
df["time"] = pd.to_datetime(df.time, utc=True)
df.to_parquet(cache)
print(f" harvested {len(df)} arrivals from {df.event.nunique() if len(df) else 0} events "
f"in {time.time() - t_start:.0f} s")
return df
def station_table(seq, picks, t0, t1):
"""Candidate stations within the radius with a pickable band, and how many
manual picks the analysts made on each inside the window."""
inv = client_for(seq["agency"]).get_stations(
latitude=seq["lat"], longitude=seq["lon"], maxradius=seq["radius"],
channel="HH?,EH?,BH?", starttime=t0, endtime=t1, level="channel")
rows = []
lo, hi = pd.Timestamp(t0.datetime, tz="UTC"), pd.Timestamp(t1.datetime, tz="UTC")
inwin = picks[(picks["time"] >= lo) & (picks["time"] <= hi) & (picks["mode"] == "manual")]
counts = inwin.groupby(["station", "phase"]).size().unstack(fill_value=0)
for net in inv:
for sta in net:
bands = {ch.code[:2]: ch.sample_rate for ch in sta}
band = next((b for b in CHANNEL_PREFERENCE if b in bands), None)
if band is None:
continue
key = f"{net.code}.{sta.code}"
rows.append(dict(station=key, band=band, rate=bands[band],
km=locations2degrees(seq["lat"], seq["lon"], sta.latitude, sta.longitude) * 111.19,
P=int(counts.loc[key, "P"]) if key in counts.index and "P" in counts else 0,
S=int(counts.loc[key, "S"]) if key in counts.index and "S" in counts else 0))
t = pd.DataFrame(rows).drop_duplicates("station")
t["picks"] = t.P + t.S
return t.sort_values(["picks", "km"], ascending=[False, True]).reset_index(drop=True)
def fetch_station(agency, key, band, t0, t1):
net, sta = key.split(".")
try:
st = client_for(agency).get_waveforms(net, sta, "*", band + "?", t0, t1)
except FDSNException:
return None
st.merge(fill_value=0)
if len(st) < 3:
return None
expected = (t1 - t0) * st[0].stats.sampling_rate
if st[0].stats.npts < 0.5 * expected:
return None
return st
def select_and_fetch(label, seq, picks, t0, t1):
"""The N_STATIONS richest-in-picks stations that actually return the window."""
table = station_table(seq, picks, t0, t1)
streams, rates, chosen = {}, {}, []
for r in table.itertuples():
if len(streams) >= N_STATIONS:
break
st = fetch_station(seq["agency"], r.station, r.band, t0, t1)
if st is None:
continue
streams[r.station] = st
rates[r.station] = st[0].stats.sampling_rate
chosen.append(r.Index)
return streams, rates, table, chosen
def reference_from(picks, stations, t0, t1):
"""Manual picks on the chosen stations inside the window, keyed (station, phase),
duplicates across events collapsed at MATCH_TOL."""
lo, hi = pd.Timestamp(t0.datetime, tz="UTC"), pd.Timestamp(t1.datetime, tz="UTC")
sel = picks[(picks["time"] >= lo) & (picks["time"] <= hi) & (picks["mode"] == "manual")
& picks.station.isin(stations)]
out = {}
for (sta, pha), g in sel.groupby(["station", "phase"]):
times, kept = sorted(UTCDateTime(t) for t in g["time"]), []
for t in times:
if not kept or t - kept[-1] > MATCH_TOL:
kept.append(t)
out[(sta, pha)] = kept
return out
def at_threshold(store, thr):
return {k: sorted(t for t, v in vals if v >= thr) for k, vals in store.items()}
def match(reference, candidate, tol=MATCH_TOL):
"""Greedy nearest match; each candidate pick is consumed at most once."""
used, residuals = set(), []
for a in reference:
best_i = best_d = None
for i, m in enumerate(candidate):
if i in used:
continue
d = m - a
if abs(d) <= tol and (best_d is None or abs(d) < abs(best_d)):
best_i, best_d = i, d
if best_i is not None:
used.add(best_i)
residuals.append(best_d)
return residuals, len(candidate) - len(used)
3. Load the weight sets¶
available = sbm.PhaseNet.list_pretrained()
models = {}
for name in WEIGHTS:
if name not in available:
print(f"{name:<16} not installed - skipping")
continue
try:
models[name] = sbm.PhaseNet.from_pretrained(name)
print(f"{name:<16} loaded")
except Exception as exc:
print(f"{name:<16} could not load ({type(exc).__name__}) - skipping")
if not models:
raise RuntimeError("no weight sets available")
names = list(models)
quakescope2026 loaded jma_wc loaded original loaded instance loaded
4. Harvest, select, fetch, pick¶
The slow part is the harvest: one request per event for GeoNet and INGV. It runs once and is cached.
results, harvested = {}, {}
for label, seq in SEQUENCES.items():
print(f"{label} (M{seq['mag']}, {seq['note']})")
t0 = seq["time"] + WINDOW_START
t1 = t0 + seq["window_min"] * 60
picks = harvest(label, seq, t0, t1)
harvested[label] = picks
if not len(picks):
print(" no arrivals returned"); continue
streams, rates, table, chosen = select_and_fetch(label, seq, picks, t0, t1)
print(f" {len(table)} candidate stations, {len(streams)} chosen by manual-pick count with data:")
for i in chosen:
r = table.loc[i]
print(f" {r.station:<10} {r.band}@{r.rate:g} Hz {r.km:5.0f} km {r.P:3d} P {r.S:3d} S")
skipped = [table.loc[i].station for i in table.index[:max(chosen) + 1] if i not in chosen] if chosen else []
if skipped:
print(f" no data in the window: {', '.join(skipped)}")
reference = reference_from(picks, list(streams), t0, t1)
n_p = sum(len(v) for (_, ph), v in reference.items() if ph == "P")
n_s = sum(len(v) for (_, ph), v in reference.items() if ph == "S")
print(f" reference: {n_p} P, {n_s} S manual picks on these stations")
store = {}
for name, model in models.items():
per = defaultdict(list)
for sta, st in streams.items():
try:
out = model.classify(st, P_threshold=DETECT_FLOOR, S_threshold=DETECT_FLOOR)
except Exception as exc:
print(f" {name} {sta}: {type(exc).__name__}"); continue
for p in out.picks:
per[(sta, p.phase)].append((p.peak_time, float(p.peak_value)))
store[name] = dict(per)
tot = sum(1 for v in per.values() for t, q in v if q >= REPORT_THRESHOLD)
print(f" {name:<16} {tot} picks at {REPORT_THRESHOLD}")
results[label] = dict(streams=streams, picks=store, reference=reference, t0=t0, t1=t1,
rates=rates, table=table, chosen=chosen)
print()
Kaikoura 2016 (M7.8, sparse permanent network; GeoNet reviews every located event)
27 candidate stations, 6 chosen by manual-pick count with data:
NZ.BSWZ EH@100 Hz 129 km 80 P 99 S
NZ.TUWZ EH@100 Hz 160 km 71 P 96 S
NZ.TCW EH@100 Hz 195 km 62 P 84 S
NZ.THZ HH@100 Hz 104 km 48 P 63 S
NZ.CMWZ EH@100 Hz 144 km 51 P 58 S
NZ.KHZ HH@100 Hz 52 km 45 P 58 S
reference: 357 P, 458 S manual picks on these stations
quakescope2026 2617 picks at 0.3
jma_wc 2844 picks at 0.3
original 3383 picks at 0.3
instance 1522 picks at 0.3 Norcia 2016 (M6.5, largest of the Amatrice-Visso-Norcia sequence; permanent plus temporary stations)
88 candidate stations, 6 chosen by manual-pick count with data:
IV.GUMA HH@100 Hz 32 km 134 P 127 S
XO.AM05 EH@100 Hz 26 km 123 P 112 S
IV.LNSS HH@100 Hz 26 km 115 P 113 S
IV.RM33 EH@100 Hz 37 km 111 P 103 S
IV.CAMP HH@100 Hz 41 km 113 P 100 S
IV.MC2 EH@100 Hz 11 km 105 P 101 S
reference: 701 P, 656 S manual picks on these stations
quakescope2026 2408 picks at 0.3
jma_wc 2537 picks at 0.3
original 3446 picks at 0.3
instance 1118 picks at 0.3 Thessaly 2021 (M6.3, doublet with an M6.0 the next day; dense Greek network with a station 5 km away)
12 candidate stations, 6 chosen by manual-pick count with data:
HL.THL HH@100 Hz 26 km 73 P 73 S
HT.TYRN HH@100 Hz 5 km 64 P 64 S
HT.LIT HH@100 Hz 46 km 62 P 55 S
HL.KZN HH@100 Hz 71 km 66 P 51 S
HT.AGG HH@100 Hz 82 km 64 P 45 S
HL.TETR HH@100 Hz 91 km 66 P 43 S
reference: 395 P, 331 S manual picks on these stations
quakescope2026 2553 picks at 0.3
jma_wc 2777 picks at 0.3
original 3877 picks at 0.3
instance 1454 picks at 0.3
4a. What the three services' QuakeML actually says¶
Inspected rather than assumed, as in the US notebook. mode is the pick's
evaluationMode; method is whatever the pick's methodID ends in;
agency is the pick's creationInfo.agencyID. Only manual picks enter
the reference.
rows = []
for label, df in harvested.items():
if not len(df):
continue
res = results.get(label)
lo, hi = pd.Timestamp(res["t0"].datetime, tz="UTC"), pd.Timestamp(res["t1"].datetime, tz="UTC")
w = df[(df["time"] >= lo) & (df["time"] <= hi)]
rows.append(dict(
sequence=label, events=w.event.nunique(), arrivals=len(w),
mode=dict(w["mode"].value_counts()), status=dict(w["status"].value_counts().head(3)),
method=dict(w["method"].value_counts().head(3)), agency=dict(w["agency"].value_counts().head(2)),
onset=dict(w["onset"].value_counts().head(3)),
uncertainty_set=f"{w.uncertainty.notna().mean():.0%}",
zero_weight=int((w.time_weight == 0).sum()),
))
prov = pd.DataFrame(rows).set_index("sequence")
with pd.option_context("display.max_colwidth", 80, "display.width", 200):
print(prov.T.to_string())
sequence Kaikoura 2016 Norcia 2016 Thessaly 2021
events 138 150 74
arrivals 7134 8873 3570
mode {'manual': 5308, 'automatic': 1826} {'manual': 8873} {'manual': 3570}
status {'None': 7134} {'None': 8873} {'None': 3570}
method {'': 5308, 'AIC': 1826} {'query?methodID=11': 8873} {'': 3570}
agency {'WEL(GNS_Primary)': 4806} {'INGV ': 8873} {'NOAreg1': 3570}
onset {'None': 7134} {'questionable': 8873} {'None': 3570}
uncertainty_set 0% 100% 0%
zero_weight 1845 617 37
5. Recall against analyst picks¶
At the shared 0.3 threshold first, because that is how everyone reads a picker. Section 7 is the fair version.
rows = []
for label, res in results.items():
for name in names:
view = at_threshold(res["picks"][name], REPORT_THRESHOLD)
for phase in ("P", "S"):
n_ref = n_hit = n_extra = 0; residuals = []
for sta in res["streams"]:
ref = res["reference"].get((sta, phase), [])
got = view.get((sta, phase), [])
r, extra = match(ref, got)
n_ref += len(ref); n_hit += len(r); n_extra += extra; residuals += r
if n_ref == 0:
continue
rows.append(dict(sequence=label, weights=name, phase=phase, analyst=n_ref, matched=n_hit,
recall=round(n_hit / n_ref, 3),
MAE=round(float(np.mean(np.abs(residuals))), 3) if residuals else np.nan,
bias=round(float(np.median(residuals)), 3) if residuals else np.nan,
extra=n_extra))
bench = pd.DataFrame(rows)
print(bench.to_string(index=False))
for phase in ("P", "S"):
piv = bench[bench.phase == phase].pivot(index="sequence", columns="weights", values="recall")
print(f"\n{phase} recall at {REPORT_THRESHOLD}")
print(piv.reindex(columns=[n for n in names if n in piv.columns]).to_string())
sequence weights phase analyst matched recall MAE bias extra Kaikoura 2016 quakescope2026 P 357 236 0.661 0.139 0.077 1390 Kaikoura 2016 quakescope2026 S 458 230 0.502 0.162 0.034 761 Kaikoura 2016 jma_wc P 357 251 0.703 0.139 0.062 1564 Kaikoura 2016 jma_wc S 458 238 0.520 0.163 -0.026 791 Kaikoura 2016 original P 357 235 0.658 0.154 0.100 1410 Kaikoura 2016 original S 458 323 0.705 0.177 0.089 1415 Kaikoura 2016 instance P 357 230 0.644 0.141 0.044 680 Kaikoura 2016 instance S 458 244 0.533 0.166 -0.039 368 Norcia 2016 quakescope2026 P 701 535 0.763 0.049 0.000 725 Norcia 2016 quakescope2026 S 656 488 0.744 0.064 0.010 660 Norcia 2016 jma_wc P 701 557 0.795 0.048 0.000 786 Norcia 2016 jma_wc S 656 510 0.777 0.069 -0.020 684 Norcia 2016 original P 701 516 0.736 0.057 0.028 1008 Norcia 2016 original S 656 595 0.907 0.085 0.050 1327 Norcia 2016 instance P 701 433 0.618 0.056 -0.010 130 Norcia 2016 instance S 656 439 0.669 0.076 -0.040 116 Thessaly 2021 quakescope2026 P 395 284 0.719 0.129 0.021 1355 Thessaly 2021 quakescope2026 S 331 129 0.390 0.129 0.034 785 Thessaly 2021 jma_wc P 395 295 0.747 0.135 0.021 1431 Thessaly 2021 jma_wc S 331 133 0.402 0.132 0.024 918 Thessaly 2021 original P 395 286 0.724 0.139 0.053 1552 Thessaly 2021 original S 331 142 0.429 0.146 0.072 1897 Thessaly 2021 instance P 395 286 0.724 0.139 0.009 646 Thessaly 2021 instance S 331 130 0.393 0.148 -0.021 392 P recall at 0.3 weights quakescope2026 jma_wc original instance sequence Kaikoura 2016 0.661 0.703 0.658 0.644 Norcia 2016 0.763 0.795 0.736 0.618 Thessaly 2021 0.719 0.747 0.724 0.724 S recall at 0.3 weights quakescope2026 jma_wc original instance sequence Kaikoura 2016 0.502 0.520 0.705 0.533 Norcia 2016 0.744 0.777 0.907 0.669 Thessaly 2021 0.390 0.402 0.429 0.393
seqs = [s for s in SEQUENCES if s in set(bench.sequence)]
fig, axes = plt.subplots(1, 2, figsize=(12, 4.2), sharey=True)
for ax, phase in zip(axes, ("P", "S")):
sub = bench[bench.phase == phase]
width = 0.8 / max(len(names), 1)
for i, name in enumerate(names):
xs, vals = [], []
for j, s in enumerate(seqs):
row = sub[(sub.sequence == s) & (sub.weights == name)]
if len(row):
xs.append(j + (i - (len(names) - 1) / 2) * width); vals.append(float(row.recall.iloc[0]))
ax.bar(xs, vals, width=width * 0.92, color=COLORS[i], label=name)
for x, v in zip(xs, vals):
ax.text(x, v + 0.015, f"{v:.2f}", ha="center", fontsize=7.5)
ax.set_xticks(range(len(seqs))); ax.set_xticklabels(seqs, fontsize=9)
ax.set_ylim(0, 1.12); ax.grid(axis="y")
ax.set_title(f"{phase} recall at the shared {REPORT_THRESHOLD} threshold", fontsize=10.5, loc="left")
axes[0].set_ylabel("recall against manual analyst picks")
axes[1].legend(frameon=False, fontsize=8.5, ncol=2, loc="upper right")
fig.tight_layout(); plt.show()
6. Timing¶
The v7 fine-tune was selected for arrival-time precision, not for recall. If that carries outside the training domain it should show here as a tighter residual on the picks both weights find. The residual is model pick minus analyst pick on matched picks at the shared threshold, pooled over the three sequences; a positive median means the model picks late.
fig, axes = plt.subplots(1, 2, figsize=(12, 3.8))
bins = np.arange(-MATCH_TOL, MATCH_TOL + 0.025, 0.025)
for ax, phase in zip(axes, ("P", "S")):
for i, name in enumerate(names):
res_all = []
for label, res in results.items():
view = at_threshold(res["picks"][name], REPORT_THRESHOLD)
for sta in res["streams"]:
r, _ = match(res["reference"].get((sta, phase), []), view.get((sta, phase), []))
res_all += r
if not res_all:
continue
res_all = np.array(res_all)
ax.hist(res_all, bins=bins, histtype="step", lw=1.8, color=COLORS[i], density=True,
label=f"{name}: median {np.median(res_all):+.3f} s, MAE {np.mean(np.abs(res_all)):.3f} s, n={len(res_all)}")
ax.axvline(0, color="#8a8a8a", lw=1); ax.set_yticks([])
ax.set_xlabel(f"{phase} pick minus analyst pick (s)")
ax.set_title(f"{phase} residuals, all three sequences", fontsize=10.5, loc="left")
ax.legend(frameon=False, fontsize=7.5, loc="upper left")
fig.tight_layout(); plt.show()
print("MAE (s) at the shared threshold, per sequence")
print(bench.pivot_table(index=["phase", "sequence"], columns="weights", values="MAE").reindex(columns=names).round(3).to_string())
MAE (s) at the shared threshold, per sequence
weights quakescope2026 jma_wc original instance
phase sequence
P Kaikoura 2016 0.139 0.139 0.154 0.141
Norcia 2016 0.049 0.048 0.057 0.056
Thessaly 2021 0.129 0.135 0.139 0.139
S Kaikoura 2016 0.162 0.163 0.177 0.166
Norcia 2016 0.064 0.069 0.085 0.076
Thessaly 2021 0.129 0.132 0.146 0.148
7. Matched pick budgets¶
A threshold is not an operating point. Each model's probabilities sit on their own scale, so at a shared 0.3 one emits more picks than another and collects recall for it. Holding the number of picks emitted fixed and asking which weight set recovers more of the analyst catalogue is the comparison that survives a change of threshold, and it is what the US comparison found erased most of the apparent ranking there.
sweep_rows = []
for label, res in results.items():
for name in names:
for thr in THRESHOLD_SWEEP:
view = at_threshold(res["picks"][name], thr)
for phase in ("P", "S"):
hit = tot = emitted = 0
for sta in res["streams"]:
ref = res["reference"].get((sta, phase), []); got = view.get((sta, phase), [])
r, _ = match(ref, got); hit += len(r); tot += len(ref); emitted += len(got)
if tot < 20:
continue
sweep_rows.append(dict(sequence=label, weights=name, phase=phase, thr=thr,
recall=hit / tot, emitted=emitted))
sweep = pd.DataFrame(sweep_rows)
def matched_budget(sweep, phase, sequence, n_points=4):
sub = sweep[(sweep.phase == phase) & (sweep.sequence == sequence)]
present = [n for n in names if n in set(sub.weights)]
if len(present) < 2:
return None
lo = max(sub[sub.weights == n].emitted.min() for n in present)
hi = min(sub[sub.weights == n].emitted.max() for n in present)
if not np.isfinite([lo, hi]).all() or hi <= lo:
ceiling = min(present, key=lambda n: sub[sub.weights == n].emitted.max())
print(f"{sequence} {phase}: no common budget - {ceiling} tops out at "
f"{int(sub[sub.weights == ceiling].emitted.max())} picks, below where the others start.")
return None
rows = []
for target in np.linspace(lo, hi, n_points):
row = {"picks_emitted": int(round(target))}
for n in present:
d = sub[sub.weights == n].sort_values("emitted")
row[n] = round(float(np.interp(target, d.emitted, d.recall)), 3)
rows.append(row)
return pd.DataFrame(rows)
budget_tables = {}
for phase in ("S", "P"):
for sequence in SEQUENCES:
tab = matched_budget(sweep, phase, sequence)
if tab is None:
continue
budget_tables[(sequence, phase)] = tab
print(f"{sequence} - {phase} recall at matched pick budgets"); print(tab.to_string(index=False)); print()
Kaikoura 2016 - S recall at matched pick budgets
picks_emitted quakescope2026 jma_wc original instance
507 0.259 0.263 0.236 0.467
771 0.377 0.382 0.359 0.609
1034 0.516 0.522 0.488 0.697
1298 0.603 0.621 0.598 0.738
Norcia 2016 - S recall at matched pick budgets
picks_emitted quakescope2026 jma_wc original instance
816 0.591 0.613 0.556 0.813
896 0.634 0.653 0.594 0.841
975 0.668 0.686 0.631 0.858
1055 0.703 0.719 0.668 0.875
Thessaly 2021 - S recall at matched pick budgets
picks_emitted quakescope2026 jma_wc original instance
708 0.358 0.350 0.302 0.416
945 0.392 0.386 0.344 0.433
1181 0.410 0.413 0.373 0.456
1418 0.428 0.432 0.397 0.474
Kaikoura 2016 - P recall at matched pick budgets
picks_emitted quakescope2026 jma_wc original instance
756 0.407 0.439 0.409 0.580
1061 0.504 0.546 0.539 0.669
1367 0.608 0.624 0.611 0.719
1672 0.667 0.681 0.660 0.739
Norcia 2016 - P recall at matched pick budgets
picks_emitted quakescope2026 jma_wc original instance
691 0.557 0.578 0.534 0.706
837 0.628 0.648 0.595 0.755
982 0.691 0.705 0.646 0.785
1128 0.736 0.751 0.680 0.807
Thessaly 2021 - P recall at matched pick budgets
picks_emitted quakescope2026 jma_wc original instance
1014 0.573 0.612 0.529 0.730
1227 0.636 0.669 0.613 0.748
1439 0.683 0.705 0.666 0.760
1652 0.720 0.736 0.693 0.762
seqs = [s for s in SEQUENCES if s in set(sweep.sequence)]
fig, axes = plt.subplots(2, len(seqs), figsize=(4.6 * len(seqs), 7.6), squeeze=False)
for row, phase in enumerate(("P", "S")):
for ax, sequence in zip(axes[row], seqs):
sub = sweep[(sweep.phase == phase) & (sweep.sequence == sequence)]
for name, color in zip(names, COLORS):
d = sub[sub.weights == name].sort_values("emitted")
if not len(d):
continue
ax.plot(d.emitted, d.recall, marker="o", ms=3.5, lw=1.8, color=color, label=name)
star = d[np.isclose(d.thr, REPORT_THRESHOLD)]
if len(star):
ax.plot(star.emitted, star.recall, marker="*", ms=13, color=color, mec="#16150f", mew=0.6, zorder=5)
ax.set_title(f"{sequence} - {phase}", fontsize=10.5, loc="left")
ax.set_xlabel(f"{phase} picks emitted"); ax.set_ylim(0, 1)
axes[row][0].set_ylabel("recall against analyst picks")
axes[0][0].legend(frameon=False, fontsize=8)
fig.suptitle(f"Recall against picks emitted; stars mark the shared {REPORT_THRESHOLD} threshold",
fontsize=10, x=0.01, ha="left")
fig.tight_layout(); plt.show()
sat_rows = []
for label, res in results.items():
for name in names:
row = {"sequence": label, "weights": name}
for thr in (DETECT_FLOOR, 0.1, REPORT_THRESHOLD):
view = at_threshold(res["picks"][name], thr)
row[f"P@{thr}"] = sum(len(v) for (_, ph), v in view.items() if ph == "P")
row[f"S@{thr}"] = sum(len(v) for (_, ph), v in view.items() if ph == "S")
sat_rows.append(row)
print("Picks emitted by threshold; the floor column is each model's ceiling")
print(pd.DataFrame(sat_rows).to_string(index=False))
Picks emitted by threshold; the floor column is each model's ceiling
sequence weights P@0.02 S@0.02 P@0.1 S@0.1 P@0.3 S@0.3
Kaikoura 2016 quakescope2026 4786 3682 2733 2072 1626 991
Kaikoura 2016 jma_wc 4855 3416 2880 2109 1815 1029
Kaikoura 2016 original 3942 3674 2445 2604 1645 1738
Kaikoura 2016 instance 2231 1644 1346 1042 910 612
Norcia 2016 quakescope2026 3618 2837 2031 1917 1260 1148
Norcia 2016 jma_wc 3667 2777 2173 1923 1343 1194
Norcia 2016 original 4074 3312 2506 2574 1524 1922
Norcia 2016 instance 1546 1321 890 878 563 555
Thessaly 2021 quakescope2026 4217 3449 2457 2077 1639 914
Thessaly 2021 jma_wc 4557 3460 2687 2287 1726 1051
Thessaly 2021 original 4135 3888 2603 2867 1838 2039
Thessaly 2021 instance 2117 1821 1382 1113 932 522
The pair that matters: quakescope2026 against jma_wc¶
One number per sequence and phase: the recall difference at a common budget, taken at the midpoint of the overlap. Positive favours the fine-tune.
rows = []
for (sequence, phase), tab in budget_tables.items():
if "quakescope2026" not in tab or "jma_wc" not in tab:
continue
mid = tab.iloc[len(tab) // 2]
at_thr = bench[(bench.sequence == sequence) & (bench.phase == phase)].set_index("weights")
rows.append(dict(sequence=sequence, phase=phase, budget=int(mid.picks_emitted),
quakescope2026=mid["quakescope2026"], jma_wc=mid["jma_wc"],
diff_at_budget=round(mid["quakescope2026"] - mid["jma_wc"], 3),
diff_at_thr=round(at_thr.loc["quakescope2026", "recall"] - at_thr.loc["jma_wc", "recall"], 3),
MAE_q26=at_thr.loc["quakescope2026", "MAE"], MAE_jma=at_thr.loc["jma_wc", "MAE"]))
pair = pd.DataFrame(rows).sort_values(["phase", "sequence"])
print(pair.to_string(index=False))
print(f"\nmean recall difference at matched budget: P {pair[pair.phase == 'P'].diff_at_budget.mean():+.3f}, "
f"S {pair[pair.phase == 'S'].diff_at_budget.mean():+.3f}")
sequence phase budget quakescope2026 jma_wc diff_at_budget diff_at_thr MAE_q26 MAE_jma Kaikoura 2016 P 1367 0.608 0.624 -0.016 -0.042 0.139 0.139 Norcia 2016 P 982 0.691 0.705 -0.014 -0.032 0.049 0.048 Thessaly 2021 P 1439 0.683 0.705 -0.022 -0.028 0.129 0.135 Kaikoura 2016 S 1034 0.516 0.522 -0.006 -0.018 0.162 0.163 Norcia 2016 S 975 0.668 0.686 -0.018 -0.033 0.064 0.069 Thessaly 2021 S 1181 0.410 0.413 -0.003 -0.012 0.129 0.132 mean recall difference at matched budget: P -0.017, S -0.009
Export for the consolidated benchmark¶
Writes the tables above to docs/benchmark/results/ so that tutorials/benchmark_summary.ipynb reads executed output.
# Export the result tables for docs/benchmark/. The consolidated benchmark
# (tutorials/benchmark_summary.ipynb) reads these files rather than numbers
# transcribed from a rendered report, so every figure there traces to an
# executed cell here. Defensive on purpose: a missing table is reported, not fatal.
import datetime as _dt, json as _json
from pathlib import Path as _Path
import seisbench as _sb
_OUT = _Path("../docs/benchmark/results/global_sequences"); _OUT.mkdir(parents=True, exist_ok=True)
_written = []
def _export(name, obj):
try:
(obj.to_csv(_OUT / f"{name}.csv", index=isinstance(obj.index, pd.MultiIndex) or obj.index.name is not None)
if hasattr(obj, "to_csv") else (_OUT / f"{name}.json").write_text(_json.dumps(obj, indent=1, default=str)))
_written.append(name)
except Exception as _exc:
print(f"{name}: not exported ({type(_exc).__name__}: {_exc})")
_export("recall_at_threshold", bench)
_export("threshold_sweep", sweep)
_export("pair_v7_vs_jma", pair)
_export("reference_provenance", prov.reset_index())
# Raw picks, so every metric can be recomputed without re-running inference.
# One row per model pick at the detection floor with its peak confidence, and
# one per reference arrival. docs/benchmark/README.md scores these with
# sb_catalog.src.benchmark_metrics; keeping them means a new metric costs a
# read rather than ten minutes of inference.
_rows = []
for _label, _res in results.items():
for _name, _store in _res["picks"].items():
for (_sta, _pha), _v in _store.items():
for _t, _q in _v:
_rows.append(dict(sequence=_label, weights=_name, station=_sta, phase=_pha,
time=pd.Timestamp(_t.datetime), conf=float(_q)))
_picks = pd.DataFrame(_rows)
_export("model_picks", _picks)
_rows = []
for _label, _res in results.items():
for (_sta, _pha), _v in _res["reference"].items():
for _t in _v:
_rows.append(dict(sequence=_label, station=_sta, phase=_pha,
time=pd.Timestamp(_t.datetime)))
_ref = pd.DataFrame(_rows)
_export("reference_picks", _ref)
print(f"exported {len(_picks):,} model picks at the {DETECT_FLOOR} floor and {len(_ref):,} reference arrivals")
_meta = dict(notebook="phasenet_global_sequences.ipynb",
executed=_dt.datetime.now(_dt.timezone.utc).isoformat(timespec="seconds"),
seisbench=_sb.__version__, weights=WEIGHTS, report_threshold=REPORT_THRESHOLD,
detect_floor=DETECT_FLOOR, match_tol_s=MATCH_TOL, window_start_s=WINDOW_START,
n_stations=N_STATIONS,
sequences={k: dict(time=str(v["time"]), mag=v["mag"], lat=v["lat"], lon=v["lon"],
agency=v["agency"], window_min=v["window_min"],
stations=sorted(results[k]["streams"]) if k in results else [])
for k, v in SEQUENCES.items()})
_export("meta", _meta)
print("wrote", ", ".join(_written), "to", _OUT)
exported 79,026 model picks at the 0.02 floor and 2,898 reference arrivals wrote recall_at_threshold, threshold_sweep, pair_v7_vs_jma, reference_provenance, model_picks, reference_picks, meta to ../docs/benchmark/results/global_sequences
8. Where the global campaign stands on these networks¶
The campaign's station table already lists these operators' networks, and it
will pick them with jma_wc. When it has, its stored picks can be scored
here exactly as the OBS report scores the obs campaign. Until then this
cell only says what is there.
BUCKET, REGION, CAMPAIGN = "quakescope-picks-2026", "us-east-2", "global"
s3 = boto3.client("s3", region_name=REGION)
_pg = s3.get_paginator("list_objects_v2")
station_tbl = pd.read_parquet(io.BytesIO(s3.get_object(Bucket=BUCKET, Key=f"{CAMPAIGN}/stations.parquet")["Body"].read()))
for label, res in results.items():
nets = sorted({s.split(".")[0] for s in res["streams"]})
t0 = res["t0"]
listed = station_tbl[station_tbl.network_code.isin(nets)]
have = 0
for net in nets:
prefix = f"{CAMPAIGN}/picks/network={net}/year={t0.year}/month={t0.month:02d}/"
have += sum(1 for page in _pg.paginate(Bucket=BUCKET, Prefix=prefix) for _ in page.get("Contents", []))
print(f"{label:<14} networks {nets}: {len(listed)} stations in the campaign table, "
f"{have} pick objects for {t0.year}-{t0.month:02d}"
+ ("" if have else " - not picked yet"))
Kaikoura 2016 networks ['NZ']: 1159 stations in the campaign table, 0 pick objects for 2016-11 - not picked yet Norcia 2016 networks ['IV', 'XO']: 0 stations in the campaign table, 0 pick objects for 2016-10 - not picked yet
Thessaly 2021 networks ['HL', 'HT']: 10 stations in the campaign table, 0 pick objects for 2021-03 - not picked yet
9. What the picks look like¶
One record section per sequence from the model with the most picks at the reporting threshold, with the analyst's picks marked above the traces so the recall numbers can be checked by eye.
def record_section(label, res, weight, span=300):
seq = SEQUENCES[label]
order = list(res["streams"])
fig, ax = plt.subplots(figsize=(11.5, 0.9 * len(order) + 1.4))
view = at_threshold(res["picks"][weight], REPORT_THRESHOLD)
for row, sta in enumerate(order):
tr = res["streams"][sta].select(component="Z")
if not tr:
continue
tr = tr[0].copy(); tr.trim(res["t0"], res["t0"] + span)
x = tr.data.astype(float); peak = np.abs(x).max()
if peak > 0:
x = x / peak * 0.38
ax.plot(tr.times(), x + row, color="#3d3d3d", lw=0.45)
ax.annotate(sta, (-span * 0.012, row + 0.1), fontsize=8, color="#52514e", ha="right")
for phase, color, ls in (("P", C_P, "-"), ("S", C_S, "--")):
for t in view.get((sta, phase), []):
dt = t - res["t0"]
if 0 <= dt <= span:
ax.plot([dt, dt], [row - 0.42, row - 0.05], color=color, lw=1.0, ls=ls, alpha=0.9)
for t in res["reference"].get((sta, phase), []):
dt = t - res["t0"]
if 0 <= dt <= span:
ax.plot([dt, dt], [row + 0.05, row + 0.42], color=color, lw=1.0, ls=ls, alpha=0.9)
handles = [plt.Line2D([], [], color=C_P, lw=1.4, label="P"), plt.Line2D([], [], color=C_S, lw=1.4, ls="--", label="S")]
ax.legend(handles=handles, frameon=False, fontsize=9, ncol=2, loc="upper right")
ax.text(0.01, 0.985, "analyst picks above each trace, model picks below", transform=ax.transAxes,
ha="left", va="top", fontsize=7.5, color="#7a7973")
ax.set_xlim(-span * 0.09, span); ax.set_ylim(-0.7, len(order) - 0.3); ax.set_yticks([])
ax.set_xlabel(f"seconds into the aftershock window ({WINDOW_START // 60} min after the mainshock)")
ax.set_title(f"{label} - M{seq['mag']} - picks from {weight} at {REPORT_THRESHOLD}", fontsize=11, loc="left")
ax.grid(axis="x")
fig.tight_layout()
return fig
for label, res in results.items():
best = max(res["picks"], key=lambda n: sum(1 for v in res["picks"][n].values() for t, q in v if q >= REPORT_THRESHOLD))
record_section(label, res, best); plt.show()
Reading the result¶
jma_wc recovers more of the analysts' picks than quakescope2026 on all
three sequences and both phases, however the comparison is framed. At the
shared 0.3 threshold the parent leads by 3–4 points on P and 1–3 on S. At
matched pick budgets, which take the threshold out of the comparison, it still
leads by 1–2 points on P and up to 2 on S, and the sign never flips. The
fine-tune was selected for timing, and there is no timing edge to see: the
mean absolute residuals agree to within 6 ms on every sequence and phase, and
both medians sit within 20 ms of zero. This is the western-states result
again, on the Marlborough faults, the Apennines and the Aegean, against three
operators' analysts. Whatever v7 learned from its training set did not travel.
instance is the most efficient picker on all three, not only in Italy.
At any common budget it recovers more analyst picks than the other three: at
Kaikōura it reaches 0.70 of the S picks with 1,034 emitted where the two
PhaseNetWC weights reach 0.52, and at Thessaly it leads P by 3–12 points along
the whole overlap. Its probabilities sit lower, so at 0.3 it emits half as
many picks as the others and looks worst in section 5, and with the threshold
on the floor it tops out below their ceilings (the saturation table). That is
the ceiling the US comparison found on Ridgecrest, but on these
regional-distance sequences the ceiling is high enough that the weight is
better per pick emitted. Which weight to run globally is therefore a
question about the budget: for a catalogue of the events an operator would
keep, instance at a low threshold; for the most complete pick set, jma_wc.
original leads S at 0.3 for the same reason it did in the US: it emits
twice the S picks. At matched budget it is the weakest of the four on every
sequence, and its S residual is the most biased, a median 55 ms late.
Thessaly's S recall is capped near 0.5 for every weight at every budget. Nothing on the sweep exceeds 0.49. A ceiling shared by four unrelated weight sets belongs to the reference or the data, not to a model: the NOA analysts pick S on stations 45–90 km out where every model finds the P and not the S. The cause is not identified here, and the number should not be quoted without that caveat.
Provenance is uneven across operators, and it shapes the reference.
GeoNet's QuakeML carries a quarter of its picks as automatic from an AIC
picker, which the reference excludes, and 1,845 arrivals with zero weight in
the location, which the reference keeps as the US notebook does. INGV marks
every pick manual, every onset questionable and every uncertainty set; NOA
every pick manual and nothing else. None publishes an analyst name. The
reference is the operator's bulletin rather than a research catalogue, and the
harvest is cached so that a re-run scores against the same picks.
The global campaign has not reached these networks. NZ, HL and HT are in
its station table (1,159 New Zealand stations, ten Greek); IV is not, because
INGV's data are not on EarthScope. Section 8 will score the stored picks once
they exist. Until then the recommendation for the global deployment rests on
this notebook and the western one, and both say the same thing: do not move
the campaign from jma_wc to quakescope2026, and consider instance at a
lower threshold where the target is the operator's catalogue rather than
everything the data hold.