From 3dffe24c251147ccc902496064ec7fc7aa2bbc0f Mon Sep 17 00:00:00 2001 From: Marius Isken Date: Mon, 15 Jan 2024 16:00:17 +0100 Subject: [PATCH] Add median magnitude and quantity field to MomentMagnitude class --- src/qseek/magnitudes/base.py | 4 +++ src/qseek/magnitudes/moment_magnitude.py | 35 ++++++++++++------------ src/qseek/models/detection.py | 25 ++++++++++++----- src/qseek/octree.py | 6 +++- 4 files changed, 45 insertions(+), 25 deletions(-) diff --git a/src/qseek/magnitudes/base.py b/src/qseek/magnitudes/base.py index b2a20063..71567758 100644 --- a/src/qseek/magnitudes/base.py +++ b/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.", diff --git a/src/qseek/magnitudes/moment_magnitude.py b/src/qseek/magnitudes/moment_magnitude.py index 36f5104e..4cdd9511 100644 --- a/src/qseek/magnitudes/moment_magnitude.py +++ b/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: diff --git a/src/qseek/models/detection.py b/src/qseek/models/detection.py index af2a09d7..f4ba512b 100644 --- a/src/qseek/models/detection.py +++ b/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( diff --git a/src/qseek/octree.py b/src/qseek/octree.py index cd069816..73566de9 100644 --- a/src/qseek/octree.py +++ b/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