moment magnitude

src/qseek/magnitudes/base.py

@@ -7,13 +7,12 @@
 from pyrocko.trace import Trace
 from typing_extensions import Self
 
-from qseek.models.detection import Receiver
 from qseek.utils import NSL
 
 if TYPE_CHECKING:
     from pyrocko.squirrel import Squirrel
 
-    from qseek.models.detection import EventDetection
+    from qseek.models.detection import EventDetection, Receiver
     from qseek.models.station import Stations
     from qseek.octree import Octree
 

src/qseek/magnitudes/local_magnitude.py

@@ -8,6 +8,7 @@
 import numpy as np
 from matplotlib.ticker import FuncFormatter
 from pydantic import Field, PositiveFloat, PrivateAttr, model_validator
+from scipy.stats import median_abs_deviation
 from typing_extensions import Self
 
 from qseek.magnitudes.base import (
@@ -44,19 +45,12 @@ class LocalMagnitude(EventMagnitude):
     station_magnitudes: list[StationLocalMagnitude] = []
     average: float = Field(
         default=0.0,
-        description="Average local magnitude.",
-    )
-    average_weighted: float = Field(
-        default=0.0,
-        description="Weighted average local magnitude.",
-    )
-    median: float = Field(
-        default=0.0,
-        description="Median local magnitude.",
+        description="The network's local magnitude, as median of"
+        " all station magnitudes.",
     )
     error: float = Field(
         default=0.0,
-        description="Average error of local magnitude.",
+        description="Average error of local magnitude, as median absolute deviation.",
     )
 
     @classmethod
@@ -80,15 +74,9 @@ async def from_model(
             self.station_magnitudes.append(station_magnitude)
 
         magnitudes = self.magnitudes
-        magnitude_errors = np.array(
-            [sta.magnitude_error for sta in self.station_magnitudes]
-        )
-        weights = 1.0 / magnitude_errors
 
-        self.average = float(np.average(magnitudes))
-        self.average_weighted = float(np.average(self.magnitudes, weights=weights))
-        self.median = float(np.median(magnitudes))
-        self.error = float(np.average(magnitudes + magnitude_errors)) - self.average
+        self.average = float(np.median(magnitudes))
+        self.error = float(median_abs_deviation(magnitudes))
 
         return self
 
@@ -127,26 +115,12 @@ def plot(self) -> None:
             capsize=1,
             ls="none",
         )
-        ax.axhline(
-            self.median,
-            color="k",
-            linestyle="--",
-            alpha=0.5,
-            label=f"Median $M_L$ {self.median:.2f}",
-        )
         ax.axhline(
             self.average,
             color="k",
             linestyle="dotted",
             alpha=0.5,
-            label=rf"Average $M_L$ {self.average:.2f} $\pm${self.error:.2f}",
-        )
-        ax.axhline(
-            self.average_weighted,
-            color="k",
-            linestyle="-",
-            alpha=0.5,
-            label=f"Weighted Average $M_L$ {self.average_weighted:.2f}",
+            label=rf"Median $M_L$ {self.average:.2f} $\pm${self.error:.2f}",
         )
         ax.set_xlabel("Distance to Hypocenter [km]")
         ax.set_ylabel("$M_L$")

src/qseek/magnitudes/moment_magnitude.py

@@ -3,7 +3,7 @@
 import itertools
 import logging
 from pathlib import Path
-from typing import TYPE_CHECKING, Any, Literal, NamedTuple
+from typing import TYPE_CHECKING, Any, Iterable, Literal, NamedTuple
 
 import numpy as np
 from pydantic import DirectoryPath, Field, PositiveFloat, PrivateAttr
@@ -20,7 +20,7 @@
     PeakAmplitudesStore,
     PeakAmplitudeStoreCache,
 )
-from qseek.utils import CACHE_DIR, NSL, ChannelSelector, ChannelSelectors
+from qseek.utils import CACHE_DIR, NSL, ChannelSelector, ChannelSelectors, Range
 
 if TYPE_CHECKING:
     from pyrocko.squirrel import Squirrel
@@ -56,7 +56,7 @@ def norm_traces(traces: list[Trace]) -> np.ndarray:
 
 
 class PeakAmplitudeDefinition(PeakAmplitudesBase):
