Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
A direct command skips the review prompt. Inspect the source before running it.
A comprehensive skill for seismological data acquisition, processing, and analysis
using the ObsPy framework. Covers FDSN waveform retrieval, instrument response
removal, P/S phase arrival picking, basic moment tensor inversion, and
seismicity visualisation.
When to Use This Skill
Use this skill when you need to:
Download seismogram waveforms from IRIS, GEOFON, ORFEUS, or any FDSN-compliant
data centre
Remove instrument response and convert raw counts to physical units (m/s, m/s², Pa)
Filter, decimate, and taper seismic traces before analysis
Automatically pick P and S phase arrivals with STA/LTA or kurtosis detectors
Estimate earthquake source parameters (origin time, location, focal mechanism)
Build seismicity catalogues and map earthquake distributions
Compute and visualise spectrograms and particle motion diagrams
This skill is not appropriate for:
Real-time continuous acquisition from hardware digitisers (use SeisComP or
Earthworm instead)
Full waveform tomography (use SPECFEM or SALVUS skill)
Distributed seismic arrays with >100 stations in a single session
Background & Key Concepts
ObsPy Data Model
ObsPy organises seismic data in three nested containers:
Class
Contains
Analogy
Stream
list of Trace objects
a multi-channel recording session
Trace
1-D NumPy array + Stats
a single channel
Stats
dict-like metadata
network, station, channel, start time, sampling rate
The SEED channel code (e.g., BHZ) encodes band (B = broad-band),
instrument (H = high-gain seismometer), and orientation (Z = vertical).
FDSN Web Services
The International Federation of Digital Seismograph Networks (FDSN) standardises
three web service endpoints:
dataselect — returns MiniSEED waveforms
station — returns StationXML inventory (instrument response)
event — returns QuakeML earthquake catalogues
ObsPy's Client class wraps all three. Major nodes: , ,
, , .
"IRIS"
"GEOFON"
"ORFEUS"
"ETH"
"NCEDC"
Instrument Response
Raw seismometer output is in digital counts. The instrument response (poles,
zeros, sensitivity) converts counts to ground motion. Removing it via spectral
division yields velocity [m/s], displacement [m], or acceleration [m/s²]
records. ObsPy reads response from StationXML and applies it with
Trace.remove_response().
Phase Picking
Seismic phase picking identifies the onset time of P (compressional) and S (shear)
waves. Classical approaches use the STA/LTA (short-term average / long-term
average) ratio: a sudden energy increase produces a ratio spike. The obspy.signal
module provides classic_sta_lta and recursive_sta_lta.
Moment Tensor
A seismic moment tensor is a 3×3 symmetric matrix describing the equivalent
force system of an earthquake. The scalar seismic moment M_0 and moment
magnitude M_w = (2/3)(log10(M_0) - 9.1) are derived from it. Full waveform
inversion (e.g., time-domain L2 misfit minimisation) fits synthetic seismograms
to observed data to recover the tensor.
Seismicity Maps
Earthquake catalogues are typically distributed as QuakeML or CSV files with
origin time, latitude, longitude, depth, and magnitude. Matplotlib with a
Cartopy or Basemap projection renders these as geographic scatter plots, colour-
coded by depth and scaled by magnitude.
FDSN data centres are open for most uses. If you access restricted data
(embargoed networks), store credentials securely:
export FDSN_USER="<your-username>"export FDSN_PASSWORD=$(cat ~/.fdsn_passwd) # read from file, never hardcode
Access in Python:
import os
from obspy.clients.fdsn import Client
user = os.getenv("FDSN_USER", "")
password = os.getenv("FDSN_PASSWORD", "")
if user:
client = Client("IRIS", user=user, password=password)
else:
client = Client("IRIS") # anonymous for open data
Core Workflow
Step 1 — Download waveforms via FDSN
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
# Connect to IRIS FDSN data centre
client = Client("IRIS")
# Define a 5-minute window around the 2011 Tohoku earthquake (Mw 9.0)
origin_time = UTCDateTime("2011-03-11T05:46:24")
t_start = origin_time - 60# 1 min before origin
t_end = origin_time + 300 - 60# 4 min after# Download broadband vertical (BHZ) from station IU.MAJO (Japan)
st = client.get_waveforms(
network="IU", station="MAJO", location="00", channel="BHZ",
starttime=t_start, endtime=t_end
)
print(st) # Stream summaryprint(st[0].stats) # Trace metadataprint(f"Sampling rate: {st[0].stats.sampling_rate} Hz")
print(f"Duration : {st[0].stats.npts / st[0].stats.sampling_rate:.1f} s")
# Save to MiniSEED for offline use
st.write("tohoku_IU_MAJO_BHZ.mseed", format="MSEED")
print("Saved tohoku_IU_MAJO_BHZ.mseed")
Step 2 — Retrieve StationXML and remove instrument response
import numpy as np
import matplotlib.pyplot as plt
from obspy.taup import TauPyModel
model = TauPyModel(model="ak135")
# Compute multiple phase arrivals at 60 degrees for a 100 km deep source
arrivals = model.get_travel_times(
source_depth_in_km=100,
distance_in_degree=60,
phase_list=["P", "pP", "PP", "S", "SS", "SKS", "ScS", "PKP"]
)
print(f"{'Phase':10s}{'Time (s)':>10s}{'Ray param':>10s}{'Purist':>15s}")
print("-" * 50)
for arr in arrivals:
print(f"{arr.name:10s}{arr.time:10.2f}{arr.ray_param:10.4f} "f"{arr.purist_name:>15s}")
# Ray path plot
arrivals_plot = model.get_ray_paths(
source_depth_in_km=100,
distance_in_degree=60,
phase_list=["P", "S", "ScS", "PKP"]
)
ax = arrivals_plot.plot_rays(plot_type="spherical", show=False,
legend=True, phase_list=["P", "S", "ScS", "PKP"])
ax.figure.savefig("ray_paths.png", dpi=150)
print("Saved ray_paths.png")
Waveform cross-correlation for relative arrival times
import numpy as np
from scipy.signal import correlate, correlation_lags
from obspy import read, Stream
defcross_correlate_picks(st_ref, st_cmp, freqmin=1.0, freqmax=10.0,
window_s=2.0, pick_sample=None):
"""Return sub-sample delay (seconds) between two aligned traces."""for tr in [st_ref[0], st_cmp[0]]:
tr.detrend("linear")
tr.taper(max_percentage=0.05)
tr.filter("bandpass", freqmin=freqmin, freqmax=freqmax,
corners=4, zerophase=True)
df = st_ref[0].stats.sampling_rate
if pick_sample isNone:
pick_sample = len(st_ref[0].data) // 2
hw = int(window_s * df / 2)
a = st_ref[0].data[pick_sample - hw : pick_sample + hw]
b = st_cmp[0].data[pick_sample - hw : pick_sample + hw]
cc = correlate(a, b, mode="full")
lags = correlation_lags(len(a), len(b), mode="full")
lag_s = lags[np.argmax(cc)] / df
cc_max = cc.max() / (np.std(a) * np.std(b) * len(a))
return lag_s, cc_max
# Demo: shift a trace by 0.3 s and recover the delayfrom obspy import Trace
import numpy as np
rng = np.random.default_rng(42)
df = 100.0
t = np.arange(0, 30, 1.0 / df)
sig = np.sin(2 * np.pi * 3.0 * t) * np.exp(-0.05 * t) + 0.1 * rng.standard_normal(len(t))
tr_ref = Trace(data=sig.copy())
tr_ref.stats.sampling_rate = df
true_delay = 0.3# seconds
shift_samples = int(true_delay * df)
tr_cmp = Trace(data=np.roll(sig, shift_samples))
tr_cmp.stats.sampling_rate = df
st_ref = Stream([tr_ref])
st_cmp = Stream([tr_cmp])
delay, cc = cross_correlate_picks(st_ref, st_cmp, freqmin=1.0, freqmax=10.0)
print(f"True delay : {true_delay:.3f} s")
print(f"Measured delay: {delay:.3f} s (CC max = {cc:.4f})")
Troubleshooting
FDSNNoDataException — no data available
# The time window or channel code may be wrong, or data was not archived.# Verify by listing available channels first:
inv = client.get_stations(
network="IU", station="MAJO", level="channel",
starttime=UTCDateTime("2011-03-11"), endtime=UTCDateTime("2011-03-12")
)
for net in inv:
for sta in net:
for cha in sta:
print(cha.code, cha.location_code, cha.start_date, cha.end_date)
Ensure the waveform and inventory cover the same time window.
Use pre_filt to suppress low-frequency integration drift, for example
pre_filt = (0.004, 0.008, 45.0, 50.0) for broadband data.
Reduce water_level from 60 to 30 dB only if the signal-to-noise ratio
is very high.
STA/LTA picks on many false triggers (noise)
# Increase the on-threshold, or apply a bandpass filter first:
st.filter("bandpass", freqmin=1.0, freqmax=20.0, corners=4, zerophase=True)
# then re-run STA/LTA
read() raises TypeError: Not a valid MiniSEED file
Some archives deliver data in SAC, GSE2, or SEISAN format. ObsPy auto-detects
most formats:
st = read("waveform.sac") # SAC
st = read("waveform.gse") # GSE2
For MSEED files with quality flags, add headonly=False, check_compression=True.
Memory error when reading a long continuous stream
Process in chunks using starttime / endtime slices:
chunk_size = 3600# 1 hour per chunk
t0 = UTCDateTime("2020-01-01")
for i inrange(24):
t_start = t0 + i * chunk_size
t_end = t_start + chunk_size
st_chunk = client.get_waveforms("IU", "ANMO", "00", "BHZ", t_start, t_end)
# ... process st_chunk ...
