From 3e7e778e8d4a0e2e3384db6f1f80b554631dce2e Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 17 Sep 2026 19:44:33 +0000 Subject: [PATCH 1/2] fix(acquire): ask for waveforms by name, so the radius applies to them too MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `max_radius_km` reached `get_stations` and could never reach the waveform call: FDSN dataselect takes no geographic parameters, so repeating the config's wildcards asked the archive for the whole network however the stations had been selected. Found on a live fetch. A Magna smoke test with `max_radius_km = 50` came back with three BH channels from a station 225 km out, absent from the archived `inventory.xml` — an artefact holding waveforms with no metadata beside them, and a config whose stated spatial restriction the download did not honour. The waveform request is now built from the inventory the station query returned: one bulk line per channel, epochs collapsed. That makes the radius apply to both halves and makes the orphan case unreachable by construction — what was asked for is exactly what has metadata. An empty station query now says so rather than sending a bulk request with no lines, which FDSN rejects with a message about syntax. `min_radius_km` had the same hole and is fixed by the same change. The four new tests fail against the previous implementation. The fake client drops `get_waveforms` entirely, so a return to the wildcard path cannot pass. Co-Authored-By: Claude Opus 5 --- src/specmod/acquire.py | 58 +++++++++++++++++++++++----- tests/test_acquire.py | 87 ++++++++++++++++++++++++++++++++++++++---- 2 files changed, 127 insertions(+), 18 deletions(-) diff --git a/src/specmod/acquire.py b/src/specmod/acquire.py index 658f8e0..71c586f 100644 --- a/src/specmod/acquire.py +++ b/src/specmod/acquire.py @@ -232,6 +232,48 @@ def _resolve_event( return resolved, catalogue +def _bulk_request( + inventory: Any, config: AcquisitionConfig, start: Any, end: Any +) -> list[tuple[Any, ...]]: + """One request line per channel the inventory holds. + + The waveforms are asked for by name, from the inventory the station query + returned, rather than by repeating the config's wildcards. FDSN dataselect + takes no geographic parameters, so a wildcard request is unbounded however + the stations were selected: ``max_radius_km`` reaches ``get_stations`` and + can never reach the waveform call. Asking channel by channel is what makes + the radius apply to both, and it is also what guarantees every trace in the + archive has metadata beside it. + + Channel epochs are collapsed. An instrument replaced mid-window appears + twice in the inventory and is one request. + """ + lines = { + ( + network.code, + station.code, + channel.location_code or "--", + channel.code, + ) + for network in inventory + for station in network + for channel in station + } + if not lines: + raise ValueError( + f"no channels matched {config.stations.network}." + f"{config.stations.station}.{config.stations.location}." + f"{config.stations.channel}" + + ( + "" + if config.stations.max_radius_km is None + else f" within {config.stations.max_radius_km} km" + ) + + f" at {config.data_centre}, so there are no waveforms to ask for." + ) + return [(*line, start, end) for line in sorted(lines)] + + def fetch( config: str | Path | AcquisitionConfig, out: str | Path, @@ -243,8 +285,11 @@ def fetch( Returns the manifest, which is also written to ``manifest.json`` beside the data. ``client`` accepts anything with the ObsPy FDSN client's - ``get_events``, ``get_stations`` and ``get_waveforms`` methods; tests pass a - fake so that no test touches the network. + ``get_events``, ``get_stations`` and ``get_waveforms_bulk`` methods; tests + pass a fake so that no test touches the network. + + The waveform request is built from the station query's answer, not from the + config's wildcards — see :func:`_bulk_request`. """ import obspy # noqa: PLC0415 @@ -276,14 +321,7 @@ def fetch( station_kwargs["minradius"] = config.stations.min_radius_km / 111.195 inventory = client.get_stations(**station_kwargs) - stream = client.get_waveforms( - network=config.stations.network, - station=config.stations.station, - location=config.stations.location, - channel=config.stations.channel, - starttime=start, - endtime=end, - ) + stream = client.get_waveforms_bulk(_bulk_request(inventory, config, start, end)) paths = EventDirectory(Path(out) / event.origin) paths.waveforms.mkdir(parents=True, exist_ok=True) diff --git a/tests/test_acquire.py b/tests/test_acquire.py index 998500c..d3e94fa 100644 --- a/tests/test_acquire.py +++ b/tests/test_acquire.py @@ -44,10 +44,13 @@ class FakeClient: hand-written double makes those readable. """ - def __init__(self, *, n_traces: int = 2, events: Any = None) -> None: + def __init__( + self, *, n_traces: int = 2, events: Any = None, inventory: Any = None + ) -> None: self.calls: dict[str, dict[str, Any]] = {} self._n_traces = n_traces self._events = events + self._inventory = inventory def get_events(self, **kwargs: Any) -> Any: self.calls["get_events"] = kwargs @@ -55,10 +58,12 @@ def get_events(self, **kwargs: Any) -> Any: def get_stations(self, **kwargs: Any) -> Any: self.calls["get_stations"] = kwargs + if self._inventory is not None: + return self._inventory return obspy.read_inventory() # ObsPy's bundled example inventory - def get_waveforms(self, **kwargs: Any) -> Any: - self.calls["get_waveforms"] = kwargs + def get_waveforms_bulk(self, bulk: Any) -> Any: + self.calls["get_waveforms_bulk"] = {"bulk": list(bulk)} stream = obspy.read() # three traces of example data return obspy.Stream(stream[: self._n_traces]) @@ -133,9 +138,75 @@ def test_the_window_is_taken_relative_to_the_origin(self, tmp_path: Path) -> Non fetch(_config(), out=tmp_path, client=client) origin = obspy.UTCDateTime(EXPLICIT["origin"]) - asked = client.calls["get_waveforms"] - assert asked["starttime"] == origin - 5.0 - assert asked["endtime"] == origin + 60.0 + bulk = client.calls["get_waveforms_bulk"]["bulk"] + assert bulk + # `UTCDateTime` is not usable in a set, so compare elementwise. + assert all(line[4] == origin - 5.0 for line in bulk) + assert all(line[5] == origin + 60.0 for line in bulk) + + def test_the_waveforms_are_asked_for_by_name_not_by_wildcard( + self, tmp_path: Path + ) -> None: + """The request that carries the radius. + + FDSN dataselect has no geographic parameters, so a wildcard waveform + request is unbounded whatever the station query asked for. A live fetch + for Magna with `max_radius_km = 50` returned three BH channels from a + station 225 km out, with no metadata in the archived inventory, because + the station pattern went to the wire verbatim. + """ + client = FakeClient() + fetch( + _config(stations=StationSpec(max_radius_km=50.0)), + out=tmp_path, + client=client, + ) + + bulk = client.calls["get_waveforms_bulk"]["bulk"] + assert bulk + assert not [line for line in bulk if "*" in line[:4] or "?" in line[:4]] + + def test_no_waveform_is_asked_for_without_metadata_beside_it( + self, tmp_path: Path + ) -> None: + """The property the bulk request gives for free. + + Whatever the station query returned is exactly what is requested, so + the archive cannot hold a trace whose station is missing from + `inventory.xml`. + """ + client = FakeClient() + fetch(_config(stations=StationSpec(station="*")), out=tmp_path, client=client) + + inventory = obspy.read_inventory() + available = { + (network.code, station.code, channel.location_code or "--", channel.code) + for network in inventory + for station in network + for channel in station + } + asked = {line[:4] for line in client.calls["get_waveforms_bulk"]["bulk"]} + assert asked == available + + def test_a_channel_with_two_epochs_is_one_request(self, tmp_path: Path) -> None: + """An instrument replaced mid-window appears twice in the inventory.""" + doubled = obspy.read_inventory() + obspy.read_inventory() + client = FakeClient(inventory=doubled) + fetch(_config(), out=tmp_path, client=client) + + bulk = client.calls["get_waveforms_bulk"]["bulk"] + assert len(bulk) == len({line[:4] for line in bulk}) + + def test_an_empty_station_query_says_so(self, tmp_path: Path) -> None: + """Better than a bulk request with no lines, which FDSN rejects with a + message about syntax rather than about the radius.""" + client = FakeClient(inventory=obspy.Inventory()) + with pytest.raises(ValueError, match="no channels matched"): + fetch( + _config(stations=StationSpec(network="LV", max_radius_km=1.0)), + out=tmp_path, + client=client, + ) def test_it_asks_for_the_response(self, tmp_path: Path) -> None: """Raw counts are only useful with the response stored beside them.""" @@ -341,7 +412,7 @@ def test_the_catalogue_is_asked_not_the_archive(self, tmp_path: Path) -> None: assert "get_events" in catalogue.calls assert "get_events" not in archive.calls # And the waveforms still come from the archive. - assert "get_waveforms" in archive.calls + assert "get_waveforms_bulk" in archive.calls def test_the_manifest_records_both(self, tmp_path: Path) -> None: manifest = fetch( @@ -360,7 +431,7 @@ def test_one_centre_is_reused_when_no_catalogue_is_named( client = FakeClient(events=_magna_catalogue()) fetch(_config(event=EventSpec(eventid="x")), out=tmp_path, client=client) assert "get_events" in client.calls - assert "get_waveforms" in client.calls + assert "get_waveforms_bulk" in client.calls def _magna_catalogue() -> _Catalogue: From 2bd4fdc1ae2f75889767ce56ad1765fa6992f38c Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 17 Sep 2026 21:59:48 +0000 Subject: [PATCH 2/2] feat(acquire): rank co-sited instruments, and cut the radius in kilometres MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two ways the thin wrapper handed the caller a trap that `MassDownloader` closes for free. `channel = "HH*,BH*,HN*,EN*"` returns every instrument at every station that has them — 507 channels for Magna's 88 stations. The waveform request is built from that answer, so all of them are downloaded, and a station with a broadband beside an accelerometer contributes the same ground motion twice: two entries in a `SpectrumSet`, weighted as two stations by anything that averages over them. `channel_priorities` and `location_priorities` rank them instead, first-match-wins per station, in ObsPy's order — channels within a location, then locations — so a list carried over from `Restrictions` selects the same instrument. The inventory is pruned to match, keeping the archive describing exactly what it holds. With no list set, the fetch warns and names the stations; `["*"]` is how a config says it wants all of them and means it. `max_radius_km` was converted with a single 111.195 km/degree. A degree of arc is not a fixed distance: measured on WGS84 across latitude and azimuth it runs 110.574 km at the equator to 111.691 km at the poles, a 1.01% spread, so one constant moves the boundary by -0.44% to +0.56% — 0.3 km on a 50 km radius, 2 km on Magna's 400 km. The station query is now sent at the loose end of that bracket, deliberately wide, and the cut made against the true epicentral distance from `gps2dist_azimuth`. The bracketing constants are asserted against the geodetics rather than trusted. Also fixed: `min_radius_km` was only sent alongside `max_radius_km`, so on its own it did nothing. A non-positive radius and an inverted annulus are now refused rather than silently selecting nothing, which is the failure mode of a metres-for-kilometres typo on a local array. Co-Authored-By: Claude Opus 5 --- datasets/magna_2020.toml | 6 + docs/REFACTOR_PLAN.md | 9 + docs/releasing-data.md | 47 +++++ src/specmod/acquire.py | 285 ++++++++++++++++++++++++-- tests/test_acquire.py | 433 ++++++++++++++++++++++++++++++++++++++- 5 files changed, 754 insertions(+), 26 deletions(-) diff --git a/datasets/magna_2020.toml b/datasets/magna_2020.toml index cd3bba8..47a322c 100644 --- a/datasets/magna_2020.toml +++ b/datasets/magna_2020.toml @@ -35,6 +35,12 @@ location = "*" # fetched too: the workflow rotates horizontals to R/T, and rotation is not a # reason to discard Z from the archive. channel = "HH*,BH*,HN*,EN*" +# Deliberately unranked. `channel_priorities` would keep one instrument per +# station — the right default for a regional dataset, and what silences the +# co-sited warning this fetch prints — but here the near-source broadbands are +# the records most likely to clip, and the accelerometer beside them is the +# reason the paper covers that distance range at all. Both are kept and the +# choice is made at processing time. max_radius_km = 400.0 [window] diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index f7bde5f..16de3a9 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -2434,6 +2434,15 @@ Two design constraints worth fixing now: reimplementation of `MassDownloader`, which is built for large restricted-data campaigns and is more machinery than this needs. + What that choice costs, recorded as it is paid: `MassDownloader` supplies + `CircularDomain`, `sanitize=True` and `channel_priorities` for free, and each + has had to be built by hand here after its absence showed up in a live fetch. + The radius reaching only the station query, an archive holding waveforms with + no metadata beside them, and 507 channels for 88 stations were all one + property of `Restrictions` each. The wrapper stays thin, but "ObsPy already + does it" is an argument for reading what ObsPy does before writing the thin + version. + TOML is suggested for the config, read by stdlib `tomllib` on 3.11+ and consistent with `pyproject.toml`. YAML is equally fine if it reads better for nested station rules; it costs one small dependency. diff --git a/docs/releasing-data.md b/docs/releasing-data.md index 18d70f8..44e5430 100644 --- a/docs/releasing-data.md +++ b/docs/releasing-data.md @@ -46,6 +46,10 @@ python -m json.tool build/magna_2020/manifest.json | head -40 - `channels` — the list *after* wildcard expansion. `HH*` does not tell you what you got; this does. A much shorter list than expected usually means a network was not open at that time, not that the request was malformed. +- `channels_available`, `channels_requested` — what the station query offered, + and what was asked for after the priority lists narrowed it. A large gap + between `channels_requested` and `channels` means the archive did not serve + what its own metadata advertises. - `data_centre`, `obspy_version`, `specmod_version`, `fetched_at` — the provenance that makes the response reproducible-ish. It cannot be fully reproducible: FDSN is not content-addressed, archives get backfilled and @@ -54,6 +58,49 @@ python -m json.tool build/magna_2020/manifest.json | head -40 If any of that is wrong, fix the config and re-fetch. The artefact is about to become immutable. +### Picking one instrument per station + +A `channel` pattern that names several instruments — `HH*,BH*,HN*,EN*` — returns +every one of them at every station that has them. Magna's 88 stations return 507 +channels that way, and a station with both a broadband and an accelerometer +contributes the same ground motion twice: two entries in a `SpectrumSet`, +weighted as two stations by anything that averages over them. + +`channel_priorities` ranks them instead. The first pattern that matches anything +at a station wins, so the list is read as "a broadband if there is one, +otherwise an accelerometer": + +```toml +[stations] +channel = "HH*,BH*,HN*,EN*" +channel_priorities = ["HH[ZNE]", "BH[ZNE]", "HN[ZNE]", "EN[ZNE]"] +location_priorities = ["", "00", "10"] +``` + +`channel` still says what the data centre may offer; the priorities say which of +what it offers to keep. A station matching no pattern at all is dropped — the +list is a restriction as well as a ranking — and the written `inventory.xml` is +pruned to match, so the archive keeps describing exactly what it holds. + +Patterns are matched case-sensitively against SEED codes, which are upper case. +`["*"]` keeps everything, which is how a config says it wants both instruments +and means it; with no list at all, `fetch` warns that it is fetching co-sited +instruments and names the stations. + +### What the radius means + +`max_radius_km` is cut against the true WGS84 distance from the epicentre, not +against the degrees FDSN takes. A degree of arc is 110.574 km at the equator and +111.691 km at the poles, so converting with a single constant moves the boundary +by up to half a percent — ±0.3 km on a 50 km radius, ±2 km on 400 km. The +station query is sent deliberately wide and the exact cut is made locally, which +costs a few extra stations in the query and makes the radius mean kilometres +wherever the event is. + +The radii are epicentral. For a local array the depth is usually the larger +term: every station of a 500 m array over a 2 km-deep event is at 2 km +hypocentral distance. + ## 2. Archive and hash ```bash diff --git a/src/specmod/acquire.py b/src/specmod/acquire.py index 71c586f..63e1dad 100644 --- a/src/specmod/acquire.py +++ b/src/specmod/acquire.py @@ -26,9 +26,12 @@ from __future__ import annotations +import fnmatch import hashlib +import itertools import json import tomllib +import warnings from dataclasses import dataclass, field, replace from datetime import UTC, datetime from pathlib import Path @@ -107,15 +110,74 @@ class StationSpec: The patterns are FDSN wildcards, so the config alone does not say what you got — which is why the manifest records the channel list after expansion. + + ``channel`` says what the data centre may offer; the priority lists say + which of what it offers to keep. A station carrying both a broadband and an + accelerometer matches ``HH*,HN*`` twice, and both records are the same + ground motion measured by two instruments — one station, not two + observations. ``channel_priorities = ["HH[ZNE]", "HN[ZNE]"]`` keeps the + broadband where there is one and the accelerometer where there is not. + + The radii are **epicentral**. For a local array the depth is usually the + larger term — every station of a 500 m array over a 2 km-deep event is at + 2 km hypocentral distance — so a small ``max_radius_km`` selects on map + distance, not on how far the waves travelled. """ network: str = "*" station: str = "*" location: str = "*" channel: str = "*" - #: Kilometres from the epicentre. ``None`` means no limit. + #: Kilometres from the epicentre, cut against the true WGS84 distance + #: rather than against the degrees FDSN takes. ``None`` means no limit. max_radius_km: float | None = None min_radius_km: float | None = None + #: Channel-code patterns in descending preference, applied per station: + #: every channel matching the first pattern that matches anything is kept + #: and the rest are dropped. ``None`` keeps every channel the query + #: returned; ``["*"]`` does too, and says so deliberately. + channel_priorities: tuple[str, ...] | None = None + #: Location codes in descending preference, ``""`` for a blank one. Applied + #: after ``channel_priorities``, in the same first-match-wins way. + location_priorities: tuple[str, ...] | None = None + + def __post_init__(self) -> None: + for name in ("max_radius_km", "min_radius_km"): + radius = getattr(self, name) + if radius is not None and radius <= 0: + raise ValueError( + f"`{name}` is {radius} km, which selects nothing. Omit it " + f"for no limit. Note the unit is kilometres: a local array " + f"is `max_radius_km = 0.5`, not 500." + ) + if ( + self.max_radius_km is not None + and self.min_radius_km is not None + and self.min_radius_km >= self.max_radius_km + ): + raise ValueError( + f"`min_radius_km` ({self.min_radius_km}) is not inside " + f"`max_radius_km` ({self.max_radius_km}), so the annulus is " + f"empty" + ) + + for name in ("channel_priorities", "location_priorities"): + patterns = getattr(self, name) + if patterns is None: + continue + if isinstance(patterns, str): + raise ValueError( + f"`{name}` is a list of patterns in descending preference, " + f"not one pattern: write [{patterns!r}] rather than " + f"{patterns!r}" + ) + patterns = tuple(str(p) for p in patterns) + if not patterns: + raise ValueError( + f"`{name}` is empty, which would select nothing. Omit it to " + f"keep every channel the station query returns." + ) + object.__setattr__(self, name, patterns) @dataclass(frozen=True) @@ -232,6 +294,187 @@ def _resolve_event( return resolved, catalogue +#: Kilometres of WGS84 surface distance per degree of great-circle arc. A +#: degree is not a fixed number of kilometres — measured across latitude and +#: azimuth it runs from 110.574 km at the equator to 111.691 km at the poles, a +#: spread of 1.01% — so these bracket it rather than approximating it. The +#: station query is widened to the loose end of the bracket and the radius is +#: then cut exactly; see :func:`_within_radius`. +_KM_PER_DEGREE_MIN = 110.5 +_KM_PER_DEGREE_MAX = 111.8 + + +def _location(channel: Any) -> str: + """A channel's location code, with a blank one as the empty string.""" + return channel.location_code or "" + + +def _pruned(inventory: Any, select: Any) -> Any: + """A copy of the inventory holding only ``select(station)`` at each station. + + Stations left with no channel are dropped, and so are networks left with no + station: the waveform request is built from the inventory, so pruning it is + both how a selection reaches the download and how the archive is kept to + describing exactly what it holds. + """ + out = inventory.copy() + for network in out: + for station in network: + station.channels = select(station) + network.stations = [s for s in network.stations if s.channels] + out.networks = [n for n in out.networks if n.stations] + return out + + +def _within_radius(inventory: Any, config: AcquisitionConfig, event: Any) -> Any: + """Drop stations outside the configured radius, measured on WGS84. + + FDSN takes a radius in degrees of arc, so a config in kilometres has to be + converted, and no conversion is exact: a degree of arc is 110.574 km at the + equator and 111.691 km at the poles. Converting with a single constant + therefore moves the boundary by up to half a percent with latitude — ±0.3 km + on a 50 km radius, ±2 km on 400 km. + + So the conversion is not relied on. The station query is widened by the + spread, which costs a few extra stations in the margin, and the cut is made + here against the true epicentral distance. + """ + from obspy.geodetics import gps2dist_azimuth # noqa: PLC0415 + + spec = config.stations + if spec.max_radius_km is None and spec.min_radius_km is None: + return inventory + + def inside(station: Any) -> bool: + metres, _, _ = gps2dist_azimuth( + event.latitude, event.longitude, station.latitude, station.longitude + ) + km = metres / 1000.0 + if spec.max_radius_km is not None and km > spec.max_radius_km: + return False + return not (spec.min_radius_km is not None and km < spec.min_radius_km) + + return _pruned(inventory, lambda s: list(s) if inside(s) else []) + + +def _channel_ids(inventory: Any) -> set[tuple[str, str, str, str]]: + """``(network, station, location, channel)`` for each channel, epochs collapsed. + + An instrument replaced mid-window appears twice in the inventory and is one + channel here, because it is one waveform request. + """ + return { + (network.code, station.code, _location(channel) or "--", channel.code) + for network in inventory + for station in network + for channel in station + } + + +def _by_priority( + channels: list[Any], code: Any, patterns: tuple[str, ...] +) -> list[Any]: + """Every channel matching the first pattern that matches anything. + + Once a pattern matches, the later ones are not considered — that is what + makes the list a ranking rather than a filter. Channels matching no pattern + at all are dropped, so a priority list is also a restriction. + + Patterns are matched case-sensitively against SEED codes, which are upper + case. + """ + for pattern in patterns: + matched = [c for c in channels if fnmatch.fnmatchcase(code(c), pattern)] + if matched: + return matched + return [] + + +def _preferred_channels(station: Any, spec: StationSpec) -> list[Any]: + """The channels at one station that survive its priority lists. + + Channel priorities rank within each location code, then location priorities + rank across what is left. That is the order ObsPy's ``MassDownloader`` + applies them in, so a priority list carried over from it selects the same + instrument here. + """ + channels = list(station) + if spec.channel_priorities is not None: + ranked: list[Any] = [] + for _, group in itertools.groupby(sorted(channels, key=_location), _location): + ranked += _by_priority( + list(group), lambda c: c.code, spec.channel_priorities + ) + channels = ranked + if spec.location_priorities is not None: + channels = _by_priority(channels, _location, spec.location_priorities) + return channels + + +def _warn_about_co_sited_instruments(inventory: Any) -> None: + """Say so when a station carries more than one instrument and nothing ranks them. + + Two instruments at one station are two recordings of the same ground + motion. Downstream they are two entries in a ``SpectrumSet``, weighted as + two stations by anything that averages over them. + """ + crowded: dict[str, list[str]] = {} + for network in inventory: + for station in network: + groups = sorted({f"{_location(c) or '--'}.{c.code[:2]}" for c in station}) + if len(groups) > 1: + crowded[f"{network.code}.{station.code}"] = groups + if not crowded: + return + + examples = "; ".join( + f"{name} has {', '.join(groups)}" + for name, groups in sorted(crowded.items())[:3] + ) + warnings.warn( + f"{len(crowded)} station(s) carry more than one instrument, and every " + f"channel of each will be fetched: {examples}. Rank them with " + f"`channel_priorities` in [stations] — e.g. " + f'["HH[ZNE]", "BH[ZNE]", "HN[ZNE]"] — or set it to ["*"] to keep them ' + f"all deliberately.", + stacklevel=3, + ) + + +def _apply_priorities(inventory: Any, config: AcquisitionConfig) -> Any: + """Prune the station query's answer to the instruments the config prefers. + + The waveform request is built from the inventory, so pruning it is what + makes the priorities reach the download. The inventory is pruned rather + than merely consulted so that the archive keeps describing exactly what it + holds — a station with no selected channel is not written at all. + """ + spec = config.stations + if spec.channel_priorities is None and spec.location_priorities is None: + _warn_about_co_sited_instruments(inventory) + return inventory + + pruned = _pruned(inventory, lambda s: _preferred_channels(s, spec)) + + available = len(_channel_ids(inventory)) + if available and not _channel_ids(pruned): + named = ", ".join( + f"{name}={list(value)}" + for name, value in ( + ("channel_priorities", spec.channel_priorities), + ("location_priorities", spec.location_priorities), + ) + if value is not None + ) + raise ValueError( + f"the station query returned {available} channel(s) and {named} " + f"matched none of them, so there are no waveforms to ask for. The " + f"patterns are matched case-sensitively against SEED codes, which " + f"are upper case." + ) + return pruned + + def _bulk_request( inventory: Any, config: AcquisitionConfig, start: Any, end: Any ) -> list[tuple[Any, ...]]: @@ -248,17 +491,7 @@ def _bulk_request( Channel epochs are collapsed. An instrument replaced mid-window appears twice in the inventory and is one request. """ - lines = { - ( - network.code, - station.code, - channel.location_code or "--", - channel.code, - ) - for network in inventory - for station in network - for channel in station - } + lines = _channel_ids(inventory) if not lines: raise ValueError( f"no channels matched {config.stations.network}." @@ -289,7 +522,10 @@ def fetch( pass a fake so that no test touches the network. The waveform request is built from the station query's answer, not from the - config's wildcards — see :func:`_bulk_request`. + config's wildcards — see :func:`_bulk_request` — after + ``[stations] channel_priorities`` has narrowed that answer to one instrument + per station. With no priorities set, a station carrying two instruments is + fetched twice and a warning says so. """ import obspy # noqa: PLC0415 @@ -313,15 +549,23 @@ def fetch( "endtime": end, "level": "response", } + # Widened in both directions, because the exact cut is made against the + # true epicentral distance below and this query must not pre-empt it. if config.stations.max_radius_km is not None: station_kwargs["latitude"] = event.latitude station_kwargs["longitude"] = event.longitude - station_kwargs["maxradius"] = config.stations.max_radius_km / 111.195 - if config.stations.min_radius_km is not None: - station_kwargs["minradius"] = config.stations.min_radius_km / 111.195 + station_kwargs["maxradius"] = config.stations.max_radius_km / _KM_PER_DEGREE_MIN + if config.stations.min_radius_km is not None: + station_kwargs["latitude"] = event.latitude + station_kwargs["longitude"] = event.longitude + station_kwargs["minradius"] = config.stations.min_radius_km / _KM_PER_DEGREE_MAX inventory = client.get_stations(**station_kwargs) - stream = client.get_waveforms_bulk(_bulk_request(inventory, config, start, end)) + available = len(_channel_ids(inventory)) + inventory = _within_radius(inventory, config, event) + inventory = _apply_priorities(inventory, config) + bulk = _bulk_request(inventory, config, start, end) + stream = client.get_waveforms_bulk(bulk) paths = EventDirectory(Path(out) / event.origin) paths.waveforms.mkdir(parents=True, exist_ok=True) @@ -364,8 +608,11 @@ def fetch( "catalogue_magnitude": event.catalogue_magnitude, "catalogue_magnitude_type": event.catalogue_magnitude_type, "eventid": spec.eventid, - # After wildcard expansion: the config alone does not say what the - # request actually returned. + # Three numbers, because they can differ: what the station query + # offered, what the priority lists narrowed it to, and what the + # archive actually served. The config alone says none of them. + "channels_available": available, + "channels_requested": len(bulk), "channels": sorted({trace.id for trace in stream}), "window": {"start": str(start), "end": str(end)}, }, diff --git a/tests/test_acquire.py b/tests/test_acquire.py index d3e94fa..a8b5915 100644 --- a/tests/test_acquire.py +++ b/tests/test_acquire.py @@ -15,7 +15,11 @@ obspy = pytest.importorskip("obspy") +from obspy.geodetics import gps2dist_azimuth, locations2degrees # noqa: E402 + from specmod.acquire import ( # noqa: E402 + _KM_PER_DEGREE_MAX, + _KM_PER_DEGREE_MIN, AcquisitionConfig, EventSpec, StationSpec, @@ -28,6 +32,14 @@ ROOT = Path(__file__).resolve().parent.parent +# Most tests here fetch without priority lists, and ObsPy's example inventory +# has four instruments at GR.FUR, so the co-sited warning fires all over the +# suite and says nothing about the test that tripped it. `TestChannelPriorities` +# asserts it with `pytest.warns`, which ignores this filter. +pytestmark = pytest.mark.filterwarnings( + "ignore:.*carry more than one instrument:UserWarning" +) + EXPLICIT = { "origin": "2019-08-26T07:30:47.000000Z", "latitude": 53.785021, @@ -155,7 +167,7 @@ def test_the_waveforms_are_asked_for_by_name_not_by_wildcard( station 225 km out, with no metadata in the archived inventory, because the station pattern went to the wire verbatim. """ - client = FakeClient() + client = FakeClient(inventory=_inventory_around_the_event()) fetch( _config(stations=StationSpec(max_radius_km=50.0)), out=tmp_path, @@ -214,11 +226,9 @@ def test_it_asks_for_the_response(self, tmp_path: Path) -> None: fetch(_config(), out=tmp_path, client=client) assert client.calls["get_stations"]["level"] == "response" - def test_a_radius_is_converted_to_degrees_about_the_epicentre( - self, tmp_path: Path - ) -> None: + def test_a_radius_is_asked_about_the_epicentre(self, tmp_path: Path) -> None: """FDSN takes degrees; a config in kilometres is the friendlier unit.""" - client = FakeClient() + client = FakeClient(inventory=_inventory_around_the_event()) fetch( _config(stations=StationSpec(max_radius_km=111.195)), out=tmp_path, @@ -226,8 +236,9 @@ def test_a_radius_is_converted_to_degrees_about_the_epicentre( ) asked = client.calls["get_stations"] - assert asked["maxradius"] == pytest.approx(1.0, rel=1e-6) + assert asked["maxradius"] == pytest.approx(1.0, rel=0.01) assert asked["latitude"] == EXPLICIT["latitude"] + assert asked["longitude"] == EXPLICIT["longitude"] def test_no_radius_means_no_geographic_filter(self, tmp_path: Path) -> None: """Sending a centre with no radius would silently narrow the request.""" @@ -244,6 +255,409 @@ def test_it_does_not_consult_a_catalogue_for_an_explicit_event( assert "get_events" not in client.calls +def _inventory_around_the_event(*, north_km: float = 10.0) -> Any: + """The example inventory, moved to sit near the test event. + + ObsPy's example stations are in Bavaria and the test event is in + Lancashire, so any radius at all excludes every one of them. Each station is + shifted onto the event's meridian, spread northwards, and the tests then + measure what they built rather than assuming a kilometres-per-degree. + """ + inventory = obspy.read_inventory() + # By code, so that BW.RJOB's three epochs stay one place. + codes = sorted({s.code for network in inventory for s in network}) + for network in inventory: + for station in network: + step = codes.index(station.code) + 1 + station.latitude = EXPLICIT["latitude"] + step * north_km / 111.0 + station.longitude = EXPLICIT["longitude"] + return inventory + + +def _epicentral_km(station: Any) -> float: + metres, _, _ = gps2dist_azimuth( + EXPLICIT["latitude"], EXPLICIT["longitude"], station.latitude, station.longitude + ) + return float(metres) / 1000.0 + + +def _two_sensors_at_one_station() -> Any: + """The example inventory, with GR.FUR's broadband and short period split. + + ``HH`` moves to location ``00`` and ``BH`` to ``10``, which is the layout + that makes the two priority lists distinguishable: channel priorities rank + within a location, location priorities rank across them. + """ + inventory = obspy.read_inventory() + for network in inventory: + for station in network: + for channel in station: + if channel.code.startswith("HH"): + channel.location_code = "00" + elif channel.code.startswith("BH"): + channel.location_code = "10" + return inventory + + +def _asked_for(client: FakeClient) -> set[tuple[str, ...]]: + return {line[:4] for line in client.calls["get_waveforms_bulk"]["bulk"]} + + +class TestRadius: + """`max_radius_km` means kilometres, at every latitude. + + FDSN takes a radius in degrees of arc, and a degree is not a fixed number of + kilometres: measured on WGS84 across latitude and azimuth it runs from + 110.574 km at the equator to 111.691 km at the poles. Converting with one + constant moves the boundary by up to half a percent — ±0.3 km on a 50 km + radius, ±2 km on 400 km — so the query is widened and the cut made here. + """ + + def test_the_cut_is_the_true_epicentral_distance(self, tmp_path: Path) -> None: + inventory = _inventory_around_the_event(north_km=10.0) + client = FakeClient(inventory=inventory) + limit = 25.0 + fetch( + _config(stations=StationSpec(max_radius_km=limit)), + out=tmp_path, + client=client, + ) + + expected = { + station.code + for network in inventory + for station in network + if _epicentral_km(station) <= limit + } + assert {line[1] for line in _asked_for(client)} == expected + + @pytest.mark.parametrize("limit", [9.9, 10.1]) + def test_a_station_either_side_of_the_boundary( + self, limit: float, tmp_path: Path + ) -> None: + """The nearest station sits at ~10 km, so the two limits straddle it.""" + inventory = _inventory_around_the_event(north_km=10.0) + client = FakeClient(inventory=inventory) + nearest = min((s for network in inventory for s in network), key=_epicentral_km) + inside = _epicentral_km(nearest) <= limit + + if not inside: + with pytest.raises(ValueError, match="no channels matched"): + fetch( + _config(stations=StationSpec(max_radius_km=limit)), + out=tmp_path, + client=client, + ) + return + + fetch( + _config(stations=StationSpec(max_radius_km=limit)), + out=tmp_path, + client=client, + ) + assert {line[1] for line in _asked_for(client)} == {nearest.code} + + def test_the_station_query_is_widened_so_it_cannot_clip( + self, tmp_path: Path + ) -> None: + """The server's degrees must not exclude what the exact cut would keep. + + Sent at the loosest kilometres-per-degree, so the margin is on the side + of asking for a few stations too many and dropping them here. + """ + client = FakeClient(inventory=_inventory_around_the_event()) + fetch( + _config(stations=StationSpec(max_radius_km=50.0)), + out=tmp_path, + client=client, + ) + + sent = client.calls["get_stations"]["maxradius"] + assert sent >= 50.0 / 110.574 # the tightest degree on the ellipsoid + assert sent <= 50.0 / 110.0 # but not so loose as to be a free-for-all + + def test_a_minimum_radius_alone_is_honoured(self, tmp_path: Path) -> None: + """It used to be sent only alongside a maximum, so on its own it did nothing.""" + inventory = _inventory_around_the_event(north_km=10.0) + client = FakeClient(inventory=inventory) + fetch( + _config(stations=StationSpec(min_radius_km=25.0)), + out=tmp_path, + client=client, + ) + + assert "minradius" in client.calls["get_stations"] + assert client.calls["get_stations"]["latitude"] == EXPLICIT["latitude"] + expected = { + station.code + for network in inventory + for station in network + if _epicentral_km(station) >= 25.0 + } + assert {line[1] for line in _asked_for(client)} == expected + + def test_a_local_array_is_cut_to_the_metre(self, tmp_path: Path) -> None: + """A sub-kilometre deployment is a radius like any other. + + Induced-seismicity and mine arrays span hundreds of metres, which is + where converting kilometres to degrees with one constant is worst in + relative terms: at 500 m the equator-to-pole spread is 5 m, and station + spacings are smaller than that. + """ + inventory = _inventory_around_the_event(north_km=0.2) + client = FakeClient(inventory=inventory) + fetch( + _config(stations=StationSpec(max_radius_km=0.5)), + out=tmp_path, + client=client, + ) + + expected = { + station.code + for network in inventory + for station in network + if _epicentral_km(station) <= 0.5 + } + assert expected # the fixture must straddle the limit, not sit inside it + assert len(expected) < 3 + assert {line[1] for line in _asked_for(client)} == expected + + @pytest.mark.parametrize("radius", [0.0, -1.0]) + def test_a_radius_that_selects_nothing_is_refused(self, radius: float) -> None: + """`max_radius_km = 0` is a typo, not a request for an empty dataset.""" + with pytest.raises(ValueError, match="selects nothing"): + StationSpec(max_radius_km=radius) + + def test_an_empty_annulus_is_refused(self) -> None: + with pytest.raises(ValueError, match="annulus is empty"): + StationSpec(min_radius_km=50.0, max_radius_km=10.0) + + def test_the_conversion_bounds_bracket_every_latitude(self) -> None: + """The two constants are claims about WGS84, so measure them. + + If a later edit tightens them towards the 111.195 they replaced, the + widened query starts clipping stations the exact cut would have kept. + """ + measured = [] + for lat in range(0, 89): + for dlat, dlon in ((1.0, 0.0), (0.0, 1.0), (0.7, 0.7)): + if lat + dlat > 89.5: + continue + degrees = locations2degrees(lat, 0.0, lat + dlat, dlon) + metres, _, _ = gps2dist_azimuth(lat, 0.0, lat + dlat, dlon) + measured.append(metres / 1000.0 / degrees) + + assert min(measured) >= _KM_PER_DEGREE_MIN + assert max(measured) <= _KM_PER_DEGREE_MAX + + +class TestChannelPriorities: + """One instrument per station, rather than every instrument at every station. + + A station query for `HH*,BH*,HN*,EN*` returns a broadband *and* an + accelerometer wherever both are installed, and the waveform request is built + from that answer, so both are downloaded and both reach the pipeline — where + they are two records of one ground motion, weighted as two stations by + anything that averages over them. + """ + + def test_the_first_matching_pattern_wins(self, tmp_path: Path) -> None: + client = FakeClient() + fetch( + _config( + stations=StationSpec(channel_priorities=["HH[ZNE]", "BH[ZNE]"]), + ), + out=tmp_path, + client=client, + ) + + codes = {line[3] for line in _asked_for(client)} + assert codes == {"HHZ", "HHN", "HHE"} + + def test_a_station_without_the_first_choice_falls_through( + self, tmp_path: Path + ) -> None: + """The point of a ranking: BW.RJOB has only EH, and is still fetched.""" + client = FakeClient() + fetch( + _config(stations=StationSpec(channel_priorities=["HH[ZNE]", "EH[ZNE]"])), + out=tmp_path, + client=client, + ) + + by_station = {} + for network, station, _, channel in _asked_for(client): + by_station.setdefault(f"{network}.{station}", set()).add(channel[:2]) + assert by_station == { + "GR.FUR": {"HH"}, + "GR.WET": {"HH"}, + "BW.RJOB": {"EH"}, + } + + def test_a_station_matching_no_pattern_is_dropped(self, tmp_path: Path) -> None: + """A priority list is a restriction as well as a ranking.""" + client = FakeClient() + fetch( + _config(stations=StationSpec(channel_priorities=["HH[ZNE]"])), + out=tmp_path, + client=client, + ) + assert "BW" not in {line[0] for line in _asked_for(client)} + + def test_the_written_inventory_holds_only_what_was_fetched( + self, tmp_path: Path + ) -> None: + """The invariant the bulk request gives, kept in the other direction. + + Pruning the inventory as well as the request is what keeps the archive + describing exactly what it holds — no metadata for an instrument that + was ranked out and never downloaded. + """ + client = FakeClient() + fetch( + _config(stations=StationSpec(channel_priorities=["HH[ZNE]"])), + out=tmp_path, + client=client, + ) + + written = obspy.read_inventory( + str(EventDirectory(tmp_path / EXPLICIT["origin"]).inventory) + ) + on_disk = { + (network.code, station.code, channel.location_code or "--", channel.code) + for network in written + for station in network + for channel in station + } + assert on_disk == _asked_for(client) + + def test_channels_are_ranked_within_each_location(self, tmp_path: Path) -> None: + """Two sensors at one station are two selections, not one. + + ObsPy's `MassDownloader` groups by location before ranking channels, so + a config carried over from it selects the same instruments. + """ + client = FakeClient(inventory=_two_sensors_at_one_station()) + fetch( + _config(stations=StationSpec(channel_priorities=["HH[ZNE]", "BH[ZNE]"])), + out=tmp_path, + client=client, + ) + + fur = {line[2:4] for line in _asked_for(client) if line[1] == "FUR"} + assert {loc for loc, _ in fur} == {"00", "10"} + + def test_location_priorities_pick_one_of_them(self, tmp_path: Path) -> None: + client = FakeClient(inventory=_two_sensors_at_one_station()) + fetch( + _config( + stations=StationSpec( + channel_priorities=["HH[ZNE]", "BH[ZNE]"], + location_priorities=["10", "00"], + ) + ), + out=tmp_path, + client=client, + ) + + fur = {line[2:4] for line in _asked_for(client) if line[1] == "FUR"} + assert {loc for loc, _ in fur} == {"10"} + assert {code[:2] for _, code in fur} == {"BH"} + + def test_a_blank_location_is_the_empty_pattern(self, tmp_path: Path) -> None: + """`""` in the config, `--` on the wire: the two conventions must meet.""" + client = FakeClient() + fetch( + _config(stations=StationSpec(location_priorities=[""])), + out=tmp_path, + client=client, + ) + assert {line[2] for line in _asked_for(client)} == {"--"} + + def test_priorities_that_match_nothing_say_so(self, tmp_path: Path) -> None: + """Otherwise this is an empty bulk request, which FDSN rejects on syntax.""" + client = FakeClient() + with pytest.raises(ValueError, match="matched none of them"): + fetch( + _config(stations=StationSpec(channel_priorities=["ZZ[ZNE]"])), + out=tmp_path, + client=client, + ) + + def test_the_patterns_are_matched_against_upper_case_codes( + self, tmp_path: Path + ) -> None: + """A lower-case pattern matches nothing, and the error says why.""" + client = FakeClient() + with pytest.raises(ValueError, match="upper case"): + fetch( + _config(stations=StationSpec(channel_priorities=["hh[zne]"])), + out=tmp_path, + client=client, + ) + + def test_an_unranked_fetch_warns_about_co_sited_instruments( + self, tmp_path: Path + ) -> None: + """The trap this exists to close, named at the moment it is sprung.""" + with pytest.warns(UserWarning, match="more than one instrument"): + fetch(_config(), out=tmp_path, client=FakeClient()) + + def test_one_instrument_per_station_is_not_warned_about( + self, tmp_path: Path, recwarn: Any + ) -> None: + only_rjob = obspy.read_inventory().select(network="BW") + fetch(_config(), out=tmp_path, client=FakeClient(inventory=only_rjob)) + assert not [w for w in recwarn if "more than one instrument" in str(w.message)] + + def test_keeping_everything_deliberately_is_not_warned_about( + self, tmp_path: Path, recwarn: Any + ) -> None: + """`["*"]` is how a config says it wants both instruments and means it.""" + client = FakeClient() + fetch( + _config(stations=StationSpec(channel_priorities=["*"])), + out=tmp_path, + client=client, + ) + + assert not [w for w in recwarn if "more than one instrument" in str(w.message)] + assert len(_asked_for(client)) == 24 + + def test_the_manifest_counts_what_was_narrowed_away(self, tmp_path: Path) -> None: + manifest = fetch( + _config(stations=StationSpec(channel_priorities=["HH[ZNE]"])), + out=tmp_path, + client=FakeClient(), + ) + + resolved = manifest["resolved"] + assert resolved["channels_available"] == 24 + assert resolved["channels_requested"] == 6 # HH at FUR and WET + assert len(resolved["channels"]) == 2 # what the fake archive serves + + @pytest.mark.parametrize("field", ["channel_priorities", "location_priorities"]) + def test_one_pattern_is_still_a_list(self, field: str) -> None: + """`channel_priorities = "HH[ZNE]"` would otherwise rank by character.""" + with pytest.raises(ValueError, match="not one pattern"): + StationSpec(**{field: "HH[ZNE]"}) + + @pytest.mark.parametrize("field", ["channel_priorities", "location_priorities"]) + def test_an_empty_list_is_refused(self, field: str) -> None: + with pytest.raises(ValueError, match="select nothing"): + StationSpec(**{field: []}) + + def test_a_config_carries_them(self) -> None: + config = AcquisitionConfig.from_dict( + { + "name": "x", + "event": EXPLICIT, + "stations": {"channel_priorities": ["HH[ZNE]", "HN[ZNE]"]}, + } + ) + assert config.stations.channel_priorities == ("HH[ZNE]", "HN[ZNE]") + + class TestManifest: def test_it_records_what_was_asked_and_what_came_back(self, tmp_path: Path) -> None: manifest = fetch(_config(), out=tmp_path, client=FakeClient()) @@ -259,7 +673,12 @@ def test_it_records_what_was_asked_and_what_came_back(self, tmp_path: Path) -> N def test_the_config_travels_with_the_data(self, tmp_path: Path) -> None: """A downloaded dataset has to be self-describing.""" source = ROOT / "datasets" / "pnr_2019.toml" - manifest = fetch(source, out=tmp_path, client=FakeClient()) + # The shipped config has a 30 km radius, and it is now cut exactly. + manifest = fetch( + source, + out=tmp_path, + client=FakeClient(inventory=_inventory_around_the_event()), + ) assert manifest["config"] == source.read_text() def test_it_is_written_beside_the_data(self, tmp_path: Path) -> None: