Skip to main content Skills Marktplatz Entdecken und erkunden Sie KI-Skills, die von der Community erstellt wurden.
Mit Codex oder Claude installieren Kopieren Sie diesen Prompt, fügen Sie ihn in Codex, Claude oder einen anderen Assistant ein und lassen Sie die Skill-Seite prüfen und installieren.
Prompt-Details anzeigen Ein direkter Befehl überspringt den Prüf-Prompt. Prüfen Sie die Quelle, bevor Sie ihn ausführen.
npx skills add https://github.com/xjtulyc/awesome-rosetta-skills --skill obspy-seismologyDer Befehl bleibt in einer Zeile. Scrollen Sie horizontal, um ihn vor dem Kopieren vollständig zu prüfen.
Sie bevorzugen eine lokale Kopie? Laden Sie die Dateien herunter, die SkillsMP derzeit vorliegen.
ZIP herunterladen Herunterladen... Mehr aus diesem Repository
Verwandte Berufe SOC
Basierend auf der SOC-Berufsklassifikation
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
