From dde2d0250c8d29a47815e47eed8570f4f614f8ba Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Fri, 7 Aug 2026 07:52:28 +0200 Subject: [PATCH 1/3] Add GeoNet (New Zealand) open-data S3 archive waveforms/miniseed/{Y}/{Y}.{DDD}/{STA}.{NET}/{Y}.{DDD}.{STA}.{LOC}-{CHA}.{NET}.D per-channel day objects in geonet-open-data (ap-southeast-2), verified live: exact-key fetch (NZ.WEL.10.HHZ, 13.2 MB), HH? wildcard discovery (three components), and station listing. NZ routes to geonet. Also corrects the NCEDC region to us-west-2 (was us-east-2; worked via S3 redirect but paid a cross-region hop on every request). Co-Authored-By: Claude Fable 5 --- seisfetch/s3.py | 70 ++++++++++++++++++++++++++++++++++++++++++------ tests/test_s3.py | 10 +++++++ 2 files changed, 72 insertions(+), 8 deletions(-) diff --git a/seisfetch/s3.py b/seisfetch/s3.py index ee2ae71..30c096f 100644 --- a/seisfetch/s3.py +++ b/seisfetch/s3.py @@ -1,7 +1,7 @@ """ S3-based backends for seismic miniSEED data. -Supports three open-data archives with different path conventions: +Supports four open-data archives with different path conventions: EarthScope s3://earthscope-geophysical-data (us-east-2) miniseed/{NET}/{YEAR}/{DOY}/{STA}.{NET}.{YEAR}.{DOY} @@ -11,14 +11,19 @@ continuous_waveforms/{YEAR}/{YEAR}_{DOY}/{NET}{STA}{LOC}{CHA}__{YEAR}{DOY}.ms One object per channel-day. - NCEDC s3://ncedc-pds (us-east-2) + NCEDC s3://ncedc-pds (us-west-2) continuous_waveforms/{NET}/{YEAR}/{YEAR}.{DOY}/{STA}.{NET}.{CHA}.{LOC}.D.{YEAR}.{DOY} One object per channel-day. + GeoNet s3://geonet-open-data (ap-southeast-2) + waveforms/miniseed/{YEAR}/{YEAR}.{DOY}/{STA}.{NET}/{YEAR}.{DOY}.{STA}.{LOC}-{CHA}.{NET}.D + One object per channel-day (New Zealand; NZ network). + The :class:`S3Router` auto-selects the right datacenter by network code. Attribution: SCEDC — Yu et al. (2021), doi:10.7909/C3WD3xH1 + GeoNet — https://www.geonet.org.nz/data/supplementary/channels (CC BY 4.0) NCEDC — doi:10.7932/NCEDC EarthScope — https://www.earthscope.org/how-to-cite/ NoisePy S3 store pattern — Jiang & Denolle (2020), doi:10.1785/0220190364 @@ -83,6 +88,21 @@ def _ncedc_key(network, station, year, doy, location="", channel="", **_): ) +def _geonet_key(network, station, year, doy, location="", channel="", **_): + """GeoNet (New Zealand): one object per channel-day. + + Layout (verified on the live bucket, 2026-08-07): + ``waveforms/miniseed/{Y}/{Y}.{DDD}/{STA}.{NET}/{Y}.{DDD}.{STA}.{LOC}-{CHA}.{NET}.D`` + e.g. ``waveforms/miniseed/2022/2022.002/WEL.NZ/2022.002.WEL.10-HHZ.NZ.D``. + GeoNet channels always carry a numeric location code (10, 20, ...). + """ + loc = location if location and location != "*" else "" + return ( + f"waveforms/miniseed/{year}/{year}.{doy:03d}/{station}.{network}/" + f"{year}.{doy:03d}.{station}.{loc}-{channel}.{network}.D" + ) + + # =========================================================================== # # Datacenter configs # =========================================================================== # @@ -103,10 +123,17 @@ def _ncedc_key(network, station, year, doy, location="", channel="", **_): }, "ncedc": { "bucket": "ncedc-pds", - "region": "us-east-2", + "region": "us-west-2", # was us-east-2: worked via redirect, but + # ncedc-pds lives in us-west-2 — direct addressing avoids the hop "key_fn": _ncedc_key, "per_channel": True, }, + "geonet": { + "bucket": "geonet-open-data", + "region": "ap-southeast-2", + "key_fn": _geonet_key, + "per_channel": True, + }, } # Network → datacenter routing (following quakescope/noisepy pattern) @@ -153,8 +180,9 @@ def route_network(network: str) -> str: """ Auto-select datacenter for a given network code. - Returns ``"scedc"``, ``"ncedc"``, or ``"earthscope"``. - SCEDC is preferred for CI; NCEDC for BK/NC; EarthScope for everything else. + Returns ``"scedc"``, ``"ncedc"``, ``"geonet"``, or ``"earthscope"``. + SCEDC is preferred for CI; NCEDC for BK/NC; GeoNet for NZ; EarthScope + for everything else. """ net = network.upper() if net == "CI" or net in _SCEDC_NETS - _NCEDC_NETS: @@ -163,6 +191,8 @@ def route_network(network: str) -> str: return "ncedc" if net in _SCEDC_NETS & _NCEDC_NETS: return "ncedc" # prefer NCEDC for shared nets (NC, NP, etc.) + if net == "NZ": + return "geonet" return "earthscope" @@ -173,13 +203,14 @@ def route_network(network: str) -> str: class S3OpenClient: """ - Anonymous S3 access to EarthScope, SCEDC, and NCEDC open-data buckets. + Anonymous S3 access to the EarthScope, SCEDC, NCEDC, and GeoNet + open-data buckets. Parameters ---------- datacenter : str or None - ``"earthscope"``, ``"scedc"``, ``"ncedc"``, or ``None`` (auto-route - by network code, default). + ``"earthscope"``, ``"scedc"``, ``"ncedc"``, ``"geonet"``, or + ``None`` (auto-route by network code, default). max_workers : int Thread pool for parallel downloads. """ @@ -300,6 +331,22 @@ def _discover_channel_keys( if location != "*" and loc != (location or ""): continue keys.append(key) + elif dc_name == "geonet": + prefix = ( + f"waveforms/miniseed/{yr}/{yr}.{doy:03d}/" + f"{station}.{network}/{yr}.{doy:03d}.{station}." + ) + for key in self._iter_keys(s3, dc["bucket"], prefix): + # {Y}.{DDD}.{STA}.{LOC}-{CHA}.{NET}.D + parts = key.rsplit("/", 1)[-1].split(".") + if len(parts) < 4 or "-" not in parts[3]: + continue + loc, cha = parts[3].split("-", 1) + if not fnmatch.fnmatch(cha, channel): + continue + if location != "*" and loc != (location or ""): + continue + keys.append(key) else: # ncedc prefix = ( f"continuous_waveforms/{network}/{yr}/{yr}.{doy:03d}/" @@ -482,6 +529,8 @@ def list_stations(self, network, year, doy, datacenter=None): prefix = f"miniseed/{network}/{year}/{doy:03d}/" elif dc_name == "scedc": prefix = f"continuous_waveforms/{year}/{year}_{doy:03d}/{network}" + elif dc_name == "geonet": + prefix = f"waveforms/miniseed/{year}/{year}.{doy:03d}/" else: prefix = f"continuous_waveforms/{network}/{year}/{year}.{doy:03d}/" stations = set() @@ -498,6 +547,11 @@ def list_stations(self, network, year, doy, datacenter=None): else fname[:5] ) stations.add(sta.rstrip("_")) + elif dc_name == "geonet": + # .../{STA}.{NET}/{Y}.{DDD}.{STA}.{LOC}-{CHA}.{NET}.D + stadir = obj["Key"].rsplit("/", 2)[-2] + if stadir.endswith(f".{network}"): + stations.add(stadir.rsplit(".", 1)[0]) else: # ncedc stations.add(fname.split(".")[0]) return sorted(stations) diff --git a/tests/test_s3.py b/tests/test_s3.py index 5e03c7d..519e759 100644 --- a/tests/test_s3.py +++ b/tests/test_s3.py @@ -8,6 +8,7 @@ from seisfetch.s3 import ( S3OpenClient, _earthscope_key, + _geonet_key, _ncedc_key, _scedc_key, route_network, @@ -35,6 +36,12 @@ def test_ncedc_no_location(self): k = _ncedc_key("NC", "JBL", 2023, 1, location="", channel="HHZ") assert "JBL.NC.HHZ..D.2023.001" in k + def test_geonet(self): + # verified against the live bucket 2026-08-07: + # waveforms/miniseed/2022/2022.002/WEL.NZ/2022.002.WEL.10-HHZ.NZ.D + k = _geonet_key("NZ", "WEL", 2022, 2, location="10", channel="HHZ") + assert k == ("waveforms/miniseed/2022/2022.002/WEL.NZ/2022.002.WEL.10-HHZ.NZ.D") + # ── Network routing ────────────────────────────────────────────────── # @@ -59,6 +66,9 @@ def test_nc_shared_prefers_ncedc(self): # NC is in both SCEDC and NCEDC; should prefer NCEDC assert route_network("NC") == "ncedc" + def test_nz_to_geonet(self): + assert route_network("NZ") == "geonet" + def test_case_insensitive(self): assert route_network("ci") == "scedc" From 018ecc84acf7c85045d5ec339104fd5b445f540c Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Fri, 7 Aug 2026 07:58:12 +0200 Subject: [PATCH 2/3] Live Earth2Studio source with mandatory physical calibration MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit SeisfetchLiveSource is time-indexed like every other Earth2Studio source (GFS, ERA5, ...): call it with timestamps and it fetches day objects on demand from EarthScope/SCEDC/NCEDC/GeoNet with a small day cache, trims sample-precise windows, and always returns physical units — gain correction by default (channel Scale from the FDSN station text service, one stdlib-parsed HTTP request via new ChannelEpoch/parse_channel_text), or full response deconvolution via contrib.response. Raw counts are never returned. The existing adapters' async fetch is now real (asyncio.to_thread) and SeismicDataFrameSource can auto-fill station coordinates. No band/RMS feature extraction by design — derived low-rate feature series are a later, separate interface. Verified live on GeoNet: NZ.WEL.10.HHZ window in m/s with real Wellington coordinates; unit tests cover parsing, epoch selection, gain mechanics, NaN gap fill, cache, and async. Co-Authored-By: Claude Fable 5 --- CHANGELOG.md | 25 ++++ seisfetch/earth2.py | 281 +++++++++++++++++++++++++++++++++++++- seisfetch/fdsn.py | 74 ++++++++++ tests/test_earth2_live.py | 137 +++++++++++++++++++ 4 files changed, 513 insertions(+), 4 deletions(-) create mode 100644 tests/test_earth2_live.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 787a751..7930446 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,31 @@ All notable changes to seisfetch are documented here. The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/); versions follow [Semantic Versioning](https://semver.org). +## Unreleased + +### Added + +- **GeoNet (New Zealand)** open-data S3 archive: `geonet-open-data` + (ap-southeast-2), per-channel day objects, wildcard discovery and + station listing; `NZ` routes there automatically. Verified live. +- **`SeisfetchLiveSource`** (`seisfetch.earth2`): a time-indexed + Earth2Studio `DataSource` that fetches on demand from the four archives + with a day-bundle cache — physical units are required (`calibrate="gain"` + by default via the FDSN station text service; `"response"` for full + deconvolution via `contrib.response`; raw counts are never returned). +- `ChannelEpoch` + `parse_channel_text` (`seisfetch.fdsn`): stdlib parsing + of `format=text&level=channel` — gain, coordinates and epochs in one + HTTP request, no StationXML. +- `SeismicDataFrameSource(auto_coords=True)` fills station lat/lon from + the FDSN station service automatically. + +### Changed + +- Earth2Studio adapters' `fetch` is now genuinely async + (`asyncio.to_thread`), so pipelines can overlap sources. +- NCEDC S3 region corrected to us-west-2 (was us-east-2; worked via + redirect but paid a cross-region hop per request). + ## 0.3.1 — 2026-08-06 Cross-station validation release: the noisepy equivalence evidence now diff --git a/seisfetch/earth2.py b/seisfetch/earth2.py index a6ec7b4..7007571 100644 --- a/seisfetch/earth2.py +++ b/seisfetch/earth2.py @@ -159,8 +159,11 @@ async def fetch( time: datetime | list[datetime] | np.ndarray, variable: str | list[str] | np.ndarray, ): - """Async fetch — just delegates to __call__.""" - return self.__call__(time, variable) + """True async fetch: runs the blocking path in a worker thread so + Earth2Studio pipelines can overlap sources.""" + import asyncio + + return await asyncio.to_thread(self.__call__, time, variable) # --------------------------------------------------------------------------- # @@ -195,6 +198,7 @@ def __init__( self, bundle_or_dataset: Any, station_coords: dict[str, tuple[float, float]] | None = None, + auto_coords: bool = False, ): try: import pandas as pd # noqa: F401 @@ -220,6 +224,31 @@ def __init__( ) self._station_coords = station_coords or {} + if auto_coords and self._bundle is not None: + # one FDSN station-text request per net.sta — fills lat/lon + # without hand-building a coordinate table (stdlib parse) + seen = set() + for tr in self._bundle.traces: + key = f"{tr.network}.{tr.station}" + if key in seen or key in self._station_coords: + continue + seen.add(key) + try: + eps = channel_metadata( + tr.network, + tr.station, + "*", + "*", + "1970-01-01", + "2100-01-01", + ) + if eps: + self._station_coords[key] = ( + eps[0].latitude, + eps[0].longitude, + ) + except Exception as e: # metadata is best-effort here + logger.warning("auto_coords failed for %s: %s", key, e) self._init_schema() def _init_schema(self): @@ -314,8 +343,10 @@ async def fetch( variable: str | list[str] | np.ndarray, fields: Any = None, ): - """Async fetch — delegates to __call__.""" - return self.__call__(time, variable, fields) + """True async fetch: runs the blocking path in a worker thread.""" + import asyncio + + return await asyncio.to_thread(self.__call__, time, variable, fields) # --------------------------------------------------------------------------- # @@ -340,3 +371,245 @@ def bundle_to_earth2(bundle, variables: list[str] | None = None): """ src = SeismicDataSource(bundle) return src + + +# --------------------------------------------------------------------------- # +# Live, time-indexed DataSource (fetch-on-call, physically calibrated) +# --------------------------------------------------------------------------- # + +_DC_TO_PROVIDER = { + "earthscope": "EARTHSCOPE", + "scedc": "SCEDC", + "ncedc": "NCEDC", + "geonet": "GEONET", +} + + +def channel_metadata(network, station, location, channel, start, end, provider=None): + """Channel epochs (gain + coordinates) via the FDSN station text service. + + One HTTP request, stdlib parsing — no StationXML, no ObsPy, no scipy. + ``provider`` defaults to the FDSN service matching the archive that + :func:`seisfetch.s3.route_network` selects for ``network``. + """ + from seisfetch.fdsn import FDSNClient, parse_channel_text + from seisfetch.s3 import route_network + + prov = provider or _DC_TO_PROVIDER[route_network(network)] + text = FDSNClient(provider=prov).get_station_text( + network=network, + station=station, + location=location if location else "--", + channel=channel, + starttime=start, + endtime=end, + level="channel", + format="text", + ) + return parse_channel_text(text) + + +class SeisfetchLiveSource: + """Time-indexed Earth2Studio ``DataSource`` backed by the cloud archives. + + Unlike :class:`SeismicDataSource` (which wraps data you already fetched), + this source is called with timestamps and fetches on demand — the shape + every other Earth2Studio source (GFS, ERA5, ...) has. Day objects are + fetched via :class:`seisfetch.s3.S3OpenClient` (auto-routed per network: + EarthScope, SCEDC, NCEDC, GeoNet) and cached; windows are trimmed + sample-precisely from the cached day bundles. + + Physical units are **required**: raw counts are never returned. + + * ``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 via + :mod:`seisfetch.contrib.response` (StationXML fetch + evalresp- + equivalent evaluator). + + Parameters + ---------- + channels : list[str] + Channel IDs as ``"NET.STA.LOC.CHA"`` (blank location: ``NET.STA..CHA``). + window_s : float + Length of the window returned per requested timestamp (default 3600). + calibrate : str + ``"gain"`` or ``"response"`` (see above). + datacenter : str, optional + Force one archive instead of per-network auto-routing. + cache_days : int + Day bundles kept in memory per channel (default 3). + """ + + def __init__( + self, + channels, + window_s: float = 3600.0, + calibrate: str = "gain", + datacenter: str | None = None, + cache_days: int = 3, + ): + if calibrate not in ("gain", "response"): + raise ValueError( + "calibrate must be 'gain' or 'response' — this source always " + "returns physical units, never raw counts" + ) + self.channels = list(channels) + self.window_s = float(window_s) + self.calibrate = calibrate + self._datacenter = datacenter + self._cache_days = cache_days + self._day_cache: dict = {} # (nslc, date) -> TraceBundle + self._meta: dict = {} # nslc -> list[ChannelEpoch] + self._resp: dict = {} # nslc -> ChannelResponse (calibrate="response") + + # -- metadata ---------------------------------------------------------- # + + def _epochs(self, nslc: str): + if nslc not in self._meta: + net, sta, loc, cha = nslc.split(".") + self._meta[nslc] = channel_metadata( + net, sta, loc, cha, "1970-01-01", "2100-01-01" + ) + if not self._meta[nslc]: + raise LookupError(f"no channel metadata for {nslc}") + return self._meta[nslc] + + def _gain(self, nslc: str, time_iso: str) -> float: + for ep in self._epochs(nslc): + if ep.covers(time_iso): + if not ep.scale: + raise LookupError( + f"{nslc}: channel epoch has no Scale (sensitivity) — " + "cannot calibrate" + ) + return ep.scale + raise LookupError(f"{nslc}: no channel epoch covers {time_iso}") + + def coords(self, nslc: str): + ep = self._epochs(nslc)[0] + return ep.latitude, ep.longitude, ep.elevation + + # -- data -------------------------------------------------------------- # + + def _day_bundle(self, nslc: str, day): + from seisfetch.convert import parse_mseed + from seisfetch.s3 import S3OpenClient + + key = (nslc, day.isoformat()) + if key not in self._day_cache: + net, sta, loc, cha = nslc.split(".") + raw = S3OpenClient(datacenter=self._datacenter).get_raw( + net, sta, day.isoformat(), location=loc, channel=cha + ) + self._day_cache[key] = parse_mseed(raw) + while len(self._day_cache) > self._cache_days * len(self.channels): + self._day_cache.pop(next(iter(self._day_cache))) + return self._day_cache[key] + + def _window(self, nslc: str, t0: datetime) -> np.ndarray: + """One calibrated window [t0, t0+window_s), NaN-padded over gaps.""" + from datetime import timedelta + + start_ns = int(t0.timestamp() * 1e9) + end_ns = start_ns + int(self.window_s * 1e9) + days = {t0.date(), (t0 + timedelta(seconds=self.window_s)).date()} + segs = [] + for day in sorted(days): + try: + cut = self._day_bundle(nslc, day).trim(start_ns, end_ns) + except Exception: + continue + segs.extend(s for s in cut.traces if s.id == nslc) + if not segs: + raise LookupError(f"{nslc}: no data in [{t0}, +{self.window_s}s)") + fs = segs[0].sampling_rate + n = int(round(self.window_s * fs)) + out = np.full(n, np.nan) + for s in segs: + i0 = int(round((s.starttime_ns - start_ns) * fs / 1e9)) + src = np.asarray(s.data, dtype=np.float64) + j0, j1 = max(i0, 0), min(i0 + s.npts, n) + if j1 > j0: + out[j0:j1] = src[j0 - i0 : j1 - i0] + t_iso = t0.isoformat() + if self.calibrate == "gain": + return out / self._gain(nslc, t_iso) + return self._deconvolve(nslc, out, fs, t_iso) + + def _deconvolve(self, nslc, x, fs, time_iso): + from seisfetch.contrib.response import ( + parse_stationxml_response, + remove_response_np, + ) + from seisfetch.fdsn import FDSNClient + from seisfetch.s3 import route_network + + if nslc not in self._resp: + net, sta, loc, cha = nslc.split(".") + xml = FDSNClient( + provider=_DC_TO_PROVIDER[route_network(net)] + ).get_station_text( + network=net, + station=sta, + location=loc if loc else "--", + channel=cha, + level="response", + format="xml", + ) + self._resp[nslc] = parse_stationxml_response( + xml.encode(), net, sta, loc, cha, time_iso + ) + mask = np.isnan(x) + filled = np.where(mask, 0.0, x) + v = remove_response_np(filled, fs, self._resp[nslc], output="VEL") + v[mask] = np.nan + return v + + # -- Earth2Studio protocol --------------------------------------------- # + + def __call__(self, time, variable=None): + import pandas as pd + import xarray as xr + + times = [time] if isinstance(time, datetime) else list(time) + variables = ( + self.channels + if variable is None + else [v.replace("_", ".") for v in np.atleast_1d(variable)] + ) + data, nsamp = [], None + for t0 in times: + row = [] + for nslc in variables: + w = self._window(nslc, t0) + nsamp = len(w) if nsamp is None else nsamp + if len(w) != nsamp: # mixed sampling rates across channels + raise ValueError( + "channels have different sampling rates; request " + "them in separate calls" + ) + row.append(w) + data.append(row) + return xr.DataArray( + np.asarray(data), + dims=["time", "variable", "sample"], + coords={ + "time": pd.to_datetime(times), + "variable": [v.replace(".", "_") for v in variables], + "sample": np.arange(nsamp), + }, + attrs={ + "units": "m/s (gain-corrected)" + if self.calibrate == "gain" + else "m/s (response-removed)", + "calibration": self.calibrate, + }, + ) + + async def fetch(self, time, variable=None): + """True async fetch: runs the blocking pipeline in a worker thread.""" + import asyncio + + return await asyncio.to_thread(self.__call__, time, variable) diff --git a/seisfetch/fdsn.py b/seisfetch/fdsn.py index 82dd43f..c75cd89 100644 --- a/seisfetch/fdsn.py +++ b/seisfetch/fdsn.py @@ -13,6 +13,7 @@ import logging import time from concurrent.futures import ThreadPoolExecutor +from dataclasses import dataclass from seisfetch.exceptions import FDSNError, FetchError from seisfetch.utils import to_epoch, to_isoformat @@ -79,6 +80,79 @@ def resolve_provider(provider: str) -> str: ) +@dataclass(frozen=True) +class ChannelEpoch: + """One row of fdsnws-station ``format=text&level=channel`` output — + the lean route to instrument gain and coordinates: one HTTP request, + stdlib parsing, no StationXML, no ObsPy. + + ``scale`` is the total sensitivity (counts per ``scale_units`` ground + motion, referenced at ``scale_frequency`` Hz). Dividing raw counts by + ``scale`` is a GAIN correction: exact at the reference frequency, flat + to the instrument's passband shape elsewhere. For full deconvolution + use :mod:`seisfetch.contrib.response`. + """ + + network: str + station: str + location: str + channel: str + latitude: float + longitude: float + elevation: float + depth: float + scale: float | None + scale_frequency: float | None + scale_units: str | None + sample_rate: float + start: str + end: str | None + + def covers(self, time_iso: str) -> bool: + t = str(time_iso) + return self.start <= t and (self.end is None or t <= self.end) + + +def parse_channel_text(text: str) -> list[ChannelEpoch]: + """Parse fdsnws-station ``format=text&level=channel`` rows. + + Column order per the FDSN spec: Network|Station|Location|Channel| + Latitude|Longitude|Elevation|Depth|Azimuth|Dip|SensorDescription| + Scale|ScaleFreq|ScaleUnits|SampleRate|StartTime|EndTime. + """ + + def _f(v): + v = v.strip() + return float(v) if v else None + + out = [] + for line in text.splitlines(): + if not line.strip() or line.startswith("#"): + continue + p = line.split("|") + if len(p) < 17: + continue + out.append( + ChannelEpoch( + network=p[0].strip(), + station=p[1].strip(), + location=p[2].strip(), + channel=p[3].strip(), + latitude=float(p[4]), + longitude=float(p[5]), + elevation=_f(p[6]) or 0.0, + depth=_f(p[7]) or 0.0, + scale=_f(p[11]), + scale_frequency=_f(p[12]), + scale_units=(p[13].strip() or None), + sample_rate=float(p[14]), + start=p[15].strip(), + end=(p[16].strip() or None), + ) + ) + return out + + def list_providers() -> dict[str, str]: return dict(PROVIDERS) diff --git a/tests/test_earth2_live.py b/tests/test_earth2_live.py new file mode 100644 index 0000000..3f4322a --- /dev/null +++ b/tests/test_earth2_live.py @@ -0,0 +1,137 @@ +"""Tests for the live Earth2Studio source: channel-text metadata parsing, +mandatory calibration, window assembly, and real async.""" + +import asyncio +from datetime import datetime, timezone + +import numpy as np +import pytest + +pytest.importorskip("xarray") +pytest.importorskip("pandas") + +from seisfetch.earth2 import SeisfetchLiveSource +from seisfetch.fdsn import ChannelEpoch, parse_channel_text +from tests.helpers import make_mseed + +TEXT = ( + "#Network|Station|Location|Channel|Latitude|Longitude|Elevation|Depth|" + "Azimuth|Dip|SensorDescription|Scale|ScaleFreq|ScaleUnits|SampleRate|" + "StartTime|EndTime\n" + "NZ|WEL|10|HHZ|-41.2865|174.768|138.0|0.0|0.0|-90.0|Broadband|" + "600000000.0|1.0|M/S|100.0|2010-01-01T00:00:00|\n" + "NZ|WEL|10|HHZ|-41.2865|174.768|138.0|0.0|0.0|-90.0|Old sensor|" + "300000000.0|1.0|M/S|100.0|2000-01-01T00:00:00|2009-12-31T23:59:59\n" +) + + +class TestChannelText: + def test_parse(self): + eps = parse_channel_text(TEXT) + assert len(eps) == 2 + ep = eps[0] + assert (ep.network, ep.station, ep.location, ep.channel) == ( + "NZ", + "WEL", + "10", + "HHZ", + ) + assert ep.scale == 600000000.0 + assert ep.latitude == pytest.approx(-41.2865) + assert ep.end is None + + def test_covers_selects_epoch(self): + eps = parse_channel_text(TEXT) + now = [e for e in eps if e.covers("2022-01-02T00:00:00")] + old = [e for e in eps if e.covers("2005-06-01T00:00:00")] + assert len(now) == 1 and now[0].scale == 600000000.0 + assert len(old) == 1 and old[0].scale == 300000000.0 + + def test_blank_scale_is_none(self): + row = TEXT.splitlines()[1].split("|") + row[11] = "" + eps = parse_channel_text("|".join(row)) + assert eps[0].scale is None + + +class TestLiveSource: + NSLC = "IU.ANMO.00.BHZ" + + def _source(self, monkeypatch, calibrate="gain", raw=None): + src = SeisfetchLiveSource([self.NSLC], window_s=5.0, calibrate=calibrate) + # inject metadata (no network in unit tests) + src._meta[self.NSLC] = [ + ChannelEpoch( + network="IU", + station="ANMO", + location="00", + channel="BHZ", + latitude=34.9, + longitude=-106.5, + elevation=1700.0, + depth=0.0, + scale=2.0, + scale_frequency=1.0, + scale_units="M/S", + sample_rate=100.0, + start="1990-01-01T00:00:00", + end=None, + ) + ] + if raw is None: + raw = make_mseed() # IU.ANMO.00.BHZ, 100 sps, 2024-01-15T00:00Z + + def fake_get_raw(self_client, net, sta, day, location="", channel=""): + return raw + + from seisfetch.s3 import S3OpenClient + + monkeypatch.setattr(S3OpenClient, "get_raw", fake_get_raw) + return src + + def test_requires_calibration(self): + with pytest.raises(ValueError, match="physical units"): + SeisfetchLiveSource(["IU.ANMO.00.BHZ"], calibrate="none") + + def test_gain_corrected_window(self, monkeypatch): + raw = make_mseed() # one payload shared by both sources + src = self._source(monkeypatch, raw=raw) + t0 = datetime(2024, 1, 15, 0, 0, 0, tzinfo=timezone.utc) + da = src(t0) + assert da.dims == ("time", "variable", "sample") + assert da.shape == (1, 1, 500) # 5 s at 100 sps + assert da.attrs["calibration"] == "gain" + # gain division: scale=2 output x2 must equal scale=1 output exactly + src2 = self._source(monkeypatch, raw=raw) + src2._meta[self.NSLC][0] = ChannelEpoch( + **{**src2._meta[self.NSLC][0].__dict__, "scale": 1.0} + ) + da2 = src2(t0) + np.testing.assert_allclose(da.values * 2.0, da2.values) + + def test_window_count_and_coords(self, monkeypatch): + src = self._source(monkeypatch) + lat, lon, elev = src.coords(self.NSLC) + assert lat == pytest.approx(34.9) + t0 = datetime(2024, 1, 15, 0, 0, 0, tzinfo=timezone.utc) + da = src([t0, t0]) + assert da.shape[0] == 2 + + def test_missing_scale_raises(self, monkeypatch): + src = self._source(monkeypatch) + src._meta[self.NSLC][0] = ChannelEpoch( + **{**src._meta[self.NSLC][0].__dict__, "scale": None} + ) + t0 = datetime(2024, 1, 15, 0, 0, 0, tzinfo=timezone.utc) + with pytest.raises(LookupError, match="Scale"): + src(t0) + + def test_async_fetch_is_real(self, monkeypatch): + src = self._source(monkeypatch) + t0 = datetime(2024, 1, 15, 0, 0, 0, tzinfo=timezone.utc) + + async def go(): + return await src.fetch(t0) + + da = asyncio.run(go()) + assert da.shape == (1, 1, 500) From ae54c9b8ad07bc3c6e046b550f40f494af08c4ec Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Fri, 7 Aug 2026 08:07:53 +0200 Subject: [PATCH 3/3] Address Copilot review on PR #5 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - ChannelEpoch.covers compares aware UTC datetimes, not ISO strings — offsets and fraction-width variants now select the right epoch (a -05:00 timestamp equal to the epoch boundary instant was the failing case); SeisfetchLiveSource passes normalized UTC ISO times downstream. - Scale=None (missing) and Scale=0 (defective) are now distinct loud failures; 0.0 is never treated as absent. - _geonet_key raises on blank location with guidance (GeoNet keys always carry a location code; wildcard discovery finds it). - DataArray units come from the channel's ScaleUnits metadata — accelerometer channels report M/S**2, mixed sets say so. - All three caches behind one lock (fetch() runs in worker threads); day fetches stay outside the lock via double-checked population. - Metadata provider honors the datacenter override instead of always routing by network, so calibration metadata and waveforms always come from the same archive. Co-Authored-By: Claude Fable 5 --- seisfetch/earth2.py | 103 +++++++++++++++++++++++++++----------- seisfetch/fdsn.py | 31 +++++++++++- seisfetch/s3.py | 9 +++- tests/test_earth2_live.py | 42 ++++++++++++++++ tests/test_s3.py | 4 ++ 5 files changed, 157 insertions(+), 32 deletions(-) diff --git a/seisfetch/earth2.py b/seisfetch/earth2.py index 7007571..180cdb2 100644 --- a/seisfetch/earth2.py +++ b/seisfetch/earth2.py @@ -16,7 +16,8 @@ from __future__ import annotations import logging -from datetime import datetime +import threading +from datetime import datetime, timezone from typing import Any import numpy as np @@ -385,6 +386,15 @@ def bundle_to_earth2(bundle, variables: list[str] | None = None): } +def _iso_utc(dt: datetime) -> str: + """Normalized UTC ISO string (no offset suffix): naive datetimes are + taken as UTC. Keeps epoch selection and StationXML parsing free of + "+00:00"-suffix ambiguity.""" + if dt.tzinfo is not None: + dt = dt.astimezone(timezone.utc).replace(tzinfo=None) + return dt.isoformat() + + def channel_metadata(network, station, location, channel, start, end, provider=None): """Channel epochs (gain + coordinates) via the FDSN station text service. @@ -463,27 +473,52 @@ def __init__( self._day_cache: dict = {} # (nslc, date) -> TraceBundle self._meta: dict = {} # nslc -> list[ChannelEpoch] self._resp: dict = {} # nslc -> ChannelResponse (calibrate="response") + self._units: dict = {} # nslc -> ScaleUnits of the last-used epoch + # one coarse lock guards all three caches: fetch() runs __call__ in + # worker threads, and unsynchronized eviction/population races + self._lock = threading.Lock() + + def _provider(self, network: str) -> str: + # metadata must come from the same archive as the waveforms: honor + # the datacenter override instead of always routing by network + from seisfetch.s3 import route_network + + return _DC_TO_PROVIDER[self._datacenter or route_network(network)] # -- metadata ---------------------------------------------------------- # def _epochs(self, nslc: str): - if nslc not in self._meta: - net, sta, loc, cha = nslc.split(".") - self._meta[nslc] = channel_metadata( - net, sta, loc, cha, "1970-01-01", "2100-01-01" - ) - if not self._meta[nslc]: - raise LookupError(f"no channel metadata for {nslc}") - return self._meta[nslc] + with self._lock: + if nslc not in self._meta: + net, sta, loc, cha = nslc.split(".") + eps = channel_metadata( + net, + sta, + loc, + cha, + "1970-01-01", + "2100-01-01", + provider=self._provider(net), + ) + if not eps: + raise LookupError(f"no channel metadata for {nslc}") + self._meta[nslc] = eps + return self._meta[nslc] def _gain(self, nslc: str, time_iso: str) -> float: for ep in self._epochs(nslc): if ep.covers(time_iso): - if not ep.scale: + if ep.scale is None: raise LookupError( f"{nslc}: channel epoch has no Scale (sensitivity) — " "cannot calibrate" ) + if ep.scale == 0.0: + raise LookupError( + f"{nslc}: channel epoch declares Scale=0 — defective " + "metadata, refusing to divide" + ) + self._units[nslc] = ep.scale_units or "unknown" return ep.scale raise LookupError(f"{nslc}: no channel epoch covers {time_iso}") @@ -498,15 +533,19 @@ def _day_bundle(self, nslc: str, day): from seisfetch.s3 import S3OpenClient key = (nslc, day.isoformat()) - if key not in self._day_cache: - net, sta, loc, cha = nslc.split(".") - raw = S3OpenClient(datacenter=self._datacenter).get_raw( - net, sta, day.isoformat(), location=loc, channel=cha - ) - self._day_cache[key] = parse_mseed(raw) + with self._lock: + if key in self._day_cache: + return self._day_cache[key] + net, sta, loc, cha = nslc.split(".") + raw = S3OpenClient(datacenter=self._datacenter).get_raw( + net, sta, day.isoformat(), location=loc, channel=cha + ) + bundle = parse_mseed(raw) + with self._lock: + self._day_cache.setdefault(key, bundle) while len(self._day_cache) > self._cache_days * len(self.channels): self._day_cache.pop(next(iter(self._day_cache))) - return self._day_cache[key] + return self._day_cache[key] def _window(self, nslc: str, t0: datetime) -> np.ndarray: """One calibrated window [t0, t0+window_s), NaN-padded over gaps.""" @@ -533,7 +572,7 @@ def _window(self, nslc: str, t0: datetime) -> np.ndarray: j0, j1 = max(i0, 0), min(i0 + s.npts, n) if j1 > j0: out[j0:j1] = src[j0 - i0 : j1 - i0] - t_iso = t0.isoformat() + t_iso = _iso_utc(t0) if self.calibrate == "gain": return out / self._gain(nslc, t_iso) return self._deconvolve(nslc, out, fs, t_iso) @@ -544,13 +583,12 @@ def _deconvolve(self, nslc, x, fs, time_iso): remove_response_np, ) from seisfetch.fdsn import FDSNClient - from seisfetch.s3 import route_network - if nslc not in self._resp: + with self._lock: + cached = nslc in self._resp + if not cached: net, sta, loc, cha = nslc.split(".") - xml = FDSNClient( - provider=_DC_TO_PROVIDER[route_network(net)] - ).get_station_text( + xml = FDSNClient(provider=self._provider(net)).get_station_text( network=net, station=sta, location=loc if loc else "--", @@ -558,9 +596,9 @@ def _deconvolve(self, nslc, x, fs, time_iso): level="response", format="xml", ) - self._resp[nslc] = parse_stationxml_response( - xml.encode(), net, sta, loc, cha, time_iso - ) + resp = parse_stationxml_response(xml.encode(), net, sta, loc, cha, time_iso) + with self._lock: + self._resp.setdefault(nslc, resp) mask = np.isnan(x) filled = np.where(mask, 0.0, x) v = remove_response_np(filled, fs, self._resp[nslc], output="VEL") @@ -601,13 +639,20 @@ def __call__(self, time, variable=None): "sample": np.arange(nsamp), }, attrs={ - "units": "m/s (gain-corrected)" - if self.calibrate == "gain" - else "m/s (response-removed)", + "units": self._units_attr(variables), "calibration": self.calibrate, }, ) + def _units_attr(self, variables) -> str: + if self.calibrate == "response": + return "m/s (response-removed)" + # gain mode: report the channels' ScaleUnits from metadata (an + # accelerometer channel is M/S**2, not m/s) + units = {self._units.get(v, "unknown") for v in variables} + label = units.pop() if len(units) == 1 else "mixed: " + ", ".join(sorted(units)) + return f"{label} (gain-corrected)" + async def fetch(self, time, variable=None): """True async fetch: runs the blocking pipeline in a worker thread.""" import asyncio diff --git a/seisfetch/fdsn.py b/seisfetch/fdsn.py index c75cd89..704c240 100644 --- a/seisfetch/fdsn.py +++ b/seisfetch/fdsn.py @@ -80,6 +80,29 @@ def resolve_provider(provider: str) -> str: ) +def _iso_to_utc(s: str): + """ISO 8601 -> aware UTC datetime (stdlib only). Handles 'Z', numeric + offsets, and fractional seconds of any width; naive times are UTC. + String comparison of ISO timestamps breaks as soon as one side carries + an offset or a different fraction width — epochs must be compared as + datetimes.""" + import re + from datetime import datetime, timezone + + s = s.strip() + if not s: + return None + if s.endswith("Z"): + s = s[:-1] + "+00:00" + m = re.match(r"^([^.]*\.)(\d+)(.*)$", s) + if m: + s = m.group(1) + m.group(2)[:6].ljust(6, "0") + m.group(3) + dt = datetime.fromisoformat(s) + if dt.tzinfo is None: + dt = dt.replace(tzinfo=timezone.utc) + return dt.astimezone(timezone.utc) + + @dataclass(frozen=True) class ChannelEpoch: """One row of fdsnws-station ``format=text&level=channel`` output — @@ -109,8 +132,12 @@ class ChannelEpoch: end: str | None def covers(self, time_iso: str) -> bool: - t = str(time_iso) - return self.start <= t and (self.end is None or t <= self.end) + t = _iso_to_utc(str(time_iso)) + if t is None: + raise ValueError(f"unparseable time {time_iso!r}") + start = _iso_to_utc(self.start) + end = _iso_to_utc(self.end) if self.end else None + return (start is None or start <= t) and (end is None or t <= end) def parse_channel_text(text: str) -> list[ChannelEpoch]: diff --git a/seisfetch/s3.py b/seisfetch/s3.py index 30c096f..6324073 100644 --- a/seisfetch/s3.py +++ b/seisfetch/s3.py @@ -94,9 +94,16 @@ def _geonet_key(network, station, year, doy, location="", channel="", **_): Layout (verified on the live bucket, 2026-08-07): ``waveforms/miniseed/{Y}/{Y}.{DDD}/{STA}.{NET}/{Y}.{DDD}.{STA}.{LOC}-{CHA}.{NET}.D`` e.g. ``waveforms/miniseed/2022/2022.002/WEL.NZ/2022.002.WEL.10-HHZ.NZ.D``. - GeoNet channels always carry a numeric location code (10, 20, ...). + GeoNet channels always carry a numeric location code (10, 20, ...), + so a blank location cannot form a valid key — use ``location="*"`` + (wildcard discovery) or pass the real code. """ loc = location if location and location != "*" else "" + if not loc: + raise ValueError( + "GeoNet keys require a location code (e.g. '10'); use " + "location='*' to discover it" + ) return ( f"waveforms/miniseed/{year}/{year}.{doy:03d}/{station}.{network}/" f"{year}.{doy:03d}.{station}.{loc}-{channel}.{network}.D" diff --git a/tests/test_earth2_live.py b/tests/test_earth2_live.py index 3f4322a..2f89fa1 100644 --- a/tests/test_earth2_live.py +++ b/tests/test_earth2_live.py @@ -47,6 +47,23 @@ def test_covers_selects_epoch(self): assert len(now) == 1 and now[0].scale == 600000000.0 assert len(old) == 1 and old[0].scale == 300000000.0 + def test_covers_is_timezone_safe(self): + # lexicographic string comparison fails on offset/fraction variants; + # these must all select the modern epoch (Copilot review, PR #5) + eps = parse_channel_text(TEXT) + modern = eps[0] + for t in ( + "2022-01-02T00:00:00+00:00", + "2022-01-02T00:00:00Z", + "2022-01-01T19:00:00-05:00", + "2022-01-02T00:00:00.0000", + ): + assert modern.covers(t), t + # 2009-12-31T19:00-05:00 == 2010-01-01T00:00Z: modern epoch, not + # the old one — exactly the case string comparison gets wrong + assert modern.covers("2009-12-31T19:00:00-05:00") + assert not eps[1].covers("2009-12-31T19:00:00-05:00") + def test_blank_scale_is_none(self): row = TEXT.splitlines()[1].split("|") row[11] = "" @@ -117,6 +134,31 @@ def test_window_count_and_coords(self, monkeypatch): da = src([t0, t0]) assert da.shape[0] == 2 + def test_zero_scale_raises(self, monkeypatch): + src = self._source(monkeypatch) + src._meta[self.NSLC][0] = ChannelEpoch( + **{**src._meta[self.NSLC][0].__dict__, "scale": 0.0} + ) + t0 = datetime(2024, 1, 15, 0, 0, 0, tzinfo=timezone.utc) + with pytest.raises(LookupError, match="Scale=0"): + src(t0) + + def test_units_from_metadata(self, monkeypatch): + # an accelerometer channel must not be labeled m/s + src = self._source(monkeypatch) + src._meta[self.NSLC][0] = ChannelEpoch( + **{**src._meta[self.NSLC][0].__dict__, "scale_units": "M/S**2"} + ) + t0 = datetime(2024, 1, 15, 0, 0, 0, tzinfo=timezone.utc) + da = src(t0) + assert da.attrs["units"] == "M/S**2 (gain-corrected)" + + def test_datacenter_override_routes_metadata(self): + src = SeisfetchLiveSource(["IU.ANMO.00.BHZ"], datacenter="geonet") + assert src._provider("IU") == "GEONET" + src2 = SeisfetchLiveSource(["IU.ANMO.00.BHZ"]) + assert src2._provider("IU") == "EARTHSCOPE" + def test_missing_scale_raises(self, monkeypatch): src = self._source(monkeypatch) src._meta[self.NSLC][0] = ChannelEpoch( diff --git a/tests/test_s3.py b/tests/test_s3.py index 519e759..dc149c3 100644 --- a/tests/test_s3.py +++ b/tests/test_s3.py @@ -36,6 +36,10 @@ def test_ncedc_no_location(self): k = _ncedc_key("NC", "JBL", 2023, 1, location="", channel="HHZ") assert "JBL.NC.HHZ..D.2023.001" in k + def test_geonet_blank_location_raises(self): + with pytest.raises(ValueError, match="location code"): + _geonet_key("NZ", "WEL", 2022, 2, location="", channel="HHZ") + def test_geonet(self): # verified against the live bucket 2026-08-07: # waveforms/miniseed/2022/2022.002/WEL.NZ/2022.002.WEL.10-HHZ.NZ.D