obspy-process-earthquake

name

obspy-process-earthquake

description

Process downloaded earthquake waveform data with ObsPy. Use when Codex needs to preprocess MiniSEED, SAC, or ObsPy-readable seismic waveforms by reading StationXML/Inventory metadata, removing instrument response, removing mean, detrending, tapering, filtering, trimming, resampling, decimating/interpolating, predicting seismic phase arrivals with TauP, cutting waveforms around predicted arrivals, converting to SAC, setting SAC headers, saving processed waveform files, or generating quick quality-control summaries.

ObsPy Earthquake Processing

Workflow

Turn the user's processing request into a reproducible pipeline before modifying data.

1. Identify inputs:
   • waveform files: MiniSEED, SAC, or another ObsPy-readable format
   • response metadata: StationXML/Inventory preferred for response removal
   • event metadata when predicting arrivals or phase-window trimming: origin time, latitude, longitude, depth
   • station coordinates when predicting arrivals, often available from StationXML
   • target output directory and output format, usually MiniSEED or SAC
2. Clarify processing choices only when missing values materially change results:
   • response output units: VEL, DISP, ACC, or DEF
   • filter type and corners, especially bandpass freqmin/freqmax
   • whether to preserve raw files and write processed copies
   • whether to trim to an event window or process full traces
   • resampling target rate or integer decimation factor
   • TauP model and phases when predicting arrivals, e.g. iasp91 with P,S
   • phase-relative window, e.g. P-30s to P+120s
3. Use the standard sequence unless the user specifies otherwise:
   • read waveform and inventory
   • merge/split and remove empty traces as needed
   • remove mean with detrend("demean")
   • remove trend with detrend("linear")
   • taper with a small Hann taper such as max_percentage=0.05
   • remove instrument response with remove_response(inventory=inv, output=..., pre_filt=..., water_level=...)
   • optionally repeat demean/detrend/taper before time-domain filtering
   • filter with filter("bandpass", freqmin=..., freqmax=..., corners=4, zerophase=True) or the requested filter
   • optionally predict arrivals with TauPyModel.get_travel_times_geo()
   • optionally trim with Stream.trim() using absolute times or arrival-relative windows
   • optionally change sampling rate with interpolate(), resample(), or decimate()
   • optionally write SAC files and set useful SAC headers (evla, evlo, evdp, stla, stlo, o, t0, t1)
   • write processed output and a processing summary
4. Never overwrite raw data unless the user explicitly asks. ObsPy processing methods mutate streams/traces in place, so use .copy() when preserving originals.

Correct API Choices

• Use Stream.remove_response() or Trace.remove_response() to deconvolve instrument response.
• Do not use Inventory.remove() for instrument correction. That method returns a copy of the inventory with selected networks/stations/channels removed; it is useful for pruning metadata, not removing instrument response from waveform samples.
• Use Stream.attach_response() only when you need to attach responses to traces before later calling remove_response().

Quick Pattern

from obspy import read, read_inventory

st_raw = read("waveforms/*.mseed")
inv = read_inventory("stations/*.xml")

st = st_raw.copy()
st.merge(method=1, fill_value="interpolate")
st.detrend("demean")
st.detrend("linear")
st.taper(max_percentage=0.05, type="hann")

st.remove_response(
    inventory=inv,
    output="VEL",
    pre_filt=(0.005, 0.01, 20.0, 25.0),
    water_level=None,
)

st.detrend("demean")
st.detrend("linear")
st.taper(max_percentage=0.05, type="hann")
st.filter("bandpass", freqmin=0.01, freqmax=1.0, corners=4, zerophase=True)
st.write("processed.mseed", format="MSEED")

Arrival prediction and phase-window trim:

from obspy import UTCDateTime
from obspy.taup import TauPyModel

origin = UTCDateTime("2026-02-22T16:57:46.698")
model = TauPyModel(model="iasp91")
arrivals = model.get_travel_times_geo(
    source_depth_in_km=629.0,
    source_latitude_in_deg=6.8241,
    source_longitude_in_deg=116.2475,
    receiver_latitude_in_deg=43.093,
    receiver_longitude_in_deg=-79.31,
    phase_list=["P", "S"],
)
p_arrival = origin + arrivals[0].time
st.trim(p_arrival - 30, p_arrival + 120)

Sampling-rate and SAC output examples:

st.interpolate(20.0, method="lanczos", a=20)
st.decimate(5)
for tr in st:
    tr.stats.sac = {"o": 0.0, "evla": 6.8241, "evlo": 116.2475, "evdp": 629.0}
    tr.write(f"{tr.id.replace('.', '_')}.sac", format="SAC")

Defaults

Use these defaults when the user asks for ordinary teleseismic earthquake processing and does not specify parameters:

• output unit: VEL for broadband velocity sensors; ACC for strong-motion acceleration channels (HN?, EN?) if the user wants acceleration
• response pre_filt: choose from the data and requested passband; keep f2 below the analysis band and f3 above it
• water level: prefer water_level=None with a sensible pre_filt, especially if converting to a non-native unit
• taper: 5 percent Hann taper
• filter: do not invent a final science band silently; if needed for teleseismic body/surface-wave inspection, propose a band such as 0.01-1.0 Hz and state it
• arrival model: iasp91 for ordinary teleseismic examples unless the user requests another Earth model
• phase list: P,S for first-pass body-wave windows; for large epicentral distances, expect Pdiff/Sdiff instead of direct P/S
• sampling-rate changes: prefer interpolate(..., method="lanczos", a=20) for arbitrary target rates; use decimate() only for integer downsampling after anti-alias protection

Bundled Resources

• scripts/process_waveforms.py: reusable command-line processor for MiniSEED plus StationXML, optional TauP arrivals, phase trimming, resampling, and SAC output.
• references/obspy-processing-guide.md: API notes, parameter choices, and documentation links. Read it before implementing a custom processing script.