uncertainty: adding uncertainty

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]
+        )

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)

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])),
+        )

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()