From 366f28e179653ffccff3fe8c87d0bc4a8923a467 Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Sun, 16 Aug 2026 17:50:03 -0700 Subject: [PATCH 1/6] Release 0.4.0: four archives, live Earth2Studio source MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Cut the Unreleased section — GeoNet, SeisfetchLiveSource, the NCEDC region fix — as 0.4.0, and bring the docs up to what the code does. README was written for three archives and the pre-live Earth2Studio adapters: no GeoNet anywhere, no SeisfetchLiveSource, NCEDC listed in us-east-2, and the notebook list stopped at 04 while 05 (response removal) and 06 (three-archive cross-correlation) went unlinked. Adds an "ObsPy-free, and what that buys" section so the footprint case and the bit-identical NoisePy equivalence are visible from the front page instead of only in benchmarks/. Also exports SeisfetchLiveSource at the top level, where the other two Earth2Studio adapters already live, and gives GeoNet an entry in THIRD_PARTY_NOTICES (which also still said NCEDC was us-east-2). Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 +++ README.md | 120 ++++++++++++++++++++++++++++++++++++++--- THIRD_PARTY_NOTICES.md | 9 +++- pyproject.toml | 2 +- seisfetch/__init__.py | 9 ++-- 5 files changed, 133 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7930446..aad4e44 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,12 @@ All notable changes to seisfetch are documented here. The format follows ## Unreleased +## 0.4.0 — 2026-08-16 + +Four-archive release: GeoNet joins SCEDC, NCEDC and EarthScope, and the +Earth2Studio integration becomes a live, time-indexed source that fetches on +demand and always returns physical units (PR #5). + ### Added - **GeoNet (New Zealand)** open-data S3 archive: `geonet-open-data` diff --git a/README.md b/README.md index a2fcde2..969be0c 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # seisfetch -Cloud-first seismic waveform access for EarthScope, SCEDC, NCEDC, and fallback FDSN services. +Cloud-first seismic waveform access for EarthScope, SCEDC, NCEDC, GeoNet, and fallback FDSN services. `seisfetch` is built around one core path: @@ -24,7 +24,7 @@ The design goal is simple: The intended acquisition order is: -1. `s3_open` for SCEDC and NCEDC open buckets +1. `s3_open` for SCEDC, NCEDC and GeoNet open buckets 2. `s3_auth` for EarthScope S3 access 3. `fdsn` only when the archive-backed path is unavailable or the network is not served from those buckets @@ -34,7 +34,8 @@ At the package level, the main interfaces are: - `SeisfetchClient.get_numpy()` -> `TraceBundle` of numpy arrays - `SeisfetchClient.get_xarray()` -> `xarray.Dataset` - `SeisfetchClient.get_waveforms()` -> ObsPy `Stream` -- `SeismicDataFrameSource` / `SeismicDataSource` -> Earth2Studio-compatible adapters +- `SeismicDataFrameSource` / `SeismicDataSource` -> Earth2Studio adapters over data you already fetched +- `SeisfetchLiveSource` -> time-indexed Earth2Studio `DataSource` that fetches on demand ## Workflow @@ -46,26 +47,31 @@ At the package level, the main interfaces are: |---|---|---| | `CI`, other SCEDC-routed networks | SCEDC open S3 | `s3_open` | | `BK`, other NCEDC-routed networks | NCEDC open S3 | `s3_open` | +| `NZ` | GeoNet open S3 | `s3_open` | | `IU`, `UW`, `TA`, other EarthScope-routed networks | EarthScope S3 | `s3_auth` | +Networks served by both SCEDC and NCEDC (`NC`, `NP`, ...) route to NCEDC. + Archive details: | Archive | Bucket | Region | Auth | |---|---|---|---| | EarthScope | `earthscope-geophysical-data` | `us-east-2` | EarthScope SDK credentials | | SCEDC | `scedc-pds` | `us-west-2` | none | -| NCEDC | `ncedc-pds` | `us-east-2` | none | +| NCEDC | `ncedc-pds` | `us-west-2` | none | +| GeoNet | `geonet-open-data` | `ap-southeast-2` | none | Notes: -- SCEDC and NCEDC are per-channel archives, so you should pass `channel=...`. +- SCEDC, NCEDC and GeoNet are per-channel archives, so you should pass `channel=...`. +- GeoNet channels always carry a numeric location code (`10`, `20`, ...), so `location=` is required there — a blank location raises rather than silently missing data. - EarthScope stores station-day miniSEED objects and currently requires authenticated access through `earthscope-sdk`. ### 2. FDSN second Use `backend="fdsn"` when: -- the desired network is not available from EarthScope / SCEDC / NCEDC S3 +- the desired network is not available from EarthScope / SCEDC / NCEDC / GeoNet S3 - the archive-backed attempt fails and you want a fallback provider - you need a non-US provider such as GEOFON, INGV, ETH, ORFEUS, etc. @@ -113,6 +119,47 @@ That makes it useful for: - Earth2Studio interoperability - data assimilation and digital twin workflows +## ObsPy-Free, and What That Buys + +ObsPy is an optional extra here, not a dependency of the data path. The +evaluation behind that choice — including whether it changes the science — is +written up in +[docs/noisepy-obspy-replacement-report.md](docs/noisepy-obspy-replacement-report.md), +with every number traceable to a committed JSON under `benchmarks/results/` or +a test in `tests/precision/`. + +| | seisfetch | obspy stack | +|---|---|---| +| Parse an 11 MB Steim2 channel-day | 21.4 ms | 37.2 ms | +| Cold `import` | 0.08 s | 0.13 s | +| Parse peak memory | 27.7 MB | 52.0 MB | +| Installed footprint | 80.4 MB | 311.4 MB | +| arm64 Linux install | wheels | needs gcc | + +The footprint line is the cloud argument: AWS Lambda caps a layer at 250 MB, +so the ObsPy stack does not fit and seisfetch does. ObsPy publishes no +linux/aarch64 wheels, so on Graviton it compiles from source; seisfetch and +pymseed install from wheels. + +### The science does not change + +Identical archive bytes were pushed through (A) `obspy.read` plus NoisePy's own +`preprocess_raw` and (B) seisfetch's parser plus the numpy/scipy ports in +`seisfetch.contrib.noisepy_adapter`, then through NoisePy's own `compute_fft` +and `correlate`. The pass criterion is bit-identity, not closeness: + +| Harness | Result | +|---|---| +| Single-station CI.PASC, EN/EZ/NZ/ZZ at 40 sps | max abs diff `0.0` | +| Cross-station SCEDC x NCEDC x EarthScope at 20 sps | max abs diff `0.0` | +| dv/v stretching | same grid cell on all pairs | + +The cross-station run matters because it exercises the Fourier-resample and +sub-sample-alignment branches inside the chain. Instrument response removal is +obspy-free end to end and agrees with ObsPy to 7.6e-16 of peak amplitude on +real data. Tables, plots and the harnesses: +[benchmarks/RESULTS.md](benchmarks/RESULTS.md). + ## Quick Start ### Open S3: SCEDC / NCEDC @@ -134,6 +181,22 @@ print(bundle.ids) arrays = bundle.to_dict() ``` +### Open S3: GeoNet (New Zealand) + +`NZ` auto-routes to the GeoNet open-data bucket. GeoNet channels carry a +numeric location code, so pass `location=`: + +```python +bundle = SeisfetchClient(backend="s3_open").get_numpy( + "NZ", + "WEL", + location="10", + channel="HHZ", + starttime="2022-01-02T00:00:00", + endtime="2022-01-02T00:10:00", +) +``` + ### EarthScope S3 (authenticated) Requires `earthscope-sdk` and an EarthScope account that has been granted the @@ -421,6 +484,7 @@ SeisfetchClient +- backend="s3_open" | +- SCEDC open bucket | +- NCEDC open bucket +| +- GeoNet open bucket | +- auto-routing by network code | +- backend="s3_auth" @@ -438,8 +502,9 @@ SeisfetchClient The package includes adapters in `seisfetch.earth2` for Earth2Studio-style usage: -- `SeismicDataSource` -- `SeismicDataFrameSource` +- `SeismicDataSource` — wraps a bundle you already fetched +- `SeismicDataFrameSource` — sparse sensor table; `auto_coords=True` fills station lat/lon from the FDSN station service +- `SeisfetchLiveSource` — time-indexed `DataSource` that fetches on demand - `bundle_to_earth2` These are intended for: @@ -455,6 +520,42 @@ Typical path: miniSEED -> numpy -> xarray / sparse dataframe -> Earth2Studio adapter ``` +### Live source + +`SeisfetchLiveSource` has the shape every other Earth2Studio source (GFS, +ERA5, ...) has: you call it with timestamps and it fetches, auto-routing per +network across all four archives and caching day bundles in memory. + +```python +from datetime import datetime +from seisfetch.earth2 import SeisfetchLiveSource + +source = SeisfetchLiveSource( + channels=["CI.PASC..BHZ", "BK.PKD.00.BHZ", "II.PFO.00.BHZ"], + window_s=3600, + calibrate="gain", +) +da = source(datetime(2022, 1, 2, 6)) # -> (time, variable, sample) DataArray +``` + +Channels are `NET.STA.LOC.CHA` strings; the returned `variable` coordinate +spells them with underscores (`CI_PASC__BHZ`), which is also what the +optional `variable=` argument accepts. All channels in one call must share a +sampling rate — request mixed rates (a 40 sps `BHZ` alongside a 100 sps +`HHZ`) in separate calls. + +Physical units are **required** — this source never returns raw counts: + +- `calibrate="gain"` (default): divide by the channel's total sensitivity from + the FDSN station service. Exact at the reference frequency, one metadata + request per channel, no extra dependencies. +- `calibrate="response"`: full spectral deconvolution through + `seisfetch.contrib.response` (StationXML fetch plus an evalresp-equivalent + evaluator). Agrees with ObsPy to 7.6e-16 of peak amplitude on real data. + +`fetch()` is genuinely async (`asyncio.to_thread`), so pipelines can overlap +this source with others. + ## Recipes A few common end-to-end tasks the package is designed for. @@ -565,6 +666,8 @@ See [notebooks/](notebooks/) for worked examples: - [02_bulk_mining.ipynb](notebooks/02_bulk_mining.ipynb) - [03_xarray_zarr_pipeline.ipynb](notebooks/03_xarray_zarr_pipeline.ipynb) - [04_earth2studio_interop.ipynb](notebooks/04_earth2studio_interop.ipynb) +- [05_response_removal.ipynb](notebooks/05_response_removal.ipynb) — instrument response removal without ObsPy +- [06_cross_correlation_three_archives.ipynb](notebooks/06_cross_correlation_three_archives.ipynb) — four stations, three archives, two months of NoisePy cross-correlations Notebook setup instructions are in [notebooks/README.md](notebooks/README.md). @@ -604,6 +707,7 @@ When using data accessed through `seisfetch`: - EarthScope: cite the network operators and NSF SAGE facility - SCEDC: doi:[10.7909/C3WD3xH1](https://doi.org/10.7909/C3WD3xH1) - NCEDC: doi:[10.7932/NCEDC](https://doi.org/10.7932/NCEDC) +- GeoNet: GNS Science, [GeoNet open data](https://www.geonet.org.nz/data/supplementary/channels) (CC BY 4.0) — cite per GeoNet's data policy - other FDSN providers: cite the underlying network/provider Software references: diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index f3ab662..0b778b5 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -106,13 +106,20 @@ authors and operators. ### NCEDC (Northern California Earthquake Data Center) - Operator: UC Berkeley Seismological Laboratory -- Bucket: `s3://ncedc-pds` (us-east-2) +- Bucket: `s3://ncedc-pds` (us-west-2) - Citation: doi:10.7932/NCEDC "Waveform data, metadata, or data products for this study were accessed through the Northern California Earthquake Data Center (NCEDC), doi:10.7932/NCEDC." - Reference: https://ncedc.org/db/cloud.html +### GeoNet (New Zealand) +- Operator: GNS Science (GeoNet programme) +- Bucket: `s3://geonet-open-data` (ap-southeast-2) +- License: CC-BY-4.0 +- Citation: cite GeoNet and the NZ network operators per GeoNet's data policy +- Reference: https://www.geonet.org.nz/data/supplementary/channels + --- ## Design Acknowledgments diff --git a/pyproject.toml b/pyproject.toml index e334a23..34ae792 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "seisfetch" -version = "0.3.1" +version = "0.4.0" description = "Fast seismic miniSEED from EarthScope, SCEDC, NCEDC, and 30+ FDSN servers. Core: numpy + boto3 + pymseed. No ObsPy required." readme = "README.md" # Compound: the package is MIT except seisfetch/contrib/obspy_ports.py, diff --git a/seisfetch/__init__.py b/seisfetch/__init__.py index 6a26721..4b0d431 100644 --- a/seisfetch/__init__.py +++ b/seisfetch/__init__.py @@ -1,5 +1,5 @@ """ -seisfetch: Fast seismic miniSEED from EarthScope, SCEDC, NCEDC, +seisfetch: Fast seismic miniSEED from EarthScope, SCEDC, NCEDC, GeoNet, and 37+ FDSN servers. Core deps: numpy + boto3 + pymseed. No ObsPy required. @@ -7,7 +7,8 @@ S3 archives: EarthScope s3://earthscope-geophysical-data (us-east-2, auth via earthscope-sdk) SCEDC s3://scedc-pds (us-west-2) - NCEDC s3://ncedc-pds (us-east-2) + NCEDC s3://ncedc-pds (us-west-2) + GeoNet s3://geonet-open-data (ap-southeast-2) Optional outputs: pandas → bundle_to_metadata_table(), write_metadata_csv() @@ -75,6 +76,7 @@ def __dir__(): # Earth2Studio adapters — lazy import (requires earth2studio + xarray) try: from seisfetch.earth2 import ( + SeisfetchLiveSource, SeismicDataFrameSource, SeismicDataSource, bundle_to_earth2, @@ -87,7 +89,7 @@ def __dir__(): __version__ = _pkg_version("seisfetch") except Exception: # not installed (e.g. vendored copy) - __version__ = "0.3.0" + __version__ = "0.4.0" __all__ = [ "SeisfetchClient", "S3OpenClient", @@ -119,5 +121,6 @@ def __dir__(): "requests_from_csv", "SeismicDataSource", "SeismicDataFrameSource", + "SeisfetchLiveSource", "bundle_to_earth2", ] From bcc78e02a99a16c43deaca55fabf354bb83788fc Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Sun, 16 Aug 2026 18:10:34 -0700 Subject: [PATCH 2/6] README: give response removal its own section MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit It was one sentence inside the obspy-free pitch, which undersells the piece of work that actually unblocked the migration — obspy has no pure-Python response evaluator at all, so contrib.response had to re-derive every stage and then match compiled evalresp. Adds the measured numbers (1.6e-10 vs evalresp, 6.6e-16 of peak vs Trace.remove_response, 1.9s vs 3.6s), both evaluation modes with the paz error envelope stated, both deconvolution styles, the loud-failure posture, and the not-implemented list. Also scopes the bit-identity claim: the 0.0 CCF result is the rm_resp=NO chain, and response removal is separate evidence. That was implicit and a reader could reasonably have merged the two. Co-Authored-By: Claude Opus 5 (1M context) --- README.md | 49 +++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 45 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 969be0c..4c3b21d 100644 --- a/README.md +++ b/README.md @@ -155,10 +155,51 @@ and `correlate`. The pass criterion is bit-identity, not closeness: | dv/v stretching | same grid cell on all pairs | The cross-station run matters because it exercises the Fourier-resample and -sub-sample-alignment branches inside the chain. Instrument response removal is -obspy-free end to end and agrees with ObsPy to 7.6e-16 of peak amplitude on -real data. Tables, plots and the harnesses: -[benchmarks/RESULTS.md](benchmarks/RESULTS.md). +sub-sample-alignment branches inside the chain. Tables, plots and the +harnesses: [benchmarks/RESULTS.md](benchmarks/RESULTS.md). + +Scope note: this bit-identity result covers the preprocessing chain at +`rm_resp=NO`. Response removal is validated separately, below. + +### Response removal without evalresp + +Response removal was the one ObsPy capability the NoisePy migration still +needed. `seisfetch.contrib.response` provides it in ~577 lines of numpy and +stdlib `xml.etree` — no evalresp C library, no ObsPy, no lxml, and no scipy in +that module. ObsPy has no pure-Python response evaluator (`remove_response` +calls compiled evalresp), so every stage was re-derived and then checked +against the compiled implementation: + +| Check | Result | +|---|---| +| `evaluate_response(mode="full")` vs compiled evalresp — both CI.PASC epochs, VEL/ACC/DISP, 1 mHz–19.9 Hz | max rel diff `1.6e-10` | +| `remove_response_np` vs `Trace.remove_response` — real 6.9M-sample Tohoku day, water_level=60, pre_filt | `6.6e-16` of peak | +| Same, CI.PASC demo hour in notebook 06 | `7.6e-16` of peak | +| Deconvolve that Tohoku day | 1.9 s vs ObsPy 3.6 s | + +Two evaluation modes: `mode="full"` is evalresp-equivalent (all stages, analog +and digital poles/zeros, FIR/Coefficients with DC normalization and the +`Decimation/CorrectionApplied` phase advance). `mode="paz"` is the SACPZ +shortcut — 0.7–1.3 % error below 4 Hz but ~23 % by 16 Hz, since the FIR +anti-alias roll-off is unmodeled; don't use it above ~Nyquist/3. + +Two deconvolution styles: `remove_response_np` ports ObsPy's water-level +method for drop-in equivalence, and `translate_resp_np` follows SeisIO.jl's +translation approach, taking stabilization from the target response's own +roll-off instead of a water level. + +Defective metadata fails loudly — zero or missing gains, degenerate +normalization references, zero-sum FIR stages, polynomial and `ResponseList` +stages all raise with the stage number named, never a silent NaN or unity +gain. Not implemented, and raising rather than approximating: IIR +`Coefficients` stages with denominators, polynomial (blockette-62) responses, +`ResponseList` stages. Metadata is StationXML only; RESP and SACPZ files are +not parsed. + +Full derivation, the conditional-A0 finding about evalresp's normalization +rule, and the SeisIO comparison: +[docs/response-removal-design.md](docs/response-removal-design.md). Tutorial: +[notebooks/05_response_removal.ipynb](notebooks/05_response_removal.ipynb). ## Quick Start From eaf867e832e7782cbeaff78e122bb000ed61db1f Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Sun, 16 Aug 2026 18:33:02 -0700 Subject: [PATCH 3/6] Evaluate pymseed 0.9.5 take_np_datasamples (EarthScope/pymseed#6) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Chad shipped take_np_datasamples() in 0.9.5 the day after the 0.9.4 keepalive landed, and asked on the issue for help testing it — he recommends it over the record_list path for the common case. Adds two profile_parse cases: the plain no-copy view (safe since 0.9.4) and the new take. Findings in docs/pymseed-issue-draft.md. It is correct — bit-identical to the record-list decode, survives GC of the trace list, one-shot and destructive (segment goes to numsamples 0). But on this channel-day it is 22.4 ms against 21.5 ms for the record-list path, so it is not a win for us: the record-list build we would be avoiding is the same pass that hands us the per-record encoding, and asking for both is 27.8 ms. The native tier is also the wrong place to judge a memory-traffic change; the cgroup tiers where the original issue lived need docker, which was unavailable. Separately confirms 0.4.0 is safe to publish: parse_mseed output is bit-identical between 0.9.4 and 0.9.5 on every fixture, so the noisepy equivalence carries over by construction, and both suites are green on 0.9.5. Pin >=0.6,<0.10 stands. Co-Authored-By: Claude Opus 5 (1M context) --- benchmarks/profile_parse.py | 33 ++++++++++++++++++++++++ docs/pymseed-issue-draft.md | 51 +++++++++++++++++++++++++++++++++++++ 2 files changed, 84 insertions(+) diff --git a/benchmarks/profile_parse.py b/benchmarks/profile_parse.py index 7da1766..18ed82f 100644 --- a/benchmarks/profile_parse.py +++ b/benchmarks/profile_parse.py @@ -89,6 +89,37 @@ def case_tracelist_unpack_copy(raw): return n +def case_tracelist_unpack_view(raw): + """unpack in C, then the no-copy view (safe since pymseed 0.9.4).""" + from pymseed import MS3TraceList + + tl = MS3TraceList.from_buffer(raw, unpack_data=True) + n = 0 + for tid in tl: + for seg in tid: + arr = seg.np_datasamples + n += arr.shape[0] + return n + + +def case_tracelist_take_np(raw): + """unpack in C, then take ownership of the buffer (pymseed >= 0.9.5). + + ``take_np_datasamples()`` transfers the decoded buffer to numpy with no + copy, so the array outlives the trace list without the 0.9.4 keepalive + pinning it. EarthScope/pymseed#6. + """ + from pymseed import MS3TraceList + + tl = MS3TraceList.from_buffer(raw, unpack_data=True) + n = 0 + for tid in tl: + for seg in tid: + arr = seg.take_np_datasamples() + n += arr.shape[0] + return n + + def case_parse_mseed(raw): from seisfetch.convert import parse_mseed @@ -149,6 +180,8 @@ def case_obspy_read(raw): "tl_reclist_only": case_tracelist_reclist_only, "tl_reclist->np": case_tracelist_reclist_np, "tl_unpack+copy": case_tracelist_unpack_copy, + "tl_unpack+view": case_tracelist_unpack_view, + "tl_unpack+take_np": case_tracelist_take_np, "parse_mseed": case_parse_mseed, "parse_mseed_flags": case_parse_mseed_flags, "parse_records_fb": case_parse_records_fallback, diff --git a/docs/pymseed-issue-draft.md b/docs/pymseed-issue-draft.md index 190812b..02e3af6 100644 --- a/docs/pymseed-issue-draft.md +++ b/docs/pymseed-issue-draft.md @@ -10,6 +10,57 @@ the copy), matching the record-list path; seisfetch tests pass on 0.9.4 (already inside the >=0.6,<0.10 pin). seisfetch stays on the record-list path for the per-record encoding metadata. Kept here for the record. +### Follow-up: pymseed 0.9.5 and `take_np_datasamples()` (2026-08-16) + +Chad released 0.9.5 on 2026-08-07 adding +`MS3TraceSeg.take_np_datasamples()` — transfers the decoded buffer to numpy +with no copy, so the array outlives the trace list without relying on the +0.9.4 keepalive — and asked on the issue for help testing it, recommending it +over `record_list=True` + `create_numpy_array_from_recordlist()` for the +common read-and-get-arrays case. Issue #6 is still open on that request. + +Evaluated on 0.9.5 (`benchmarks/profile_parse.py`, cases `tl_unpack+view` +and `tl_unpack+take_np`; same 11.3 MB CI.PASC.00.BHZ Steim2 channel-day, +min of 7, native macOS arm64): + +| case | 0.9.4 (ms) | 0.9.5 (ms) | +|---|---|---| +| `record_list=True` → `create_numpy_array_from_recordlist()` | 21.0 | 21.5 | +| `unpack_data=True` → `np_datasamples` (no copy, safe since 0.9.4) | 23.2 | 22.8 | +| `unpack_data=True` → `take_np_datasamples()` | n/a | **22.4** | +| `unpack_data=True` → `np_datasamples.copy()` (old idiom) | 25.0 | 26.8 | +| `seisfetch.parse_mseed` (record-list path) | 20.9 | 21.4 | + +Findings: + +- **Correct.** `take_np_datasamples()` output is bit-identical to the + record-list decode on the channel-day, survives `del tracelist` + GC, and + is `int32` as expected. Semantics are one-shot and destructive: after + taking, the segment reports `numsamples == 0` and a second take returns an + empty array. +- **Not faster than the record-list path here.** 22.4 ms vs 21.5 ms natively + — the two are within noise of each other, with the record-list path + marginally ahead. Chad's "huge win" framing is against the *copy* idiom and + the record-list *setup cost*; on this file the record-list build is only + ~4.5 ms and it is amortized by giving us what we need anyway. +- **seisfetch stays on the record-list path.** Not for speed but for + metadata: we want the per-record encoding, which `record_list=True` + provides in the same pass. Reaching the same place via `take_np` would mean + `unpack_data=True, record_list=True` (27.8 ms — slower than either) or + giving up the encoding. +- **Not yet tested under cgroup limits**, which is the tier the original + issue was about. Docker was unavailable for this pass; the native numbers + are the weakest place to judge a memory-traffic change. + +Release safety, checked separately: `parse_mseed` output is **bit-identical +between 0.9.4 and 0.9.5** across every fixture (Steim2 day, float32, float64, +int16, 3-segment gap, overlap) — sample SHA-256, npts, dtype, rate and +integer-ns start times all equal. The full unit suite (270 passed) and the +precision suite (53 passed) are green on 0.9.5. The noisepy equivalence +results therefore carry over by construction: path B of both harnesses starts +from `parse_mseed`, whose output does not change. The `>=0.6,<0.10` pin +stands. + --- ## Title From 0069b340fda2d4219ebf8a9782fd9ea88b1e37aa Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Sun, 16 Aug 2026 18:43:55 -0700 Subject: [PATCH 4/6] Run the pymseed 0.9.4/0.9.5 A/B under cgroup limits MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The native pass could not separate take_np_datasamples() from the record-list path, and native is the wrong tier for a memory-traffic question anyway. A first container pass was unusable — at 0.5 cpu the means ran 2-4x the mins and tl_unpack+view swung 40 to 88 ms between runs of identical code. Adds benchmarks/ab_pymseed_paths.py: only the four cases that matter, shuffled and interleaved each round so quota refill and host drift cancel instead of biasing whichever case runs first. Two containers x 15 reps per version per tier. That resolves it. take_np_datasamples() is indistinguishable from the plain 0.9.4 safe view (27.6 vs 27.7 ms fargate, 28.6 vs 28.3 lambda) — which makes sense, since the keepalive already made the view zero-copy, so the new call changes ownership rather than data movement. Its value is lifetime semantics, not throughput. The record-list path stays fastest in both tiers and leads by ~25% there against ~1 ms natively. Caveats recorded: single trace/segment file, and Docker Desktop on macOS is cgroup-in-VM rather than bare Linux. Co-Authored-By: Claude Opus 5 (1M context) --- benchmarks/ab_pymseed_paths.py | 75 ++++++++++++++++++++++++++++++++++ docs/pymseed-issue-draft.md | 39 ++++++++++++++++-- 2 files changed, 110 insertions(+), 4 deletions(-) create mode 100644 benchmarks/ab_pymseed_paths.py diff --git a/benchmarks/ab_pymseed_paths.py b/benchmarks/ab_pymseed_paths.py new file mode 100644 index 0000000..f943693 --- /dev/null +++ b/benchmarks/ab_pymseed_paths.py @@ -0,0 +1,75 @@ +"""Focused A/B of pymseed decode paths under cgroup limits. + +Only the four cases that matter, interleaved in a shuffled order each round +so host drift and CPU-quota refill cancel instead of biasing whichever case +runs first. +""" + +import gc +import json +import random +import statistics as st +import sys +import time +from pathlib import Path + +import pymseed +from pymseed import MS3TraceList + +RAW = (Path(__file__).resolve().parent.parent / "tests" / "bench.mseed").read_bytes() +REPS = int(sys.argv[1]) if len(sys.argv) > 1 else 15 + + +def reclist_np(raw): + tl = MS3TraceList.from_buffer(raw, unpack_data=False, record_list=True) + return sum( + seg.create_numpy_array_from_recordlist().shape[0] for tid in tl for seg in tid + ) + + +def view(raw): + tl = MS3TraceList.from_buffer(raw, unpack_data=True) + return sum(seg.np_datasamples.shape[0] for tid in tl for seg in tid) + + +def take_np(raw): + tl = MS3TraceList.from_buffer(raw, unpack_data=True) + return sum(seg.take_np_datasamples().shape[0] for tid in tl for seg in tid) + + +def copy(raw): + tl = MS3TraceList.from_buffer(raw, unpack_data=True) + return sum(seg.np_datasamples.copy().shape[0] for tid in tl for seg in tid) + + +CASES = {"reclist->np": reclist_np, "view": view, "take_np": take_np, "copy": copy} +if not hasattr( + next(iter(next(iter(MS3TraceList.from_buffer(RAW, unpack_data=True))))), + "take_np_datasamples", +): + del CASES["take_np"] + +for fn in CASES.values(): # warm up every case before timing any + fn(RAW) + gc.collect() + +times = {k: [] for k in CASES} +order = list(CASES) +rng = random.Random(20260816) +for _ in range(REPS): + rng.shuffle(order) + for name in order: + gc.collect() + t0 = time.perf_counter() + CASES[name](RAW) + times[name].append((time.perf_counter() - t0) * 1e3) + +out = {"pymseed": pymseed.__version__, "reps": REPS, "cases": {}} +print(f"pymseed {pymseed.__version__} reps={REPS}") +print(f"{'case':<14}{'min':>8}{'p25':>8}{'median':>8}") +for k, v in times.items(): + v.sort() + rec = {"min": v[0], "p25": v[len(v) // 4], "median": st.median(v)} + out["cases"][k] = rec + print(f"{k:<14}{rec['min']:>8.1f}{rec['p25']:>8.1f}{rec['median']:>8.1f}") +print("JSON " + json.dumps(out)) diff --git a/docs/pymseed-issue-draft.md b/docs/pymseed-issue-draft.md index 02e3af6..2a01a24 100644 --- a/docs/pymseed-issue-draft.md +++ b/docs/pymseed-issue-draft.md @@ -42,15 +42,46 @@ Findings: — the two are within noise of each other, with the record-list path marginally ahead. Chad's "huge win" framing is against the *copy* idiom and the record-list *setup cost*; on this file the record-list build is only - ~4.5 ms and it is amortized by giving us what we need anyway. + ~4.5 ms and it is amortized by giving us what we need anyway. The container + tiers below widen this gap rather than closing it. - **seisfetch stays on the record-list path.** Not for speed but for metadata: we want the per-record encoding, which `record_list=True` provides in the same pass. Reaching the same place via `take_np` would mean `unpack_data=True, record_list=True` (27.8 ms — slower than either) or giving up the encoding. -- **Not yet tested under cgroup limits**, which is the tier the original - issue was about. Docker was unavailable for this pass; the native numbers - are the weakest place to judge a memory-traffic change. +Under cgroup limits (`benchmarks/ab_pymseed_paths.py`, the four cases +interleaved in shuffled order, 2 containers x 15 reps per version per tier, +same channel-day; Docker Desktop linux/arm64 on M1). Minimum ms — medians are +2-3x higher at the Lambda tier from CPU-quota throttling and are reported in +the script output: + +| case | fargate 2cpu/4g, 0.9.4 | 0.9.5 | lambda 0.5cpu/512m, 0.9.4 | 0.9.5 | +|---|---|---|---|---| +| `record_list` → `create_numpy_array_from_recordlist()` | **21.7** | **21.6** | **21.5** | **21.2** | +| `unpack_data=True` → `np_datasamples` (safe view) | 28.1 | 27.7 | 28.1 | 28.3 | +| `unpack_data=True` → `take_np_datasamples()` | n/a | 27.6 | n/a | 28.6 | +| `unpack_data=True` → `np_datasamples.copy()` | 35.0 | 35.2 | 36.4 | 47.6 | + +This sharpens the native result rather than overturning it: + +- **`take_np_datasamples()` is indistinguishable from the plain 0.9.4 safe + view** — 27.6 vs 27.7 ms at the Fargate tier, 28.6 vs 28.3 at Lambda. + Expected in hindsight: since the 0.9.4 keepalive the view is already + zero-copy, so `take_np` changes *ownership*, not data movement. Its value + is lifetime semantics — the array no longer pins the trace list, so the + rest of that memory can be released earlier — not throughput. +- **The record-list path is the fastest in every tier**, and by a wider + margin than native: ~6-7 ms (about 25 %) ahead of both unpack-then-take + paths at both tiers, where natively it led by ~1 ms. +- The copy idiom remains the worst path and is the only one that degrades + with tighter limits, reproducing the original issue's finding. + +Caveats on these numbers: one file with a single trace/segment, so the +record-list advantage may not generalize to many-segment or many-channel +objects; and Docker Desktop on macOS runs a Linux VM, so this is +cgroup-in-VM rather than bare Linux — the same property the original issue's +measurements had, which makes them comparable to each other but not a +substitute for real Fargate/Lambda. Release safety, checked separately: `parse_mseed` output is **bit-identical between 0.9.4 and 0.9.5** across every fixture (Steim2 day, float32, float64, From d4a48a58e625a571b49cf841d2436b6b07d7a0fa Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Sun, 16 Aug 2026 18:53:28 -0700 Subject: [PATCH 5/6] A/B the decode paths on real station-day objects, not just a channel-day MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The single-segment 11 MB channel-day is the friendliest case for the unpack paths, so it could not answer whether the record-list lead holds on what EarthScope actually serves. ab_pymseed_paths.py now takes a glob and times the whole set per round. Re-run over II.PFO station-days (44 traceids / ~50 segments each): the record-list lead grows from ~1.3x to ~2.1x — 216 vs 447 ms on 194 segments where it was 21.6 vs 27.7 on one. That fits the mechanism, since mstl3_unpack_recordlist decodes into one numpy-owned allocation per segment while unpack_data=True wraps libmseed's buffers and pays per-segment overhead. So the path seisfetch already uses wins by the widest margin exactly where our heaviest objects live. One counterpoint recorded honestly: take_np did beat the plain view at lambda-1g with 97 segments (341.6 vs 384.2), consistent across min/p25/median but a single run, so noted as suggestive only. Co-Authored-By: Claude Opus 5 (1M context) --- benchmarks/ab_pymseed_paths.py | 48 +++++++++++++++++++++++++++++---- docs/pymseed-issue-draft.md | 49 +++++++++++++++++++++++++++++----- 2 files changed, 86 insertions(+), 11 deletions(-) diff --git a/benchmarks/ab_pymseed_paths.py b/benchmarks/ab_pymseed_paths.py index f943693..cbc6a73 100644 --- a/benchmarks/ab_pymseed_paths.py +++ b/benchmarks/ab_pymseed_paths.py @@ -3,6 +3,14 @@ Only the four cases that matter, interleaved in a shuffled order each round so host drift and CPU-quota refill cancel instead of biasing whichever case runs first. + +Usage: + python benchmarks/ab_pymseed_paths.py [REPS] [GLOB] + +With no GLOB the single-segment 11 MB Steim2 channel-day in tests/ is used. +Pass a glob (e.g. '/data/earthscope.*.ms') to time the whole set per round +instead — station-day objects carry tens of traceids and segments, and the +per-segment cost structure differs from one big contiguous segment. """ import gc @@ -16,8 +24,17 @@ import pymseed from pymseed import MS3TraceList -RAW = (Path(__file__).resolve().parent.parent / "tests" / "bench.mseed").read_bytes() REPS = int(sys.argv[1]) if len(sys.argv) > 1 else 15 +if len(sys.argv) > 2: + import glob as _glob + + PATHS = sorted(Path(p) for p in _glob.glob(sys.argv[2])) + if not PATHS: + sys.exit(f"no files matched {sys.argv[2]!r}") +else: + PATHS = [Path(__file__).resolve().parent.parent / "tests" / "bench.mseed"] +BUFFERS = [p.read_bytes() for p in PATHS] # read once; this measures decode +RAW = BUFFERS[0] def reclist_np(raw): @@ -49,8 +66,14 @@ def copy(raw): ): del CASES["take_np"] + +def run_case(fn): + """One timed round = the whole file set, so many-segment objects count.""" + return sum(fn(b) for b in BUFFERS) + + for fn in CASES.values(): # warm up every case before timing any - fn(RAW) + run_case(fn) gc.collect() times = {k: [] for k in CASES} @@ -61,11 +84,26 @@ def copy(raw): for name in order: gc.collect() t0 = time.perf_counter() - CASES[name](RAW) + run_case(CASES[name]) times[name].append((time.perf_counter() - t0) * 1e3) -out = {"pymseed": pymseed.__version__, "reps": REPS, "cases": {}} -print(f"pymseed {pymseed.__version__} reps={REPS}") +nseg = sum( + len(seg_list) + for b in BUFFERS + for seg_list in [[s for t in MS3TraceList.from_buffer(b) for s in t]] +) +out = { + "pymseed": pymseed.__version__, + "reps": REPS, + "files": len(BUFFERS), + "mb": round(sum(len(b) for b in BUFFERS) / 1e6, 1), + "segments": nseg, + "cases": {}, +} +print( + f"pymseed {pymseed.__version__} reps={REPS} files={len(BUFFERS)} " + f"({out['mb']} MB, {nseg} segments)" +) print(f"{'case':<14}{'min':>8}{'p25':>8}{'median':>8}") for k, v in times.items(): v.sort() diff --git a/docs/pymseed-issue-draft.md b/docs/pymseed-issue-draft.md index 2a01a24..5de5950 100644 --- a/docs/pymseed-issue-draft.md +++ b/docs/pymseed-issue-draft.md @@ -76,12 +76,49 @@ This sharpens the native result rather than overturning it: - The copy idiom remains the worst path and is the only one that degrades with tighter limits, reproducing the original issue's finding. -Caveats on these numbers: one file with a single trace/segment, so the -record-list advantage may not generalize to many-segment or many-channel -objects; and Docker Desktop on macOS runs a Linux VM, so this is -cgroup-in-VM rather than bare Linux — the same property the original issue's -measurements had, which makes them comparable to each other but not a -substitute for real Fargate/Lambda. +### Does it change with realistic station-day objects? Yes — it doubles + +The tests above use one 11 MB channel-day that merges to a single +trace/segment, which is the friendliest possible case for the unpack paths. +Re-run over real EarthScope station-day objects from the xcorr cache +(`earthscope.II.PFO.all.*`, 44 traceids / ~50 segments each), which is what +the archive actually serves: + +| case | fargate 2cpu/4g, 4 files (106 MB, 194 seg) | lambda-1g, 2 files (52 MB, 97 seg) | +|---|---|---| +| | 0.9.4 / 0.9.5 | 0.9.4 / 0.9.5 | +| `record_list` → `create_numpy_array_from_recordlist()` | **216.9 / 215.7** | **155.9 / 151.6** | +| `unpack_data=True` → `np_datasamples` (safe view) | 451.2 / 447.4 | 391.4 / 384.2 | +| `unpack_data=True` → `take_np_datasamples()` | n/a / 450.6 | n/a / 341.6 | +| `unpack_data=True` → `np_datasamples.copy()` | 483.1 / 481.1 | 407.6 / 409.3 | + +Minimum ms per round, where a round decodes every file in the set. The +fargate column is 2 containers x 7 reps per version; the lambda-1g column is +a single container x 7 reps. + +- **The record-list lead grows from ~1.3x to ~2.1x.** On the single-segment + file it was 21.6 vs 27.7 ms; on 194-segment station-days it is 216 vs 447. + The advantage scales with segment count, which fits the mechanism: + `mstl3_unpack_recordlist` decodes straight into one numpy-owned allocation + per segment, while `unpack_data=True` decodes into libmseed's own buffers + and then wraps each segment, paying per-segment overhead that the + record-list path never incurs. +- **This is the case that matters for us.** EarthScope objects are + station-days with tens of channels; the single-segment channel-day is the + SCEDC/NCEDC shape. So the path seisfetch already uses wins by the largest + margin exactly where our heaviest objects live. +- **One place `take_np` does appear to help**: at lambda-1g with 97 segments + it came in at 341.6 ms against 384.2 for the plain view, consistent across + min/p25/median. Plausible mechanism — releasing the trace list earlier + relieves memory pressure when many segments are live at once. Single run, + so treat as suggestive rather than established. It is still ~2.2x slower + than the record-list path. + +Caveats on all of the above: Docker Desktop on macOS runs a Linux VM, so this +is cgroup-in-VM rather than bare Linux — the same property the original +issue's measurements had, which makes them comparable to each other but not a +substitute for real Fargate/Lambda. Files are read into memory before timing, +so these measure decode only, not S3 transfer. Release safety, checked separately: `parse_mseed` output is **bit-identical between 0.9.4 and 0.9.5** across every fixture (Steim2 day, float32, float64, From f0ddecf84e23af8630add8f611939dd306696cfe Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Sun, 16 Aug 2026 19:00:57 -0700 Subject: [PATCH 6/6] Address Copilot review on PR #6 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit __all__ listed the four Earth2Studio adapter names unconditionally while they are only bound when the optional seisfetch.earth2 import succeeds, so `from seisfetch import *` raised AttributeError on a minimal install. Reproduced by forcing that import to fail — three of the four names predate this branch, and adding SeisfetchLiveSource widened it. __all__ now gains those names only when the import succeeded, which also keeps __all__ an honest description of what the install actually provides. The import needs an explicit noqa: F401 now that ruff can no longer see the re-export through a dynamic __all__. profile_parse's take_np case called take_np_datasamples() unguarded. run_table catches per-case failures, but run_rss does not, so `--rss tl_unpack+take_np` on an older pymseed printed a bare AttributeError rather than the version reason. Now raises with the requirement and the installed version named; table mode still degrades to one FAILED row, with a better message than before. Co-Authored-By: Claude Opus 5 (1M context) --- benchmarks/profile_parse.py | 11 +++++++++++ seisfetch/__init__.py | 19 +++++++++++++------ 2 files changed, 24 insertions(+), 6 deletions(-) diff --git a/benchmarks/profile_parse.py b/benchmarks/profile_parse.py index 18ed82f..ea5e9ae 100644 --- a/benchmarks/profile_parse.py +++ b/benchmarks/profile_parse.py @@ -109,7 +109,18 @@ def case_tracelist_take_np(raw): copy, so the array outlives the trace list without the 0.9.4 keepalive pinning it. EarthScope/pymseed#6. """ + import pymseed from pymseed import MS3TraceList + from pymseed.mstracelist import MS3TraceSeg + + # run_table catches per-case exceptions, but run_rss does not — so + # `--rss tl_unpack+take_np` on an older pymseed would otherwise print an + # unexplained AttributeError instead of the version reason. + if not hasattr(MS3TraceSeg, "take_np_datasamples"): + raise RuntimeError( + f"take_np_datasamples() requires pymseed >= 0.9.5 " + f"(installed: {pymseed.__version__})" + ) tl = MS3TraceList.from_buffer(raw, unpack_data=True) n = 0 diff --git a/seisfetch/__init__.py b/seisfetch/__init__.py index 4b0d431..e2d53e0 100644 --- a/seisfetch/__init__.py +++ b/seisfetch/__init__.py @@ -74,15 +74,22 @@ def __dir__(): # Earth2Studio adapters — lazy import (requires earth2studio + xarray) +_EARTH2_EXPORTS = [ + "SeismicDataSource", + "SeismicDataFrameSource", + "SeisfetchLiveSource", + "bundle_to_earth2", +] try: - from seisfetch.earth2 import ( + # re-exported via the dynamic __all__ below, which ruff cannot follow + from seisfetch.earth2 import ( # noqa: F401 SeisfetchLiveSource, SeismicDataFrameSource, SeismicDataSource, bundle_to_earth2, ) except ImportError: # earth2studio / xarray not installed - pass + _EARTH2_EXPORTS = [] try: from importlib.metadata import version as _pkg_version @@ -119,8 +126,8 @@ def __dir__(): "fetch_bulk_numpy", "requests_from_list", "requests_from_csv", - "SeismicDataSource", - "SeismicDataFrameSource", - "SeisfetchLiveSource", - "bundle_to_earth2", ] +# Earth2Studio adapter names only when the optional import succeeded: they are +# bound eagerly above, and a name in __all__ that is not resolvable makes +# `from seisfetch import *` raise AttributeError on a minimal install. +__all__ += _EARTH2_EXPORTS