PhaseNet smoke test — 2019 Ridgecrest¶
Runs a PhaseNet weight set over five SCSN stations around the 2019 Ridgecrest sequence and checks the picks against the waveforms.
The question this answers: do the picks land on real arrivals, in the right order, with moveout that matches distance? It is a sanity check you can run on a laptop in about two minutes, not a benchmark.
Two events are used, deliberately:
| Event | Why |
|---|---|
| M7.1 mainshock, 2019-07-06 03:19:53 UTC | Strong signal everywhere; good for the record section |
| M4.6 aftershock, 2019-07-06 08:32:58 UTC | Short, impulsive source — the honest test of S picking |
A magnitude 7 ruptures for tens of seconds, so at these distances its S arrival is buried inside ongoing rupture radiation and pickers routinely miss it. That is expected behaviour rather than a broken setup, which is why the S−P timing check below uses the moderate aftershock instead.
Data comes from the SCEDC public S3 bucket (scedc-pds) over anonymous
access — no AWS credentials required.
import io
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.geodetics import gps2dist_azimuth
from s3fs import S3FileSystem
%matplotlib inline
1. Configuration¶
# --- Events (USGS catalog) --------------------------------------------------
# Both fall on 2019 day-of-year 187, so they share the same SCEDC day files.
MAINSHOCK = dict(
name="M7.1 mainshock",
time=UTCDateTime("2019-07-06T03:19:53.040000Z"),
lat=35.770, lon=-117.599, depth_km=8.0,
)
AFTERSHOCK = dict(
name="M4.6 aftershock",
time=UTCDateTime("2019-07-06T08:32:57.550000Z"),
lat=35.639, lon=-117.491, depth_km=3.1,
)
YEAR, DOY = 2019, 187
# --- Stations: CI (SCSN), all verified present in SCEDC on this day --------
STATIONS = [
("CLC", 35.8157, -117.5975),
("TOW2", 35.8086, -117.7649),
("SRT", 35.6923, -117.7505),
("WRC2", 35.9479, -117.6504),
("JRC2", 35.9825, -117.8089),
]
NETWORK, CHANNEL = "CI", "HH" # broadband, 100 Hz
P_THRESHOLD = S_THRESHOLD = 0.3
VP, VS = 6.0, 3.5 # crustal averages for reference curves
# Categorical colors, CVD-validated. Phase is also encoded by linestyle and
# text label, so identity never depends on color alone.
C_P, C_S = "#2a78d6", "#eb6834"
def distance_km(event, lat, lon):
return gps2dist_azimuth(event["lat"], event["lon"], lat, lon)[0] / 1000.0
2. Data access¶
SCEDC stores one whole-day miniSEED file per channel (~20 MB), so each event costs roughly 300 MB of download and a minute or so on a warm connection. Streams are trimmed to the event window immediately.
def scedc_key(sta, comp, net=NETWORK, cha=CHANNEL, year=YEAR, doy=DOY, loc=""):
"""Build an SCEDC S3 object key.
The SCEDC layout differs from NCEDC - the network is not a directory, the
day folder uses an underscore, and the station is padded to five characters:
scedc-pds/continuous_waveforms/<year>/<year>_<doy>/
<net><sta:_<5><cha><comp><loc:_<3><year><doy>.ms
e.g. CICLC__HHZ___2019187.ms
"""
base = f"{net}{sta.ljust(5, '_')}{cha}{comp}{loc.ljust(3, '_')}{year}{doy:03d}.ms"
return f"scedc-pds/continuous_waveforms/{year}/{year}_{doy:03d}/{base}"
def load_event(fs, event, pre=30, post=120):
"""Fetch every station for one event, trimmed to its time window."""
out = {}
for sta, lat, lon 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 FileNotFoundError:
print(f" {sta} {comp}: not found")
except Exception as exc:
print(f" {sta} {comp}: {type(exc).__name__}")
if len(st) == 0:
print(f" {sta}: no data, skipping")
continue
st.merge(fill_value=0) # close gaps into one trace
st.trim(event["time"] - pre, event["time"] + post)
if len(st):
out[sta] = st
return out
def pick_all(model, streams):
picks = {}
for sta, st in streams.items():
out = model.classify(st, P_threshold=P_THRESHOLD, S_threshold=S_THRESHOLD)
picks[sta] = list(out.picks)
return picks
3. Load the picker¶
original is the published Zhu & Beroza (2019) PhaseNet. If the
QuakeScope 2026 weights (quakescope2026) have been converted and
installed locally they are used instead — see
sb_catalog/models/phasenet/README.md for that step.
WEIGHT = "quakescope2026" if "quakescope2026" in sbm.PhaseNet.list_pretrained() else "original"
if WEIGHT == "original":
print("quakescope2026 not installed - using published 'original' weights")
print(" (convert with sb_catalog/models/phasenet/convert_checkpoint.py)")
model = sbm.PhaseNet.from_pretrained(WEIGHT)
print(f"\nWeight set : {WEIGHT}")
print(f"Input : {model.in_samples} samples @ {model.sampling_rate:g} Hz")
Weight set : quakescope2026 Input : 3001 samples @ 100 Hz
4. Mainshock — picks on the waveforms¶
Each pick should sit on the onset of visible energy: P on the first break, clearest on the vertical; S on the later, larger arrival, clearest on the horizontals.
fs = S3FileSystem(anon=True)
print(f"Fetching {MAINSHOCK['name']} ...")
main_streams = load_event(fs, MAINSHOCK)
main_picks = pick_all(model, main_streams)
print(f"\n{len(main_streams)}/{len(STATIONS)} stations loaded\n")
for sta, pk in main_picks.items():
n_p = sum(1 for p in pk if p.phase == 'P')
n_s = sum(1 for p in pk if p.phase == 'S')
print(f" {sta:<6} {len(pk):3d} picks (P={n_p}, S={n_s})")
Fetching M7.1 mainshock ...
5/5 stations loaded CLC 6 picks (P=4, S=2) TOW2 6 picks (P=5, S=1) SRT 4 picks (P=3, S=1) WRC2 9 picks (P=5, S=4) JRC2 8 picks (P=5, S=3)
def plot_station(sta, st, sta_picks, dist, event, window=(-5, 40)):
fig, axes = plt.subplots(3, 1, figsize=(12, 5.5), sharex=True)
for ax, comp in zip(axes, "ZNE"):
tr = st.select(component=comp)
ax.set_ylabel(f"{CHANNEL}{comp}", fontsize=9)
ax.grid(alpha=0.25, lw=0.5)
ax.tick_params(labelsize=9)
if not tr:
continue
tr = tr[0]
ax.plot(tr.times(reftime=event["time"]), tr.data, color="#3d3d3d", lw=0.6)
for p in sta_picks:
dt = p.peak_time - event["time"]
if not (window[0] <= dt <= window[1]):
continue
is_p = p.phase == "P"
ax.axvline(dt, color=C_P if is_p else C_S, lw=1.6,
ls="-" if is_p else "--", alpha=0.9)
for p in sta_picks:
dt = p.peak_time - event["time"]
if window[0] <= dt <= window[1]:
axes[0].annotate(p.phase, (dt, 1.03), xycoords=("data", "axes fraction"),
ha="center", fontsize=9,
color=C_P if p.phase == "P" else C_S)
handles = [plt.Line2D([], [], color=C_P, lw=1.6, ls="-", label="P"),
plt.Line2D([], [], color=C_S, lw=1.6, ls="--", label="S")]
axes[0].legend(handles=handles, loc="upper right", fontsize=9,
frameon=False, ncol=2)
axes[0].set_title(f"{NETWORK}.{sta} - {dist:.1f} km - {event['name']}",
fontsize=11, loc="left")
axes[-1].set_xlim(*window)
axes[-1].set_xlabel("Seconds after origin time")
fig.tight_layout()
return fig
for sta, st in main_streams.items():
lat, lon = next((la, lo) for c, la, lo in STATIONS if c == sta)
plot_station(sta, st, main_picks[sta],
distance_km(MAINSHOCK, lat, lon), MAINSHOCK)
plt.show()
Record section¶
Traces ordered by distance. Real arrivals form a coherent moveout, so picks should follow the reference curves rather than scattering.
def record_section(streams, picks, event, window=(-5, 45), gain=3.0):
fig, ax = plt.subplots(figsize=(11, 6))
dmax = 0
for sta, st in streams.items():
lat, lon = next((la, lo) for c, la, lo in STATIONS if c == sta)
d = distance_km(event, lat, lon)
dmax = max(dmax, d)
tr = st.select(component="Z")
if not tr:
continue
tr = tr[0]
x = tr.data.astype(float)
peak = np.abs(x).max()
if peak > 0:
x = x / peak * gain
ax.plot(tr.times(reftime=event["time"]), x + d, color="#3d3d3d", lw=0.5)
ax.annotate(sta, (window[0] + 0.5, d + 0.8), fontsize=9, color="#52514e")
for p in picks[sta]:
dt = p.peak_time - event["time"]
if not (window[0] <= dt <= window[1]):
continue
is_p = p.phase == "P"
ax.plot([dt, dt], [d - gain * 0.8, d + gain * 0.8],
color=C_P if is_p else C_S, lw=1.8,
ls="-" if is_p else "--")
dd = np.linspace(0, dmax * 1.1, 60)
hyp = np.hypot(dd, event["depth_km"])
ax.plot(hyp / VP, dd, color=C_P, lw=1, alpha=0.45)
ax.plot(hyp / VS, dd, color=C_S, lw=1, alpha=0.45, ls="--")
handles = [
plt.Line2D([], [], color=C_P, lw=1.8, ls="-", label="P pick"),
plt.Line2D([], [], color=C_S, lw=1.8, ls="--", label="S pick"),
plt.Line2D([], [], color="#8a8a8a", lw=1,
label=f"predicted (Vp {VP}, Vs {VS} km/s)"),
]
ax.legend(handles=handles, loc="lower right", fontsize=9, frameon=False)
ax.set_xlim(*window)
ax.set_ylim(-gain, dmax * 1.15)
ax.set_xlabel("Seconds after origin time")
ax.set_ylabel("Epicentral distance (km)")
ax.set_title(f"Record section - vertical component - {event['name']}",
fontsize=11, loc="left")
ax.grid(alpha=0.25, lw=0.5)
fig.tight_layout()
return fig
record_section(main_streams, main_picks, MAINSHOCK)
plt.show()
5. Aftershock — S−P timing check¶
The moderate aftershock has a short, impulsive source, so both phases should be picked cleanly. Observed S−P is compared against the interval implied by hypocentral distance.
An S is only paired with a P when it falls within a plausible interval of it — otherwise the next S in an active aftershock sequence can easily belong to a different earthquake.
print(f"Fetching {AFTERSHOCK['name']} ...")
after_streams = load_event(fs, AFTERSHOCK, pre=30, post=90)
after_picks = pick_all(model, after_streams)
print(f"\n{len(after_streams)}/{len(STATIONS)} stations loaded")
Fetching M4.6 aftershock ...
5/5 stations loaded
def first_arrivals(sta_picks, event, dist, max_ratio=2.5, pad=3.0):
"""First P after origin, and the first S that could belong to the same event.
An S is only accepted if it lands within a physically plausible interval of
the P. Without that guard the "next S after the P" can easily belong to a
later aftershock, which is a real hazard in this sequence.
"""
hyp = np.hypot(dist, event["depth_km"])
pred_sp = hyp / VS - hyp / VP
ps = [p for p in sta_picks if p.phase == "P" and p.peak_time >= event["time"]]
if not ps:
return None, None, pred_sp
p0 = min(ps, key=lambda p: p.peak_time)
limit = pred_sp * max_ratio + pad
ss = [s for s in sta_picks
if s.phase == "S"
and 0 < (s.peak_time - p0.peak_time) <= limit]
s0 = min(ss, key=lambda s: s.peak_time) if ss else None
return p0, s0, pred_sp
rows = []
for sta, lat, lon in STATIONS:
if sta not in after_picks:
continue
dist = distance_km(AFTERSHOCK, lat, lon)
p0, s0, pred = first_arrivals(after_picks[sta], AFTERSHOCK, dist)
obs = (s0.peak_time - p0.peak_time) if (p0 and s0) else np.nan
rows.append(dict(
station=sta,
dist_km=round(dist, 1),
P=round(p0.peak_time - AFTERSHOCK['time'], 2) if p0 else np.nan,
S=round(s0.peak_time - AFTERSHOCK['time'], 2) if s0 else np.nan,
obs_SP=round(obs, 2) if p0 and s0 else np.nan,
pred_SP=round(pred, 2),
diff=round(obs - pred, 2) if p0 and s0 else np.nan,
))
df = pd.DataFrame(rows)
print(df.to_string(index=False))
paired = df.dropna(subset=['obs_SP'])
print(f"\nStations with a P/S pair : {len(paired)}/{len(df)}")
if len(paired):
print(f"S-P within 2 s of predicted : {int(paired['diff'].abs().lt(2).sum())}/{len(paired)}")
print(f"P precedes S everywhere : {bool((paired['S'] > paired['P']).all())}")
station dist_km P S obs_SP pred_SP diff
CLC 21.8 3.92 6.99 3.07 2.63 0.44
TOW2 31.1 5.97 10.81 4.84 3.72 1.12
SRT 24.2 4.73 8.80 4.07 2.91 1.16
WRC2 37.2 0.39 1.39 1.00 4.44 -3.44
JRC2 47.7 8.43 14.57 6.14 5.69 0.45
Stations with a P/S pair : 5/5
S-P within 2 s of predicted : 4/5
P precedes S everywhere : True
for sta, st in after_streams.items():
lat, lon = next((la, lo) for c, la, lo in STATIONS if c == sta)
plot_station(sta, st, after_picks[sta],
distance_km(AFTERSHOCK, lat, lon), AFTERSHOCK, window=(-5, 30))
plt.show()
How to read the result¶
The setup is healthy when picks sit on visible onsets, P precedes S at every station, observed S−P grows with distance and tracks the prediction to within a second or two, and picks in the record section follow the moveout curves.
Observed S−P running slightly longer than predicted is normal — a single constant-velocity layer is a crude approximation, and the real Vp/Vs in this crust is a little higher than the 1.71 assumed here.
Worth investigating is S before P, S−P that does not scale with distance (often a component-mapping or resampling problem), dense picks in the pre-event noise, or a station that produces nothing while its neighbours at similar distance work.
Two things are expected rather than alarming: missing S for the mainshock at close range, and extra picks throughout — this window sits inside one of the most active aftershock sequences on record, so much of that extra energy is real earthquakes.