Phase-picker benchmark — Ridgecrest aftershocks¶
Measures how much of the analyst-picked catalog each PhaseNet weight set recovers, on a dense half hour of the 2019 Ridgecrest aftershock sequence. The question driving it: do the fine-tuned weights recover more S arrivals than the stock models?
S is the interesting phase. It is the harder pick, it is what constrains depth and location, and it is where the fine-tuning trade-offs show up.
Why this window¶
Twenty-one M3.2–4.3 aftershocks in thirty minutes, several separated by only a few seconds, so their waveforms overlap and a picker cannot treat each event as an isolated snippet. The mainshock itself is deliberately excluded — a magnitude 7 ruptures for tens of seconds and buries its own S arrival, which tests something other than picking skill.
Ground truth, and what it is not¶
SCEDC analyst picks are the reference. They are authoritative for what they contain, but they are not exhaustive — in a sequence this dense no analyst picks every arrival on every station. So:
- Recall is meaningful. Of the arrivals a human marked, how many did the model find?
- Precision is not computable here. A model pick with no analyst counterpart may be a false positive or a real arrival the analyst skipped. This notebook reports those as extra detections and never calls them false positives.
Separating those two cases needs either exhaustive re-picking or an association step, neither of which belongs in a benchmark this size.
import io
from collections import defaultdict
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 s3fs import S3FileSystem
%matplotlib inline
1. Configuration¶
# --- Benchmark window: 21 aftershocks, M3.2-4.3, several overlapping ---
T0 = UTCDateTime("2019-07-06T05:00:00")
T1 = T0 + 1800 # 30 minutes
YEAR, DOY = 2019, 187
REGION = dict(minlatitude=35.5, maxlatitude=36.1,
minlongitude=-117.9, maxlongitude=-117.3)
MIN_MAG = 3.2
# Stations carrying the most analyst picks in this window (network CI)
STATIONS = ["CLC", "WCS2", "MPM", "WBS", "JRC2", "WOR", "TEH", "LRL"]
NETWORK, CHANNEL = "CI", "HH"
# Weight sets to compare. Missing ones are skipped with a message.
WEIGHTS = [
"quakescope2026", # v7 fine-tune, this project
"jma_wc", # the parent v7 was distilled from
"original", # Zhu & Beroza (2019)
"instance", # stock baseline
]
P_THRESHOLD = S_THRESHOLD = 0.3
MATCH_TOL = 0.5 # seconds; a pick counts as recovered within this
# Categorical colors, CVD-validated. Every chart also carries direct value
# labels, which the palette requires and which makes color non-essential.
COLORS = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100"]
2. Analyst picks from SCEDC¶
Each event carries arrivals that reference picks; the arrival supplies
the phase label and the pick supplies the time and station.
client = Client("SCEDC")
catalog = client.get_events(
starttime=T0, endtime=T1, minmagnitude=MIN_MAG,
includearrivals=True, **REGION,
)
print(f"{len(catalog)} events, M>={MIN_MAG}, in {(T1 - T0) / 60:.0f} minutes")
# analyst[(station, phase)] -> sorted list of UTCDateTime
analyst = defaultdict(list)
origins = []
for ev in catalog:
o = ev.preferred_origin() or ev.origins[0]
mag = ev.preferred_magnitude() or ev.magnitudes[0]
origins.append((o.time, mag.mag))
by_id = {p.resource_id.id: p for p in ev.picks}
for arr in o.arrivals:
pick = by_id.get(arr.pick_id.id)
if pick is None or not arr.phase:
continue
phase = arr.phase[0].upper() # Pg/Pn -> P, Sg/Sn -> S
if phase not in ("P", "S"):
continue
wid = pick.waveform_id
if wid.network_code != NETWORK or wid.station_code not in STATIONS:
continue
analyst[(wid.station_code, phase)].append(pick.time)
for key in analyst:
analyst[key] = sorted(analyst[key])
origins.sort()
n_p = sum(len(v) for (_, ph), v in analyst.items() if ph == 'P')
n_s = sum(len(v) for (_, ph), v in analyst.items() if ph == 'S')
print(f"analyst picks on the {len(STATIONS)} benchmark stations: {n_p} P, {n_s} S")
21 events, M>=3.2, in 30 minutes analyst picks on the 8 benchmark stations: 125 P, 144 S
How much do these events overlap?¶
Gaps between consecutive origins, against the ~20–30 s an S coda lasts at these distances. Anything below that means two events share a window.
gaps = np.diff([t for t, _ in origins])
print(f"median gap between origins : {np.median(gaps):.0f} s")
print(f"gaps under 30 s : {(gaps < 30).sum()} of {len(gaps)}")
print(f"shortest gap : {gaps.min():.0f} s")
fig, ax = plt.subplots(figsize=(11, 2.2))
for t, m in origins:
ax.vlines(t - T0, 0, m, color='#3d3d3d', lw=1.4)
ax.set_xlim(0, T1 - T0)
ax.set_ylim(3.0, 4.6)
ax.set_xlabel('Seconds into the benchmark window')
ax.set_ylabel('Magnitude')
ax.set_title(f'{len(origins)} aftershocks, M>={MIN_MAG}', fontsize=11, loc='left')
ax.grid(alpha=0.25, lw=0.5)
fig.tight_layout()
plt.show()
median gap between origins : 56 s gaps under 30 s : 5 of 20 shortest gap : 4 s
3. Waveforms¶
One 30-minute window per station from the SCEDC public bucket. Whole-day files are ~20 MB per channel, so this fetches roughly 500 MB the first time and takes a couple of minutes.
def scedc_key(sta, comp, loc=""):
base = (f"{NETWORK}{sta.ljust(5, '_')}{CHANNEL}{comp}"
f"{loc.ljust(3, '_')}{YEAR}{DOY:03d}.ms")
return f"scedc-pds/continuous_waveforms/{YEAR}/{YEAR}_{DOY:03d}/{base}"
fs = S3FileSystem(anon=True)
streams = {}
for sta in STATIONS:
st = obspy.Stream()
for comp in "ZNE":
try:
with fs.open(scedc_key(sta, comp)) as fh:
st += obspy.read(io.BytesIO(fh.read()))
except Exception as exc:
print(f" {sta} {comp}: {type(exc).__name__}")
if len(st) == 0:
print(f" {sta}: no data")
continue
st.merge(fill_value=0)
st.trim(T0, T1)
streams[sta] = st
print(f" {sta:<6} {len(st)} traces, {st[0].stats.npts} samples")
print(f"\n{len(streams)}/{len(STATIONS)} stations loaded")
CLC 3 traces, 180001 samples
WCS2 3 traces, 180001 samples
MPM 3 traces, 180001 samples
WBS 3 traces, 180001 samples
JRC2 3 traces, 180001 samples
WOR 3 traces, 180001 samples
TEH 3 traces, 180001 samples
LRL 3 traces, 180001 samples 8/8 stations loaded
4. Run each weight set¶
models = {}
available = sbm.PhaseNet.list_pretrained()
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")
# picks[model][(station, phase)] -> sorted list of UTCDateTime
picks = {name: defaultdict(list) for name in models}
for name, model in models.items():
for sta, st in streams.items():
out = model.classify(st, P_threshold=P_THRESHOLD, S_threshold=S_THRESHOLD)
for p in out.picks:
picks[name][(sta, p.phase)].append(p.peak_time)
for key in picks[name]:
picks[name][key] = sorted(picks[name][key])
tot_p = sum(len(v) for (_, ph), v in picks[name].items() if ph == 'P')
tot_s = sum(len(v) for (_, ph), v in picks[name].items() if ph == 'S')
print(f"{name:<16} {tot_p:4d} P, {tot_s:4d} S")
quakescope2026 loaded jma_wc loaded original loaded instance loaded
quakescope2026 530 P, 435 S
jma_wc 519 P, 448 S
original 629 P, 710 S
instance 255 P, 239 S
5. Match model picks to analyst picks¶
Greedy nearest-neighbour within the tolerance, each model pick consumed at most once, so two analyst picks cannot both claim the same detection.
def match(analyst_times, model_times, tol=MATCH_TOL):
"""Return (residuals of matched picks, count of unmatched model picks)."""
used = set()
residuals = []
for a in analyst_times:
best_i, best_d = None, None
for i, m in enumerate(model_times):
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(model_times) - len(used)
rows, residual_store = [], {}
for name in models:
for phase in ("P", "S"):
n_analyst = n_matched = n_extra = 0
res_all = []
for sta in streams:
a = analyst.get((sta, phase), [])
m = picks[name].get((sta, phase), [])
res, extra = match(a, m)
n_analyst += len(a)
n_matched += len(res)
n_extra += extra
res_all += res
residual_store[(name, phase)] = res_all
rows.append(dict(
weights=name, phase=phase,
analyst=n_analyst, matched=n_matched,
recall=round(n_matched / n_analyst, 3) if n_analyst else np.nan,
MAE=round(float(np.mean(np.abs(res_all))), 3) if res_all else np.nan,
bias=round(float(np.median(res_all)), 3) if res_all else np.nan,
extra=n_extra,
))
bench = pd.DataFrame(rows)
print(bench.to_string(index=False))
print(f"\nmatch tolerance: +/-{MATCH_TOL} s")
print("recall = analyst picks recovered / analyst picks")
print("bias = median (model - analyst); positive means the model picks late")
print("extra = model picks with no analyst counterpart - NOT false positives,")
print(" the analyst catalog is not exhaustive in a sequence this dense")
weights phase analyst matched recall MAE bias extra
quakescope2026 P 125 116 0.928 0.034 -0.012 414
quakescope2026 S 144 114 0.792 0.045 -0.012 321
jma_wc P 125 117 0.936 0.030 -0.012 402
jma_wc S 144 116 0.806 0.053 -0.032 332
original P 125 109 0.872 0.037 0.018 520
original S 144 125 0.868 0.056 0.028 585
instance P 125 98 0.784 0.036 -0.022 157
instance S 144 105 0.729 0.080 -0.062 134
match tolerance: +/-0.5 s
recall = analyst picks recovered / analyst picks
bias = median (model - analyst); positive means the model picks late
extra = model picks with no analyst counterpart - NOT false positives,
the analyst catalog is not exhaustive in a sequence this dense
S recall, the headline number¶
Bars carry their values directly, so the comparison does not depend on reading colour.
fig, axes = plt.subplots(1, 2, figsize=(11, 4), sharey=True)
names = list(models)
for ax, phase in zip(axes, ('S', 'P')):
sub = bench[bench.phase == phase].set_index('weights').loc[names]
bars = ax.bar(range(len(names)), sub['recall'],
color=COLORS[:len(names)], width=0.62)
for i, (v, m) in enumerate(zip(sub['recall'], sub['matched'])):
ax.text(i, v + 0.015, f'{v:.2f}', ha='center', fontsize=10)
ax.text(i, 0.02, f'{m}/{sub["analyst"].iloc[i]}', ha='center',
fontsize=8, color='#52514e')
ax.set_xticks(range(len(names)))
ax.set_xticklabels(names, rotation=20, ha='right', fontsize=9)
ax.set_title(f'{phase} recall', fontsize=11, loc='left')
ax.grid(alpha=0.25, lw=0.5, axis='y')
ax.set_axisbelow(True)
axes[0].set_ylabel('Fraction of analyst picks recovered')
axes[0].set_ylim(0, 1.05)
fig.tight_layout()
plt.show()
Timing of the picks that were recovered¶
Residual is model minus analyst, so negative means the model picks early. Recall and timing trade against each other, and a model that recovers more arrivals by loosening its threshold usually pays for it here.
fig, axes = plt.subplots(1, 2, figsize=(11, 3.6), sharey=True)
bins = np.linspace(-MATCH_TOL, MATCH_TOL, 25)
for ax, phase in zip(axes, ('S', 'P')):
for name, color in zip(names, COLORS):
res = residual_store[(name, phase)]
if not res:
continue
ax.hist(res, bins=bins, histtype='step', lw=1.8,
color=color, label=f'{name} (n={len(res)})')
ax.axvline(0, color='#8a8a8a', lw=1)
ax.set_title(f'{phase} residuals', fontsize=11, loc='left')
ax.set_xlabel('model - analyst (s)')
ax.grid(alpha=0.25, lw=0.5)
ax.legend(fontsize=8, frameon=False)
axes[0].set_ylabel('picks')
fig.tight_layout()
plt.show()
Per-station S recall¶
A model that wins on aggregate but collapses at particular stations is a different proposition from one that is uniformly decent.
rows = []
for sta in streams:
row = {'station': sta, 'analyst_S': len(analyst.get((sta, 'S'), []))}
for name in names:
a = analyst.get((sta, 'S'), [])
m = picks[name].get((sta, 'S'), [])
res, _ = match(a, m)
row[name] = round(len(res) / len(a), 2) if a else np.nan
rows.append(row)
per_sta = pd.DataFrame(rows)
print(per_sta.to_string(index=False))
print("\nvalues are S recall per station")
station analyst_S quakescope2026 jma_wc original instance
CLC 20 0.85 0.85 0.85 0.60
WCS2 18 0.83 0.83 0.89 0.83
MPM 18 0.72 0.78 0.78 0.72
WBS 18 0.72 0.78 0.89 0.72
JRC2 15 0.87 0.87 0.93 0.80
WOR 19 0.74 0.79 0.84 0.79
TEH 18 0.89 0.83 0.89 0.67
LRL 18 0.72 0.72 0.89 0.72
values are S recall per station
6. Sensitivity to the matching tolerance¶
A single tolerance can flatter whichever model happens to sit just inside it, so the ranking should hold as the window widens.
tols = [0.25, 0.5, 1.0, 2.0]
rows = []
for tol in tols:
row = {'tolerance_s': tol}
for name in names:
tot = hit = 0
for sta in streams:
a = analyst.get((sta, 'S'), [])
m = picks[name].get((sta, 'S'), [])
res, _ = match(a, m, tol=tol)
tot += len(a)
hit += len(res)
row[name] = round(hit / tot, 3) if tot else np.nan
rows.append(row)
sens = pd.DataFrame(rows)
print(sens.to_string(index=False))
print("\nS recall as a function of match tolerance")
tolerance_s quakescope2026 jma_wc original instance
0.25 0.792 0.799 0.868 0.715
0.50 0.792 0.806 0.868 0.729
1.00 0.792 0.806 0.868 0.729
2.00 0.792 0.806 0.868 0.729
S recall as a function of match tolerance
7. Reading the result¶
S recall answers the question this benchmark was built for. Compare
quakescope2026 against instance for the stock-versus-fine-tuned
question, and against jma_wc for what the fine-tuning itself changed,
since that is the model it was distilled from.
Extra detections are not errors. A model with high recall and many extra picks may be finding the aftershocks the analysts did not have time to pick — this sequence had hundreds of events per hour. Deciding that needs association: real arrivals line up across stations into locatable events, noise does not. That is the natural follow-on, and QuakeScope already runs PyOcto for it.
One window, one region, one magnitude band. Nothing here generalises to teleseismic distances, other networks, or quiet periods. It is a targeted measurement of S recovery under dense overlapping aftershocks, which is a regime that matters for catalog completeness and one where pickers differ sharply.
For the cross-domain benchmark these weights were selected on — and the
caveats attached to that selection — see
docs/phasenet_v7_model_description.md.