Skip to main content Skills Marketplace Discover and explore AI skills built by the community.
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.
Show prompt details A direct command skips the review prompt. Inspect the source before running it.
npx skills add https://github.com/xjtulyc/awesome-rosetta-skills --skill obspy-seismologyThe command stays on one line. Scroll horizontally to inspect it before copying.
Prefer a local copy? Download the files currently available to SkillsMP.
Download Zip Downloading... name obspy-seismology description Use this Skill for seismological analysis with ObsPy: waveform download, filtering, P/S phase picking, moment magnitude, and spectral analysis.
tags ["earth-science","seismology","obspy","waveform","earthquake"] version 1.0.0 authors [{"name":"Rosetta Skills Contributors","github":"@xjtulyc"}] license MIT platforms ["claude-code","codex","gemini-cli","cursor"] dependencies {"python":["obspy>=1.4","numpy>=1.24","scipy>=1.11","matplotlib>=3.7"]} last_updated 2026-03-17 status stable
Seismological Analysis with ObsPy
One-line summary : Download, process, and analyze seismic waveforms with ObsPy: bandpass filtering, P/S phase picking, spectral analysis, focal mechanisms, and moment magnitude estimation.
When to Use This Skill
When downloading seismic waveforms from IRIS/FDSN data centers
When filtering and deconvolving instrument response from seismograms
When picking P and S wave arrivals automatically or manually
When computing seismic spectra and corner frequencies
When estimating moment magnitude from seismic records
When computing traveltimes using 1D Earth models (iasp91, PREM)
Trigger keywords : ObsPy, seismology, seismogram, waveform, earthquake, P-wave, S-wave, FDSN, IRIS, bandpass filter, instrument response, moment magnitude, seismic, travel time
Background & Key Concepts
Seismic Wave Types
P-waves : Compressional (primary) waves, fastest, arrive firstS-waves : Shear waves, ~60% of P velocity, arrive secondSurface waves : Rayleigh and Love waves, slowest, largest amplitude
Instrument Response Removal
Raw data in counts → ground motion (m, m/s, or m/s²) by deconvolving the instrument response:
$$
U(\omega) = \frac{X(\omega)}{I(\omega)}
$$
where $X(\omega)$ is the raw spectrum and $I(\omega)$ is the instrument response (poles/zeros + sensitivity).
Moment Magnitude
$$
M_w = \frac{2}{3}\log_{10}(M_0) - 10.7
$$
where $M_0$ is the seismic moment (N·m), estimated from the plateau of the displacement spectrum.
Environment Setup
Install Dependencies
pip install obspy>=1.4 numpy>=1.24 scipy>=1.11 matplotlib>=3.7
Verify Installation
import obspy
from obspy import read_events
print (f"ObsPy version: {obspy.__version__} " )
from obspy import Trace, Stream
import numpy as np
t = np.linspace(0
10
1000
2
5
3
100.0
"XX"
"TEST"
print
f"Test trace: {tr} "
Core Workflow
Step 1: Download Waveforms and Remove Instrument Response from obspy import UTCDateTime
from obspy.clients.fdsn import Client
import matplotlib.pyplot as plt
import numpy as np
client = Client("IRIS" )
event_time = UTCDateTime("2011-03-11T05:46:24" )
origin_lat, origin_lon = 38.297 , 142.373
network, station = "IU" , "MAJO"
starttime = event_time - 60
endtime = event_time + 600
print (f"Downloading waveforms for {network} .{station} ..." )
try :
st = client.get_waveforms(
network=network, station=station, location="00" ,
channel="BH?" ,
starttime=starttime, endtime=endtime,
)
print (f"Downloaded: {st} " )
inv = client.get_stations(
network=network, station=station,
starttime=starttime, endtime=endtime,
level="response" ,
)
print (f"Inventory: {inv} " )
except Exception as e:
print (f"Download failed: {e} " )
print ("Creating synthetic data for demonstration..." )
from obspy import Trace, Stream, Inventory
import numpy as np
dt = 0.01 ; N = 66000
t = np.arange(N) * dt
p_wave = np.zeros(N); s_wave = np.zeros(N)
p_idx = int (100 /dt)
s_idx = int (180 /dt)
p_wave[p_idx:p_idx+500 ] = np.sin(2 *np.pi*1.5 *t[:500 ]) * np.exp(-t[:500 ]/10 )
s_wave[s_idx:s_idx+800 ] = np.sin(2 *np.pi*0.8 *t[:800 ]) * np.exp(-t[:800 ]/15 )
noise = np.random.randn(N) * 0.02
st = Stream()
for ch, amp in [("BHZ" , 1.0 ), ("BHN" , 0.7 ), ("BHE" , 0.5 )]:
tr = Trace(data=(p_wave + s_wave + noise) * amp)
tr.stats.network = network; tr.stats.station = station
tr.stats.channel = ch; tr.stats.sampling_rate = 100.0
tr.stats.starttime = starttime
st.append(tr)
inv = None
st_raw = st.copy()
st.detrend('demean' )
st.detrend('linear' )
st.taper(max_percentage=0.05 , type ='cosine' )
if inv is not None :
st.remove_response(inventory=inv, output='VEL' ,
pre_filt=(0.005 , 0.01 , 40 , 45 ))
y_unit = "Velocity (m/s)"
else :
st.filter ('bandpass' , freqmin=0.01 , freqmax=10.0 )
y_unit = "Counts (filtered)"
fig, axes = plt.subplots(3 , 1 , figsize=(14 , 8 ), sharex=True )
t_axis = np.arange(len (st[0 ].data)) * st[0 ].stats.delta
t0 = float (event_time - starttime)
for ax, tr in zip (axes, st):
ax.plot(t_axis, tr.data, 'k-' , linewidth=0.6 )
ax.axvline(t0, color='red' , linestyle='--' , linewidth=1.5 , label='Origin time' )
ax.set_ylabel(f"{tr.stats.channel} \n{y_unit} " , fontsize=9 )
ax.grid(True , alpha=0.3 )
ax.legend(loc='upper right' , fontsize=8 )
axes[-1 ].set_xlabel("Time after start (s)" )
plt.suptitle(f"Seismogram: {network} .{station} | 2011 Tohoku M9.0" )
plt.tight_layout()
plt.savefig("seismogram.png" , dpi=150 )
plt.show()
Step 2: Phase Picking and Travel Time Calculation from obspy import UTCDateTime, Stream
from obspy.taup import TauPyModel
import numpy as np
import matplotlib.pyplot as plt
model = TauPyModel(model="iasp91" )
source_depth_km = 25.0
dist_deg = 15.0
arrivals = model.get_travel_times(
source_depth_in_km=source_depth_km,
distance_in_degree=dist_deg,
phase_list=["P" , "pP" , "PP" , "S" , "SS" , "Rayleigh" , "Love" ],
)
print (f"Theoretical travel times (Δ={dist_deg} °, h={source_depth_km} km):" )
print (f"{'Phase' :12s} {'Time (s)' :12s} {'Slowness' :12s} " )
print ("-" * 40 )
for arr in arrivals:
print (f"{arr.name:12s} {arr.time:12.2 f} {arr.ray_param_sec_degree:12.4 f} " )
from obspy.signal.trigger import classic_sta_lta, trigger_onset
np.random.seed(42 )
sps = 100.0
N = int (300 * sps)
t_arr = np.arange(N) / sps
data = np.random.randn(N) * 0.5
p_onset = int (100 * sps)
data[p_onset:p_onset+int (50 *sps)] += (
np.sin(2 *np.pi*3 *t_arr[:int (50 *sps)]) * np.exp(-t_arr[:int (50 *sps)]/8 ) * 10
)
s_onset = int (180 * sps)
data[s_onset:s_onset+int (80 *sps)] += (
np.sin(2 *np.pi*1.5 *t_arr[:int (80 *sps)]) * np.exp(-t_arr[:int (80 *sps)]/15 ) * 15
)
cft = classic_sta_lta(data, int (0.5 *sps), int (10 *sps))
on_off = trigger_onset(cft, thres1=3.5 , thres2=0.5 )
print (f"\nSTA/LTA detected {len (on_off)} trigger(s):" )
for on, off in on_off:
print (f" ON at t={on/sps:.2 f} s, OFF at t={off/sps:.2 f} s" )
fig, axes = plt.subplots(2 , 1 , figsize=(14 , 6 ), sharex=True )
axes[0 ].plot(t_arr, data, 'k-' , linewidth=0.5 )
for on, _ in on_off:
axes[0 ].axvline(on/sps, color='red' , linewidth=1.5 , linestyle='--' )
axes[0 ].set_ylabel("Amplitude" ); axes[0 ].set_title("Seismogram + Phase Picks" )
axes[0 ].grid(True , alpha=0.3 )
axes[1 ].plot(t_arr, cft, 'b-' , linewidth=1 )
axes[1 ].axhline(3.5 , color='red' , linestyle='--' , label='ON threshold' )
axes[1 ].axhline(0.5 , color='orange' , linestyle='--' , label='OFF threshold' )
axes[1 ].set_ylabel("STA/LTA ratio" ); axes[1 ].set_xlabel("Time (s)" )
axes[1 ].set_title("Characteristic Function (STA/LTA)" )
axes[1 ].legend(); axes[1 ].grid(True , alpha=0.3 )
plt.tight_layout()
plt.savefig("phase_picking.png" , dpi=150 )
plt.show()
Step 3: Spectral Analysis and Magnitude Estimation import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch
from scipy.optimize import curve_fit
np.random.seed(0 )
sps = 100.0
N = int (60 * sps)
t = np.arange(N) / sps
t_r = 2.0
signal = np.zeros(N)
onset = int (5 * sps)
pulse_t = t[:N-onset]
signal[onset:] = pulse_t * np.exp(-pulse_t / t_r)
signal += np.random.randn(N) * 1e-8
freq = np.fft.rfftfreq(N, d=1 /sps)
spec = np.abs (np.fft.rfft(signal)) / sps
freq_fit = freq[(freq > 0.05 ) & (freq < 20 )]
spec_fit = spec[(freq > 0.05 ) & (freq < 20 )]
def brune_model (f, omega0, fc ):
return omega0 / (1 + (f/fc)**2 )
try :
popt, _ = curve_fit(brune_model, freq_fit, spec_fit,
p0=[spec_fit.max (), 1.0 ], maxfev=5000 )
omega0_hat, fc_hat = popt
print (f"Brune fit: Ω₀ = {omega0_hat:.4 e} , fc = {fc_hat:.3 f} Hz" )
rho = 2700.0
v_p = 6000.0
R = 100e3
F_p = 0.52
M0 = 4 * np.pi * rho * v_p**3 * R * omega0_hat / (2 * F_p)
Mw = (2 /3 ) * np.log10(M0) - 10.7
print (f"Seismic moment M₀ = {M0:.3 e} N·m" )
print (f"Moment magnitude Mw = {Mw:.2 f} " )
except Exception as e:
print (f"Fit failed: {e} " )
omega0_hat, fc_hat = spec_fit.max (), 1.0
fig, axes = plt.subplots(1 , 2 , figsize=(13 , 5 ))
axes[0 ].plot(t, signal * 1e6 , 'k-' , linewidth=0.8 )
axes[0 ].set_xlabel("Time (s)" ); axes[0 ].set_ylabel("Displacement (μm)" )
axes[0 ].set_title("Ground Displacement Seismogram" ); axes[0 ].grid(True , alpha=0.3 )
f_smooth = np.logspace(np.log10(0.05 ), np.log10(20 ), 200 )
axes[1 ].loglog(freq_fit, spec_fit, 'b.' , markersize=3 , alpha=0.5 , label='Data' )
axes[1 ].loglog(f_smooth, brune_model(f_smooth, omega0_hat, fc_hat),
'r-' , linewidth=2.5 , label=f'Brune fit (fc={fc_hat:.2 f} Hz)' )
axes[1 ].axvline(fc_hat, color='gray' , linestyle='--' , linewidth=1 , alpha=0.7 )
axes[1 ].set_xlabel("Frequency (Hz)" ); axes[1 ].set_ylabel("Spectral amplitude" )
axes[1 ].set_title("Displacement Spectrum — Brune Source Model" )
axes[1 ].legend(); axes[1 ].grid(True , which='both' , alpha=0.3 )
plt.tight_layout()
plt.savefig("seismic_spectrum.png" , dpi=150 )
plt.show()
Advanced Usage
Receiver Function Computation from obspy import Stream, Trace
import numpy as np
from scipy.signal import correlate, fftconvolve
def compute_receiver_function (Z_data, R_data, dt, water_level=0.01 , f_gauss=2.5 ):
"""
Estimate P-to-S receiver function via iterative time-domain deconvolution.
Parameters
----------
Z_data : array — vertical component (reference)
R_data : array — radial component (response)
dt : float — sample interval (s)
water_level : float — damping for spectral division
f_gauss : float — Gaussian filter width (Hz)
Returns
-------
rf : array — receiver function time series
t : array — time axis (s)
"""
N = len (Z_data)
freq = np.fft.rfftfreq(N, d=dt)
gauss = np.exp(-(freq/(2 *f_gauss))**2 )
Z_f = np.fft.rfft(Z_data)
R_f = np.fft.rfft(R_data)
denom = np.abs (Z_f)**2 + water_level * np.max (np.abs (Z_f)**2 )
rf_f = R_f * np.conj(Z_f) / denom * gauss
rf = np.fft.irfft(rf_f, n=N)
t = np.arange(N) * dt - N*dt/2
return np.roll(rf, N//2 ), t
np.random.seed(42 )
dt = 0.025
N = int (100 /dt)
t_syn = np.arange(N) * dt - 50
from scipy.signal import ricker
Z = ricker(N, 10 ) + np.random.randn(N) * 0.01
ps_delay = int (5 / dt)
R = np.roll(ricker(N, 10 ), ps_delay) * 0.4 + np.random.randn(N) * 0.01
rf, t_rf = compute_receiver_function(Z, R, dt)
import matplotlib.pyplot as plt
fig, axes = plt.subplots(3 , 1 , figsize=(10 , 7 ), sharex=False )
axes[0 ].plot(t_syn, Z, 'k-' ); axes[0 ].set_title("Vertical (Z)" ); axes[0 ].grid(alpha=0.3 )
axes[1 ].plot(t_syn, R, 'b-' ); axes[1 ].set_title("Radial (R)" ); axes[1 ].grid(alpha=0.3 )
mask = (t_rf > -5 ) & (t_rf < 20 )
axes[2 ].plot(t_rf[mask], rf[mask], 'r-' , linewidth=1.5 )
axes[2 ].axvline(0 , color='gray' , linestyle='--' ); axes[2 ].axhline(0 , color='gray' , linestyle='-' , linewidth=0.5 )
axes[2 ].set_title("Receiver Function" ); axes[2 ].set_xlabel("Time (s)" ); axes[2 ].grid(alpha=0.3 )
plt.tight_layout(); plt.savefig("receiver_function.png" , dpi=150 ); plt.show()
Troubleshooting
Error: No data available from FDSN client Cause : Station or channel not available for the requested time window.
from obspy.clients.fdsn import Client
client = Client("IRIS" )
inv = client.get_stations(network="IU" , station="MAJO" , level="channel" ,
starttime=UTCDateTime("2011-01-01" ),
endtime=UTCDateTime("2012-01-01" ))
print (inv)
Error: NonLinearLSQError in Brune fit Cause : Poor initial guess for curve_fit.
p0 = [spec_fit[0 ], freq_fit[np.argmax(np.diff(np.log(spec_fit+1e-30 )) < -0.1 )]]
popt, _ = curve_fit(brune_model, freq_fit, spec_fit, p0=p0, maxfev=10000 )
Version Compatibility Package Tested versions Notes obspy 1.4.x Stable API; 1.4 added FDSN mass downloader scipy 1.11, 1.12 curve_fit behavior stablenumpy 1.24, 1.26 No issues
External Resources
Official Documentation
Key Papers
Beyreuther, M. et al. (2010). ObsPy: A Python Toolbox for Seismology . Seismological Research Letters, 81(3), 530–533.
Examples
Example 1: Earthquake Catalog Download from obspy import UTCDateTime
from obspy.clients.fdsn import Client
client = Client("USGS" )
start = UTCDateTime("2023-01-01" )
end = UTCDateTime("2023-12-31" )
catalog = client.get_events(
starttime=start, endtime=end,
minmagnitude=6.0 ,
orderby="magnitude" ,
)
print (f"Found {len (catalog)} events (M≥6) in 2023:" )
for event in catalog[:5 ]:
mag = event.magnitudes[0 ].mag
orig = event.origins[0 ]
print (f" M{mag:.1 f} {orig.time.date} "
f"lat={orig.latitude:.2 f} , lon={orig.longitude:.2 f} , "
f"depth={orig.depth/1000 :.1 f} km" )
Example 2: Waveform Filtering and Spectrogram import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import spectrogram
np.random.seed(1 )
sps = 100.0
N = int (120 * sps)
t = np.arange(N) / sps
data = np.random.randn(N) * 0.1
p_idx = int (30 * sps)
data[p_idx:p_idx+int (20 *sps)] += np.sin(2 *np.pi*5 *(t[:int (20 *sps)])) * np.exp(-t[:int (20 *sps)]/5 )
s_idx = int (60 * sps)
data[s_idx:s_idx+int (40 *sps)] += np.sin(2 *np.pi*2 *(t[:int (40 *sps)])) * np.exp(-t[:int (40 *sps)]/8 ) * 2
f_spec, t_spec, Sxx = spectrogram(data, fs=sps, nperseg=256 , noverlap=200 )
fig, axes = plt.subplots(2 , 1 , figsize=(12 , 6 ), sharex=False )
axes[0 ].plot(t, data, 'k-' , linewidth=0.5 )
axes[0 ].axvline(30 , color='r' , linestyle='--' , label='P arrival' ); axes[0 ].axvline(60 , color='b' , linestyle='--' , label='S arrival' )
axes[0 ].set_ylabel("Amplitude" ); axes[0 ].set_title("Seismogram" ); axes[0 ].legend(); axes[0 ].grid(alpha=0.3 )
axes[1 ].pcolormesh(t_spec, f_spec[f_spec<20 ], 10 *np.log10(Sxx[f_spec<20 ]+1e-20 ), shading='gouraud' , cmap='viridis' )
axes[1 ].set_ylabel("Frequency (Hz)" ); axes[1 ].set_xlabel("Time (s)" ); axes[1 ].set_title("Spectrogram" )
plt.tight_layout(); plt.savefig("spectrogram_seismic.png" , dpi=150 ); plt.show()
Last updated: 2026-03-17 | Maintainer: @xjtulyc
Issues: GitHub Issues
More from this repository
Related occupations SOC
Based on SOC occupation classification
