Offshore pickers on ocean-bottom seismometers¶
QuakeScope will process ocean-bottom data, and ocean-bottom data is not land data with more noise. The instrument sits under a water column that reverberates, often on soft sediment, with horizontals that are arbitrarily oriented because nobody was there to point them north. Most stations also carry a hydrophone, a channel that has no land equivalent at all.
SeisBench ships three ocean-bottom pickers, and the differences between them are the interesting part:
| Model | Backbone | Components | Hydrophone |
|---|---|---|---|
PickBlue(base="phasenet") |
PhaseNet | Z12H |
yes |
PickBlue(base="eqtransformer") |
EQTransformer | Z12H |
yes |
OBSTransformer (Niksejel & Zhang, 2024) |
EQTransformer | ZNE |
no |
PickBlue is a constructor rather than a model — it returns the obs
weights on whichever backbone you ask for. OBSTransformer is the useful
control: trained on ocean-bottom data but taking only three components,
which separates trained offshore from uses the fourth channel.
Two land models quantify what running them offshore costs. All five see identical windows on three deployments in contrasting settings.
| Deployment | Network | When | Setting |
|---|---|---|---|
| Cascadia Initiative | 7D | 2012–13 | Subduction margin, offshore northern California |
| AACSE | XO | 2018–19 | Alaska Peninsula subduction, shelf to trench |
| Blanco | X9 | 2012–13 | Oceanic transform, strike-slip in young crust |
How sections 5–7 score detection, and why sections 8–11 do it differently¶
Regional catalogs do not pick temporary OBS deployments, so there are no analyst arrivals at these stations to score against. Instead each catalog event gets a predicted P arrival from iasp91, and a model counts as having detected it if it places a P pick within a tolerance of that time.
That prediction is the weak link and the tolerance has to absorb it: iasp91 has no water layer and no sediment column, both of which delay the true arrival relative to the model. Treat the numbers as relative between weight sets on identical data, not as absolute detection rates.
Since this was first written, published picks have turned up for two of
the three settings. Sections 8–11 use them: the obs campaign's stored
picks are scored against analyst-checked arrivals on the AACSE array and
against the University of Washington's real-time catalogue at Axial
Seamount, and the five models above are re-scored on AACSE windows against
analyst picks instead of a prediction. Blanco stays as it is - the
published picks there are available on request, not for download.
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 obspy.geodetics import gps2dist_azimuth, locations2degrees
from obspy.taup import TauPyModel
%matplotlib inline
1. Configuration¶
# Three deployments, chosen for contrasting settings and instrumentation.
# Stations were selected by confirming all four components actually return
# data, not from metadata - and the hydrophone sampling rate differs by an
# order of magnitude between experiments, which section 6 tests the effect of.
EXPERIMENTS = {
"Cascadia (7D)": dict(
network="7D", start=UTCDateTime("2012-09-01"), end=UTCDateTime("2013-05-01"),
lat=40.3, lon=-124.8, min_mag=3.5, max_radius_deg=2.0,
stations=[("FS07B", "HH", "HDH"), ("FS06B", "BH", "BDH"), ("FS05B", "BH", "BDH")],
note="Cascadia subduction margin, offshore northern California",
),
"AACSE (XO)": dict(
network="XO", start=UTCDateTime("2018-08-01"), end=UTCDateTime("2019-06-01"),
lat=54.8, lon=-155.5, min_mag=3.5, max_radius_deg=2.0,
stations=[("LA23", "HH", "EDH"), ("LA25", "HH", "EDH"), ("LD36", "HH", "HDH")],
note="Alaska Peninsula subduction zone, shelf and trench",
),
"Blanco (X9)": dict(
network="X9", start=UTCDateTime("2012-10-01"), end=UTCDateTime("2013-08-01"),
lat=43.1, lon=-126.4, min_mag=3.5, max_radius_deg=2.0,
stations=[("BB060", "HH", "BDH"), ("BB090", "HH", "BDH"), ("BB070", "HH", "BDH")],
note="Blanco oceanic transform - strike-slip in young oceanic crust",
),
}
# SeisBench ships three ocean-bottom pickers and they differ in ways that
# matter here.
#
# PickBlue is a constructor, not a model: PickBlue(base=...) returns the
# "obs" weights on either a PhaseNet or an EQTransformer backbone. Both take
# four components, Z12H, the last being the hydrophone.
#
# OBSTransformer (Niksejel & Zhang, 2024) is also OBS-trained but takes only
# three components - no hydrophone. It is the control that separates
# "trained on ocean-bottom data" from "uses the fourth channel".
#
# The land models quantify what running them offshore costs.
MODEL_SPECS = {
"pickblue_phasenet": dict(build=lambda: sbm.PickBlue(base="phasenet"),
hydrophone=True),
"pickblue_eqt": dict(build=lambda: sbm.PickBlue(base="eqtransformer"),
hydrophone=True),
"obstransformer": dict(build=lambda: sbm.OBSTransformer.from_pretrained("obst2024"),
hydrophone=False),
"quakescope2026": dict(build=lambda: sbm.PhaseNet.from_pretrained("quakescope2026"),
hydrophone=False),
"original": dict(build=lambda: sbm.PhaseNet.from_pretrained("original"),
hydrophone=False),
}
# Detection is scored against a predicted arrival rather than an analyst pick,
# because regional catalogs do not pick temporary OBS stations. iasp91 has no
# water layer and no sediments, so the prediction is systematically early at
# these sites and the tolerance has to be generous.
ABLATION_MODEL = "pickblue_phasenet" # the 4-component model used in section 6
GALLERY_N = 10 # windows shown per deployment in section 7
TOLERANCE = 10.0 # seconds around the predicted P
PRE, POST = 60, 120 # window around the predicted arrival
DETECT_FLOOR = 0.02 # run once here, threshold offline
REPORT_THRESHOLD = 0.3
COLORS = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100"]
C_P, C_S = "#2a78d6", "#eb6834"
2. Data access and arrival prediction¶
These are temporary networks with no public S3 bucket, so everything comes over EarthScope FDSN.
_es = Client("EARTHSCOPE", timeout=300)
_usgs = Client("USGS", timeout=300)
_taup = TauPyModel(model="iasp91")
_station_coords = {}
def station_coords(net, sta, t0, t1):
"""Coordinates for an OBS station, cached."""
key = (net, sta)
if key not in _station_coords:
inv = _es.get_stations(network=net, station=sta, level="station",
starttime=t0, endtime=t1)
s = inv[0][0]
_station_coords[key] = (s.latitude, s.longitude)
return _station_coords[key]
def predicted_p(origin, slat, slon):
"""iasp91 P arrival at the station. None if no ray reaches it."""
deg = locations2degrees(slat, slon, origin.latitude, origin.longitude)
depth = max((origin.depth or 0) / 1000.0, 0.0)
arrivals = _taup.get_travel_times(source_depth_in_km=depth,
distance_in_degree=deg, phase_list=["p", "P"])
if not arrivals:
return None, deg
return origin.time + arrivals[0].time, deg
def fetch_obs(net, sta, band, hydro, t0, t1):
"""Four-component OBS stream: Z, 1, 2 and the hydrophone."""
try:
st = _es.get_waveforms(net, sta, "*", f"{band}?,{hydro}", t0, t1)
except Exception:
return None
st.merge(fill_value=0)
comps = {tr.stats.channel[-1] for tr in st}
if not {"Z", "1", "2"} <= comps:
return None
return st
def as_three_component(st):
"""Land models expect ZNE. OBS horizontals are 1 and 2, arbitrarily
oriented; renaming them is the usual convention and costs nothing here
because neither model uses absolute orientation."""
out = obspy.Stream()
for tr in st:
if tr.stats.channel[-1] == "H": # drop the hydrophone
continue
tr = tr.copy()
if tr.stats.channel[-1] == "1":
tr.stats.channel = tr.stats.channel[:-1] + "N"
elif tr.stats.channel[-1] == "2":
tr.stats.channel = tr.stats.channel[:-1] + "E"
out += tr
return out
3. Load the weight sets¶
models, needs_h = {}, {}
for name, spec in MODEL_SPECS.items():
try:
m = spec['build']()
except Exception as exc:
print(f"{name:<20} could not load ({type(exc).__name__}) - skipping")
continue
models[name] = m
needs_h[name] = spec['hydrophone']
print(f"{name:<20} {type(m).__name__:<15} comps={m.component_order:<6} "
f"in_samples={m.in_samples:<5} hydrophone={needs_h[name]}")
names = list(models)
pickblue_phasenet PhaseNet comps=Z12H in_samples=3001 hydrophone=True pickblue_eqt EQTransformer comps=Z12H in_samples=6000 hydrophone=True obstransformer OBSTransformer comps=ZNE in_samples=6000 hydrophone=False quakescope2026 PhaseNet comps=ZNE in_samples=3001 hydrophone=False original PhaseNet comps=ENZ in_samples=3001 hydrophone=False
4. Run every model over every deployment¶
One short window per catalogued event per station. The slow part is the network, not the inference.
records = []
examples, gallery = {}, {}
for label, exp in EXPERIMENTS.items():
print(f"{label} - {exp['note']}")
try:
cat = _usgs.get_events(starttime=exp['start'], endtime=exp['end'],
latitude=exp['lat'], longitude=exp['lon'],
maxradius=exp['max_radius_deg'],
minmagnitude=exp['min_mag'])
except Exception as exc:
print(f" catalog failed: {type(exc).__name__}")
continue
print(f" {len(cat)} catalogued events M>={exp['min_mag']}")
for sta, band, hydro in exp['stations']:
try:
slat, slon = station_coords(exp['network'], sta, exp['start'], exp['end'])
except Exception as exc:
print(f" {sta}: metadata failed ({type(exc).__name__})")
continue
n_win = 0
for ev in cat:
origin = ev.preferred_origin() or (ev.origins[0] if ev.origins else None)
if origin is None:
continue
tp, deg = predicted_p(origin, slat, slon)
if tp is None:
continue
st = fetch_obs(exp['network'], sta, band, hydro, tp - PRE, tp + POST)
if st is None:
continue
n_win += 1
st3 = as_three_component(st)
mag = ev.preferred_magnitude() or ev.magnitudes[0]
has_h = any(tr.stats.channel[-1] == 'H' for tr in st)
win_picks = {}
for name, model in models.items():
use = st if needs_h[name] else st3
try:
out = model.classify(use, P_threshold=DETECT_FLOOR,
S_threshold=DETECT_FLOOR)
except Exception:
continue
best = None
for p in out.picks:
if p.phase != 'P':
continue
dt = p.peak_time - tp
if abs(dt) <= TOLERANCE and (best is None or
float(p.peak_value) > best[1]):
best = (dt, float(p.peak_value))
if best:
win_picks[name] = best
records.append(dict(
experiment=label, station=sta, weights=name,
mag=round(mag.mag, 1), dist_deg=round(deg, 3),
has_hydrophone=has_h,
dt=round(best[0], 2) if best else np.nan,
conf=round(best[1], 3) if best else 0.0,
n_picks=len(out.picks),
))
if name == ABLATION_MODEL and best and best[1] >= REPORT_THRESHOLD:
examples.setdefault(label, []).append(
(st, tp, sta, mag.mag, best))
per_sta = sum(1 for g in gallery.get(label, []) if g['sta'] == sta)
room = int(np.ceil(GALLERY_N / max(len(exp['stations']), 1)))
if per_sta < room and len(gallery.get(label, [])) < GALLERY_N + room:
keep = st.slice(tp - 30, tp + 70).copy()
gallery.setdefault(label, []).append(dict(
stream=keep, tp=tp, sta=sta, mag=mag.mag,
deg=deg, picks=dict(win_picks)))
print(f" {sta:<6} {n_win} windows with data")
print()
det = pd.DataFrame(records)
print(f"{len(det)} (event, station, model) rows")
Cascadia (7D) - Cascadia subduction margin, offshore northern California
19 catalogued events M>=3.5
FS07B 19 windows with data
FS06B 19 windows with data
FS05B 19 windows with data AACSE (XO) - Alaska Peninsula subduction zone, shelf and trench
12 catalogued events M>=3.5
LA23 12 windows with data
LA25 12 windows with data
LD36 12 windows with data Blanco (X9) - Blanco oceanic transform - strike-slip in young oceanic crust
15 catalogued events M>=3.5
BB060 15 windows with data
BB090 15 windows with data
BB070 15 windows with data 690 (event, station, model) rows
5. Detection rate¶
A detection is a P pick within the tolerance of the predicted arrival, counted at the reporting threshold.
hit = det[det.conf >= REPORT_THRESHOLD]
rows = []
for label in EXPERIMENTS:
sub_all = det[det.experiment == label]
if not len(sub_all):
continue
for name in names:
total = len(sub_all[sub_all.weights == name])
found = len(hit[(hit.experiment == label) & (hit.weights == name)])
if not total:
continue
res = hit[(hit.experiment == label) & (hit.weights == name)]['dt']
rows.append(dict(experiment=label, weights=name,
windows=total, detected=found,
rate=round(found / total, 3),
median_dt=round(float(res.median()), 2) if len(res) else np.nan))
summary = pd.DataFrame(rows)
print(summary.to_string(index=False))
print(f'\ndetection = a P pick within {TOLERANCE:g} s of the iasp91 prediction, '
f'at confidence >= {REPORT_THRESHOLD}')
print('median_dt = median offset from the prediction. What matters is that it '
'is small and similar across models, which says they lock onto the same '
'arrival and the offset belongs to the prediction rather than the picker.')
experiment weights windows detected rate median_dt Cascadia (7D) pickblue_phasenet 57 42 0.737 -1.22 Cascadia (7D) pickblue_eqt 57 45 0.789 -1.26 Cascadia (7D) obstransformer 57 35 0.614 -1.10 Cascadia (7D) quakescope2026 57 41 0.719 -0.93 Cascadia (7D) original 57 37 0.649 -0.93 AACSE (XO) pickblue_phasenet 36 28 0.778 -1.49 AACSE (XO) pickblue_eqt 36 26 0.722 -1.59 AACSE (XO) obstransformer 36 32 0.889 -1.69 AACSE (XO) quakescope2026 36 28 0.778 -1.43 AACSE (XO) original 36 24 0.667 -1.32 Blanco (X9) pickblue_phasenet 45 24 0.533 -3.42 Blanco (X9) pickblue_eqt 45 24 0.533 -3.39 Blanco (X9) obstransformer 45 23 0.511 -3.11 Blanco (X9) quakescope2026 45 24 0.533 -2.83 Blanco (X9) original 45 22 0.489 -3.31 detection = a P pick within 10 s of the iasp91 prediction, at confidence >= 0.3 median_dt = median offset from the prediction. What matters is that it is small and similar across models, which says they lock onto the same arrival and the offset belongs to the prediction rather than the picker.
if len(summary):
labels = [l for l in EXPERIMENTS if l in set(summary.experiment)]
fig, ax = plt.subplots(figsize=(10.5, 4.3))
width = 0.8 / max(len(names), 1)
for i, name in enumerate(names):
xs, vals = [], []
for j, lab in enumerate(labels):
r = summary[(summary.experiment == lab) & (summary.weights == name)]
if len(r):
xs.append(j + (i - (len(names) - 1) / 2) * width)
vals.append(float(r['rate'].iloc[0]))
ax.bar(xs, vals, width=width * 0.9, color=COLORS[i % len(COLORS)], label=name)
for x, v in zip(xs, vals):
ax.text(x, v + 0.015, f'{v:.2f}', ha='center', fontsize=8)
ax.set_xticks(range(len(labels)))
ax.set_xticklabels(labels, fontsize=9)
ax.set_ylabel('fraction of catalogued events detected')
ax.set_ylim(0, 1.05)
ax.grid(alpha=0.25, lw=0.5, axis='y')
ax.set_axisbelow(True)
ax.legend(frameon=False, fontsize=9, ncol=len(names))
ax.set_title('Detection of catalogued events on ocean-bottom stations',
fontsize=11, loc='left')
fig.tight_layout()
plt.show()
Where the picks land relative to the prediction¶
What matters is that the cluster is tight and in the same place for every model — that indicates they are all locking onto the same physical arrival, and that the offset is a property of the prediction rather than of any weight set.
The offset direction is not a clean diagnostic here. A water column and sediments delay the true arrival relative to iasp91, but offshore catalog locations and depths are themselves poorly constrained — many are fixed rather than solved — and that error enters the prediction directly and in either direction. A consistent offset of a second or two says the prediction is biased, not that the picks are wrong.
fig, ax = plt.subplots(figsize=(10.5, 3.8))
bins = np.linspace(-TOLERANCE, TOLERANCE, 33)
for name, color in zip(names, COLORS):
d = hit[hit.weights == name]['dt'].dropna()
if not len(d):
continue
ax.hist(d, bins=bins, histtype='step', lw=1.8, color=color,
label=f'{name} (n={len(d)})')
ax.axvline(0, color='#8a8a8a', lw=1)
ax.set_xlabel('pick minus iasp91 prediction (s)')
ax.set_ylabel('detections')
ax.set_title('Offset from the predicted arrival', fontsize=11, loc='left')
ax.grid(alpha=0.25, lw=0.5)
ax.legend(frameon=False, fontsize=9)
fig.tight_layout()
plt.show()
6. Does the hydrophone earn its place?¶
obs is the only model here that can use the fourth channel. The
hydrophone is also the channel whose sampling rate varies most between
these deployments — 100 Hz at AACSE, 40 Hz at Cascadia, 10 Hz at
Blanco — while the model resamples everything to 100 Hz. An upsampled
10 Hz trace carries very little of the band a P onset lives in.
Re-running obs on the same windows with the hydrophone withheld
separates the value of the channel from the value of the training.
if ABLATION_MODEL in models:
rows = []
for label, entries in examples.items():
for st, tp, sta, mag, best in entries:
with_h = st
without_h = obspy.Stream([tr for tr in st if tr.stats.channel[-1] != 'H'])
entry = {'experiment': label, 'station': sta, 'M': mag}
for tag, stream in (('with_H', with_h), ('without_H', without_h)):
out = models[ABLATION_MODEL].classify(stream, P_threshold=DETECT_FLOOR,
S_threshold=DETECT_FLOOR)
cand = [float(p.peak_value) for p in out.picks
if p.phase == 'P' and abs(p.peak_time - tp) <= TOLERANCE]
entry[tag] = round(max(cand), 3) if cand else 0.0
entry['hydrophone_rate'] = next(
(tr.stats.sampling_rate for tr in st if tr.stats.channel[-1] == 'H'), np.nan)
rows.append(entry)
abl = pd.DataFrame(rows)
if len(abl):
abl['delta'] = (abl.with_H - abl.without_H).round(3)
print(f'{len(abl)} detected windows re-run with the hydrophone withheld\n')
print(abl.groupby(['experiment', 'hydrophone_rate'])[['with_H', 'without_H', 'delta']]
.agg(['mean', 'count']).round(3).to_string())
print(f'\nmean change in P confidence when the hydrophone is included: '
f'{abl.delta.mean():+.4f}')
worse = int((abl.delta < -0.01).sum()); better = int((abl.delta > 0.01).sum())
print(f'windows where it helped by >0.01: {better}; hurt by >0.01: {worse}; '
f'unchanged: {len(abl) - better - worse}')
else:
print('no detections captured for the ablation')
94 detected windows re-run with the hydrophone withheld
with_H without_H delta
mean count mean count mean count
experiment hydrophone_rate
AACSE (XO) 100.0 0.807 28 0.866 28 -0.059 28
Blanco (X9) 10.0 0.855 24 0.876 24 -0.022 24
Cascadia (7D) 40.0 0.772 30 0.693 30 0.079 30
125.0 0.699 12 0.714 12 -0.015 12
mean change in P confidence when the hydrophone is included: +0.0002
windows where it helped by >0.01: 32; hurt by >0.01: 35; unchanged: 27
7. Scanning the records¶
Ten windows per deployment, spread across its stations rather than drawn from one, vertical component, with every model's pick overlaid so disagreements are visible at a glance. The dotted line is the iasp91 prediction, which per section 5 carries the catalog's location error and should be read as an approximate marker rather than truth.
What to look for: picks clustering on a visible onset is the model working; picks spread across seconds of an emergent arrival is the hard case that ocean-bottom noise creates; and a pick with nothing visible under it is worth following up.
def gallery_figure(label, entries, window=(-25, 60)):
n = len(entries)
fig, axes = plt.subplots(n, 1, figsize=(11.5, 1.35 * n + 1.0), sharex=True)
axes = np.atleast_1d(axes)
for ax, e in zip(axes, entries):
tr = e['stream'].select(component='Z')
if not tr:
continue
tr = tr[0]
t = tr.times(reftime=e['tp'])
x = tr.data.astype(float)
peak = np.abs(x).max()
if peak > 0:
x = x / peak
ax.plot(t, x, color='#3d3d3d', lw=0.5)
ax.axvline(0, color='#8a8a8a', lw=1.1, ls=':')
for i, name in enumerate(names):
if name not in e['picks']:
continue
dt, conf = e['picks'][name]
ax.axvline(dt, color=COLORS[i % len(COLORS)], lw=1.5,
alpha=0.85 if conf >= REPORT_THRESHOLD else 0.35)
found = len(e['picks'])
ax.set_ylabel(f"{e['sta']}\nM{e['mag']:.1f} {e['deg']*111:.0f}km",
fontsize=7.5, rotation=0, ha='right', va='center', labelpad=32)
ax.set_yticks([])
ax.grid(alpha=0.18, lw=0.4, axis='x')
ax.text(0.995, 0.82, f'{found}/{len(names)} models',
transform=ax.transAxes, ha='right', fontsize=7, color='#7a7973')
handles = [plt.Line2D([], [], color=COLORS[i % len(COLORS)], lw=1.6, label=n)
for i, n in enumerate(names)]
handles.append(plt.Line2D([], [], color='#8a8a8a', lw=1.1, ls=':',
label='iasp91 prediction'))
axes[0].legend(handles=handles, fontsize=7.5, frameon=False, ncol=3,
loc='lower left', bbox_to_anchor=(0, 1.05))
axes[-1].set_xlim(*window)
axes[-1].set_xlabel('seconds from the predicted P')
fig.suptitle(label, fontsize=11, x=0.01, ha='left', y=0.998)
fig.tight_layout(rect=[0, 0, 1, 0.985])
return fig
for label in EXPERIMENTS:
entries = gallery.get(label) or []
if not entries:
continue
gallery_figure(label, entries[:GALLERY_N])
plt.show()
One record in full¶
All four components for a single detected event per deployment, including the hydrophone, at its native sampling rate.
def plot_obs(label, entry, window=(-30, 90)):
st, tp, sta, mag, best = entry
order = ['Z', '1', '2', 'H']
traces = [next((tr for tr in st if tr.stats.channel[-1] == c), None) for c in order]
traces = [t for t in traces if t is not None]
fig, axes = plt.subplots(len(traces), 1, figsize=(11, 1.5 * len(traces) + 1.1),
sharex=True)
axes = np.atleast_1d(axes)
for ax, tr in zip(axes, traces):
ax.plot(tr.times(reftime=tp), tr.data, color='#3d3d3d', lw=0.5)
ax.axvline(0, color='#8a8a8a', lw=1.2, ls=':')
ax.axvline(best[0], color=C_P, lw=1.6)
ax.set_ylabel(f'{tr.stats.channel}\n{tr.stats.sampling_rate:g} Hz', fontsize=8)
ax.grid(alpha=0.22, lw=0.5)
ax.tick_params(labelsize=8)
axes[-1].set_xlim(*window)
axes[-1].set_xlabel('seconds from the iasp91 predicted P '
'(dotted); solid line is the pick')
axes[0].set_title(f'{label} - {sta} - M{mag} - '
f'pick {best[0]:+.1f} s at confidence {best[1]:.2f}',
fontsize=10, loc='left')
fig.tight_layout()
return fig
for label in EXPERIMENTS:
entries = examples.get(label) or []
if entries:
plot_obs(label, max(entries, key=lambda e: e[4][1]))
plt.show()
8. The campaign's own picks against published arrivals¶
Everything above scores models on windows fetched for the purpose. The
obs campaign has now run: 6,231 shards over 23 networks, written as
Parquet on S3, using one of the weight sets from section 3. That output is
what QuakeScope will hand to people, so it is what needs validating, and
against something better than a travel-time prediction.
Two published references are obtainable today. Both are compared the way
western_pick_validation.ipynb compares
the western campaign against a re-pick: match on station, phase and
time within a stated tolerance; report what the campaign holds separately
from what it found; and never turn the reference's silence into a false
positive. A reference catalogue is not exhaustive, and an OBS deployment's
archive is not complete, so an unmatched campaign pick is unscored, not
wrong.
| Reference | What it is | Where it comes from |
|---|---|---|
| AACSE 2018: Barcheck (2023), Ocean-bottom P and S arrival waveform dataset from the Alaska Amphibious Community Seismic Experiment, 2018–19, Cornell eCommons, doi:10.7298/01da-ka24, CC BY 4.0 | Analyst-checked P and S arrival times from the Alaska Earthquake Center catalogue (Ruppert, Barcheck & Abers 2023, SRL 94, doi:10.1785/0220220226), one row per (event, OBS station) for every event within 350 km. Each pick carries a manual/automatic status. A companion land dataset (doi:10.7298/q2fq-9688) gives the same for the 30 land stations of the same network. |
Downloaded here from the repository's public API. Only the metadata tables are needed - the picks are in them - so the 2 GB of waveforms stay where they are. |
Axial Seamount 2015–2025: the University of Washington near-real-time catalogue (Wilcock et al. 2016, Science 354, doi:10.1126/science.aah5563), phase file ph2dtInputCatalog.dat at axial.ocean.washington.edu |
Automatic P and S picks at the seven OOI cabled stations, as travel times from each located event, with a pick weight. Roughly 320,000 events. | Downloaded here. It is regenerated hourly, so counts drift slightly between runs. |
Axial Seamount 2014–2021: Wang et al. (2024) ML-DD catalogue Axial.MLDD.v202112.2 at axialdd.ldeo.columbia.edu |
144,329 relocated hypocentres. Events only, no picks; used for an event-level detection check. | Downloaded here. |
The AACSE comparison is the one that matters most: the reference is analyst truth, both sides are picks, and it covers 65 instruments across three instrument types. Axial is a picker-against-picker comparison - the UW picks are automatic - on a very different kind of seismicity, and it is kept because that difference is informative.
What the campaign ran¶
Read from the campaign's own run records, never typed in.
import datetime, hashlib, io, json, ssl, tarfile, time, urllib.error, urllib.request
from concurrent.futures import ThreadPoolExecutor
from pathlib import Path
import boto3
BUCKET, REGION, CAMPAIGN = "quakescope-picks-2026", "us-east-2", "obs"
CACHE = Path("obs_benchmark_cache") # gitignored. Delete it to re-download everything.
CACHE.mkdir(exist_ok=True)
s3 = boto3.client("s3", region_name=REGION)
_pg = s3.get_paginator("list_objects_v2")
TOL_ANALYST = 1.0 # s, against AACSE analyst arrivals (sections 9 and 10)
TOL_AXIAL = 0.5 # s, against Axial picks - sources sit within a few km of the array
PALETTE = COLORS + ["#8e5bd6"] # a fifth slot, so five models never share a colour
runs, token = [], None
while len(runs) < 300: # a sample of the ~6,500 run records
kw = dict(Bucket=BUCKET, Prefix=f"{CAMPAIGN}/runs/", MaxKeys=150)
if token:
kw["ContinuationToken"] = token
page = s3.list_objects_v2(**kw)
for o in page.get("Contents", []):
runs.append(json.loads(s3.get_object(Bucket=BUCKET, Key=o["Key"])["Body"].read()))
token = page.get("NextContinuationToken")
if not token:
break
cfg = pd.DataFrame(runs)
fields = ["model", "weight", "p_threshold", "s_threshold", "components_loaded",
"seisbench_version", "weight_version"]
print(f"{len(cfg)} run records sampled")
for f in fields:
print(f" {f:20s} {cfg[f].value_counts().to_dict()}")
assert all(cfg[f].nunique() == 1 for f in fields), "the obs campaign was not run under one configuration"
PROD = {f: cfg[f].iloc[0] for f in fields}
# The campaign's weight file and the one PickBlue(base="phasenet") loaded in section 3.
repo_w = Path("../sb_catalog/models/v3/phasenet/obs.pt.v1")
cache_w = Path.home() / ".seisbench/models/v3/phasenet/obs.pt.v1"
md5 = {p.name + " (" + tag + ")": hashlib.md5(p.read_bytes()).hexdigest()[:12]
for tag, p in (("shipped in the image", repo_w), ("seisbench cache", cache_w)) if p.exists()}
print(f"\nThe campaign is PhaseNet/{PROD['weight']} at P={PROD['p_threshold']}, S={PROD['s_threshold']}, "
f"components {PROD['components_loaded']}, seisbench {PROD['seisbench_version']}.")
print("Weight file md5:", md5)
print("`obs` is the PickBlue PhaseNet weight (component order Z12H). The campaign loads ZNE12 and no\n"
"hydrophone, so SeisBench zero-fills the H slot - section 6 measured that channel's effect at +0.0002.")
300 run records sampled
model {'PhaseNet': 300}
weight {'obs': 300}
p_threshold {'0.2': 300}
s_threshold {'0.2': 300}
components_loaded {'ZNE12': 300}
seisbench_version {'0.12.5': 300}
weight_version {'1': 300}
The campaign is PhaseNet/obs at P=0.2, S=0.2, components ZNE12, seisbench 0.12.5.
Weight file md5: {'obs.pt.v1 (shipped in the image)': '21e57268fcd9', 'obs.pt.v1 (seisbench cache)': '21e57268fcd9'}
`obs` is the PickBlue PhaseNet weight (component order Z12H). The campaign loads ZNE12 and no
hydrophone, so SeisBench zero-fills the H slot - section 6 measured that channel's effect at +0.0002.
Reference data¶
Fetched into obs_benchmark_cache/ beside this notebook and reused on later
runs. The Barcheck tables are the SeisBench metadata format: one row per
waveform, with the P and S arrival as a sample index into a trace whose
start time is given to the microsecond.
ECOMMONS = "https://ecommons.cornell.edu/server/api/core/bitstreams/{}/content"
REFERENCES = {
"barcheck_obs_readme": (ECOMMONS.format("77528b48-2397-41c4-8bfc-558a29faa28d"), "AACSE_2018_ReadMe.txt"),
"barcheck_obs_2018": (ECOMMONS.format("08969cf9-30ff-4067-8ceb-8f732696d2ff"), "archive_metadata2018.tar.gz"),
"barcheck_land": (ECOMMONS.format("20b7a01a-c576-4b44-b8c1-a592ee4d2980"), "archive_AACSE_metadata_land.tar.gz"),
"uw_axial_phases": ("http://axial.ocean.washington.edu/ph2dtInputCatalog.dat", "ph2dtInputCatalog.dat"),
"ldeo_axial_readme": ("https://axialdd.ldeo.columbia.edu/catalog/README", "Axial.MLDD.README"),
"ldeo_axial_catalogue": ("https://axialdd.ldeo.columbia.edu/catalog/Axial.MLDD.v202112.2", "Axial.MLDD.v202112.2"),
}
def fetch(url, name):
dest = CACHE / name
if dest.exists() and dest.stat().st_size > 0:
return dest
try:
with urllib.request.urlopen(url, timeout=900) as r:
dest.write_bytes(r.read())
except urllib.error.URLError as exc:
if "ldeo.columbia.edu" not in url or "CERTIFICATE" not in str(exc).upper():
raise
# The LDEO host serves a chain Python cannot verify. Public catalogue, so accept it.
with urllib.request.urlopen(url, timeout=900, context=ssl._create_unverified_context()) as r:
dest.write_bytes(r.read())
return dest
files = {}
for key, (url, name) in REFERENCES.items():
t0 = time.time()
files[key] = fetch(url, name)
print(f" {key:22s} {files[key].stat().st_size / 1e6:8.1f} MB {time.time() - t0:5.1f} s {name}")
for key, sub in (("barcheck_obs_2018", "barcheck_obs_2018"), ("barcheck_land", "barcheck_land")):
out = CACHE / sub
if not out.exists():
with tarfile.open(files[key]) as tf:
tf.extractall(out)
print("\nBarcheck readme, first lines:")
print("\n".join(files["barcheck_obs_readme"].read_text().splitlines()[:3]))
print("\nLDEO readme, overview:")
print("\n".join(l for l in files["ldeo_axial_readme"].read_text().splitlines()[18:23]))
barcheck_obs_readme 0.0 MB 0.0 s AACSE_2018_ReadMe.txt barcheck_obs_2018 0.7 MB 0.0 s archive_metadata2018.tar.gz barcheck_land 4.0 MB 0.0 s archive_AACSE_metadata_land.tar.gz uw_axial_phases 102.2 MB 0.0 s ph2dtInputCatalog.dat ldeo_axial_readme 0.0 MB 0.0 s Axial.MLDD.README ldeo_axial_catalogue 14.7 MB 0.0 s Axial.MLDD.v202112.2 Barcheck readme, first lines: Title: Dataset: Ocean-bottom P and S arrival waveform dataset from the Alaska Amphibious Community Seismic Experiment, 2018 LDEO readme, overview: OVERVIEW: This file includes high-precision hypocenter locations for 144,329 earthquakes recorded between November 4, 2014 and December 31, 2021 by the Ocean Observatories Initiatives (OOI) cabled OBS array. They
The campaign's picks¶
Every pick object under obs/picks/network=XO/year=2018/ and
obs/picks/network=OO/, read once from S3 and cached as Parquet. XO 2018
is the only AACSE year the campaign holds - 2019 is recorded complete with
nothing read, see
27_obs_literature_benchmark.md -
and OO runs from the start of the cabled array in 2014 to the campaign's
end of 2025.
def campaign_picks(prefix, name):
dest = CACHE / name
if dest.exists():
return pd.read_parquet(dest)
keys = [o["Key"] for page in _pg.paginate(Bucket=BUCKET, Prefix=prefix)
for o in page.get("Contents", [])]
def read(k):
return pd.read_parquet(io.BytesIO(s3.get_object(Bucket=BUCKET, Key=k)["Body"].read()))
with ThreadPoolExecutor(16) as ex:
frames = list(ex.map(read, keys))
df = pd.concat(frames, ignore_index=True) if frames else pd.DataFrame()
df.to_parquet(dest)
return df
t0 = time.time()
camp_xo = campaign_picks(f"{CAMPAIGN}/picks/network=XO/year=2018/", "campaign_XO_2018.parquet")
camp_oo = campaign_picks(f"{CAMPAIGN}/picks/network=OO/", "campaign_OO.parquet")
print(f"read in {time.time() - t0:.0f} s\n")
for label, df in (("XO 2018", camp_xo), ("OO 2014-2025", camp_oo)):
print(f"{label:14s} {len(df):>10,} picks {df.tid.nunique():3d} stations "
f"bands {df.cha.value_counts().to_dict()} P/S {df.pha.value_counts().to_dict()} "
f"{df.peak.min():%Y-%m-%d} .. {df.peak.max():%Y-%m-%d}")
print(f"\nconfidence floor in the data: P {camp_xo[camp_xo.pha == 'P'].conf.min():.3f}, "
f"S {camp_xo[camp_xo.pha == 'S'].conf.min():.3f} (the campaign wrote everything at or above its thresholds)")
read in 1 s
XO 2018 4,750,741 picks 97 stations bands {'HH': 4750741} P/S {'P': 2887289, 'S': 1863452} 2018-05-11 .. 2018-12-31
OO 2014-2025 4,045,124 picks 13 stations bands {'EH': 2671822, 'HH': 1373302} P/S {'P': 2993075, 'S': 1052049} 2014-11-04 .. 2026-01-01
confidence floor in the data: P 0.200, S 0.200 (the campaign wrote everything at or above its thresholds)
9. AACSE 2018: sixty-five ocean-bottom stations against analyst arrivals¶
The reference is every P and S arrival the Alaska Earthquake Center's
analysts kept for an OBS station, for events within 350 km, May–December
2018. Picks flagged automatic in that catalogue were machine-made and not
individually reviewed, so the headline is manual picks only; the
automatic ones are reported alongside because they are a second machine's
opinion of the same onsets.
Two numbers are kept apart on purpose:
- coverage - does the campaign hold any pick for that station on that day? If not, the arrival is not scored. A station-day with no picks at all is one the campaign did not read, not one where the picker found nothing; section 9.5 traces those days back to their shards.
- recall - on covered station-days, the fraction of analyst arrivals with a campaign pick of the same phase within 1 s.
One second is deliberately strict. The campaign wrote its picks at millisecond resolution and the analyst times are given to the microsecond, so the tolerance is set by physics rather than bookkeeping: for these paths a second is well under the S–P time and comparable to the width of an emergent OBS onset.
def analyst_arrivals(ref):
"""One row per analyst arrival: station, phase, absolute time, status, and the
event it belongs to. The Barcheck tables give the arrival as a sample index."""
ref = ref.copy()
ref["t0"] = pd.to_datetime(ref.trace_start_time, utc=True).dt.tz_localize(None)
ref["tid"] = "XO." + ref.station_code.str.strip() + "."
ref["dist_km"] = locations2degrees(ref.station_latitude, ref.station_longitude,
ref.source_latitude_deg, ref.source_longitude_deg) * 111.19
keep = ["tid", "pha", "t", "status", "source_id", "source_magnitude", "source_depth_km",
"trace_region", "trace_instrument_type", "dist_km"]
out = []
for pha in ("P", "S"):
s = ref[ref[f"trace_{pha}_arrival_sample"].notna()].copy()
s["t"] = s.t0 + pd.to_timedelta(s[f"trace_{pha}_arrival_sample"] / s.trace_sampling_rate_hz, unit="s")
s["pha"], s["status"] = pha, s[f"trace_{pha}_status"]
out.append(s[keep])
return pd.concat(out, ignore_index=True)
def _sec(ts):
return pd.to_datetime(ts).values.astype("datetime64[ns]").astype("int64") / 1e9
def nearest(ref_t, cand_t, cand_conf):
"""For each reference time, the nearest candidate (candidates sorted): offset and confidence."""
if len(cand_t) == 0:
return np.full(len(ref_t), np.nan), np.full(len(ref_t), np.nan)
i = np.searchsorted(cand_t, ref_t)
lo, hi = np.clip(i - 1, 0, len(cand_t) - 1), np.clip(i, 0, len(cand_t) - 1)
j = np.where(np.abs(cand_t[hi] - ref_t) < np.abs(cand_t[lo] - ref_t), hi, lo)
return cand_t[j] - ref_t, cand_conf[j]
def match_to_campaign(arr, picks, tol):
"""dt = campaign pick minus reference pick, in seconds, for the nearest campaign pick of
the same phase at the same station; `covered` = the campaign holds picks on that station-day."""
arr = arr.copy()
arr["day"] = arr.t.dt.floor("D")
held = picks.assign(day=picks.peak.dt.floor("D")).groupby(["tid", "day"]).size()
arr["covered"] = pd.MultiIndex.from_arrays([arr.tid, arr.day]).isin(held.index)
arr["dt"], arr["conf"] = np.nan, np.nan
groups = {k: (_sec(g.peak), g.conf.values) for k, g in picks.sort_values("peak").groupby(["tid", "pha"])}
for (tid, pha), g in arr.groupby(["tid", "pha"]):
if (tid, pha) in groups:
dt, conf = nearest(_sec(g.t), *groups[(tid, pha)])
arr.loc[g.index, "dt"], arr.loc[g.index, "conf"] = dt, conf
arr["hit"] = arr.dt.abs() <= tol
return arr
ref_obs = pd.concat([pd.read_csv(f) for f in sorted((CACHE / "barcheck_obs_2018").glob("*.csv"))],
ignore_index=True)
ref_land = pd.concat([pd.read_csv(f) for f in sorted((CACHE / "barcheck_land").glob("metadata2018*.csv"))],
ignore_index=True)
ref_land = ref_land[ref_land.station_network_code == "XO"]
print(f"Barcheck OBS 2018: {len(ref_obs):,} station-events, {ref_obs.source_id.nunique():,} events, "
f"{ref_obs.station_code.nunique()} stations, M {ref_obs.source_magnitude.min():.1f}-{ref_obs.source_magnitude.max():.1f}")
print(f"Barcheck land 2018: {len(ref_land):,} station-events on the {ref_land.station_code.nunique()} XO land stations")
t0 = time.time()
arr_obs = match_to_campaign(analyst_arrivals(ref_obs), camp_xo, TOL_ANALYST)
arr_land = match_to_campaign(analyst_arrivals(ref_land), camp_xo, TOL_ANALYST)
print(f"matched in {time.time() - t0:.1f} s\n")
for label, arr in (("OBS", arr_obs), ("land", arr_land)):
print(f"{label}: {len(arr):,} analyst arrivals; {arr.covered.mean():.1%} fall on station-days the campaign holds picks for "
f"({(~arr.covered).sum():,} do not, on {arr[~arr.covered].groupby(['tid', 'day']).ngroups} station-days)")
cov = arr_obs[arr_obs.covered]
print(f"\nRecall within {TOL_ANALYST:g} s on covered station-days, OBS:")
print(cov.groupby(["pha", "status"]).hit.agg(recall="mean", n="size").round(3).to_string())
man = cov[cov.status == "manual"]
man_land = arr_land[arr_land.covered & (arr_land.status == "manual")]
print("\nManual picks, by tolerance (OBS / land):")
tols = (0.25, 0.5, 1.0, 2.0)
tab = pd.DataFrame({f"{t:g} s": [(g.dt.abs() <= t).mean() for _, g in man.groupby("pha")] for t in tols},
index=[f"OBS {p}" for p in ("P", "S")])
tab = pd.concat([tab, pd.DataFrame({f"{t:g} s": [(g.dt.abs() <= t).mean() for _, g in man_land.groupby("pha")] for t in tols},
index=[f"land {p}" for p in ("P", "S")])])
print(tab.round(3).to_string())
hits = man[man.hit]
print("\nResidual, campaign minus analyst, on hits (s):")
print(hits.groupby("pha").dt.agg(median="median", q25=lambda x: x.quantile(.25), q75=lambda x: x.quantile(.75), n="size").round(3).to_string())
Barcheck OBS 2018: 20,828 station-events, 2,598 events, 65 stations, M 0.6-6.1 Barcheck land 2018: 20,710 station-events on the 30 XO land stations
matched in 6.5 s
OBS: 30,618 analyst arrivals; 95.3% fall on station-days the campaign holds picks for (1,436 do not, on 350 station-days)
land: 31,448 analyst arrivals; 96.2% fall on station-days the campaign holds picks for (1,201 do not, on 169 station-days)
Recall within 1 s on covered station-days, OBS:
recall n
pha status
P automatic 0.980 5413
manual 0.856 9150
S automatic 0.910 4047
manual 0.863 10572
Manual picks, by tolerance (OBS / land):
0.25 s 0.5 s 1 s 2 s
OBS P 0.730 0.805 0.856 0.882
OBS S 0.589 0.765 0.863 0.904
land P 0.785 0.871 0.906 0.916
land S 0.657 0.797 0.852 0.867
Residual, campaign minus analyst, on hits (s):
median q25 q75 n
pha
P 0.018 -0.050 0.091 7834
S 0.041 -0.091 0.209 9126
9.1 Where the campaign's picks land relative to the analyst's¶
The land stations run through the same obs weight, so they are the control
for "is this an ocean-bottom problem or a picker problem". If the two
distributions sit in the same place with the same width, the weight set is
timing onsets the same way on both.
fig, axes = plt.subplots(1, 2, figsize=(10.5, 3.6), sharey=False)
bins = np.arange(-TOL_ANALYST, TOL_ANALYST + 0.05, 0.05)
for ax, pha in zip(axes, ("P", "S")):
for arr, color, label in ((man, PALETTE[0], "ocean-bottom"), (man_land, PALETTE[1], "land, same weight")):
d = arr[(arr.pha == pha) & arr.hit].dt
ax.hist(d, bins=bins, histtype="step", lw=1.8, color=color, density=True,
label=f"{label} (median {d.median():+.3f} s, n={len(d):,})")
ax.axvline(0, color="#8a8a8a", lw=1)
ax.set_title(f"{pha}: campaign pick minus analyst pick", fontsize=10.5, loc="left")
ax.set_xlabel("seconds"); ax.set_yticks([])
ax.grid(alpha=0.25, lw=0.5); ax.set_axisbelow(True)
ax.legend(frameon=False, fontsize=8, loc="upper left")
fig.tight_layout(); plt.show()
9.2 Recall by instrument, setting, magnitude and distance¶
Three instrument types went into the water: shielded shallow-shelf instruments, and two makes of unshielded deep-water OBS (WHOI and LDEO). The catalogue's own regions are the shelf, the shelf-slope, and the outer rise seaward of the trench.
INSTR = {"obs_shallowshelf_shielded": "shelf, shielded", "obs_deep_whoi_unshielded": "deep WHOI",
"obs_deep_ldeo_unshielded": "deep LDEO"}
REGION = {"shelf_obs": "shelf", "shelf_slope": "shelf slope", "outer_rise": "outer rise"}
m2 = man.assign(instrument=man.trace_instrument_type.map(INSTR), region=man.trace_region.map(REGION),
magnitude=pd.cut(man.source_magnitude, [0, 2, 2.5, 3, 3.5, 4, 8],
labels=["<2", "2-2.5", "2.5-3", "3-3.5", "3.5-4", ">4"]),
distance=pd.cut(man.dist_km, [0, 50, 100, 150, 200, 350],
labels=["<50 km", "50-100", "100-150", "150-200", "200-350"]))
def breakdown(ax, col, title):
g = m2.groupby([col, "pha"], observed=True).hit.agg(["mean", "size"]).unstack("pha")
x = np.arange(len(g)); w = 0.38
for k, (pha, color) in enumerate((("P", C_P), ("S", C_S))):
vals = g[("mean", pha)].values
ax.bar(x + (k - 0.5) * w, vals, width=w * 0.92, color=color, label=pha)
for xi, v, n in zip(x + (k - 0.5) * w, vals, g[("size", pha)].values):
ax.text(xi, v + 0.012, f"{v:.2f}", ha="center", fontsize=7.2, color="#3d3d3d")
ax.text(xi, 0.02, f"n={n:,}", ha="center", fontsize=6, color="white", rotation=90, va="bottom")
ax.set_xticks(x); ax.set_xticklabels([str(i) for i in g.index], fontsize=8.5)
ax.set_ylim(0, 1.06); ax.set_title(title, fontsize=10, loc="left")
ax.grid(alpha=0.25, lw=0.5, axis="y"); ax.set_axisbelow(True)
fig, axes = plt.subplots(2, 2, figsize=(11, 7))
breakdown(axes[0, 0], "instrument", "by instrument type")
breakdown(axes[0, 1], "region", "by setting")
breakdown(axes[1, 0], "magnitude", "by magnitude")
breakdown(axes[1, 1], "distance", "by epicentral distance")
axes[0, 0].set_ylabel(f"recall within {TOL_ANALYST:g} s"); axes[1, 0].set_ylabel(f"recall within {TOL_ANALYST:g} s")
fig.legend(handles=[plt.Rectangle((0, 0), 1, 1, color=C_P, label="P"), plt.Rectangle((0, 0), 1, 1, color=C_S, label="S")],
frameon=False, fontsize=9, ncol=2, loc="upper right", bbox_to_anchor=(0.99, 0.995))
fig.suptitle("AACSE 2018, manual picks on covered station-days", fontsize=11, x=0.01, ha="left")
fig.tight_layout(rect=[0, 0, 1, 0.97]); plt.show()
9.3 What the confidence floor buys¶
The campaign kept everything at or above 0.2 for both phases. Raising the threshold offline is free, and this curve says what it costs: recall of analyst picks as a function of the confidence the campaign attached to the matching pick. The two vertical lines are the campaign's floor and the 0.3 that the rest of this notebook reports at.
thr = np.arange(0.2, 0.96, 0.05)
fig, ax = plt.subplots(figsize=(8.5, 3.8))
for arr, ls, label in ((man, "-", "ocean-bottom"), (man_land, "--", "land")):
for pha, color in (("P", C_P), ("S", C_S)):
g = arr[arr.pha == pha]
rec = [((g.dt.abs() <= TOL_ANALYST) & (g.conf >= t)).mean() for t in thr]
ax.plot(thr, rec, ls=ls, lw=2, color=color, label=f"{pha}, {label} ({rec[6]:.2f} at 0.5)")
for t, txt in ((float(PROD["p_threshold"]), "campaign floor"), (REPORT_THRESHOLD, "reported at")):
ax.axvline(t, color="#8a8a8a", lw=1, ls=":"); ax.text(t + 0.005, 0.03, txt, fontsize=7.5, color="#7a7973")
ax.set_xlabel("confidence threshold applied to the campaign's picks")
ax.set_ylabel(f"recall of manual picks within {TOL_ANALYST:g} s"); ax.set_ylim(0, 1)
ax.grid(alpha=0.25, lw=0.5); ax.set_axisbelow(True); ax.legend(frameon=False, fontsize=8.5, ncol=2)
ax.set_title("Recall against the analyst as the threshold rises", fontsize=10.5, loc="left")
fig.tight_layout(); plt.show()
9.4 Station by station¶
The spread between stations is the more useful number than the mean, because it separates a picker that is uniformly mediocre from one that is good almost everywhere and broken somewhere specific. For any station well below the pack, the month-by-month median offset between the campaign's nearest pick and the analyst's tells whether the picks are missing or mis-timed: a picker failure gives offsets scattered within a few seconds, whereas an archive whose clock correction differs from the one the analysts worked with gives offsets that move by minutes or hours from month to month.
per = (man[man.pha == "P"].groupby("tid")
.agg(recall=("hit", "mean"), n=("hit", "size"), instrument=("trace_instrument_type", "first"))
.sort_values("recall"))
per["instrument"] = per.instrument.map(INSTR)
print(f"manual P recall per station: median {per.recall.median():.3f}, "
f"quartiles {per.recall.quantile(.25):.3f}-{per.recall.quantile(.75):.3f}, "
f"{(per.recall >= 0.8).sum()} of {len(per)} stations at or above 0.8")
fig, ax = plt.subplots(figsize=(8, 0.16 * len(per) + 1.2))
icol = {v: PALETTE[i] for i, v in enumerate(INSTR.values())}
ax.barh(np.arange(len(per)), per.recall, color=[icol[i] for i in per.instrument], height=0.7)
ax.set_yticks(np.arange(len(per))); ax.set_yticklabels([t.split(".")[1] for t in per.index], fontsize=6.5)
ax.set_xlim(0, 1); ax.set_xlabel(f"manual P recall within {TOL_ANALYST:g} s"); ax.grid(alpha=0.25, lw=0.5, axis="x")
ax.set_axisbelow(True)
ax.legend(handles=[plt.Rectangle((0, 0), 1, 1, color=c, label=k) for k, c in icol.items()],
frameon=False, fontsize=8, loc="lower right")
ax.set_title("AACSE 2018, one bar per OBS station", fontsize=10.5, loc="left")
fig.tight_layout(); plt.show()
weak = per[per.recall < 0.7]
print(f"\n{len(weak)} stations below 0.7. Month-by-month median offset of the nearest campaign P pick (s):")
for tid in weak.index:
a = man[(man.tid == tid) & (man.pha == "P")]
monthly = a.groupby(a.t.dt.month).dt.median()
within = {f"{t:g} s": f"{(a.dt.abs() <= t).mean():.2f}" for t in (1, 30, 600)}
npicks = (camp_xo.tid == tid).sum()
print(f" {tid:9s} recall {weak.loc[tid, 'recall']:.2f} over {weak.loc[tid, 'n']:3d} picks; "
f"campaign wrote {npicks:,} picks; within {within}")
print(f" {', '.join(f'{datetime.date(1900, m, 1):%b} {v:+.0f}' for m, v in monthly.items())}")
manual P recall per station: median 0.919, quartiles 0.858-0.949, 56 of 64 stations at or above 0.8
6 stations below 0.7. Month-by-month median offset of the nearest campaign P pick (s):
XO.WD46. recall 0.00 over 326 picks; campaign wrote 2,292 picks; within {'1 s': '0.00', '30 s': '0.01', '600 s': '0.08'}
Jul -173, Aug -2866, Sep -2532, Oct +4954, Nov +5623, Dec -1934
XO.WD47. recall 0.01 over 81 picks; campaign wrote 2,264 picks; within {'1 s': '0.01', '30 s': '0.01', '600 s': '0.02'}
Jul -25, Aug -1072, Sep -43697, Oct +10474, Nov -2455, Dec -19574
XO.WS72. recall 0.34 over 339 picks; campaign wrote 20,909 picks; within {'1 s': '0.34', '30 s': '0.37', '600 s': '0.50'}
Jul +0, Aug +0, Sep +102, Oct +419, Nov -6366, Dec -3650
XO.LA26. recall 0.56 over 39 picks; campaign wrote 5,956 picks; within {'1 s': '0.56', '30 s': '0.62', '600 s': '0.72'}
May +0, Jun +0, Jul -0, Aug +0, Sep +0
XO.LA23. recall 0.66 over 53 picks; campaign wrote 50,834 picks; within {'1 s': '0.66', '30 s': '0.66', '600 s': '0.68'}
May -1726, Jun -0, Jul +0, Aug -0, Sep -0, Dec -0
XO.LT01. recall 0.67 over 151 picks; campaign wrote 161,244 picks; within {'1 s': '0.67', '30 s': '0.70', '600 s': '0.97'}
May +0, Jun -0, Jul +0, Aug +0, Sep +0, Oct +0, Nov -0, Dec -0
9.5 The station-days the campaign does not hold¶
An analyst arrival on a station-day with no campaign pick at all is not a
miss - a whole day of OBS data never produces zero picks at a 0.2 floor - it is
a day the campaign did not read. Each such day belongs to a shard, and each
shard left a completion record with the number of station-days it wrote
(picks_record) against the number it was planned for.
shards_file = CACHE / "obs_shards.jsonl"
if not shards_file.exists():
s3.download_file(BUCKET, f"{CAMPAIGN}/shards.jsonl", str(shards_file))
shards = [json.loads(l) for l in shards_file.open()]
def _yd(s):
y, d = s.split(".")
return datetime.date(int(y), 1, 1) + datetime.timedelta(days=int(d) - 1)
by_station = {}
for sh in shards:
for t in sh["stations"]:
by_station.setdefault(t, []).append(sh)
gaps = pd.concat([arr_obs[~arr_obs.covered], arr_land[~arr_land.covered]]).groupby(["tid", "day"]).size().reset_index()
gaps["day"] = gaps.day.dt.date
gaps["shard"] = [next((s["shard_id"] for s in by_station.get(r.tid, []) if _yd(s["start"]) <= r.day < _yd(s["end"])), None)
for r in gaps.itertuples()]
print(f"{len(gaps)} uncovered station-days across OBS and land, {gaps.shard.notna().sum()} inside a planned shard, "
f"{gaps.shard.nunique()} shards")
recs = {}
for sid in gaps.shard.dropna().unique():
try:
recs[sid] = json.loads(s3.get_object(Bucket=BUCKET, Key=f"{CAMPAIGN}/complete/{sid}.json")["Body"].read())
except s3.exceptions.NoSuchKey:
recs[sid] = {"picks_record": np.nan, "station_days": np.nan, "seconds": np.nan}
tab = pd.DataFrame([dict(shard=k, planned=r.get("station_days"), written=r.get("picks_record"),
seconds=r.get("seconds"), uncovered_here=(gaps.shard == k).sum())
for k, r in recs.items()]).sort_values("uncovered_here", ascending=False)
tab["written_frac"] = (tab.written / tab.planned).round(2)
print(tab.to_string(index=False))
worst = tab.iloc[0].shard
sh = next(s for s in shards if s["shard_id"] == worst)
days = pd.date_range(_yd(sh["start"]), _yd(sh["end"]) - datetime.timedelta(days=1))
held = set(zip(camp_xo.tid, camp_xo.peak.dt.floor("D")))
grid = pd.DataFrame([[(tid, d) in held for d in days] for tid in sh["stations"]],
index=sh["stations"], columns=[f"{d:%m-%d}" for d in days])
print(f"\nshard {worst}: {len(sh['stations'])} stations x {len(days)} days. Stations with picks, per day:")
print(grid.sum(axis=0).to_string())
519 uncovered station-days across OBS and land, 519 inside a planned shard, 23 shards
shard planned written seconds uncovered_here written_frac 2018239-2018259-0cae12e6b282 380 116 1460.7 86 0.31 2018179-2018199-fe9201992769 380 100 1461.1 71 0.26 2018279-2018299-7cfbaa476119 380 220 3483.5 68 0.58 2018359-2019014-be62b51256aa 580 83 1353.9 57 0.14 2018239-2018259-69b557f9b5da 580 340 5003.2 53 0.59 2018199-2018219-a34a2aaa9907 380 300 3653.5 33 0.79 2018219-2018239-a12d8e3d1d88 580 20 317.9 29 0.03 2018299-2018319-bb658359bfb5 580 20 211.3 17 0.03 2018299-2018319-8a78748e7f1a 160 40 493.8 16 0.25 2018339-2018359-744ce5cc5167 140 120 1791.8 13 0.86 2018299-2018319-e90c9487a633 140 120 2033.2 11 0.86 2018319-2018339-5cf7f0be5c8c 140 120 1661.5 11 0.86 2018279-2018299-894ab4397d9a 140 120 2058.2 9 0.86 2018339-2018359-b5a7fb8ff313 140 120 1211.4 9 0.86 2018299-2018319-3971f13e7dda 380 240 2867.8 9 0.63 2018259-2018279-e1875a586b2b 140 120 2106.8 7 0.86 2018239-2018259-0bd6cccdf239 140 113 2068.8 6 0.81 2018359-2019014-c1e9321013a6 140 42 630.3 4 0.30 2018279-2018299-98be32a218f0 140 120 1789.2 3 0.86 2018319-2018339-73dd30eeeb5e 140 120 1678.3 3 0.86 2018339-2018359-2b2d59c14ded 580 580 6712.0 2 1.00 2018279-2018299-7b063bd2e8bd 580 580 8101.3 1 1.00 2018359-2019014-bd47425ffd2b 140 42 466.6 1 0.30
shard 2018239-2018259-0cae12e6b282: 19 stations x 20 days. Stations with picks, per day: 08-27 16 08-28 16 08-29 16 08-30 16 08-31 16 09-01 0 09-02 0 09-03 0 09-04 0 09-05 0 09-06 0 09-07 4 09-08 19 09-09 19 09-10 19 09-11 19 09-12 19 09-13 19 09-14 19 09-15 19
9.6 What the campaign has that the analyst does not¶
The other direction cannot be scored, but it can be sized. On the same station-days, the campaign wrote far more P picks than the analysts kept - which is expected, since the catalogue holds events large enough to locate across the network and an OBS also records everything small and local. The ratio is reported here so nobody reads the recall above as a precision.
cp = camp_xo[camp_xo.pha == "P"].assign(day=camp_xo.peak.dt.floor("D"))
per_day = cp.groupby(["tid", "day"]).agg(at_02=("conf", "size"), at_05=("conf", lambda c: (c >= 0.5).sum()),
at_07=("conf", lambda c: (c >= 0.7).sum()))
analyst_days = arr_obs[arr_obs.pha == "P"].groupby(["tid", "day"]).size().rename("analyst")
j = analyst_days.to_frame().join(per_day, how="inner")
print(f"{len(j):,} OBS station-days with at least one analyst P pick. Median P picks that day:")
print(f" analyst {j.analyst.median():.0f} campaign at >=0.2: {j.at_02.median():.0f} "
f">=0.5: {j.at_05.median():.0f} >=0.7: {j.at_07.median():.0f}")
print(f" campaign-to-analyst ratio at the floor: {j.at_02.sum() / j.analyst.sum():.0f}x; at 0.7: {j.at_07.sum() / j.analyst.sum():.1f}x")
5,874 OBS station-days with at least one analyst P pick. Median P picks that day: analyst 2 campaign at >=0.2: 90 >=0.5: 11 >=0.7: 4 campaign-to-analyst ratio at the floor: 69x; at 0.7: 3.2x
10. The five models re-scored against analyst picks¶
Section 5 scored every model against an iasp91 prediction with a 10 s tolerance, because nothing better existed for OBS stations. For AACSE it now does. The same windows-around-an-arrival procedure runs again on a random sample of the analyst's manual P picks - four components fetched live from EarthScope, every model in section 3 on identical data - and a detection is now a pick within 1 s of the analyst, not 10 s of a prediction.
Two things this adds over section 9. The campaign only ever ran one weight
set, so this is the only place the other four are held to the analyst; and
because the campaign's stored pick for each window is also to hand, the
stored pick and a live pickblue_phasenet pick on the same 3-minute window
can be compared directly: a reproduction check across a different data
path, the same idea as the western validation, on a window rather than a
full day of context.
import warnings
warnings.filterwarnings("ignore", message="Discarding nonzero nanoseconds")
N_RESCORE = 400
rescore_file = CACHE / "rescoring_aacse_2018.parquet"
# The hydrophone code per station, from the inventory - asking for ?DH also returns 1 Hz
# hydrophones, which the OBS model's 0.5 Hz high-pass cannot be applied to.
inv = _es.get_stations(network="XO", starttime=UTCDateTime("2018-05-01"), endtime=UTCDateTime("2019-01-01"),
level="channel")
hydro = {}
for s in inv[0]:
codes = {c.code for c in s if c.code.endswith("DH") and c.code[0] in "HEB"}
if codes:
hydro[s.code] = sorted(codes, key=lambda c: "HEB".index(c[0]))[0]
pool = ref_obs[ref_obs.trace_P_status == "manual"].copy()
pool["tP"] = pd.to_datetime(pool.trace_start_time, utc=True) + pd.to_timedelta(pool.trace_P_arrival_sample / pool.trace_sampling_rate_hz, unit="s")
pool["tS"] = pd.to_datetime(pool.trace_start_time, utc=True) + pd.to_timedelta(pool.trace_S_arrival_sample / pool.trace_sampling_rate_hz, unit="s")
pool.loc[pool.trace_S_status != "manual", "tS"] = pd.NaT
sample = pool.sample(min(N_RESCORE, len(pool)), random_state=2026)
stored = {k: (_sec(g.peak), g.conf.values) for k, g in camp_xo.sort_values("peak").groupby(["tid", "pha"])}
def best_pick(picks, pha, tref):
cand = [(p.peak_time - tref, float(p.peak_value)) for p in picks if p.phase == pha]
if not cand:
return np.nan, 0.0
return min(cand, key=lambda c: abs(c[0]))
if rescore_file.exists():
rescore = pd.read_parquet(rescore_file)
print(f"{len(rescore)} rows from cache")
else:
rows, t_start = [], time.time()
for i, r in enumerate(sample.itertuples(), 1):
tp = UTCDateTime(r.tP.to_pydatetime())
ts = UTCDateTime(r.tS.to_pydatetime()) if pd.notna(r.tS) else None
chans = "HH?" + (f",{hydro[r.station_code]}" if r.station_code in hydro else "")
try:
st = _es.get_waveforms("XO", r.station_code, "*", chans, tp - PRE, tp + POST)
except Exception:
continue
st.merge(fill_value=0)
if not {"Z", "1", "2"} <= {tr.stats.channel[-1] for tr in st}:
continue
st3 = as_three_component(st)
base = dict(source_id=r.source_id, station=r.station_code, mag=r.source_magnitude,
instrument=INSTR.get(r.trace_instrument_type, r.trace_instrument_type),
has_h=any(tr.stats.channel[-1] == "H" for tr in st), has_s=ts is not None)
for name, model in models.items():
try:
out = model.classify(st if needs_h[name] else st3, P_threshold=DETECT_FLOOR, S_threshold=DETECT_FLOOR)
except Exception:
continue
dtp, cp_ = best_pick(out.picks, "P", tp)
dts, cs_ = best_pick(out.picks, "S", ts) if ts else (np.nan, np.nan)
rows.append(dict(base, weights=name, P_dt=dtp, P_conf=cp_, S_dt=dts, S_conf=cs_))
tid = f"XO.{r.station_code}."
rec = dict(base, weights="campaign (stored)")
for pha, tref in (("P", tp), ("S", ts)):
if tref is None or (tid, pha) not in stored:
rec[f"{pha}_dt"], rec[f"{pha}_conf"] = np.nan, np.nan
continue
dt, conf = nearest(np.array([tref.timestamp]), *stored[(tid, pha)])
rec[f"{pha}_dt"], rec[f"{pha}_conf"] = float(dt[0]), float(conf[0])
rows.append(rec)
if i % 50 == 0:
print(f" {i}/{len(sample)} windows, {time.time() - t_start:.0f} s", flush=True)
rescore = pd.DataFrame(rows)
rescore.to_parquet(rescore_file)
print(f"{rescore.source_id.nunique()} windows scored in {time.time() - t_start:.0f} s")
2400 rows from cache
def recall_table(df, thr, tol=TOL_ANALYST):
rows = []
for name, g in df.groupby("weights", sort=False, observed=True):
p = g
s = g[g.has_s]
rows.append(dict(weights=name, windows=len(g),
P=((p.P_dt.abs() <= tol) & (p.P_conf >= thr)).mean(),
S=((s.S_dt.abs() <= tol) & (s.S_conf >= thr)).mean() if len(s) else np.nan,
P_median_dt=p[(p.P_dt.abs() <= tol) & (p.P_conf >= thr)].P_dt.median()))
return pd.DataFrame(rows).set_index("weights")
order = list(names) + ["campaign (stored)"]
rescore["weights"] = pd.Categorical(rescore.weights, order)
rescore = rescore.sort_values("weights")
floor = float(PROD["p_threshold"])
t_report = recall_table(rescore, REPORT_THRESHOLD)
t_floor = recall_table(rescore, floor)
print(f"Recall of analyst picks within {TOL_ANALYST:g} s, {rescore.source_id.nunique()} windows, "
f"{int(rescore.groupby('source_id').has_s.first().sum())} of them with a manual S\n")
print(pd.concat({f"conf >= {REPORT_THRESHOLD}": t_report[["P", "S"]], f"conf >= {floor:g}": t_floor[["P", "S"]],
"median P dt (s)": t_report[["P_median_dt"]]}, axis=1).round(3).to_string())
print("\nThe stored campaign row is scored at its own floor only - it holds nothing below it.")
fig, axes = plt.subplots(1, 2, figsize=(11, 3.9), gridspec_kw=dict(width_ratios=[1.35, 1]))
ax = axes[0]
x = np.arange(len(order)); w = 0.38
for k, (pha, color) in enumerate((("P", C_P), ("S", C_S))):
vals = t_floor[pha].values
ax.bar(x + (k - 0.5) * w, vals, width=w * 0.92, color=color, label=pha)
for xi, v in zip(x + (k - 0.5) * w, vals):
ax.text(xi, v + 0.012, f"{v:.2f}", ha="center", fontsize=7.5, color="#3d3d3d")
ax.set_xticks(x); ax.set_xticklabels(order, fontsize=8, rotation=15, ha="right")
ax.set_ylim(0, 1.2); ax.set_ylabel(f"recall within {TOL_ANALYST:g} s at conf >= {floor:g}")
ax.grid(alpha=0.25, lw=0.5, axis="y"); ax.set_axisbelow(True); ax.legend(frameon=False, fontsize=9, loc="upper right", ncol=2)
ax.set_title("Against the analyst, same windows, every model", fontsize=10.5, loc="left")
ax = axes[1]
bins = np.arange(-1, 1.0001, 0.05)
for i, name in enumerate(order):
d = rescore[(rescore.weights == name) & (rescore.P_conf >= floor)].P_dt
d = d[d.abs() <= 1]
ax.hist(d, bins=bins, histtype="step", lw=1.6, color=PALETTE[i % len(PALETTE)] if i < len(PALETTE) else "#3d3d3d",
label=f"{name} ({d.median():+.2f})")
ax.axvline(0, color="#8a8a8a", lw=1); ax.set_yticks([])
ax.set_xlabel("P pick minus analyst P (s)"); ax.grid(alpha=0.25, lw=0.5); ax.set_axisbelow(True)
ax.legend(frameon=False, fontsize=7, title="median", title_fontsize=7)
ax.set_title("Where each model's P lands", fontsize=10.5, loc="left")
fig.tight_layout(); plt.show()
# stored campaign pick against a live pickblue_phasenet pick on the same window
live = rescore[rescore.weights == "pickblue_phasenet"].set_index(["source_id", "station"])
st_ = rescore[rescore.weights == "campaign (stored)"].set_index(["source_id", "station"])
both = live.join(st_, lsuffix="_live", rsuffix="_stored", how="inner")
both = both[(both.P_conf_live >= floor) & (both.P_conf_stored >= floor)]
agree = (both.P_dt_live - both.P_dt_stored).abs()
print(f"\nStored campaign pick vs live pickblue_phasenet on the same window, both at conf >= {floor:g}: "
f"{len(both)} windows; |difference| median {agree.median():.3f} s, within 0.1 s {(agree <= 0.1).mean():.1%}, "
f"within 1 s {(agree <= 1).mean():.1%}")
print("The two are not expected to be identical: the campaign picked a full day with SeisBench's overlapping windows "
"and the live run picked 3 minutes, so the window boundaries and the normalisation differ.")
Recall of analyst picks within 1 s, 358 windows, 184 of them with a manual S
conf >= 0.3 conf >= 0.2 median P dt (s)
P S P S P_median_dt
weights
pickblue_phasenet 0.695 0.651 0.770 0.737 0.03
pickblue_eqt 0.652 0.569 0.735 0.675 0.04
obstransformer 0.752 0.842 0.805 0.866 0.07
quakescope2026 0.638 0.292 0.742 0.459 0.06
original 0.448 0.478 0.488 0.560 0.10
campaign (stored) 0.802 0.703 0.825 0.809 0.01
The stored campaign row is scored at its own floor only - it holds nothing below it.
Stored campaign pick vs live pickblue_phasenet on the same window, both at conf >= 0.2: 331 windows; |difference| median 0.170 s, within 0.1 s 39.6%, within 1 s 86.4% The two are not expected to be identical: the campaign picked a full day with SeisBench's overlapping windows and the live run picked 3 minutes, so the window boundaries and the normalisation differ.
11. Axial Seamount: eleven years against the University of Washington catalogue¶
The campaign holds OO picks from the start of the cabled array in
November 2014 through 2025, on seven stations around the caldera (five on
short-period EH, two on broadband HH, all sampled at 200 Hz and
downsampled to the weights' 100 Hz by the campaign) and six more at
Southern Hydrate Ridge that the reference does not cover. The reference is
the phase file behind the UW near-real-time catalogue: automatic P and S
picks with weights, for events located by HYPOINVERSE.
This is not the same test as section 9. The reference is a machine, and the
seismicity is nothing like Alaska's: sources are inside or under a caldera a
few kilometres across, magnitudes are mostly below 0.5, P arrives within a
second of the origin, and S follows it by half a second to a second. That
is outside the regional-OBS regime the obs weights were trained on. The
tolerance is 0.5 s, and everything is reported per station so the two
instrument bands can be read separately.
def parse_ph2dt(path):
"""HypoDD phase format as served by the UW catalogue: '# YR MO DY HR MN SEC LAT LON DEP MAG
EH EZ RMS ID' event headers, then 'STA TT WGHT PHA' lines with travel time from the origin."""
ev, ph, origin, eid = [], [], None, None
with open(path) as f:
for line in f:
p = line.split()
if line.startswith("#"):
y, mo, d, h, mi = map(int, p[1:6])
origin = datetime.datetime(y, mo, d, h, mi) + datetime.timedelta(seconds=float(p[6]))
eid = int(p[14])
ev.append((eid, origin, float(p[7]), float(p[8]), float(p[9]), float(p[10])))
elif len(p) >= 4:
ph.append((eid, p[0], p[3], float(p[1]), float(p[2])))
ev = pd.DataFrame(ev, columns=["eid", "origin", "lat", "lon", "dep", "mag"])
ph = pd.DataFrame(ph, columns=["eid", "sta", "pha", "tt", "wt"]).merge(ev[["eid", "origin", "mag"]], on="eid")
ph["t"] = ph.origin + pd.to_timedelta(ph.tt, unit="s")
ph["tid"] = "OO." + ph.sta + "."
return ev, ph.drop(columns=["origin"])
uw_ev_file, uw_ph_file = CACHE / "uw_axial_events.parquet", CACHE / "uw_axial_picks.parquet"
if uw_ph_file.exists():
uw_ev, uw_ph = pd.read_parquet(uw_ev_file), pd.read_parquet(uw_ph_file)
else:
uw_ev, uw_ph = parse_ph2dt(files["uw_axial_phases"])
uw_ev.to_parquet(uw_ev_file); uw_ph.to_parquet(uw_ph_file)
print(f"UW catalogue as fetched: {len(uw_ev):,} events, {len(uw_ph):,} picks, "
f"{uw_ev.origin.min():%Y-%m-%d} .. {uw_ev.origin.max():%Y-%m-%d}; magnitude median {uw_ev.mag.median():.1f}")
print(f"UW travel times: P median {uw_ph[uw_ph.pha == 'P'].tt.median():.2f} s, S median {uw_ph[uw_ph.pha == 'S'].tt.median():.2f} s after the origin")
lo, hi = uw_ph.t.min().floor("D"), pd.Timestamp("2026-01-01")
ax_camp = camp_oo[(camp_oo.peak >= lo) & (camp_oo.peak < hi) & camp_oo.tid.str.startswith("OO.AX")]
ax_ref = uw_ph[(uw_ph.t >= lo) & (uw_ph.t < hi)][["tid", "pha", "t", "wt", "mag", "eid"]]
t0 = time.time()
arr_ax = match_to_campaign(ax_ref, ax_camp, TOL_AXIAL)
arr_ax["cha"] = arr_ax.tid.map(ax_camp.groupby("tid").cha.first())
print(f"matched {len(arr_ax):,} UW picks in {time.time() - t0:.0f} s; coverage {arr_ax.covered.mean():.1%}\n")
ndays = (hi - lo).days
rate = pd.DataFrame({
"band": ax_camp.groupby("tid").cha.first(),
"campaign P/day": ax_camp[ax_camp.pha == "P"].groupby("tid").size() / ndays,
"campaign P/day >=0.5": ax_camp[(ax_camp.pha == "P") & (ax_camp.conf >= 0.5)].groupby("tid").size() / ndays,
"campaign S/day": ax_camp[ax_camp.pha == "S"].groupby("tid").size() / ndays,
"UW P/day": ax_ref[ax_ref.pha == "P"].groupby("tid").size() / ndays,
"UW P/day wt>=0.5": ax_ref[(ax_ref.pha == "P") & (ax_ref.wt >= 0.5)].groupby("tid").size() / ndays,
"UW S/day": ax_ref[ax_ref.pha == "S"].groupby("tid").size() / ndays,
}).round(1)
print("Pick rates, both sides:"); print(rate.to_string())
covx = arr_ax[arr_ax.covered]
print(f"\nRecall of UW picks within {TOL_AXIAL:g} s, per station (covered station-days):")
print(covx.groupby(["tid", "cha", "pha"]).hit.agg(recall="mean", n="size").round(3).unstack("pha").to_string())
print("\nBy UW pick weight, P:")
wp = covx[covx.pha == "P"]
print(wp.groupby(pd.cut(wp.wt, [-0.01, 0.05, 0.3, 0.6, 1.0], labels=["~0", "0.05-0.3", "0.3-0.6", "0.6-1"]), observed=True)
.hit.agg(recall="mean", n="size").round(3).to_string())
print("\nBy tolerance, P and S, per band:")
print(pd.DataFrame({f"{t:g} s": covx.groupby(["cha", "pha"]).dt.apply(lambda d: (d.abs() <= t).mean()) for t in (0.1, 0.25, 0.5, 1.0)}).round(3).to_string())
print("\nResidual on hits, campaign minus UW (s):")
print(covx[covx.hit].groupby(["cha", "pha"]).dt.agg(median="median", q25=lambda x: x.quantile(.25), q75=lambda x: x.quantile(.75), n="size").round(3).to_string())
# The other direction: how many of the campaign's own P picks does the UW system also have?
rev = []
for tid, g in ax_camp[ax_camp.pha == "P"].groupby("tid"):
u = np.sort(_sec(ax_ref[(ax_ref.tid == tid) & (ax_ref.pha == "P")].t))
dt, _ = nearest(_sec(g.sort_values("peak").peak), u, np.zeros(len(u)))
rev.append(g.sort_values("peak").assign(hit=np.abs(dt) <= TOL_AXIAL))
rev = pd.concat(rev)
print("\nCampaign P picks with a UW P within 0.5 s, by the campaign's confidence:")
print(rev.groupby(["cha", pd.cut(rev.conf, [0.2, 0.3, 0.5, 0.7, 1.0])], observed=True).hit.agg(agree="mean", n="size").round(3).to_string())
UW catalogue as fetched: 319,076 events, 3,468,987 picks, 2015-01-22 .. 2026-09-07; magnitude median 0.1
UW travel times: P median 0.58 s, S median 1.19 s after the origin
matched 3,188,218 UW picks in 3 s; coverage 90.3%
Pick rates, both sides:
band campaign P/day campaign P/day >=0.5 campaign S/day UW P/day UW P/day wt>=0.5 UW S/day
tid
OO.AXAS1. EH 64.4 11.5 16.7 56.4 34.7 54.7
OO.AXAS2. EH 43.4 7.4 47.2 47.6 26.1 41.1
OO.AXBA1. HH 21.8 2.7 2.1 NaN NaN NaN
OO.AXCC1. HH 87.9 15.6 26.7 53.9 38.4 55.5
OO.AXEC1. EH 61.8 10.0 34.5 65.4 47.1 71.2
OO.AXEC2. HH 36.5 7.7 4.9 71.4 57.9 71.8
OO.AXEC3. EH 66.9 10.9 22.6 68.4 50.9 72.5
OO.AXID1. EH 56.2 7.9 34.2 34.4 18.5 33.4
Recall of UW picks within 0.5 s, per station (covered station-days):
recall n
pha P S P S
tid cha
OO.AXAS1. EH 0.291 0.013 213759 207824
OO.AXAS2. EH 0.207 0.046 180703 156651
OO.AXCC1. HH 0.342 0.014 193977 206532
OO.AXEC1. EH 0.174 0.053 247264 269959
OO.AXEC2. HH 0.255 0.006 203643 205433
OO.AXEC3. EH 0.192 0.057 259198 274872
OO.AXID1. EH 0.170 0.017 130032 128209
By UW pick weight, P:
recall n
wt
~0 0.105 138205
0.05-0.3 0.088 188076
0.3-0.6 0.169 275250
0.6-1 0.309 827045
By tolerance, P and S, per band:
0.1 s 0.25 s 0.5 s 1 s
cha pha
EH P 0.148 0.197 0.208 0.213
S 0.025 0.035 0.041 0.050
HH P 0.198 0.282 0.298 0.304
S 0.002 0.005 0.010 0.022
Residual on hits, campaign minus UW (s):
median q25 q75 n
cha pha
EH P -0.060 -0.103 -0.030 214541
S -0.070 -0.151 -0.020 42146
HH P -0.070 -0.115 -0.035 118391
S -0.215 -0.359 -0.096 4151
Campaign P picks with a UW P within 0.5 s, by the campaign's confidence:
agree n
cha conf
EH (0.2, 0.3] 0.097 528355
(0.3, 0.5] 0.181 451636
(0.5, 0.7] 0.331 102322
(0.7, 1.0] 0.534 88078
HH (0.2, 0.3] 0.096 274414
(0.3, 0.5] 0.210 205564
(0.5, 0.7] 0.389 59461
(0.7, 1.0] 0.566 44673
good = covx[(covx.pha == "P") & (covx.wt >= 0.5)]
fig, axes = plt.subplots(1, 2, figsize=(11, 3.8))
ax = axes[0]
mbins = [-3, -0.5, 0, 0.5, 1, 1.5, 2, 5]
mlab = ["<-0.5", "-0.5-0", "0-0.5", "0.5-1", "1-1.5", "1.5-2", ">2"]
for cha, color in (("EH", PALETTE[0]), ("HH", PALETTE[1])):
g = good[good.cha == cha].groupby(pd.cut(good[good.cha == cha].mag, mbins, labels=mlab), observed=True).hit.agg(["mean", "size"])
ax.plot(range(len(g)), g["mean"], marker="o", ms=5, lw=2, color=color, label=f"{cha} band")
for i, (v, n) in enumerate(zip(g["mean"], g["size"])):
ax.text(i, v + (0.03 if cha == "HH" else -0.06), f"{n:,}", fontsize=6.5, ha="center", color="#7a7973")
ax.set_xticks(range(len(mlab))); ax.set_xticklabels(mlab, fontsize=8.5); ax.set_xlabel("UW magnitude")
ax.set_ylim(0, 1); ax.set_ylabel(f"recall of UW P picks (weight >= 0.5) within {TOL_AXIAL:g} s")
ax.grid(alpha=0.25, lw=0.5); ax.set_axisbelow(True); ax.legend(frameon=False, fontsize=9)
ax.set_title("Recall rises with magnitude, then falls again", fontsize=10.5, loc="left")
ax = axes[1]
for cha, color in (("EH", PALETTE[0]), ("HH", PALETTE[1])):
g = good[good.cha == cha].groupby(good[good.cha == cha].t.dt.year).hit.agg(["mean", "size"])
ax.plot(g.index, g["mean"], marker="o", ms=5, lw=2, color=color, label=f"{cha} band")
ax.set_ylim(0, 1); ax.set_xlabel("year"); ax.grid(alpha=0.25, lw=0.5); ax.set_axisbelow(True)
ax.legend(frameon=False, fontsize=9)
ax.set_title("Year by year (2015 is the eruption year)", fontsize=10.5, loc="left")
fig.tight_layout(); plt.show()
11.1 Event-level detection against the ML-DD catalogue¶
The LDEO catalogue carries no picks, so it is scored as events: an event counts as detected if the campaign has a P pick on at least three of the seven stations in the three seconds after its origin. The same test on random origin times gives the chance rate.
rows = []
for line in files["ldeo_axial_catalogue"].open():
p = line.split()
if len(p) == 15 and p[0].isdigit() and len(p[0]) == 4:
y, mo, d, h, mi = map(int, p[:5])
rows.append((datetime.datetime(y, mo, d, h, mi) + datetime.timedelta(seconds=float(p[5])),
float(p[6]), float(p[7]), float(p[8]), float(p[13]), int(p[14])))
ldeo = pd.DataFrame(rows, columns=["origin", "lat", "lon", "dep", "mag", "id"])
print(f"LDEO ML-DD: {len(ldeo):,} events, {ldeo.origin.min():%Y-%m-%d} .. {ldeo.origin.max():%Y-%m-%d}, "
f"magnitude median {ldeo.mag.median():.1f}")
o = _sec(ldeo.origin)
oo_p = camp_oo[(camp_oo.pha == "P") & camp_oo.tid.str.startswith("OO.AX")]
def stations_within(origins, picks, after=(0.05, 3.0)):
c = np.zeros(len(origins), int)
for _, g in picks.groupby("tid"):
ct = np.sort(_sec(g.peak))
c += np.searchsorted(ct, origins + after[1]) > np.searchsorted(ct, origins + after[0])
return c
ldeo["n_sta"] = stations_within(o, oo_p)
ldeo["n_sta_05"] = stations_within(o, oo_p[oo_p.conf >= 0.5])
rng = np.random.default_rng(0)
chance = stations_within(rng.uniform(o.min(), o.max(), 20000), oo_p)
print(f"detected (P on >=3 stations within 3 s): {(ldeo.n_sta >= 3).mean():.1%}; with conf >= 0.5: {(ldeo.n_sta_05 >= 3).mean():.1%}; "
f"on >=1 station: {(ldeo.n_sta >= 1).mean():.1%}. Chance at random times: >=3 {(chance >= 3).mean():.2%}, >=1 {(chance >= 1).mean():.2%}")
print(ldeo.assign(magnitude=pd.cut(ldeo.mag, mbins, labels=mlab)).groupby("magnitude", observed=True)
.agg(detected=("n_sta", lambda x: (x >= 3).mean()), detected_conf05=("n_sta_05", lambda x: (x >= 3).mean()), n=("n_sta", "size")).round(3).to_string())
LDEO ML-DD: 144,329 events, 2014-11-04 .. 2021-12-30, magnitude median 0.1
detected (P on >=3 stations within 3 s): 18.0%; with conf >= 0.5: 7.3%; on >=1 station: 41.3%. Chance at random times: >=3 0.11%, >=1 1.07%
detected detected_conf05 n
magnitude
<-0.5 0.007 0.000 284
-0.5-0 0.032 0.001 53735
0-0.5 0.217 0.058 68991
0.5-1 0.500 0.345 13659
1-1.5 0.398 0.306 4588
1.5-2 0.211 0.145 1948
>2 0.474 0.286 325
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/obs_offshore"); _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("detection_vs_iasp91", summary)
_aacse = []
for _label, _arr in (("OBS", arr_obs), ("land", arr_land)):
_m = _arr[_arr.covered & (_arr.status == "manual")]
for _pha, _g in _m.groupby("pha"):
_aacse.append(dict(stations=_label, pha=_pha, n=len(_g), coverage=float(_arr[_arr.pha == _pha].covered.mean()),
**{f"recall_{t:g}s": float((_g.dt.abs() <= t).mean()) for t in (0.25, 0.5, 1.0, 2.0)},
median_residual_s=float(_g[_g.dt.abs() <= TOL_ANALYST].dt.median())))
_export("aacse_campaign_vs_analyst", pd.DataFrame(_aacse))
_export("aacse_per_station_P", per.reset_index())
_export("aacse_rescoring", pd.concat({f"conf>={REPORT_THRESHOLD}": t_report[["P", "S"]],
f"conf>={floor:g}": t_floor[["P", "S"]],
"median_P_dt_s": t_report[["P_median_dt"]]}, axis=1).reset_index())
_export("axial_pick_rates", rate.reset_index())
_export("axial_recall_per_station", covx.groupby(["tid", "cha", "pha"]).hit.agg(recall="mean", n="size").reset_index())
_export("axial_event_detection_by_magnitude",
ldeo.assign(magnitude=pd.cut(ldeo.mag, mbins, labels=mlab)).groupby("magnitude", observed=True)
.agg(detected=("n_sta", lambda x: (x >= 3).mean()), detected_conf05=("n_sta_05", lambda x: (x >= 3).mean()),
n=("n_sta", "size")).reset_index())
_meta = dict(notebook="phasenet_obs_offshore_benchmark.ipynb",
executed=_dt.datetime.now(_dt.timezone.utc).isoformat(timespec="seconds"),
seisbench=_sb.__version__, campaign=CAMPAIGN, production=PROD,
tol_analyst_s=TOL_ANALYST, tol_axial_s=TOL_AXIAL, report_threshold=REPORT_THRESHOLD,
experiments={k: dict(network=v["network"], start=str(v["start"]), end=str(v["end"]),
stations=[s[0] for s in v["stations"]]) for k, v in EXPERIMENTS.items()})
_export("meta", _meta)
print("wrote", ", ".join(_written), "to", _OUT)
wrote detection_vs_iasp91, aacse_campaign_vs_analyst, aacse_per_station_P, aacse_rescoring, axial_pick_rates, axial_recall_per_station, axial_event_detection_by_magnitude, meta to ../docs/benchmark/results/obs_offshore
Reading the result¶
On the AACSE array the campaign's picks agree with the analysts', and the
disagreement that exists is in the data, not the picker. Across 65 OBS
stations and 2018, 86% of the manual P arrivals and 86% of the manual S
arrivals have a campaign pick within 1 s; 73% and 59% within a quarter
second. The residual is centred: median +0.02 s for P and +0.04 s for S,
with quartiles inside ±0.1 s for P. The 30 land stations of the same
network, picked with the same obs weight, come out at 91% and 85% with the
same residual shape, so ocean-bottom recording costs a few points on P and
nothing on S. The two makes of deep-water instrument differ more than sea
and land do: the WHOI instruments are ten points below the LDEO ones on both
phases. Recall is flat in magnitude and distance out to 350 km, which is the
range the reference covers.
Two stations have zero recall, and neither is a picker failure. WD46 and WD47 wrote a handful of picks a day where their neighbours wrote hundreds, and the offset between the nearest campaign pick and the analyst's wanders by hours from month to month. WS72 is fine until August and then drifts by minutes. Of the rest, only LT01 is a picker shortfall: its offsets are centred and it still misses a third of the analyst's P picks, on a station that wrote 160,000 picks in the year. That is the signature of an archive whose clock correction is not the one the analysts worked with - the Barcheck readme notes that drift was removed "for all ocean-bottom seismometers that were still running and able to get a clock lock upon recovery" and refers the rest to a supplement. The campaign read the archive as it is. Anyone using these picks for location should treat those three stations' 2018 timing as unverified.
Five percent of the analyst arrivals fall on station-days the campaign does not hold, in runs of consecutive days inside shards that reported themselves complete having written fewer station-days than they were planned for - the worst of them 31%. That is the silent-skip pattern already documented for the western campaign, now measured on this one: the day-level coverage a completion record implies is not the coverage the Parquet has.
Against analyst picks rather than a prediction, the ordering of section 5
holds and the gaps widen. On 358 AACSE windows scored at the campaign's
0.2 floor, obstransformer recovers 81% of the analyst's P and 87% of the
S, the two PickBlue weights 74–77% of P and 68–74% of S, quakescope2026
matches them on P but finds under half the S, and original finds under
half of either. Training on ocean-bottom data is worth 25–30 points of P
recall here, far more than the 10 s tolerance of section 5 could show. The
stored campaign pick beats a live pickblue_phasenet run on the same
3-minute window, 83% against 77% on P, and the two land within 1 s of each
other on 86% of the windows where both exist but within 0.1 s on only 40%.
A whole day of context, with SeisBench's overlapping windows, is not the
same picker as a 3-minute cut, which is why section 9 scores the stored
picks and not a re-pick.
Axial Seamount is where agreement breaks down, and the reason is the
seismicity, not the archive. The campaign and the UW real-time system
each write 40–90 P picks a day per station, but only 20–35% of the UW picks
have a campaign pick within half a second, and only half of the campaign's
most confident picks have a UW counterpart. Recall climbs to about 0.5 for
M 0.5–1.5 and collapses below M 0 - which is most of the catalogue. S is
essentially unmatched: S–P here is 0.5–0.7 s, and the obs weights, trained
on regional OBS records at 100 Hz from data the campaign downsampled from
200 Hz, do not resolve it. At the event level the campaign sees 18% of the
ML-DD catalogue's 144,000 events on three or more stations, 50% of those
between M 0.5 and 1, against a chance rate of 0.1%. The campaign catalogue
at Axial is real but much shallower than the dedicated systems produce, and
it should not be used to say anything about the sub-M0 seismicity there.
The comparison in sections 5–7 remains relative, not absolute. Those sections score against iasp91 with a 10 s tolerance; section 10 shows what that was hiding by scoring the same models against analyst picks at 1 s.
The land models are being used out of domain on purpose. That is the
measurement: QuakeScope will encounter ocean-bottom data, and the question
is what running a land-trained picker over it costs. Renaming the
horizontals from 1/2 to N/E is the usual convention and is
harmless here, since none of these models uses absolute orientation.
Thresholds still belong to the weight set. Section 9.3 gives the curve
for obs: raising the campaign's floor from 0.2 to 0.5 keeps 76% of the
analyst's P picks but only 52% of the S picks. obs ships different
defaults from the rest, P 0.2 and S 0.1, which is itself a signal that its
probabilities are not on the same scale - and the S curve says the campaign's
S=0.2 is already above the weight's own default.
Blanco is still open. The GJI comparison of EQTransformer, PickBlue and OBSTransformer on the 2012–13 deployment puts its picks behind a request to the corresponding author, and the campaign holds only 2012 of it. Both halves of that arm are a wait, not a computation.