Reading the QuakeScope pick catalogue¶
Machine-learning phase picks for the western United States, as Parquet on S3. No credentials, no database, no account. The catalogue is public-read and every cell below runs anonymously, on Colab or a laptop, in about five minutes.
This notebook:
- installs what it needs,
- looks at the station table,
- mirrors one month of one network and reads it,
- plots what is in it,
- fetches the original waveforms from FDSN with ObsPy and draws the picks on top,
- reads the run record that says which model produced them,
- shows how to read at scale without waiting an afternoon.
| bucket | s3://quakescope-picks-2026, us-east-2, public-read |
| catalogue used here | western: PhaseNet original, 1986 to 2026, 24,008 station-locations, about 1.4 billion picks |
| layout | <campaign>/picks/network=<NET>/year=<YYYY>/month=<MM>/*.parquet |
| provenance | <campaign>/runs/<rid>.json |
| reference | docs/data_access.md |
For a region, a period, coverage checks and bulk export, continue with
download_the_catalogue.ipynb.
1. Install¶
Safe to re-run; skip it if you already have these.
# Colab already has pandas/numpy/matplotlib. These are the rest; installed only if missing.
import importlib.util
missing = [m for m in ("pyarrow", "s3fs", "obspy", "duckdb") if importlib.util.find_spec(m) is None]
if missing:
%pip install -q {" ".join(missing)}
print("ready" + (f" (installed {', '.join(missing)})" if missing else ""))
ready
import time
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import s3fs
BUCKET = "quakescope-picks-2026"
CAMPAIGN = "western"
ANON = {"anon": True} # public-read: no keys, no config, no account
fs = s3fs.S3FileSystem(anon=True)
MIRROR = Path("quakescope_mirror") # where partitions land when we copy them down
MIRROR.mkdir(exist_ok=True)
C_P, C_S = "#2a78d6", "#eb6834"
2. The station table¶
stations.parquet is the table the campaign was planned from: every
station-location with its coordinates, the bands the archive lists for it, its
operating window and (for the western campaign) the state it sits in. It is
the index for choosing what to read.
stations = pd.read_parquet(f"s3://{BUCKET}/{CAMPAIGN}/stations.parquet", storage_options=ANON)
print(f"{len(stations):,} station-locations, {stations.network_code.nunique()} networks")
print(stations.state.value_counts(dropna=False).to_string(), "\n")
stations.head()
24,008 station-locations, 119 networks state California 10455 Washington 7718 Wyoming 2798 Nevada 1445 Oregon 1281 Idaho 251 None 60
| id | network_code | station_code | location_code | channels | latitude | longitude | elevation | start_date | end_date | state | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1A.SIA. | 1A | SIA | CN | 44.215005 | -122.254682 | 485.280039 | 2023.136 | 2023.138 | Oregon | |
| 1 | 1B.1. | 1B | 1 | DP | 35.688107 | -120.400028 | 435.400000 | 2018.207 | 2018.304 | California | |
| 2 | 1B.10. | 1B | 10 | DP | 35.689528 | -120.401427 | 424.800000 | 2018.207 | 2018.304 | California | |
| 3 | 1B.11. | 1B | 11 | DP | 35.689002 | -120.402311 | 422.000000 | 2018.207 | 2018.304 | California | |
| 4 | 1B.12. | 1B | 12 | DP | 35.688172 | -120.402680 | 422.000000 | 2018.207 | 2018.304 | California |
fig, ax = plt.subplots(figsize=(7.5, 6.5))
for state, g in stations.groupby(stations.state.fillna("unset")):
ax.scatter(g.longitude, g.latitude, s=3, alpha=0.6, label=f"{state} ({len(g):,})")
ax.set_xlabel("longitude"); ax.set_ylabel("latitude"); ax.set_aspect(1 / np.cos(np.radians(40)))
ax.legend(fontsize=8, markerscale=4, frameon=False)
ax.set_title("western campaign: station-locations by state", loc="left", fontsize=11)
ax.spines[["top", "right"]].set_visible(False)
fig.tight_layout()
3. Read some picks¶
The catalogue is Hive-partitioned by network / year / month, and the fast way to read a partition is to copy it down first, then read it locally: the same month takes about ten seconds to mirror and under a second to read, against a minute and a half read straight from S3 file by file.
Start with the Southern California network CI in July 2019, the month of the
Ridgecrest sequence: 271 files, 149 MB.
NET, YEAR, MONTH = "CI", 2019, 7
part = f"picks/network={NET}/year={YEAR}/month={MONTH:02d}/"
t0 = time.time()
fs.get(f"{BUCKET}/{CAMPAIGN}/{part}", str(MIRROR / CAMPAIGN / part), recursive=True)
n_files = sum(1 for _ in (MIRROR / CAMPAIGN / part).glob("*.parquet"))
print(f"mirrored {n_files} files in {time.time() - t0:.1f} s")
t0 = time.time()
picks = pd.read_parquet(MIRROR / CAMPAIGN / "picks") # partition keys become columns
print(f"read {len(picks):,} picks in {time.time() - t0:.1f} s")
picks.head()
mirrored 273 files in 54.9 s
read 5,528,787 picks in 0.7 s
| tid | cha | pha | start | peak | end | conf | amp | amp_vel | rid | network | year | month | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | CI.SGL. | HH | P | 2019-07-05 00:13:18.198 | 2019-07-05 00:13:18.358 | 2019-07-05 00:13:18.588 | 0.545521 | 0.000045 | 6.442633e-07 | d345adb7-865c-40a6-9c93-894f94952b2c | CI | 2019 | 7 |
| 1 | CI.SGL. | HH | S | 2019-07-05 00:13:19.828 | 2019-07-05 00:13:19.998 | 2019-07-05 00:13:20.228 | 0.824467 | 0.000045 | 6.442633e-07 | d345adb7-865c-40a6-9c93-894f94952b2c | CI | 2019 | 7 |
| 2 | CI.SGL. | HH | P | 2019-07-05 00:18:55.038 | 2019-07-05 00:18:55.258 | 2019-07-05 00:18:55.428 | 0.470023 | NaN | 4.306997e-07 | d345adb7-865c-40a6-9c93-894f94952b2c | CI | 2019 | 7 |
| 3 | CI.SGL. | HH | P | 2019-07-05 00:19:05.888 | 2019-07-05 00:19:06.028 | 2019-07-05 00:19:06.268 | 0.328247 | NaN | 2.477334e-06 | d345adb7-865c-40a6-9c93-894f94952b2c | CI | 2019 | 7 |
| 4 | CI.SGL. | HH | P | 2019-07-05 00:33:08.228 | 2019-07-05 00:33:08.388 | 2019-07-05 00:33:08.858 | 0.361804 | NaN | 3.270497e-06 | d345adb7-865c-40a6-9c93-894f94952b2c | CI | 2019 | 7 |
Equivalent from a shell, if you have the AWS CLI:
aws s3 sync --no-sign-request s3://quakescope-picks-2026/western/picks/network=CI/year=2019/month=07/ \
quakescope_mirror/western/picks/network=CI/year=2019/month=07/
What the columns mean¶
| column | meaning |
|---|---|
tid |
trace id, NET.STA.LOC; the location code may be empty (CI.CLC.) |
cha |
band + instrument code the station-day was picked on (HH, EH, BH, HN, ...), one per station-day |
pha |
P or S |
start, peak, end |
the probability window; peak is the arrival time, millisecond resolution |
conf |
model confidence, 0 to 1. Everything at or above 0.2 is stored |
amp |
Wood-Anderson displacement, metres, for local magnitude. Only measured for conf >= 0.5, NaN otherwise |
amp_vel |
raw peak amplitude in counts, no response removed: a detection-strength proxy, not a physical unit |
rid |
run id, resolves to runs/<rid>.json for the model and thresholds |
network, year, month |
the partition keys, read from the path |
conf is a detection score, not a probability of correctness. The stored
floor of 0.2 is deliberately permissive so you can choose your own threshold;
most analyses want considerably more than that.
print(picks["pha"].value_counts().to_string(), "\n")
print(picks["conf"].describe().round(3).to_string(), "\n")
print(f"stations : {picks['tid'].nunique():,}")
print(f"bands : {picks['cha'].value_counts().to_dict()}")
print(f"amp measured : {picks['amp'].notna().mean():.1%} of picks, "
f"{picks.loc[picks.conf >= 0.5, 'amp'].notna().mean():.1%} of those with conf >= 0.5")
print(f"amp_vel set : {picks['amp_vel'].notna().mean():.1%}")
pha S 2882228 P 2646559
count 5528787.000 mean 0.530 std 0.253 min 0.200 25% 0.298 50% 0.471 75% 0.761 max 0.991
stations : 305
bands : {'HH': 5405025, 'HN': 89830, 'EH': 33932}
amp measured : 46.4% of picks, 99.1% of those with conf >= 0.5
amp_vel set : 99.3%
4. What the month looks like¶
import matplotlib.dates as mdates
fig, ax = plt.subplots(1, 3, figsize=(14, 3.6))
daily = picks.groupby([picks["peak"].dt.floor("D"), "pha"]).size().unstack(fill_value=0)
for ph, c in (("P", C_P), ("S", C_S)):
ax[0].plot(daily.index, daily[ph], lw=1.2, color=c, label=ph)
ax[0].legend(title="phase", frameon=False)
ax[0].set_title("picks per day", loc="left"); ax[0].set_ylabel("picks")
ax[0].xaxis.set_major_formatter(mdates.DateFormatter("%b %d"))
ax[0].xaxis.set_major_locator(mdates.DayLocator(interval=7))
for ph, c in (("P", C_P), ("S", C_S)):
ax[1].hist(picks.loc[picks.pha == ph, "conf"], bins=60, histtype="step", lw=1.3, color=c, label=ph)
ax[1].set_yscale("log"); ax[1].set_title("confidence", loc="left"); ax[1].legend(frameon=False)
amp = picks["amp"].dropna()
amp = amp[amp > 0]
ax[2].hist(np.log10(amp), bins=60, color="#555")
ax[2].set_title("log10 Wood-Anderson amplitude (m), conf >= 0.5", loc="left")
for a in ax:
a.spines[["top", "right"]].set_visible(False)
fig.tight_layout()
The daily curve is the useful sanity check: P and S should track each other, and a day where one collapses is usually a data outage rather than a quiet Earth. The step on July 4 and 6 is the Ridgecrest foreshock and mainshock.
5. Fetch the waveforms and check the picks¶
Picks are only as good as they look on the record. The catalogue stores no waveforms, so fetch them from the SCEDC FDSN service with ObsPy: the same data the picker read, from the archive of record.
Take the busiest station-minute of confident picks on the mainshock day, 2019-07-06: an aftershock cluster minutes after the M7.1.
import obspy
from obspy.clients.fdsn import Client
day = picks[(picks.conf >= 0.5) & (picks.peak.dt.floor("D") == pd.Timestamp("2019-07-06"))].copy()
day["minute"] = day["peak"].dt.floor("min")
tid, minute = day.groupby(["tid", "minute"]).size().idxmax()
net, sta, loc = tid.split(".")
band = picks.loc[picks.tid == tid, "cha"].iloc[0]
print(f"{tid} {minute} band {band}")
t0 = obspy.UTCDateTime(minute) - 10
st = Client("SCEDC", timeout=60).get_waveforms(
network=net, station=sta, location=loc or "--", channel=f"{band}?",
starttime=t0, endtime=t0 + 80,
)
st.merge(fill_value="interpolate").detrend("demean").filter("bandpass", freqmin=2, freqmax=20)
print(st)
CI.CCC. 2019-07-06 16:43:00 band HH
3 Trace(s) in Stream: CI.CCC..HHE | 2019-07-06T16:42:49.998300Z - 2019-07-06T16:44:09.998300Z | 100.0 Hz, 8001 samples CI.CCC..HHN | 2019-07-06T16:42:49.998300Z - 2019-07-06T16:44:09.998300Z | 100.0 Hz, 8001 samples CI.CCC..HHZ | 2019-07-06T16:42:49.998300Z - 2019-07-06T16:44:09.998300Z | 100.0 Hz, 8001 samples
window = picks[(picks.tid == tid) &
(picks.peak >= minute - pd.Timedelta(seconds=10)) &
(picks.peak < minute + pd.Timedelta(seconds=70))]
fig, axes = plt.subplots(len(st), 1, figsize=(12, 6), sharex=True)
for ax, tr in zip(np.atleast_1d(axes), st):
t = tr.times() + (tr.stats.starttime - t0)
ax.plot(t, tr.data / np.abs(tr.data).max(), lw=0.45, color="#333")
for _, p in window.iterrows():
ax.axvline(obspy.UTCDateTime(p["peak"]) - t0, lw=1.3, alpha=0.85,
color=C_P if p["pha"] == "P" else C_S)
ax.set_ylabel(tr.stats.channel)
ax.spines[["top", "right"]].set_visible(False)
axes[-1].set_xlabel(f"seconds after {t0.isoformat()}")
axes[0].set_title(f"{tid}: P in blue, S in orange ({len(window)} picks)", loc="left")
fig.tight_layout()
Zoom in on a single pick to judge it properly. Change i to step through them;
low-conf picks are where the interesting disagreements live.
i = 0
p = window.sort_values("peak").iloc[i]
tp = obspy.UTCDateTime(p["peak"])
tr = st.select(component="Z")[0].slice(tp - 3, tp + 7)
plt.figure(figsize=(11, 2.8))
plt.plot(tr.times() - 3, tr.data, lw=0.7, color="#333")
plt.axvline(0, color=C_P if p["pha"] == "P" else C_S, lw=1.5)
plt.title(f"{tid} {p['pha']} conf={p['conf']:.3f} {p['peak']}", loc="left")
plt.xlabel("seconds from pick"); plt.gca().spines[["top", "right"]].set_visible(False)
plt.tight_layout()
6. Provenance¶
Every pick carries a run id. The run record says exactly which model,
weights, thresholds and library version produced it. A campaign runs under one
configuration, but every worker start creates a new run id, so there are
thousands of records that differ only in run_id and created: quote the
configuration, not the id.
import json, urllib.request
rid = picks["rid"].iloc[0]
url = f"https://{BUCKET}.s3.us-east-2.amazonaws.com/{CAMPAIGN}/runs/{rid}.json"
run = json.load(urllib.request.urlopen(url))
print(json.dumps(run, indent=2))
# and the whole month ran under one configuration:
print(f"\n{picks.rid.nunique()} run ids in this month")
{
"run_id": "d345adb7-865c-40a6-9c93-894f94952b2c",
"created": "2026-09-04T01:40:26.051924+00:00",
"model": "PhaseNet",
"weight": "original",
"p_threshold": "0.2",
"s_threshold": "0.2",
"components_loaded": "ZNE12",
"seisbench_version": "0.12.5",
"weight_version": "2"
}
60 run ids in this month
7. Reading it at scale¶
One month of one network is 5.5 M picks and reads in under a second once
mirrored. The catalogue will not. western is about 1.4 billion picks in
417,000 Parquet objects (42 GB). At that size the thing that hurts is not
bytes, it is files: reading a Parquet dataset means opening every file's
footer, and that cost is paid per object however small your query.
Measured on 2026-09-16, anonymous, from a laptop:
| operation | objects touched | time |
|---|---|---|
| mirror one month partition (271 files, 149 MB) | 271 | 9 to 12 s |
| read the mirrored month with pandas | 271 | 0.8 s |
read the same month straight from S3 with pandas + filters= |
271 | 86 s |
DuckDB count(*), partition in the path |
271 | 7 s |
DuckDB count(*), picks/**/*.parquet with WHERE network= AND year= AND month= |
276,829 (417,000 today) | not back after 7 min |
Three rules, in the order they matter.
1. Put the partition in the path¶
network, year and month come from the path. A glob over picks/**
lists every object in the catalogue before any WHERE clause can prune it.
Naming the partition in the path lists 271 objects instead of 417,000.
import duckdb
con = duckdb.connect()
con.sql("INSTALL httpfs; LOAD httpfs; SET s3_region='us-east-2';")
con.sql("SET s3_access_key_id=''; SET s3_secret_access_key='';") # anonymous
PART = f"s3://{BUCKET}/{CAMPAIGN}/picks/network={NET}/year={YEAR}/month={MONTH:02d}/*.parquet"
t0 = time.time()
n = con.sql(f"SELECT count(*) FROM read_parquet('{PART}', hive_partitioning=1)").fetchone()[0]
print(f"{n:,} picks in {time.time() - t0:.1f} s, straight from S3, partition in the path")
5,528,787 picks in 11.1 s, straight from S3, partition in the path
2. For anything bigger than a month, mirror first¶
aws s3 sync --no-sign-request and s3fs.get(recursive=True) both fetch
concurrently, at 15 to 20 MB/s from here. A network-year is about 0.7 GB; the
whole of western is 42 GB and takes about 90 minutes on a fast link. Every
query after that is local, and a glob is fine on disk.
t0 = time.time()
busiest = con.sql(f'''
SELECT tid,
count(*) AS picks,
count(*) FILTER (pha = 'S') AS s_picks,
round(avg(conf), 3) AS mean_conf
FROM read_parquet('{MIRROR / CAMPAIGN}/picks/**/*.parquet', hive_partitioning=1)
GROUP BY tid
ORDER BY picks DESC
LIMIT 10
''').df()
print(f"{time.time() - t0:.2f} s against the local mirror")
busiest
0.08 s against the local mirror
| tid | picks | s_picks | mean_conf | |
|---|---|---|---|---|
| 0 | CI.WMF. | 243137 | 123669 | 0.673 |
| 1 | CI.WRC2. | 225225 | 115686 | 0.666 |
| 2 | CI.CLC. | 224689 | 113076 | 0.718 |
| 3 | CI.JRC2. | 204043 | 104194 | 0.699 |
| 4 | CI.MPM. | 188881 | 97789 | 0.699 |
| 5 | CI.TOW2. | 180075 | 90572 | 0.658 |
| 6 | CI.SRT. | 138907 | 66672 | 0.675 |
| 7 | CI.CCC. | 120306 | 60693 | 0.666 |
| 8 | CI.WHF. | 112802 | 60510 | 0.660 |
| 9 | CI.WBM. | 85723 | 45836 | 0.687 |
3. Aggregate in the engine, select only the columns you need¶
.df() materialises everything into memory. A count, a histogram or a
per-station summary should come back already reduced. Parquet is columnar, so
leaving amp and amp_vel out of a query that wants arrival times halves the
bytes read; SELECT * reads them whether you use them or not.
And the one that is not yours to fix: file size¶
A shard writes its own Parquet objects, so a partition accumulates hundreds of files of about 120 KB. Compaction into ~1 MB files is planned and will change object names, not content. Until then the rules above are what make the catalogue fast.
8. Going further¶
A region and a period, not a network-month. The download notebook selects stations from the station table, mirrors the matching partitions, checks which station-days were actually processed, and exports for an associator.
Other catalogues. obs is the ocean-bottom stations on the obs
(PickBlue) weight, 1993 to 2026; global is every network EarthScope serves
on jma_wc, 3.6 % done. Both share this layout and schema. Progress: the
campaign dashboard.
Choosing a threshold. Everything above conf 0.2 is stored, which is far
more permissive than most uses want. Raise it and check the effect on your own
target events rather than adopting a number from elsewhere; the benchmark
reports show that the same threshold means different things for different
weight sets.
Associating into events. These are per-station picks, not locations. A phase associator (PyOcto, GaMMA) turns them into events.
Amplitudes. amp is Wood-Anderson displacement in metres, ready for local
magnitude, and exists only for picks with conf >= 0.5. amp_vel is
pre-response counts, useful mainly for QC.
Catalogue and code: SeisSCOPED/QuakeScope.
Reference for the layout, the access rules and the known limits:
docs/data_access.md.