Make source-to-site distance a configurable registry, and specify the magnitude calculation

src/specmod/config/__init__.py

@@ -16,6 +16,7 @@
     AcquireConfig,
     Config,
     FittingConfig,
+    GeometryConfig,
     ModelConfig,
     SmoothingConfig,
     SnrConfig,
@@ -29,6 +30,7 @@
     "AcquireConfig",
     "Config",
     "FittingConfig",
+    "GeometryConfig",
     "ModelConfig",
     "Provenance",
     "ResolvedConfig",

src/specmod/config/sections.py

@@ -66,7 +66,6 @@ class WindowsConfig:
     #: run used s=3.4; 2.9 is the shipped default and is kept as such.
     p_velocity: float = 5.9
     s_velocity: float = 2.9
-    distance_metric: Literal["repi", "rhyp"] = "repi"
 
     #: Used when an S pick is missing: s_time = p_time + emergency_ratio * (p - o).
     emergency_ratio: float = 1.7
@@ -239,7 +238,7 @@ class FittingConfig:
     #: How stations are weighted into the event value. The published choice is
     #: inverse hypocentral distance: the nearer station has less path, so less
     #: of its falloff can be attenuation. See ``specmod.staged.WEIGHT_MODELS``.
-    event_weighting: str = "inverse_hypocentral_distance"
+    event_weighting: str = "inverse_distance"
 
     #: Which channels contribute to the event value, as shell globs matched
     #: against the trace id and each of its SEED components — so ``"AQ07"``
@@ -271,12 +270,45 @@ class VizConfig:
     plot_columns: int = 3
 
 
+@dataclass(frozen=True, slots=True)
+class GeometryConfig:
+    """Source-to-site geometry.
+
+    Its own section because more than one stage needs it. Distance feeds the
+    ensemble weighting of the two-stage fit (:mod:`specmod.staged`) and the
+    geometric spreading a moment calculation corrects for, and a setting two
+    consumers each keep their own copy of is how the two come to disagree.
+
+    It lived in ``[windows]`` until there was a second reader, which was the
+    wrong home even then: cutting a window does not depend on how distance is
+    measured.
+    """
+
+    #: Which distance, resolved through :data:`specmod.distance.DISTANCE_MEASURES`.
+    #:
+    #: ``repi`` is the default and is the honest one wherever sensor depths are
+    #: not known. ``rhyp`` is built from the source depth and the station
+    #: *elevation*, so it assumes every sensor sits at the surface — for a
+    #: borehole deployment that is wrong by the burial depth, and nothing in
+    #: the metadata says so.
+    #:
+    #: The choice is not a detail at short range: on the PNR data the nearest
+    #: station is 0.89 km epicentral against 2.30 km hypocentral, a factor of
+    #: 2.57, while the farthest agree to 1.00 — so anything weighted by inverse
+    #: distance is most sensitive to it exactly where it matters most.
+    #:
+    #: ``rrup`` and ``rjb`` are registered and raise: both need a rupture
+    #: surface, and for a point source they degenerate to ``rhyp`` and ``repi``.
+    distance_measure: str = "repi"
+
+
 @dataclass(frozen=True, slots=True)
 class Config:
     """The whole resolved configuration."""
 
     acquire: AcquireConfig = field(default_factory=AcquireConfig)
     windows: WindowsConfig = field(default_factory=WindowsConfig)
+    geometry: GeometryConfig = field(default_factory=GeometryConfig)
     transform: TransformConfig = field(default_factory=TransformConfig)
     smoothing: SmoothingConfig = field(default_factory=SmoothingConfig)
     snr: SnrConfig = field(default_factory=SnrConfig)

src/specmod/distance.py

