Add median magnitude and quantity field to MomentMagnitude class

src/qseek/magnitudes/base.py

@@ -26,6 +26,10 @@ class EventMagnitude(BaseModel):
         default=0.0,
         description="Average local magnitude.",
     )
+    median: float = Field(
+        default=0.0,
+        description="Average local magnitude.",
+    )
     error: float = Field(
         default=0.0,
         description="Average error of local magnitude.",

src/qseek/magnitudes/moment_magnitude.py

@@ -21,7 +21,14 @@
     PeakAmplitudesStore,
     PeakAmplitudeStoreCache,
 )
-from qseek.utils import CACHE_DIR, NSL, ChannelSelector, ChannelSelectors, Range
+from qseek.utils import (
+    CACHE_DIR,
+    NSL,
+    ChannelSelector,
+    ChannelSelectors,
+    MeasurementUnit,
+    Range,
+)
 
 if TYPE_CHECKING:
     from pyrocko.squirrel import Squirrel
@@ -120,6 +127,7 @@ def filter_receivers_by_range(
 
 
 class StationMomentMagnitude(NamedTuple):
+    quantity: MeasurementUnit
     distance_epi: float
     magnitude: float
     error: float
@@ -130,23 +138,14 @@ class StationMomentMagnitude(NamedTuple):
 class MomentMagnitude(EventMagnitude):
     magnitude: Literal["MomentMagnitude"] = "MomentMagnitude"
 
-    average: float = Field(
-        default=0.0,
-        description="Average moment magnitude.",
-    )
-    error: float = Field(
-        default=0.0,
-        description="Average error of moment magnitude.",
-    )
-    std: float = Field(
-        default=0.0,
-        description="Standard deviation of moment magnitude.",
-    )
-
     stations_magnitudes: list[StationMomentMagnitude] = Field(
         default_factory=list,
         description="The station moment magnitudes.",
     )
+    quantity: MeasurementUnit = Field(
+        default="displacement",
+        description="The quantity of the traces.",
+    )
 
     @property
     def m0(self) -> float:
@@ -209,6 +208,7 @@ async def add_traces(
                 error_lower = error_upper
 
             station_magnitude = StationMomentMagnitude(
+                quantity=store.quantity,
                 distance_epi=station.distance_epi,
                 magnitude=magnitude,
                 error=(error_upper + abs(error_lower)) / 2,
@@ -222,7 +222,8 @@ async def add_traces(
         magnitudes = np.array([sta.magnitude for sta in self.stations_magnitudes])
         median = np.median(magnitudes)
 
-        self.average = float(median)
+        self.median = float(median)
+        self.average = float(np.mean(magnitudes))
         self.error = float(np.median(np.abs(magnitudes - median)))  # MAD
 
 
@@ -281,7 +282,6 @@ async def add_magnitude(
         event: EventDetection,
     ) -> None:
         moment_magnitude = MomentMagnitude()
-
         receivers = list(event.receivers)
         for store, definition in zip(self._stores, self.models, strict=True):
             store_receivers = definition.filter_receivers_by_nsl(receivers)
@@ -308,9 +308,10 @@ async def add_magnitude(
                 demean=True,
                 seconds_fade=self.padding_seconds,
                 cut_off_fade=False,
+                remove_clipped=True,
             )
             if not traces:
-                return
+                continue
 
             for tr in traces:
                 if store.frequency_range.min != 0.0:

src/qseek/models/detection.py

@@ -22,7 +22,6 @@
 from pyrocko import io
 from pyrocko.gui import marker
 from pyrocko.model import Event, dump_events
-from pyrocko.squirrel.error import NotAvailable
 from rich.table import Table
 from typing_extensions import Self
 
@@ -320,15 +319,27 @@ async def get_waveforms_restituted(
         )
         traces = filter_clipped_traces(traces) if remove_clipped else traces
 
+        tmin = min(tr.tmin for tr in traces)
+        tmax = max(tr.tmax for tr in traces)
+        responses = await asyncio.to_thread(
+            squirrel.get_responses,
+            tmin=tmin,
+            tmax=tmax,
+            codes=[tr.nslc_id for tr in traces],
+        )
+
+        def get_response(tr: Trace) -> Any:
+            for response in responses:
+                if response.codes[:4] == tr.nslc_id:
+                    return response
+            raise ValueError(f"cannot find response for {tr.nslc_id}")
+
         restituted_traces = []
         for tr in traces:
             try:
-                response = squirrel.get_response(
-                    tmin=tr.tmin,
-                    tmax=tr.tmax,
-                    codes=[tr.nslc_id],
-                )
-            except NotAvailable:
+                response = get_response(tr)
+            except ValueError:
+                logger.debug("cannot find response for %s", ".".join(tr.nslc_id))
                 continue
 
             restituted_traces.append(

src/qseek/octree.py

@@ -66,9 +66,11 @@ class Node:
     north: float
     depth: float
     size: float
+    level: int
     semblance: float = 0.0
 
     tree: Octree | None = None
+    parent: Node | None = None
     children: tuple[Node, ...] = ()
 
     _hash: bytes | None = None
@@ -93,6 +95,8 @@ def split(self) -> tuple[Node, ...]:
                     depth=self.depth + depth * half_size / 2,
                     size=half_size,
                     tree=self.tree,
+                    parent=self,
+                    level=self.level + 1,
                 )
                 for east in (-1, 1)
                 for north in (-1, 1)
@@ -311,7 +315,7 @@ def get_root_nodes(self, length: float) -> list[Node]:
         depth_nodes = np.arange(ext_depth // ln) * ln + ln / 2 + self.depth_bounds.min
 
         return [
-            Node(east=east, north=north, depth=depth, size=ln, tree=self)
+            Node(east=east, north=north, depth=depth, size=ln, tree=self, level=0)
             for east in east_nodes
             for north in north_nodes
             for depth in depth_nodes