Seismology Data Guide

from obspy.clients.fdsn import Client
from obspy import UTCDateTime

client = Client("IRIS")

# Fetch earthquake catalog for a region and time window
catalog = client.get_events(
    starttime=UTCDateTime("2024-01-01"),
    endtime=UTCDateTime("2024-12-31"),
    minmagnitude=5.0,
    maxmagnitude=9.0,
    minlatitude=30.0, maxlatitude=45.0,
    minlongitude=125.0, maxlongitude=150.0,
    orderby="magnitude",
)

print(f"Found {len(catalog)} events")
for event in catalog[:5]:
    origin = event.preferred_origin()
    mag = event.preferred_magnitude()
    print(f"  M{mag.mag:.1f} {origin.time} "
          f"({origin.latitude:.2f}, {origin.longitude:.2f}) "
          f"depth={origin.depth/1000:.1f} km")

from obspy import UTCDateTime
from obspy.clients.fdsn import Client

client = Client("IRIS")

# Download waveform data for a specific event
t = UTCDateTime("2024-01-01T07:10:00")
st = client.get_waveforms(
    network="IU", station="ANMO", location="00", channel="BHZ",
    starttime=t, endtime=t + 600,  # 10 minutes of data
)

# Standard preprocessing pipeline
st.detrend("demean")        # Remove mean
st.detrend("linear")        # Remove linear trend
st.taper(max_percentage=0.05, type="cosine")  # Taper edges
st.filter("bandpass", freqmin=0.01, freqmax=5.0, corners=4)

# Remove instrument response to get ground velocity (m/s)
inv = client.get_stations(
    network="IU", station="ANMO", location="00", channel="BHZ",
    starttime=t, endtime=t + 600, level="response",
)
st.remove_response(inventory=inv, output="VEL", pre_filt=[0.005, 0.01, 8, 10])

import numpy as np
from scipy.signal import welch

def compute_psd(trace, nperseg=256):
    """
    Compute power spectral density of a seismic trace.
    Returns frequencies (Hz) and PSD (dB relative to 1 (m/s)^2/Hz).
    """
    freqs, psd = welch(
        trace.data,
        fs=trace.stats.sampling_rate,
        nperseg=nperseg,
        noverlap=nperseg // 2,
    )
    psd_db = 10 * np.log10(psd + 1e-30)
    return freqs, psd_db

from obspy.signal.trigger import recursive_sta_lta, trigger_onset

def pick_arrivals(trace, sta_seconds=1.0, lta_seconds=30.0,
                  threshold_on=3.5, threshold_off=1.0):
    """
    STA/LTA trigger for P-wave arrival detection.
    sta_seconds: short-term average window
    lta_seconds: long-term average window
    Returns list of (on_sample, off_sample) trigger windows.
    """
    df = trace.stats.sampling_rate
    cft = recursive_sta_lta(
        trace.data,
        int(sta_seconds * df),
        int(lta_seconds * df),
    )
    triggers = trigger_onset(cft, threshold_on, threshold_off)
    return triggers, cft

# Simplified grid search earthquake location
def grid_search_locate(stations, arrival_times, velocity_model,
                       lat_range, lon_range, depth_range, grid_spacing):
    """
    Brute-force grid search for earthquake location.
    Minimizes sum of squared travel-time residuals.
    """
    best_misfit = float("inf")
    best_location = None

    for lat in np.arange(*lat_range, grid_spacing):
        for lon in np.arange(*lon_range, grid_spacing):
            for depth in np.arange(*depth_range, grid_spacing):
                residuals = []
                for sta, obs_time in zip(stations, arrival_times):
                    dist = geodetic_distance(lat, lon, sta.lat, sta.lon)
                    pred_time = velocity_model.get_travel_time(dist, depth)
                    residuals.append((obs_time - pred_time) ** 2)
                misfit = sum(residuals)
                if misfit < best_misfit:
                    best_misfit = misfit
                    best_location = (lat, lon, depth)

    return best_location, best_misfit

def moment_magnitude(seismic_moment_nm: float) -> float:
    """
    Compute moment magnitude from seismic moment (in Newton-meters).
    Mw = (2/3) * log10(M0) - 6.07  (Hanks and Kanamori, 1979)
    """
    return (2.0 / 3.0) * np.log10(seismic_moment_nm) - 6.07