@@ -0,0 +1,190 @@
+"""Source-to-site distance, as a registry rather than a stat name.
+
+Which distance you mean is a modelling choice, and at short range it is not a
+small one. On the PNR data the nearest station is **0.89 km epicentral against
+2.30 km hypocentral** — a factor of 2.57 — while the farthest agree to 1.00.
+Anything weighted by inverse distance, or corrected for geometric spreading,
+therefore depends on the choice most strongly at exactly the station that
+matters most.
+
+Two are implemented here because they are the two a point source supports.
+Both read a value :func:`specmod.preprocess.set_stream_distance` has already
+computed:
+
+``repi``
+    Epicentral. Horizontal distance from the epicentre.
+``rhyp``
+    Hypocentral. Slant distance from the hypocentre.
+
+**Epicentral is the honest choice when sensor depths are unknown**, and that
+is more often than it sounds. ``rhyp`` is built from the source depth and the
+station *elevation*, which silently assumes every sensor sits at the surface.
+For a borehole deployment that is wrong by the burial depth, and nothing in
+the metadata announces it — the PNR inventory records channel ``depth`` as
+``123456.0``, a placeholder, so on that dataset ``rhyp`` is an assumption
+wearing a measurement's name.
+
+Finite-fault measures
+---------------------
+``Rrup`` (closest distance to the rupture surface) and ``Rjb`` (Joyner-Boore,
+closest horizontal distance to the surface projection of the rupture) are the
+measures ground-motion work generally wants, and they are **not implemented**
+— deliberately, rather than by omission.
+
+Both need a rupture *surface*: strike, dip, length, width and a hypocentre
+position on it. SpecMod carries a point source, so there is nothing to compute
+them from, and a version that quietly degenerated to ``rhyp`` and ``repi``
+would be worse than an error — those are exactly what `Rrup` and `Rjb` reduce
+to for a point source, so the substitution would be invisible in the output
+and wrong for any event large enough to warrant asking.
+
+They are registered all the same, raising with what they would need. A name
+that resolves to a clear failure is a better extension point than a name that
+does not resolve at all, and it puts the requirement where someone adding
+finite-fault support will read it.
+
+The registry is the same shape as :data:`specmod.transforms.ESTIMATORS`,
+:data:`specmod.core.noise.NOISE_MODELS` and
+:data:`specmod.staged.WEIGHT_MODELS`, so a study names a distance the way it
+names anything else and the choice travels with the resolved configuration.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING, Any, Protocol, runtime_checkable
+
+import numpy as np
+
+from .config import load_config
+
+if TYPE_CHECKING:  # pragma: no cover
+    from collections.abc import Sequence
+
+    from numpy.typing import NDArray
+
+__all__ = [
+    "DISTANCE_MEASURES",
+    "DistanceMeasure",
+    "Epicentral",
+    "FiniteFaultDistance",
+    "Hypocentral",
+    "get_distance_measure",
+    "resolve_distance_measure",
+]
+
+
+@runtime_checkable
+class DistanceMeasure(Protocol):
+    """One distance per channel, in kilometres."""
+
+    name: str
+
+    def distances(
+        self, spectra: Any, ids: Sequence[str]
+    ) -> NDArray[np.float64]: ...  # pragma: no cover
+
+
+@dataclass(frozen=True, slots=True)
+class _FromMeta:
+    """A distance already computed onto the trace metadata."""
+
+    key: str
+    name: str
+
+    def distances(self, spectra: Any, ids: Sequence[str]) -> NDArray[np.float64]:
+        out = np.empty(len(ids), dtype=np.float64)
+        for i, id in enumerate(ids):
+            meta = spectra[id].signal.meta
+            if self.key not in meta:
+                raise ValueError(
+                    f"{id} carries no {self.key!r}, so its {self.name} distance "
+                    f"is unknown. Set the geometry with "
+                    f"specmod.preprocess.set_stream_distance."
+                )
+            value = float(meta[self.key])
+            if value <= 0:
+                raise ValueError(
+                    f"{id} has {self.key}={value}, which is not a distance"
+                )
+            out[i] = value
+        return out
+
+
+@dataclass(frozen=True, slots=True)
+class Epicentral(_FromMeta):
+    key: str = "repi"
+    name: str = "epicentral"
+
+
+@dataclass(frozen=True, slots=True)
+class Hypocentral(_FromMeta):
+    key: str = "rhyp"
+    name: str = "hypocentral"
+
+
+@dataclass(frozen=True, slots=True)
+class FiniteFaultDistance:
+    """``Rrup`` and ``Rjb``: registered, and not implemented.
+
+    Raising here rather than omitting the name is the point. For a point source
+    these degenerate exactly to hypocentral and epicentral, so an
+    implementation that silently fell back would produce plausible numbers that
+    are wrong for any event big enough to justify asking for them.
+    """
+
+    name: str
+    needs: str
+
+    def distances(self, spectra: Any, ids: Sequence[str]) -> NDArray[np.float64]:
+        raise NotImplementedError(
+            f"{self.name} is not implemented. It needs {self.needs}, and "
+            f"SpecMod carries a point source — there is no rupture surface to "
+            f"measure from. For a point source {self.name} degenerates to "
+            f"{'hypocentral' if self.name == 'rrup' else 'epicentral'}; name "
+            f"that instead if it is what you mean, rather than getting it by "
+            f"accident."
+        )
+
+
+#: Registered distance measures, resolved by name from configuration.
+DISTANCE_MEASURES: dict[str, Any] = {
+    "repi": Epicentral,
+    "rhyp": Hypocentral,
+    "rrup": lambda: FiniteFaultDistance(
+        name="rrup", needs="a rupture surface — strike, dip, length and width"
+    ),
+    "rjb": lambda: FiniteFaultDistance(
+        name="rjb",
+        needs="the surface projection of a rupture — strike, dip, length and width",
+    ),
+}
+
+
+def get_distance_measure(name: str) -> DistanceMeasure:
+    """Resolve a registered measure by name."""
+    try:
+        factory = DISTANCE_MEASURES[name]
+    except KeyError:
+        raise ValueError(
+            f"Unknown distance measure {name!r}. "
+            f"Available: {sorted(DISTANCE_MEASURES)}."
+        ) from None
+    measure: DistanceMeasure = factory()
+    return measure
+
+
+def resolve_distance_measure(
+    measure: str | DistanceMeasure | None = None,
+) -> DistanceMeasure:
+    """A measure from a name, an instance, or the configuration.
+
+    ``None`` takes ``[geometry] distance_measure``, which is the project-wide
+    choice. It lived in ``[windows]`` and had no reader at all until this
+    module; cutting a window does not depend on how distance is measured.
+    """
+    if measure is None:
+        measure = str(load_config().config.geometry.distance_measure)
+    if isinstance(measure, str):
+        return get_distance_measure(measure)
+    return measure

src/specmod/staged.py

@@ -92,6 +92,7 @@
 import numpy as np
 
 from .config import load_config
+from .distance import DistanceMeasure, resolve_distance_measure
 from .fitting import FitSpectra
 
 if TYPE_CHECKING:  # pragma: no cover
@@ -134,35 +135,24 @@ class InverseDistance:
     modelling choice rather than a derivation, which is why this is a registry
     and not a hardcoded expression.
 
-    ``metric`` names a trace stat — ``rhyp`` for hypocentral, ``repi`` for
-    epicentral. A station missing it is a hard error rather than a silent
-    weight of zero: a geometry that was never set is a broken run, not a
-    station that should quietly stop contributing.
+    **Which distance is itself a choice**, and at short range not a small one:
+    see :mod:`specmod.distance`. ``measure=None`` takes the project-wide
+    setting, so a study that has decided on epicentral does not have to say so
+    again here.
     """
 