-    nsl_id: str | None = Field(
+    nsl_id: list[str] | None = Field(
         default=None,
         pattern=r"^[A-Z0-9]{2}\.[A-Z0-9]{0,5}?\.[A-Z0-9]{0,2}?$",
         description="Network, station, location id.",
@@ -65,8 +65,12 @@ class PeakAmplitudeDefinition(PeakAmplitudesBase):
         default="absolute",
         description="The peak amplitude to use.",
     )
+    station_epicentral_range: Range = Field(
+        default=Range(min=1 * KM, max=100 * KM),
+        description="The epicentral distance range of the stations.",
+    )
 
-    def filter_receivers(self, receivers: list[Receiver]) -> list[Receiver]:
+    def filter_receiver_by_nsl(self, receivers: Iterable[Receiver]) -> set[Receiver]:
         """
         Filters the list of receivers based on the NSL ID.
 
@@ -77,9 +81,36 @@ def filter_receivers(self, receivers: list[Receiver]) -> list[Receiver]:
             list[Receiver]: The filtered list of receivers.
         """
         if self.nsl_id is None:
-            return receivers
-        nsl = NSL.parse(self.nsl_id)
-        return [rcv for rcv in receivers if rcv.nsl.match(nsl)]
+            return set(receivers)
+
+        matched_receivers = []
+        for nsl in self.nsl_id:
+            nsl = NSL.parse(nsl)
+            matched_receivers.append([rcv for rcv in receivers if rcv.nsl.match(nsl)])
+
+        return set(itertools.chain.from_iterable(matched_receivers))
+
+    def filter_receiver_by_range(
+        self,
+        receivers: Iterable[Receiver],
+        event: EventDetection,
+    ) -> set[Receiver]:
+        """
+        Filters the list of receivers based on the distance range.
+
+        Args:
+            receivers (Iterable[Receiver]): The list of receivers to filter.
+            event (EventDetection): The event detection object.
+
+        Returns:
+            set[Receiver]: The filtered set of receivers.
+        """
+        receivers = [
+            rcv
+            for rcv in receivers
+            if self.station_epicentral_range.inside(rcv.surface_distance_to(event))
+        ]
+        return set(receivers)
 
 
 class StationMomentMagnitude(NamedTuple):
@@ -235,14 +266,14 @@ def model_post_init(self, __context: Any) -> None:
 
     async def prepare(self, octree: Octree, stations: Stations) -> None:
         logger.info("Preparing peak amplitude stores...")
-        root_nodes = octree.get_root_nodes(octree.size_initial)
-        depths = np.array([node.as_location().effective_depth for node in root_nodes])
+        octree_depth = octree.location.effective_depth
 
         for store in self._stores:
-            for depth in depths:
-                if store.has_depth(depth):
-                    continue
-                await store.fill_source_depth(depth)
+            await store.fill_source_depth_range(
+                depth_min=octree.depth_bounds.min + octree_depth,
+                depth_max=octree.depth_bounds.max + octree_depth,
+                depth_delta=octree.size_initial,
+            )
 
     async def add_magnitude(
         self,
@@ -253,19 +284,20 @@ async def add_magnitude(
 
         receivers = list(event.receivers)
         for store, definition in zip(self._stores, self.models, strict=True):
-            store_receivers = definition.filter_receivers(receivers)
+            store_receivers = definition.filter_receiver_by_nsl(receivers)
+            if not store_receivers:
+                continue
             for rcv in store_receivers:
                 receivers.remove(rcv)
-
+            store_receivers = definition.filter_receiver_by_range(
+                store_receivers, event
+            )
             if not store_receivers:
+                logger.debug("No receivers in range for peak amplitude %s", store.id)
+                continue
+            if not store.source_depth_range.inside(event.effective_depth):
+                logger.info("Event depth outside of store depth range.")
                 continue
-
-            if definition.max_distance is not None:
-                store_receivers = [
-                    rcv
-                    for rcv in store_receivers
-                    if rcv.surface_distance_to(event) <= definition.max_distance
-                ]
 
             traces = await event.receivers.get_waveforms_restituted(
                 squirrel,