pnsn-phase-detection

# PNSN Phase Detection and Earthquake Monitoring This skill wraps the skill-local `pnsn/` project stored at `seismo_skill/skills/pnsn_phase_detection/pnsn/` for automatic seismic phase picking and event monitoring. Use it when a task involves: - detecting or picking Pg, Sg, Pn, Sn, P, or S arrivals; - running deep-learning pickers on continuous 3-component waveforms; - scanning SAC, MSEED, SEED, or miniseed waveform directories; - converting picks into earthquake events with FastLink, REAL, or GaMMA; - building an earthquake monitoring workflow from waveform data. Project root assumptions: - Run commands from the SAGE repo root. - The pnsn code and model files are managed inside this skill folder, not at the repository root. - Reusable instructions, scripts, and workflows are managed inside this skill folder. Do not depend on temporary demo folders such as `publish_mini`; they may be cleaned or omitted from uploads. - If `seismo_skill/skills/pnsn_phase_detection/pnsn/picker.py` is missing, instruct the user to run `git clone https://github.com/cangyeone/pnsn.git seismo_skill/skills/pnsn_phase_detection/pnsn` from the SAGE repo root. - Default model files live at `seismo_skill/skills/pnsn_phase_detection/pnsn/pickers/`. - Picker configuration lives at `seismo_skill/skills/pnsn_phase_detection/pnsn/config/picker.py`. - Define `PNSN_ROOT = Path("seismo_skill/skills/pnsn_phase_detection/pnsn")` in generated code and build all model/script paths from it. Preferred Python API: ```python from seismo_skill.skills.pnsn_phase_detection.pnsn import PNSNPicker picker = PNSNPicker() # uses pnsn/pickers/pnsn.v3.jit by default picks = picker.pick_stream(obspy_stream, incomplete="skip") ``` Use `pick_stream` for uploaded or already-read three-component waveforms in chat, science-analysis, and parameter-optimization tasks. Use `pick_directory` for batch monitoring directories. Do not silently fall back to STA/LTA. If the PNSN picker is unavailable, print a clear `[SAGE_TEST] PNSN unavailable: ...` diagnostic and stop unless the user explicitly requested STA/LTA or another classical trigger method. Default generated code pattern: ```python from seismo_skill.skills.pnsn_phase_detection.pnsn import PNSNPicker if not PNSNPicker.is_available(): raise RuntimeError( "PNSNPicker unavailable. Check seismo_skill/skills/pnsn_phase_detection/pnsn/pickers/" ) picker = PNSNPicker() picks = picker.pick_stream(st, incomplete="skip") if not picks: raise RuntimeError("PNSNPicker returned no picks for this stream") print(f"[SAGE_TEST] PNSN picks: {len(picks)}") ``` ## SeismicX-Cont / Continuous HDF5 Mode Use this mode when the task mentions `publish_mini`, SeismicX-Cont, continuous HDF5 waveforms, picker JSONL, label JSON, recall/precision, or annotation plots. The demo folder may be absent; the durable skill resources are: - `references/seismicx_cont_picker.md`: data layout, JSONL pick schema, metrics, and benchmark rules. - `references/annotation_plotting.md`: waveform label/auto-pick plotting rules. - `scripts/plot_picks_and_labels.py`: reusable annotation plotting script. - `workflows/seismicx_continuous_dataset_creation.md`: build a continuous waveform dataset from user data. - `workflows/picker_benchmark_and_annotation_plots.md`: evaluate a provided picker and generate recall/precision figures. - `workflows/continuous_detection_and_association.md`: run continuous detection and associate picks into events. Large SeismicX-Cont waveform data are not stored in this skill. Point users to: ```bash modelscope download --dataset cangyeone/SeismicX-Cont --local_dir /path/to/SeismicX-Cont ``` or [https://www.modelscope.cn/datasets/cangyeone/SeismicX-Cont](https://www.modelscope.cn/datasets/cangyeone/SeismicX-Cont). If a project contains `scripts/run_picker_to_jsonl.py`, use it for continuous HDF5 picking because it normally handles the project dataloader, resume, and JSONL output. If it is missing, write custom code using the contracts in `references/seismicx_cont_picker.md`. Example annotation plot command: ```bash python seismo_skill/skills/pnsn_phase_detection/scripts/plot_picks_and_labels.py \ --project-root /path/to/project \ --h5-input "data/hdf5/*.h5" \ --label-json data/label/annotations_mini_two_hours.json \ --auto-jsonl data/picks/pnsn.v3.diff.phase.jsonl \ --outdir "$SAGE_OUTDIR/annotation_plots" \ --max-panels 12 \ --window-seconds 180 ``` Science-analysis and parameter-optimization agents should save picker metrics, annotation plots, and manifests under the current project output directory and then cite them as evidence in reports or papers. ## Autonomous Programming Mode This skill is not limited to calling `seismo_skill/skills/pnsn_phase_detection/pnsn/picker.py`. When the user asks to "自己编程", "写代码实现拾取", "检测这几条波形", "画出拾取结果", "调试检测流程", or when the waveform organization does not match `seismo_skill/skills/pnsn_phase_detection/pnsn/config/picker.py`, write a custom Python program instead of only giving a CLI command. Use custom code for: - one station or a small set of files; - uploaded SAC/MSEED files with known paths; - non-standard filenames or component names; - tasks requiring waveform plots, pick tables, SNR/AMP statistics, or intermediate diagnostics; - debugging and mini-tests; - integrating picks into a larger SAGE analysis pipeline. Important fallback rule: - If `pnsn/picker.py` and `pnsn/pickers/pnsn.v3.jit` exist, prefer the PNSN TorchScript picker over a naive STA/LTA script. - Do not write a classical STA/LTA fallback for a generic "pick phases" request. Use STA/LTA only when the user explicitly asks for STA/LTA/classical triggering, or when the code is clearly labeled as a diagnostic comparison. - If you must write a classical STA/LTA fallback, never use the first `trigger_onset` window as the final pick. Print candidate trigger windows, ignore edge triggers near the record start, compare candidates with waveform energy/SNR, and choose plausible P/S arrivals. A common failure is picking a taper/filter transient within the first few seconds while the real event is much later in the trace. ## Plot Existing Pick Results When the user asks to draw or overlay an existing pick result, such as “把这个拾取结果绘制到波形上”, this is a visualization task, not a new phase-picking task. Do not re-run STA/LTA just because a file named `picks_table.csv` is missing. Search the current execution directory, `SAGE_OUTDIR`, and the authorized waveform/data directories for recent pick outputs whose names look like `*pick*.csv`, `*pick*.txt`, `pnsn_picks.csv`, `sage_picks_*.txt`, `phase_picks.*`, or `picks.*`. Accept these schemas: - CSV columns including `phase`, `time_abs`/`absolute_time`/`time`, or `relative_time_s`/`time_rel_s`; - PNSN text outputs with comment headers and comma rows: ```text # path/to/waveform/file phase_name,relative_time_s,confidence,absolute_time,SNR,AMP,station,extra ``` Skip comment lines, parse comma rows, and preserve both relative and absolute times when available. Only use a generic `data.csv` if it contains explicit phase-pick columns; otherwise reject it as unrelated data. When picks are produced in the same script by `PNSNPicker.pick_stream(st)`, pass those pick dictionaries directly into `plot_stream(st, picks=picks, ...)`. `plot_stream` accepts the PNSN keys `time_abs` and `time_rel_s`, so do not discard or re-filter valid PNSN picks into an empty plotting list. Do not convert fresh PNSN picks with legacy text-file fields such as `phase_name`/`absolute_time`; those names are only for parsing existing PNSN text output files. Custom picking code should implement this workflow: 1. Discover waveform files. 2. Group files into 3-component station sets. 3. Read with ObsPy. 4. Align traces to a common time window. 5. Detrend, taper, bandpass, and resample to 100 Hz. 6. Stack data as `[n_samples, 3]`. 7. Load a TorchScript picker with `torch.jit.load`. 8. Run inference under `torch.no_grad()`. 9. Convert sample indices to relative and absolute times. 10. Save a CSV/text pick table. 11. Generate a waveform figure with vertical pick markers. 12. Print `[SAGE_TEST]` and file paths. Prefer writing all outputs to `SAGE_OUTDIR` when available: ```python import os from pathlib import Path OUTDIR = Path(os.environ.get("SAGE_OUTDIR", "outputs/pnsn_phase_detection")) OUTDIR.mkdir(parents=True, exist_ok=True) PNSN_ROOT = Path("seismo_skill/skills/pnsn_phase_detection/pnsn") if not (PNSN_ROOT / "picker.py").exists(): raise FileNotFoundError( "Missing skill-local pnsn/. Run: git clone https://github.com/cangyeone/pnsn.git " "seismo_skill/skills/pnsn_phase_detection/pnsn" ) ``` ## Component Discovery Pattern When writing custom code, do not assume exact filenames unless the user gives them. Search common waveform suffixes and group by station/time key. ```python from pathlib import Path from collections import defaultdict import re def discover_waveforms(root): root = Path(root).expanduser() suffixes = {".sac", ".SAC", ".mseed", ".MSEED", ".seed", ".SEED", ".miniseed"} return [p for p in root.rglob("*") if p.is_file() and p.suffix in suffixes] def component_of(path): name = path.name.upper() # Handles BHE/BHN/BHZ, HHE/HHN/HHZ, E/N/Z endings. for comp in ["BHE", "BHN", "BHZ", "HHE", "HHN", "HHZ", "EHE", "EHN", "EHZ"]: if comp in name: return comp[-1] m = re.search(r"([ENZ])(?:\.[^.]+)?$", name) return m.group(1) if m else None def station_key(path): parts = path.name.split(".") # Default SAGE/PNSN style: NET.STA.LOC.CHN... if len(parts) >= 4: return ".".join(parts[:3]) return path.parent.name def group_three_components(files): groups = defaultdict(dict) for p in files: comp = component_of(p) if comp: groups[station_key(p)][comp] = p return {k: v for k, v in groups.items() if {"E", "N", "Z"} <= set(v)} ``` If no 3-component group is found, print a clear diagnostic with example filenames and recommend updating grouping logic or `seismo_skill/skills/pnsn_phase_detection/pnsn/config/picker.py`. ## Robust Custom Picker Template Use this as the default pattern when coding the picking workflow yourself. Adapt `DATA_ROOT` when the user provides a path. ```python import os from pathlib import Path import numpy as np import pandas as pd import torch import obspy import matplotlib.pyplot as plt OUTDIR = Path(os.environ.get("SAGE_OUTDIR", "outputs/pnsn_phase_detection")) OUTDIR.mkdir(parents=True, exist_ok=True) DATA_ROOT = Path("/path/to/waveforms") PNSN_ROOT = Path("seismo_skill/skills/pnsn_phase_detection/pnsn") if not (PNSN_ROOT / "picker.py").exists(): raise FileNotFoundError( "Missing skill-local pnsn/. Run: git clone https://github.com/cangyeone/pnsn.git " "seismo_skill/skills/pnsn_phase_detection/pnsn" ) MODEL_PATH = PNSN_ROOT / "pickers" / "pnsn.v3.jit" DEVICE = torch.device("cuda:0" if torch.cuda.is_available() else "cpu") PHASE_NAMES = {0: "Pg", 1: "Sg", 2: "Pn", 3: "Sn"} def prepare_stream(paths, target_rate=100.0): st = obspy.Stream() for p in paths: st += obspy.read(str(p)) st.merge(method=1, fill_value="interpolate") st.detrend("demean") st.detrend("linear") st.taper(0.02) st.filter("bandpass", freqmin=1.0, freqmax=10.0, corners=4, zerophase=True) if abs(st[0].stats.sampling_rate - target_rate) > 1e-6: st.resample(target_rate) start = max(tr.stats.starttime for tr in st) end = min(tr.stats.endtime for tr in st) if end <= start: raise ValueError("No overlapping time window among three components") st.trim(start, end, pad=False) return st def stream_to_array(st): comp_map = {} for tr in st: ch = tr.stats.channel.upper() if ch.endswith("E"): comp_map["E"] = tr elif ch.endswith("N"): comp_map["N"] = tr elif ch.endswith("Z"): comp_map["Z"] = tr missing = {"E", "N", "Z"} - set(comp_map) if missing: raise ValueError(f"Missing components: {sorted(missing)}") n = min(len(comp_map[c].data) for c in ["E", "N", "Z"]) x = np.stack([comp_map[c].data[:n] for c in ["E", "N", "Z"]], axis=1) return x.astype(np.float32), comp_map["Z"].stats.starttime, comp_map["Z"].stats.sampling_rate def run_picker(x, model): with torch.no_grad(): picks = model(torch.tensor(x, dtype=torch.float32, device=DEVICE)).cpu().numpy() rows = [] for phase_type, sample, confidence in picks: rows.append({ "phase": PHASE_NAMES.get(int(phase_type), str(int(phase_type))), "sample": int(sample), "relative_time_s": float(sample) / 100.0, "confidence": float(confidence), }) return rows model = torch.jit.load(str(MODEL_PATH), map_location=DEVICE).to(DEVICE).eval() files = discover_waveforms(DATA_ROOT) groups = group_three_components(files) print(f"[INFO] found_files={len(files)} three_component_groups={len(groups)}") all_rows = [] for station, comps in sorted(groups.items()): try: paths = [comps["E"], comps["N"], comps["Z"]] st = prepare_stream(paths) x, starttime, fs = stream_to_array(st) rows = run_picker(x, model) for r in rows: r["station"] = station r["absolute_time"] = str(starttime + r["relative_time_s"]) r["model"] = str(MODEL_PATH) all_rows.extend(rows) if rows: fig_path = OUTDIR / f"{station.replace('.', '_')}_picks.png"