-    metric: str = "rhyp"
+    #: ``None`` means "whatever the configuration says". Distance is needed by
+    #: geometric spreading as well as by weighting, so the choice belongs in
+    #: one place rather than being restated per consumer.
+    measure: str | DistanceMeasure | None = None
     name: str = "inverse_distance"
 
     def weights(
         self, table: Any, spectra: Any, ids: Sequence[str]
     ) -> NDArray[np.float64]:
-        out = np.empty(len(ids), dtype=np.float64)
-        for i, id in enumerate(ids):
-            meta = spectra[id].signal.meta
-            if self.metric not in meta:
-                raise ValueError(
-                    f"{id} has no {self.metric!r}, so it cannot be weighted by "
-                    f"distance. Set the geometry with "
-                    f"specmod.preprocess.set_stream_distance, or choose a "
-                    f"weighting that does not need it."
-                )
-            distance = float(meta[self.metric])
-            if distance <= 0:
-                raise ValueError(
-                    f"{id} has {self.metric}={distance}, which is not a distance"
-                )
-            out[i] = 1.0 / distance
-        return out
+        distances = resolve_distance_measure(self.measure).distances(spectra, ids)
+        weights: NDArray[np.float64] = 1.0 / distances
+        return weights
 
 
 @dataclass(frozen=True, slots=True)
@@ -216,8 +206,13 @@ def weights(
 
 #: Registered weightings, resolved by name from ``[fitting] event_weighting``.
 WEIGHT_MODELS: dict[str, Any] = {
-    "inverse_hypocentral_distance": lambda: InverseDistance(metric="rhyp"),
-    "inverse_epicentral_distance": lambda: InverseDistance(metric="repi"),
+    # Follows the configured distance measure, so a project-wide choice is
+    # honoured in one place. The shipped default.
+    "inverse_distance": InverseDistance,
+    # And explicit spellings, for a study that wants to say which it used
+    # regardless of what the rest of the configuration says.
+    "inverse_hypocentral_distance": lambda: InverseDistance(measure="rhyp"),
+    "inverse_epicentral_distance": lambda: InverseDistance(measure="repi"),
     "uniform": Uniform,
     "inverse_variance": InverseVariance,
 }

studies/magna_2020_paper.toml

@@ -47,7 +47,6 @@ remove_response = false
 # Group velocities from Pechmann et al. (2007).
 p_velocity = 5.9
 s_velocity = 3.4
-distance_metric = "rhyp"
 # 20 s window opening at 80% of the elapsed Pg-Sg time, on the transverse.
 s_start_ratio = 0.8
 s_length = 20.0
@@ -129,3 +128,9 @@ motion = "velocity"
 [fitting]
 # "We use Powell's minimization technique (Powell, 1964; Press et al., 1997)."
 method = "powell"
+
+[geometry]
+# The published run used hypocentral distance. Kept as it was: this is the
+# record of what was done, not a recommendation. `repi` is the shipped default
+# because it does not assume sensor depths that are often unknown.
+distance_measure = "rhyp"