upd

src/qseek/apps/qseek.py

@@ -218,20 +218,28 @@ async def run() -> None:
             search.data_provider.prepare(search.stations)
 
             async def extract() -> None:
+                for magnitude in search.magnitudes:
+                    await magnitude.prepare(search.octree, search.stations)
+
                 iterator = asyncio.as_completed(
                     tuple(
                         search.add_magnitude_and_features(detection)
                         for detection in search._detections
                     )
                 )
+
                 for result in track(
                     iterator,
                     description="Extracting features",
                     total=search._detections.n_detections,
                 ):
-                    detection = await result
-                    await detection.save(update=True)
+                    event = await result
+                    if event.magnitudes:
+                        for mag in event.magnitudes:
+                            print(f"{mag.magnitude} {mag.average:.2f}±{mag.error:.2f}")
+                        print("--")
 
+                await search._detections.save()
                 await search._detections.export_detections(
                     jitter_location=search.octree.smallest_node_size()
                 )

src/qseek/magnitudes/__init__.py

@@ -3,7 +3,10 @@
 from pydantic import Field
 
 # Has to be imported to register as subclass
-from qseek.magnitudes import local_magnitude  # noqa: F401
+from qseek.magnitudes import (
+    local_magnitude,  # noqa: F401
+    moment_magnitude,  # noqa: F401
+)
 from qseek.magnitudes.base import EventMagnitude, EventMagnitudeCalculator
 
 EventMagnitudeType = Annotated[

src/qseek/magnitudes/base.py

@@ -16,6 +16,8 @@
     from qseek.models.station import Stations
     from qseek.octree import Octree
 
+KM = 1e3
+
 
 class EventMagnitude(BaseModel):
     magnitude: Literal["EventMagnitude"] = "EventMagnitude"
@@ -48,6 +50,49 @@ def csv_row(self) -> dict[str, float]:
             "error": self.error,
         }
 
+    def plot(self) -> None:
+        import matplotlib.pyplot as plt
+        from matplotlib.ticker import FuncFormatter
+
+        station_distances_hypo = np.array(
+            [sta.distance_hypo for sta in self.station_magnitudes]
+        )
+
+        fig = plt.figure()
+        ax = fig.gca()
+        ax.errorbar(
+            station_distances_hypo,
+            self.magnitudes,
+            yerr=[sta.magnitude_error for sta in self.station_magnitudes],
+            marker="o",
+            mec="k",
+            mfc="k",
+            ms=2,
+            ecolor=(0.0, 0.0, 0.0, 0.1),
+            capsize=1,
+            ls="none",
+        )
+        ax.axhline(
+            self.average,
+            color="k",
+            linestyle="dotted",
+            alpha=0.5,
+            label=rf"Median $M_L$ {self.average:.2f} $\pm${self.error:.2f}",
+        )
+        ax.set_xlabel("Distance to Hypocenter [km]")
+        ax.set_ylabel("$M_L$")
+        ax.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: x / KM))
+        ax.grid(alpha=0.3)
+        ax.legend(title=f"Estimator: {self.model}", loc="lower right")
+        ax.text(
+            0.05,
+            0.05,
+            f"{self.n_observations} Stations",
+            transform=ax.transAxes,
+            alpha=0.5,
+        )
+        plt.show()
+
 
 class EventMagnitudeCalculator(BaseModel):
     magnitude: Literal["MagnitudeCalculator"] = "MagnitudeCalculator"
@@ -93,7 +138,7 @@ async def prepare(
         Raises:
             NotImplementedError: This method must be implemented by subclasses.
         """
-        raise NotImplementedError
+        ...
 
 
 class StationAmplitudes(NamedTuple):

src/qseek/magnitudes/local_magnitude.py

@@ -8,7 +8,6 @@
 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 (
@@ -73,10 +72,9 @@ async def from_model(
                 continue
             self.station_magnitudes.append(station_magnitude)
 
-        magnitudes = self.magnitudes
-
-        self.average = float(np.median(magnitudes))
-        self.error = float(median_abs_deviation(magnitudes))
+        median = np.median(self.magnitudes)
+        self.average = float(median)
+        self.error = float(np.median(np.abs(self.magnitudes - median)))  # MAD
 
         return self
 
@@ -85,7 +83,7 @@ def magnitudes(self) -> np.ndarray:
         return np.array([sta.magnitude for sta in self.station_magnitudes])
 
     @property
-    def n_observations(self) -> int:
+    def n_stations(self) -> int:
         return len(self.station_magnitudes)
 
     def csv_row(self) -> dict[str, float]:
@@ -130,7 +128,7 @@ def plot(self) -> None:
         ax.text(
             0.05,
             0.05,
-            f"{self.n_observations} Stations",
+            f"{self.n_stations} Stations",
             transform=ax.transAxes,
             alpha=0.5,
         )
@@ -140,9 +138,18 @@ def plot(self) -> None:
 class LocalMagnitudeExtractor(EventMagnitudeCalculator):
     magnitude: Literal["LocalMagnitude"] = "LocalMagnitude"
 
-    seconds_before: PositiveFloat = 10.0
-    seconds_after: PositiveFloat = 10.0
-    padding_seconds: PositiveFloat = 10.0
+    seconds_before: PositiveFloat = Field(
+        default=10.0,
+        description="Seconds before first phase arrival to extract.",
+    )
+    seconds_after: PositiveFloat = Field(
+        default=10.0,
+        description="Seconds after last phase arrival to extract.",
+    )
+    padding_seconds: PositiveFloat = Field(
+        default=10.0,
+        description="Seconds padding before and after the extraction window.",
+    )
     model: ModelName = Field(
         default="iaspei-southern-california",
         description="The estimator to use for calculating the local magnitude.",

src/qseek/magnitudes/local_magnitude_model.py

@@ -188,7 +188,7 @@ def get_station_magnitude(
         return StationLocalMagnitude(
             station_nsl=sta.station_nsl,
             magnitude=magnitude,
-            magnitude_error=(magnitude_error_upper - magnitude_error_lower) / 2,
+            magnitude_error=(magnitude_error_upper + abs(magnitude_error_lower)) / 2,
             peak_amp=sta.peak,
             distance_epi=sta.distance_epi,
             distance_hypo=sta.distance_hypo,