distance weights: using FWHM scale, improved parameter naming

Author: Marius Isken

pyproject.toml

@@ -3,7 +3,7 @@ requires = [
     "wheel",
     "setuptools>=80",
     "numpy>=2",
-    "setuptools_scm[toml]>=8",
+    "setuptools_scm>=8",
     "simde-py==0.8.4rc1",
 ]
 build-backend = "setuptools.build_meta"
@@ -32,8 +32,8 @@ keywords = [
 dependencies = [
     "numpy>=2.3",
     "scipy>=1.16",
-    "pyrocko>=2024.12",
-    "seisbench>=0.9",
+    "pyrocko>=2025.12",
+    "seisbench>=0.10",
     "pydantic>=2.12",
     "aiohttp>=3.8",
     "aiohttp_cors>=0.7.0",
@@ -134,6 +134,8 @@ addopts = "--cov=qseek --cov-report html"
 [tool.cibuildwheel.macos]
 before-build = "brew install libomp"
 
+[tool.setuptools_scm]
+
 [dependency-groups]
 dev = [
     "pre-commit>=4.5",

src/qseek/distance_weights.py

@@ -21,18 +21,6 @@
 logger = logging.getLogger(__name__)
 
 
-def weights_gaussian(
-    distances: np.ndarray,
-    radius: float,
-    exponent: float = 2.0,
-    normalize: bool = True,
-) -> np.ndarray:
-    weights = np.exp(-((distances / radius) ** exponent))
-    if normalize:
-        weights /= weights.max()
-    return weights
-
-
 def weights_exponential(
     distances: np.ndarray,
     radius: float,
@@ -46,10 +34,9 @@ def weights_exponential(
     return weights
 
 
-def weights_gaussian_min_stations(
+def weights_gaussian(
     distances: np.ndarray,
-    radius: float,
-    exponent: float = 2.0,
+    distance_taper: float,
     required_stations: int = 4,
     waterlevel: float = 0.0,
 ) -> np.ndarray:
@@ -60,40 +47,35 @@ def weights_gaussian_min_stations(
     sorted_distances = np.sort(distances, axis=1)
     threshold_distance = sorted_distances[:, required_stations - 1, np.newaxis]
 
-    weights = np.exp(-(((distances - threshold_distance) / radius) ** exponent))
+    weights = np.exp(
+        -(((distances - threshold_distance) ** 2) / (2 * (distance_taper / 2) ** 2))
+    )
     weights[distances <= threshold_distance] = 1.0
     if waterlevel > 0.0:
         weights = (1 - waterlevel) * weights + waterlevel
     return weights
 
 
 class DistanceWeights(BaseModel):
-    radius_meters: PositiveFloat | Literal["mean_interstation"] = Field(
+    distance_taper: PositiveFloat | Literal["mean_interstation"] = Field(
         default="mean_interstation",
-        description="Cutoff distance for the spatial decay function in meters."
-        " 'mean_interstation' uses the mean interstation distance for the radius."
-        " Default is 'mean_interstation'.",
+        description="Taper distance for Gaussian the distance weighting function"
+        " in meters. 'mean_interstation' uses twice the mean interstation distance for"
+        " the radius. Default is 'mean_interstation'.",
     )
-    min_required_stations: PositiveInt = Field(
+    required_closest_stations: PositiveInt = Field(
         default=4,
-        description="Minimum number of stations to assign full weight in the"
-        " exponential decay function. Default is 4.",
-    )
-    exponent: float = Field(
-        default=2.0,
-        description="Exponent of the spatial decay function. For 'gaussian' decay an"
-        " exponent of 0.5 is recommended. Default is 2.",
-        ge=0.0,
+        description="Number of stations to assign full weight in the"
+        " spatial weighting function, only more distant stations are tapered with a"
+        " Gaussian decay. This ensures that the closest _N_ stations have an equal and"
+        " the highest contribution to the detection and localization. Default is 4.",
     )
     waterlevel: float = Field(
         default=0.0,
         ge=0.0,
         le=1.0,
-        description="Waterlevel for the exponential decay function. Default is 0.0.",
-    )
-    normalize: bool = Field(
-        default=True,
-        description="Normalize the weights to the range [0, 1]. Default is True.",
+        description="Stations outside the taper are lifted by this fraction. "
+        "Default is 0.0.",
     )
 
     _node_lut: ArrayLRUCache[bytes] = PrivateAttr()
@@ -110,14 +92,17 @@ def get_distances(self, nodes: Sequence[Node]) -> np.ndarray:
         )
 
     def prepare(self, stations: StationInventory, octree: Octree) -> None:
-        logger.info("preparing distance weights")
-
-        if self.radius_meters == "mean_interstation":
-            self.radius_meters = stations.mean_interstation_distance()
+        if self.distance_taper == "mean_interstation":
+            self.distance_taper = 2 * stations.mean_interstation_distance()
             logger.info(
-                "using mean interstation distance as radius: %g m",
-                self.radius_meters,
+                "using 2x mean interstation distance as distance taper: %g m",
+                self.distance_taper,
             )
+        logger.info(
+            "distance weighting uses %d closest stations and a taper of %g m",
+            self.required_closest_stations,
+            self.distance_taper,
+        )
 
         self._stations = StationList.from_inventory(stations)
         self._node_lut = ArrayLRUCache(name="distance_weights", short_name="DW")
@@ -144,11 +129,10 @@ def get_node_weights(self, node: Node, stations: list[Station]) -> np.ndarray:
             self.fill_lut([node])
             return self.get_node_weights(node, stations)
 
-        return weights_gaussian_min_stations(
+        return weights_gaussian(
             distances,
-            required_stations=self.min_required_stations,
-            radius=self.radius_meters,
-            exponent=self.exponent,
+            required_stations=self.required_closest_stations,
+            distance_taper=self.distance_taper,
             waterlevel=self.waterlevel,
         )
 
@@ -163,11 +147,10 @@ async def get_weights(
 
         try:
             distances = [node_lut[node.hash][station_indices] for node in nodes]
-            return weights_gaussian_min_stations(
+            return weights_gaussian(
                 np.array(distances),
-                required_stations=self.min_required_stations,
-                radius=self.radius_meters,
-                exponent=self.exponent,
+                required_stations=self.required_closest_stations,
+                distance_taper=self.distance_taper,
                 waterlevel=self.waterlevel,
             )
         except KeyError:

src/qseek/models/station.py

@@ -41,9 +41,21 @@ def __get_pydantic_core_schema__(
 
 
 class Station(Location):
-    network: str = Field(..., max_length=2)
-    station: str = Field(..., max_length=5)
-    location: str = Field(default="", max_length=2)
+    network: str = Field(
+        ...,
+        description="Seismic network code",
+        max_length=2,
+    )
+    station: str = Field(
+        ...,
+        description="Station code",
+        max_length=5,
+    )
+    location: str = Field(
+        default="",
+        description="Location code",
+        max_length=2,
+    )
 
     @classmethod
     def from_pyrocko_station(cls, station: PyrockoStation) -> Station:

src/qseek/search.py

@@ -173,7 +173,8 @@ def add_processed_batch(
 
     def log(self) -> None:
         log_str = (
-            f"{self.batch_count + 1}/{self.batch_count_total or '?'} {self.batch_time}"
+            f"{self.batch_count + 1}/{self.batch_count_total or '?'} "
+            f"{str(self.batch_time)[:22]}"
         )
         logger.info(
             "%s%% processed - batch %s in %s",