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..4c3b21d 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,88 @@ 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. 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 ### Open S3: SCEDC / NCEDC @@ -134,6 +222,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 +525,7 @@ SeisfetchClient +- backend="s3_open" | +- SCEDC open bucket | +- NCEDC open bucket +| +- GeoNet open bucket | +- auto-routing by network code | +- backend="s3_auth" @@ -438,8 +543,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 +561,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 +707,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 +748,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/benchmarks/ab_pymseed_paths.py b/benchmarks/ab_pymseed_paths.py new file mode 100644 index 0000000..cbc6a73 --- /dev/null +++ b/benchmarks/ab_pymseed_paths.py @@ -0,0 +1,113 @@ +"""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. + +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 +import json +import random +import statistics as st +import sys +import time +from pathlib import Path + +import pymseed +from pymseed import MS3TraceList + +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): + 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"] + + +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 + run_case(fn) + 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() + run_case(CASES[name]) + times[name].append((time.perf_counter() - t0) * 1e3) + +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() + 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/benchmarks/profile_parse.py b/benchmarks/profile_parse.py index 7da1766..ea5e9ae 100644 --- a/benchmarks/profile_parse.py +++ b/benchmarks/profile_parse.py @@ -89,6 +89,48 @@ 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. + """ + 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 + 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 +191,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..5de5950 100644 --- a/docs/pymseed-issue-draft.md +++ b/docs/pymseed-issue-draft.md @@ -10,6 +10,125 @@ 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. 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. +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. + +### 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, +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 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..e2d53e0 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() @@ -73,21 +74,29 @@ 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 __version__ = _pkg_version("seisfetch") except Exception: # not installed (e.g. vendored copy) - __version__ = "0.3.0" + __version__ = "0.4.0" __all__ = [ "SeisfetchClient", "S3OpenClient", @@ -117,7 +126,8 @@ def __dir__(): "fetch_bulk_numpy", "requests_from_list", "requests_from_csv", - "SeismicDataSource", - "SeismicDataFrameSource", - "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