QuakeXNet and the picker on Alaska events

QuakeXNet sorts waveforms into earthquake / explosion / noise / surface event. It was trained in the Pacific Northwest. This notebook asks a narrow question: does it still separate those classes in Alaska, and what does the phase picker do on the same windows?

Alaska is a genuine out-of-region test with all three interesting classes occurring naturally and catalogued by the Alaska Earthquake Center:

Catalog type	QuakeXNet class	Where
earthquake	eq	throughout
quarry blast, explosion	px	a working quarry north of Fairbanks
ice quake, landslide	su	Columbia Glacier, Prince William Sound

Two stations carry the test so that source–receiver distances stay comparable across classes: AK.POKR near the Fairbanks quarry and AK.GLI near Columbia Glacier. Each also has local earthquakes within a similar distance range, so eq is not being judged at a different scale from the others.

Waveforms come from EarthScope FDSN, which needs no credentials.

This is a transfer test, not a validation. The catalog type is what an analyst assigned; agreement means the model reproduces that judgement out of region. Counts here are small, so read the confusion matrix as a direction to investigate rather than a performance figure.

Read section 7 before trusting any number here

The single biggest control on the result is not the region — it is where the arrival sits inside the 100 s window. Training only ever placed it 4–20 s from the start (randint(-20, -4) in the training repo's helper_functions.py), and the model learned that convention. Violate it and agreement on this same set of events falls from 78% to 16%, with no change to the data at all.

Any evaluation of this classifier that does not state its window convention is uninterpretable — including the first version of this notebook, which centred the arrival and concluded the model did not transfer.

import sys
from collections import Counter, 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

sys.path.insert(0, '..')          # QuakeXNet lives in the repo, not SeisBench
from sb_catalog.src.classifier import QuakeXNet

%matplotlib inline

1. Configuration

# Stations anchoring the test, with their coordinates
SITES = {
    "POKR": dict(network="AK", lat=65.1180, lon=-147.4307,
                 note="near the Fairbanks quarry"),
    "GLI":  dict(network="AK", lat=61.0000, lon=-146.7000,
                 note="near Columbia Glacier"),
}

# Catalog event type -> expected QuakeXNet label
TYPE_TO_CLASS = {
    "earthquake":   "eq",
    "quarry blast": "px",
    "explosion":    "px",
    "ice quake":    "su",
    "landslide":    "su",
}

SEARCH_START = UTCDateTime("2025-01-01")
SEARCH_END = UTCDateTime("2026-01-01")
MAX_RADIUS_DEG = 0.8              # ~90 km around each station
N_PER_GROUP = 6                   # events per (station, type)

CHANNEL = "BH?"                   # broadband, 50 Hz at these sites
VP = 6.0                          # km/s, for the expected P arrival

# QuakeXNet consumes 5000 samples at 50 Hz = exactly 100 s.
#
# WHERE the arrival sits in that window matters a great deal. Training cut
# windows starting 4-20 s before the P, so the model has only ever seen
# arrivals near the start. Section 7 measures how far that convention can be
# violated before the classifier stops working.
WINDOW_SECONDS = 100.0
WINDOW_LEAD = 15.0     # pre-arrival seconds; inside the training band
LEAD_SWEEP = [5, 10, 15, 20, 30, 40, 50, 60, 75]

PICKER_WEIGHT = "quakescope2026"  # falls back below if not installed
P_THRESHOLD = S_THRESHOLD = 0.3

CLASSES = ["eq", "px", "no", "su"]
CLASS_NAME = {"eq": "earthquake", "px": "explosion",
              "no": "noise", "su": "surface event"}
COLORS = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100"]
C_P, C_S = "#2a78d6", "#eb6834"

2. Build the event list

Events are drawn per station and per catalog type, so each class is represented at a comparable distance range rather than whatever the catalog happens to offer.

usgs = Client("USGS")

events = []
for sta, meta in SITES.items():
    for etype in TYPE_TO_CLASS:
        try:
            cat = usgs.get_events(
                starttime=SEARCH_START, endtime=SEARCH_END,
                latitude=meta['lat'], longitude=meta['lon'],
                maxradius=MAX_RADIUS_DEG, eventtype=etype,
                minmagnitude=1.0, limit=N_PER_GROUP,
            )
        except Exception:
            continue
        for ev in cat:
            o = ev.preferred_origin() or ev.origins[0]
            mag = ev.preferred_magnitude() or ev.magnitudes[0]
            dist = gps2dist_azimuth(meta['lat'], meta['lon'],
                                    o.latitude, o.longitude)[0] / 1000
            events.append(dict(
                station=sta, network=meta['network'],
                etype=etype, expected=TYPE_TO_CLASS[etype],
                time=o.time, mag=round(mag.mag, 1),
                dist_km=round(dist, 1),
                depth_km=round((o.depth or 0) / 1000, 1),
            ))

ev_df = pd.DataFrame(events)
print(ev_df.groupby(['station', 'etype']).size().to_string())
print(f'\ntotal candidate events: {len(ev_df)}')

station  etype       
GLI      earthquake      6
         ice quake       6
         landslide       2
POKR     earthquake      6
         explosion       6
         quarry blast    6

total candidate events: 32

3. Load the models

QuakeXNet comes from this repository rather than the SeisBench model zoo, so its weights must be in the SeisBench cache — the Dockerfile copies them there, and locally you can do the same from sb_catalog/models/quakexnet/.

classifier = QuakeXNet.from_pretrained('base')
classifier.eval()
print(f"QuakeXNet: {classifier.in_samples} samples @ "
      f"{classifier.sampling_rate:g} Hz, classes {classifier.labels}")

available = sbm.PhaseNet.list_pretrained()
weight = PICKER_WEIGHT if PICKER_WEIGHT in available else 'original'
if weight != PICKER_WEIGHT:
    print(f"{PICKER_WEIGHT} not installed - picking with '{weight}'")
picker = sbm.PhaseNet.from_pretrained(weight)
print(f"picker: {weight}")

QuakeXNet: 5000 samples @ 50 Hz, classes ['eq', 'px', 'no', 'su']

picker: quakescope2026

4. Fetch, classify, pick

Each event is fetched generously and then cropped to a window that begins WINDOW_LEAD seconds before the expected P — the convention the model was trained with. The picker runs on the identical crop, and the full fetch is kept so section 7 can re-cut it without downloading again.

fdsn = Client('EARTHSCOPE')

results, waveforms = [], {}
for i, ev in ev_df.iterrows():
    arrival = ev['time'] + ev['dist_km'] / VP
    try:
        # generous fetch; the analysis window is cropped from it below
        st = fdsn.get_waveforms(ev['network'], ev['station'], '*', CHANNEL,
                                arrival - 90, arrival + 110)
    except Exception as exc:
        results.append({**ev, 'predicted': None, 'note': type(exc).__name__})
        continue
    st.merge(fill_value=0)
    if len(st) < 3:
        results.append({**ev, 'predicted': None, 'note': f'{len(st)} traces'})
        continue

    # classification on a crop matching the training convention
    win = st.slice(arrival - WINDOW_LEAD, arrival - WINDOW_LEAD + WINDOW_SECONDS)
    try:
        ann = classifier.annotate(win)
        probs = {c: float(ann.select(channel=f'*{c}')[0].data.max())
                 for c in CLASSES}
        predicted = max(probs, key=probs.get)
    except Exception as exc:
        results.append({**ev, 'predicted': None, 'note': type(exc).__name__})
        continue

    # picking on the same samples
    try:
        out = picker.classify(win, P_threshold=P_THRESHOLD, S_threshold=S_THRESHOLD)
        n_p = sum(1 for p in out.picks if p.phase == 'P')
        n_s = sum(1 for p in out.picks if p.phase == 'S')
        picks = list(out.picks)
    except Exception:
        n_p = n_s = 0
        picks = []

    key = f"{ev['station']}_{i}"
    waveforms[key] = (st, picks, arrival)   # full fetch, for the sweep
    results.append({**ev, 'predicted': predicted,
                    **{f'p_{c}': round(probs[c], 3) for c in CLASSES},
                    'n_P': n_p, 'n_S': n_s, 'key': key, 'note': ''})

res = pd.DataFrame(results)
ok = res[res['predicted'].notna()]
print(f'{len(ok)}/{len(res)} events processed')
if len(ok) < len(res):
    print('\nskipped:')
    print(res[res['predicted'].isna()][['station', 'etype', 'time', 'note']]
          .to_string(index=False))

32/32 events processed

5. Did the classes survive the move to Alaska?

cm = pd.crosstab(ok['expected'], ok['predicted'],
                 rownames=['catalog'], colnames=['QuakeXNet'])
cm = cm.reindex(index=['eq', 'px', 'su'], columns=CLASSES, fill_value=0)
print('confusion matrix (rows: catalog type, columns: predicted)')
print(cm.to_string())

print('\nper-class agreement with the catalog:')
for cls in ['eq', 'px', 'su']:
    sub = ok[ok['expected'] == cls]
    if len(sub) == 0:
        continue
    hit = (sub['predicted'] == cls).sum()
    print(f'  {CLASS_NAME[cls]:<14} {hit:2d}/{len(sub):2d}  ({hit / len(sub):.0%})')
overall = (ok['predicted'] == ok['expected']).mean()
print(f'\noverall agreement: {overall:.0%} on {len(ok)} events')

confusion matrix (rows: catalog type, columns: predicted)
QuakeXNet  eq  px  no  su
catalog                  
eq          9   2   0   1
px          1  11   0   0
su          0   1   2   5

per-class agreement with the catalog:
  earthquake      9/12  (75%)
  explosion      11/12  (92%)
  surface event   5/ 8  (62%)

overall agreement: 78% on 32 events

fig, ax = plt.subplots(figsize=(6.2, 4))
im = ax.imshow(cm.values, cmap='Blues', vmin=0)
ax.set_xticks(range(len(CLASSES)))
ax.set_xticklabels([f'{c}\n{CLASS_NAME[c]}' for c in CLASSES], fontsize=9)
ax.set_yticks(range(len(cm.index)))
ax.set_yticklabels([f'{c} ({CLASS_NAME[c]})' for c in cm.index], fontsize=9)
for r in range(cm.shape[0]):
    for c in range(cm.shape[1]):
        v = cm.values[r, c]
        ax.text(c, r, str(v), ha='center', va='center', fontsize=11,
                color='white' if v > cm.values.max() * 0.6 else '#0b0b0b')
ax.set_xlabel('QuakeXNet prediction')
ax.set_ylabel('Catalog type')
ax.set_title('Class agreement, Alaska', fontsize=11, loc='left')
fig.tight_layout()
plt.show()

Confidence behind each decision

A wrong call made at 0.9 is a different problem from one made at 0.35 — the first says the model is confidently out of domain, the second that it is merely undecided.

show = ok[['station', 'etype', 'mag', 'dist_km', 'expected', 'predicted']
          + [f'p_{c}' for c in CLASSES] + ['n_P', 'n_S']].copy()
show['correct'] = np.where(show['expected'] == show['predicted'], 'yes', '')
print(show.to_string(index=False))

station        etype  mag  dist_km expected predicted  p_eq  p_px  p_no  p_su  n_P  n_S correct
   POKR   earthquake  1.9     88.2       eq        eq 0.943 0.051 0.000 0.005    1    1     yes
   POKR   earthquake  1.6     14.6       eq        px 0.010 0.953 0.000 0.037    1    1        
   POKR   earthquake  1.7     75.3       eq        eq 0.909 0.065 0.002 0.024    0    1     yes
   POKR   earthquake  1.6     74.7       eq        eq 0.952 0.043 0.001 0.004    1    1     yes
   POKR   earthquake  1.8     13.9       eq        px 0.119 0.756 0.002 0.124    1    1        
   POKR   earthquake  1.6     74.1       eq        eq 0.962 0.030 0.002 0.006    0    0     yes
   POKR quarry blast  1.6     15.6       px        px 0.197 0.708 0.001 0.094    1    1     yes
   POKR quarry blast  1.6     16.2       px        px 0.039 0.853 0.001 0.107    1    0     yes
   POKR quarry blast  1.4     15.9       px        px 0.116 0.808 0.001 0.076    1    1     yes
   POKR quarry blast  1.7     16.2       px        eq 0.688 0.171 0.002 0.139    1    0        
   POKR quarry blast  1.3     16.8       px        px 0.021 0.902 0.000 0.076    2    2     yes
   POKR quarry blast  1.5     15.2       px        px 0.014 0.802 0.001 0.184    1    0     yes
   POKR    explosion  1.3     14.9       px        px 0.433 0.550 0.001 0.016    1    0     yes
   POKR    explosion  1.4     16.9       px        px 0.029 0.945 0.000 0.026    1    1     yes
   POKR    explosion  1.5     14.4       px        px 0.046 0.879 0.001 0.075    1    1     yes
   POKR    explosion  1.4     15.8       px        px 0.307 0.399 0.002 0.292    1    1     yes
   POKR    explosion  1.2     15.0       px        px 0.079 0.741 0.002 0.179    1    1     yes
   POKR    explosion  1.2     15.4       px        px 0.023 0.886 0.001 0.091    1    1     yes
    GLI   earthquake  2.2     75.1       eq        eq 0.851 0.101 0.028 0.021    0    2     yes
    GLI   earthquake  1.8     37.8       eq        eq 0.962 0.016 0.001 0.021    1    1     yes
    GLI   earthquake  1.9     73.1       eq        eq 0.880 0.070 0.042 0.008    2    0     yes
    GLI   earthquake  2.0     61.9       eq        eq 0.962 0.030 0.001 0.008    1    1     yes
    GLI   earthquake  1.5     49.1       eq        su 0.002 0.218 0.002 0.778    1    0        
    GLI   earthquake  1.6     41.8       eq        eq 0.865 0.118 0.008 0.009    1    1     yes
    GLI    ice quake  1.9     46.2       su        su 0.003 0.154 0.010 0.832    2    1     yes
    GLI    ice quake  1.7     67.3       su        no 0.029 0.085 0.749 0.137    2    0        
    GLI    ice quake  1.7     63.1       su        su 0.007 0.439 0.004 0.550    0    0     yes
    GLI    ice quake  2.1     23.1       su        no 0.238 0.157 0.519 0.086    1    1        
    GLI    ice quake  1.6     20.9       su        px 0.004 0.649 0.005 0.342    0    0        
    GLI    ice quake  1.6     46.2       su        su 0.002 0.296 0.001 0.701    0    0     yes
    GLI    landslide  1.6     60.4       su        su 0.001 0.008 0.002 0.988    0    0     yes
    GLI    landslide  1.4     53.6       su        su 0.001 0.005 0.003 0.991    0    0     yes

fig, ax = plt.subplots(figsize=(10, 4))
groups = ['eq', 'px', 'su']
for gi, cls in enumerate(groups):
    sub = ok[ok['expected'] == cls]
    if not len(sub):
        continue
    conf = [row[f"p_{row['predicted']}"] for _, row in sub.iterrows()]
    right = sub['predicted'] == cls
    x = np.full(len(sub), gi, dtype=float) + np.linspace(-0.22, 0.22, len(sub))
    ax.scatter(x[right.values], np.array(conf)[right.values],
               marker='o', s=52, color=COLORS[0], label='matches catalog' if gi == 0 else None)
    ax.scatter(x[~right.values], np.array(conf)[~right.values],
               marker='X', s=62, color=COLORS[1], label='differs' if gi == 0 else None)
ax.set_xticks(range(len(groups)))
ax.set_xticklabels([f'{c}\n{CLASS_NAME[c]}' for c in groups])
ax.set_ylim(0, 1.02)
ax.set_ylabel('probability of the predicted class')
ax.set_title('Confidence, split by whether the call matched the catalog',
             fontsize=11, loc='left')
ax.grid(alpha=0.25, lw=0.5, axis='y')
ax.set_axisbelow(True)
ax.legend(frameon=False, fontsize=9)
fig.tight_layout()
plt.show()

6. Waveforms, with the picks and the call

One example per catalog type. The vertical component is shown with the picker's arrivals overlaid, so the classification can be judged against what the trace actually looks like.

def plot_event(row, pad=(-30, 70)):
    st, picks, arrival = waveforms[row['key']]
    tr = st.select(component='Z')
    if not tr:
        return
    tr = tr[0]
    fig, ax = plt.subplots(figsize=(11, 2.9))
    t = tr.times(reftime=row['time'])
    ax.plot(t, tr.data, color='#3d3d3d', lw=0.6)
    for p in picks:
        is_p = p.phase == 'P'
        ax.axvline(p.peak_time - row['time'], color=C_P if is_p else C_S,
                   lw=1.6, ls='-' if is_p else '--', alpha=0.9)
    ax.axvline(0, color='#8a8a8a', lw=1)
    verdict = 'matches catalog' if row['predicted'] == row['expected'] else 'differs'
    ax.set_title(
        f"{row['network']}.{row['station']}  |  catalog: {row['etype']}  |  "
        f"QuakeXNet: {row['predicted']} ({row[f'p_' + row['predicted']]:.2f}) - {verdict}  |  "
        f"M{row['mag']}, {row['dist_km']} km",
        fontsize=10, loc='left')
    ax.set_xlim(*pad)
    ax.set_xlabel('seconds after origin')
    ax.grid(alpha=0.25, lw=0.5)
    handles = [plt.Line2D([], [], color=C_P, lw=1.6, label='P'),
               plt.Line2D([], [], color=C_S, lw=1.6, ls='--', label='S')]
    ax.legend(handles=handles, fontsize=8, frameon=False, loc='upper right', ncol=2)
    fig.tight_layout()
    plt.show()


for cls in ['eq', 'px', 'su']:
    sub = ok[ok['expected'] == cls]
    if len(sub):
        plot_event(sub.iloc[0])

7. What the picker did

Picks on non-earthquake sources are informative in both directions. An explosion has a genuine P and often a weak S, so picks there are correct and would enter a catalog as a locatable event. A surface event has no sharp body-wave arrival, so confident picks on one indicate the picker is responding to something other than an impulsive arrival.

pick_summary = ok.groupby('expected')[['n_P', 'n_S']].agg(['mean', 'max'])
print('picks per 100 s window, by catalog type')
print(pick_summary.round(2).to_string())

print('\nevents with no picks at all:')
none_picked = ok[(ok['n_P'] == 0) & (ok['n_S'] == 0)]
if len(none_picked):
    print(none_picked[['station', 'etype', 'mag', 'dist_km', 'predicted']]
          .to_string(index=False))
else:
    print('  none - every window produced at least one pick')

picks per 100 s window, by catalog type
           n_P       n_S    
          mean max  mean max
expected                    
eq        0.83   2  0.83   2
px        1.08   2  0.75   2
su        0.62   2  0.25   1

events with no picks at all:
station      etype  mag  dist_km predicted
   POKR earthquake  1.6     74.1        eq
    GLI  ice quake  1.7     63.1        su
    GLI  ice quake  1.6     20.9        px
    GLI  ice quake  1.6     46.2        su
    GLI  landslide  1.6     60.4        su
    GLI  landslide  1.4     53.6        su

7. How much does window placement matter?

Every event above is re-classified across a range of lead times, reusing the waveforms already fetched. Nothing changes except where the window is cut.

This is the diagnostic that separates the model does not transfer from the windows were cut wrong — and the two look identical in a confusion matrix.

sweep_rows = []
for _, row in ok.iterrows():
    st, _, arrival = waveforms[row['key']]
    for lead in LEAD_SWEEP:
        w = st.slice(arrival - lead, arrival - lead + WINDOW_SECONDS)
        if len(w) < 3 or w[0].stats.npts < 0.98 * WINDOW_SECONDS * 50:
            continue
        try:
            ann = classifier.annotate(w)
            probs = {c: float(ann.select(channel=f'*_{c}')[0].data.max())
                     for c in CLASSES}
        except Exception:
            continue
        sweep_rows.append(dict(expected=row['expected'], lead=lead,
                               predicted=max(probs, key=probs.get)))

sweep = pd.DataFrame(sweep_rows)
tab = []
for lead in LEAD_SWEEP:
    sub = sweep[sweep.lead == lead]
    if not len(sub):
        continue
    entry = {'lead_s': lead,
             'overall': round((sub.predicted == sub.expected).mean(), 2)}
    for cls in ['eq', 'px', 'su']:
        c = sub[sub.expected == cls]
        entry[cls] = round((c.predicted == cls).mean(), 2) if len(c) else np.nan
    entry['n'] = len(sub)
    tab.append(entry)
sweep_tab = pd.DataFrame(tab)
print(sweep_tab.to_string(index=False))
print('\nlead_s = seconds of data before the arrival; training used 4-20 s')

 lead_s  overall   eq   px   su  n
      5     0.56 0.83 0.25 0.62 32
     10     0.69 0.75 0.67 0.62 32
     15     0.78 0.75 0.92 0.62 32
     20     0.66 0.75 0.75 0.38 32
     30     0.72 0.75 0.83 0.50 32
     40     0.66 0.67 0.75 0.50 32
     50     0.38 0.42 0.17 0.62 32
     60     0.28 0.25 0.08 0.62 32
     75     0.16 0.08 0.00 0.50 32

lead_s = seconds of data before the arrival; training used 4-20 s

fig, ax = plt.subplots(figsize=(9.5, 4.2))
ax.axvspan(4, 20, color='#2a78d6', alpha=0.10)
ax.annotate('training convention\n(4-20 s before P)', xy=(12, 0.06),
            ha='center', fontsize=9, color='#52514e')
for cls, color in zip(['eq', 'px', 'su'], COLORS):
    ax.plot(sweep_tab['lead_s'], sweep_tab[cls], marker='o', lw=2,
            color=color, label=f'{cls} ({CLASS_NAME[cls]})')
ax.plot(sweep_tab['lead_s'], sweep_tab['overall'], marker='s', lw=2.4,
        color='#3d3d3d', ls='--', label='overall')
for _, r in sweep_tab.iterrows():
    ax.annotate(f"{r['overall']:.0%}", (r['lead_s'], r['overall']),
                textcoords='offset points', xytext=(0, 9),
                ha='center', fontsize=8, color='#3d3d3d')
ax.set_xlabel('seconds of data before the arrival')
ax.set_ylabel('agreement with the catalog')
ax.set_ylim(0, 1.08)
ax.set_title('Classifier agreement against window placement',
             fontsize=11, loc='left')
ax.grid(alpha=0.25, lw=0.5)
ax.legend(frameon=False, fontsize=9)
fig.tight_layout()
plt.show()

What this means for the production pipeline

QuakeScope classifies continuous data by sliding a window with a stride of 2500 samples — 50 s at this sampling rate. An event therefore lands at an arbitrary offset inside whichever windows contain it, and the curve above says the answer depends on that offset.

It also explains something that otherwise looks like noise: one earthquake spanning several consecutive windows can come back eq, then px, then su, because each window presents the same arrival at a different position.

The pipeline already knows where the arrivals are — it picked them. Cutting classification windows at a fixed lead before each P, instead of on a blind slide, costs nothing and puts every event in the band the model was trained for.

Reading the result

Does the model separate the classes in Alaska? At a window placement matching training, yes — well above chance, and the confusions that remain are structured rather than random. The learned representation transfers; what is fragile is the interface around it.

Which class is weakest? Surface events — and they are also the least sensitive to placement, both consistent with an emergent signal that offers no sharp onset to key on. Explosions are the opposite: the most placement- sensitive class, and the strongest once the window is cut correctly.

Which direction do the confusions run? Explosions read as earthquakes inflates the earthquake catalog; earthquakes read as explosions suppresses it. These are not interchangeable errors, so the confusion matrix matters more than the headline number.

Two limits worth holding onto: the sample is small — a handful of events per class at two stations — so single misfires move the percentages a lot; and the catalog type is an analyst's judgement rather than ground truth, which for quarry blasts often rests on time of day and a known source location as much as on the waveform.

If the goal is a classifier that works across regions without this much care about window cutting, the levers are at training time: augment over arrival position so the model stops depending on it, and add labelled examples from outside the Pacific Northwest.