From c3efe8b7c8198d37edaf1f301aa9d23c48348efa Mon Sep 17 00:00:00 2001 From: miili Date: Fri, 15 Dec 2023 00:05:36 +0100 Subject: [PATCH] uncertainty: adding uncertainty --- src/qseek/corrections/simple.py | 32 +++++++++++++ src/qseek/event_uncertainty.py | 47 ------------------- src/qseek/models/detection.py | 2 + src/qseek/models/detection_uncertainty.py | 57 +++++++++++++++++++++++ src/qseek/search.py | 6 +++ 5 files changed, 97 insertions(+), 47 deletions(-) create mode 100644 src/qseek/corrections/simple.py delete mode 100644 src/qseek/event_uncertainty.py create mode 100644 src/qseek/models/detection_uncertainty.py diff --git a/src/qseek/corrections/simple.py b/src/qseek/corrections/simple.py new file mode 100644 index 00000000..f8766c0d --- /dev/null +++ b/src/qseek/corrections/simple.py @@ -0,0 +1,32 @@ +from typing import Iterable, Literal + +import numpy as np + +from qseek.corrections.base import NSL, StationCorrections +from qseek.utils import PhaseDescription + + +class SimpleCorrections(StationCorrections): + corrections: Literal["SimpleCorrections"] = "SimpleCorrections" + + stations: dict[NSL, dict[PhaseDescription, float]] = {} + + @property + def n_stations(self) -> int: + return len(self.stations) + + def get_delay(self, station_nsl: NSL, phase: PhaseDescription) -> float: + if station_nsl not in self.stations: + return 0.0 + if phase not in self.stations[station_nsl]: + return 0.0 + return self.stations[station_nsl][phase] + + def get_delays( + self, + station_nsls: Iterable[NSL], + phase: PhaseDescription, + ) -> np.ndarray: + return np.array( + [self.get_delay(station_nsl, phase) for station_nsl in station_nsls] + ) diff --git a/src/qseek/event_uncertainty.py b/src/qseek/event_uncertainty.py deleted file mode 100644 index 4979f78a..00000000 --- a/src/qseek/event_uncertainty.py +++ /dev/null @@ -1,47 +0,0 @@ -from typing import TYPE_CHECKING, Literal - -import numpy as np -from pydantic import BaseModel -from typing_extensions import Self - -if TYPE_CHECKING: - from qseek.octree import Node, Octree - - -# Equivalent to one standard deviation -THRESHOLD = 1.0 / np.sqrt(np.e) - - -class EventUncertainty(BaseModel): - measure: Literal["standard_deviation"] = "standard_deviation" - - east_uncertainties: tuple[float, float] - north_uncertainties: tuple[float, float] - depth_uncertainties: tuple[float, float] - - @classmethod - def from_event( - cls, source_node: Node, octree: Octree, width: float = THRESHOLD - ) -> Self: - """ - Calculate the uncertainty of an event detection. - - Args: - event: The event detection to calculate the uncertainty for. - octree: The octree to use for the calculation. - - Returns: - The calculated uncertainty. - """ - nodes = octree.get_nodes(semblance_threshold=width) - vicinity_coords = np.array( - [(node.east, node.north, node.depth) for node in nodes] - ) - eastings = vicinity_coords[:, 0] - source_node.east - northings = vicinity_coords[:, 1] - source_node.north - depths = vicinity_coords[:, 2] - source_node.depth - return cls( - east_uncertainties=(float(np.min(eastings)), float(np.max(eastings))), - north_uncertainties=(float(np.min(northings)), float(np.max(northings))), - depth_uncertainties=(float(np.min(depths)), float(np.max(depths))), - ) diff --git a/src/qseek/models/detection.py b/src/qseek/models/detection.py index 83d2b086..97b0be7e 100644 --- a/src/qseek/models/detection.py +++ b/src/qseek/models/detection.py @@ -27,6 +27,7 @@ from qseek.console import console from qseek.features import EventFeaturesTypes, ReceiverFeaturesTypes from qseek.images.images import ImageFunctionPick +from qseek.models.detection_uncertainty import DetectionUncertainty from qseek.models.location import Location from qseek.models.station import Station, Stations from qseek.stats import Stats @@ -346,6 +347,7 @@ class EventDetection(Location): description="Number of stations in the detection.", ) + uncertainty: DetectionUncertainty | None = None features: EventFeatures = EventFeatures() _receivers: EventReceivers | None = PrivateAttr(None) diff --git a/src/qseek/models/detection_uncertainty.py b/src/qseek/models/detection_uncertainty.py new file mode 100644 index 00000000..f08ce7ff --- /dev/null +++ b/src/qseek/models/detection_uncertainty.py @@ -0,0 +1,57 @@ +from typing import TYPE_CHECKING + +import numpy as np +from pydantic import BaseModel, Field +from typing_extensions import Self + +if TYPE_CHECKING: + from qseek.octree import Node, Octree + + +# Equivalent to one standard deviation +THRESHOLD = 1.0 / np.sqrt(np.e) + + +class DetectionUncertainty(BaseModel): + east_uncertainties: tuple[float, float] = Field( + ..., + description="Uncertainty in east direction in [m].", + ) + north_uncertainties: tuple[float, float] = Field( + ..., + description="Uncertainty in north direction in [m].", + ) + depth_uncertainties: tuple[float, float] = Field( + ..., + description="Uncertainty in depth in [m].", + ) + + @classmethod + def from_event( + cls, source_node: Node, octree: Octree, width: float = THRESHOLD + ) -> Self: + """ + Calculate the uncertainty of an event detection. + + Args: + event: The event detection to calculate the uncertainty for. + octree: The octree to use for the calculation. + + Returns: + The calculated uncertainty. + """ + nodes = octree.get_nodes(semblance_threshold=width) + vicinity_coords = np.array( + [(node.east, node.north, node.depth) for node in nodes] + ) + relative_node_offsets = vicinity_coords - np.array( + [source_node.east, source_node.north, source_node.depth] + ) + min_offsets = np.min(relative_node_offsets, axis=0) + max_offsets = np.max(relative_node_offsets, axis=0) + + return cls( + east_uncertainties=(float(min_offsets[0]), float(max_offsets[0])), + north_uncertainties=(float(min_offsets[1]), float(max_offsets[1])), + depth_uncertainties=(float(min_offsets[2]), float(max_offsets[2])), + ) diff --git a/src/qseek/search.py b/src/qseek/search.py index 9dd3c8d8..ea18666b 100644 --- a/src/qseek/search.py +++ b/src/qseek/search.py @@ -25,6 +25,7 @@ from qseek.images.images import ImageFunctions, WaveformImages from qseek.models import Stations from qseek.models.detection import EventDetection, EventDetections, PhaseDetection +from qseek.models.detection_uncertainty import DetectionUncertainty from qseek.models.semblance import Semblance from qseek.octree import NodeSplitError, Octree from qseek.signals import Signal @@ -684,6 +685,11 @@ async def search( phase_arrivals=phase_detections, ) + detection.uncertainty = DetectionUncertainty.from_event( + source_node=source_node, + octree=octree, + ) + detections.append(detection) return detections, semblance.get_trace()