Skip to main content Instalar con Codex o Claude Copia este prompt, pégalo en Codex, Claude u otro asistente, y deja que revise la página de la skill y la instale por ti.
Un comando directo omite el prompt de revisión. Revisa el origen antes de ejecutarlo.
npx skills add https://github.com/xjtulyc/awesome-rosetta-skills --skill obspy-seismologyEl comando permanece en una sola línea. Desplázate horizontalmente para revisarlo antes de copiarlo.
¿Prefieres una copia local? Descarga los archivos que SkillsMP tiene disponibles ahora.
Ocupaciones relacionadasSOC
Basado en la clasificación ocupacional SOC
| 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 first
- S-waves: Shear waves, ~60% of P velocity, arrive second
- Surface 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(, , )
data = np.sin( * np.pi * * t) * np.exp(-t/)
tr = Trace(data=data)
tr.stats.sampling_rate =
tr.stats.network = ; tr.stats.station =
()
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.2f} {arr.ray_param_sec_degree:12.4f}")
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:.2f}s, OFF at t={off/sps:.2f}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:.4e}, fc = {fc_hat:.3f} 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:.3e} N·m")
print(f"Moment magnitude Mw = {Mw:.2f}")
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:.2f} 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 stable |
| numpy | 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:.1f} {orig.time.date} "
f"lat={orig.latitude:.2f}, lon={orig.longitude:.2f}, "
f"depth={orig.depth/1000:.1f} 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
