Seismogram

Python Code

"""Seismology — Seismogram (ObsPy TSPAIR sample data)."""


import io
import numpy as np
import pandas as pd
import requests
import plotly.graph_objects as go

URL = "https://raw.githubusercontent.com/obspy/obspy/master/obspy/io/ascii/tests/data/tspair.ascii"


def synthetic_seismogram():
    """Generate a realistic synthetic seismogram with P and S wave arrivals."""
    rng = np.random.default_rng(42)
    dt = 0.01  # 100 Hz
    duration = 60.0
    t = np.arange(0, duration, dt)
    n = len(t)

    # Noise floor
    signal = rng.normal(0, 15, n)

    # P-wave arrival at ~8 s
    t_P = 8.0
    idx_P = int(t_P / dt)
    t_rel_P = t[idx_P:] - t_P
    env_P = np.exp(-t_rel_P / 2.0) * np.sin(2 * np.pi * 3.0 * t_rel_P)
    amp_P = 300 * np.exp(-t_rel_P / 3.0) * env_P
    signal[idx_P:] += amp_P

    # S-wave arrival at ~18 s (stronger)
    t_S = 18.0
    idx_S = int(t_S / dt)
    t_rel_S = t[idx_S:] - t_S
    env_S = np.exp(-t_rel_S / 5.0) * np.sin(2 * np.pi * 1.5 * t_rel_S)
    amp_S = 1200 * np.exp(-t_rel_S / 8.0) * env_S
    signal[idx_S:] += amp_S

    return t, signal, t_P, t_S


def parse_tspair(text):
    """Parse ObsPy TSPAIR ASCII format into time and amplitude arrays."""
    lines = text.strip().splitlines()
    times = []
    amplitudes = []
    for line in lines:
        line = line.strip()
        if not line:
            continue
        # Skip header lines (TIMESERIES ...)
        if line.upper().startswith("TIMESERIES") or line.upper().startswith("#"):
            continue
        parts = line.split()
        if len(parts) >= 2:
            try:
                amp = float(parts[-1])
                amplitudes.append(amp)
                times.append(len(times))  # index-based; convert to seconds later
            except ValueError:
                continue
    return np.array(times), np.array(amplitudes)


def generate():
    print("fetching seismogram data …")
    used_synthetic = False
    t_s = None
    amplitude = None
    t_P, t_S = None, None

    try:
        r = requests.get(URL, verify=False, timeout=15)
        r.raise_for_status()
        t_idx, amplitude = parse_tspair(r.text)
        if len(amplitude) < 100:
            raise ValueError(f"too few samples: {len(amplitude)}")

        # Try to get sample rate from header
        header_line = ""
        for line in r.text.splitlines():
            if line.upper().startswith("TIMESERIES"):
                header_line = line
                break
        # TSPAIR header: TIMESERIES ..., 7201 samples, 100 sps, ...
        import re
        m = re.search(r"(\d+(?:\.\d+)?)\s+sps", header_line, re.IGNORECASE)
        sps = float(m.group(1)) if m else 100.0
        t_s = t_idx / sps
        print(f"  parsed {len(amplitude)} samples at {sps} sps, duration {t_s[-1]:.1f} s")

    except Exception as e:
        print(f"  fetch failed ({e}); using synthetic seismogram")
        used_synthetic = True
        t_s, amplitude, t_P, t_S = synthetic_seismogram()

    fig = go.Figure()
    fig.add_trace(go.Scatter(
        x=t_s,
        y=amplitude,
        mode="lines",
        line=dict(color=VIOLET, width=0.8),
        name="Seismogram",
        hovertemplate="t=%{x:.2f} s<br>amp=%{y:.1f}<extra></extra>",
    ))

    annotations = []
    shapes = []
    if t_P is not None:
        shapes.append(dict(
            type="line",
            x0=t_P, x1=t_P,
            y0=0, y1=1,
            yref="paper",
            line=dict(color=TEAL, width=1.5, dash="dash"),
        ))
        annotations.append(dict(
            x=t_P, y=0.92, xref="x", yref="paper",
            text="P-wave",
            showarrow=False,
            font=dict(color=TEAL, size=12),
            xanchor="left",
        ))
    if t_S is not None:
        shapes.append(dict(
            type="line",
            x0=t_S, x1=t_S,
            y0=0, y1=1,
            yref="paper",
            line=dict(color="#DA5597", width=1.5, dash="dash"),
        ))
        annotations.append(dict(
            x=t_S, y=0.92, xref="x", yref="paper",
            text="S-wave",
            showarrow=False,
            font=dict(color="#DA5597", size=12),
            xanchor="left",
        ))

    fig.update_layout(
        xaxis=dict(title="Time (s)"),
        yaxis=dict(title="Amplitude (counts)"),
        annotations=annotations,
        shapes=shapes,
        margin=dict(t=50, b=60, l=80, r=60),
        showlegend=False,
    )

    return fig


if __name__ == "__main__":
    generate()

Made with Plotly