Skip to content

Commit

Permalink
Fix NSL parsing and update StationAmplitudes class
Browse files Browse the repository at this point in the history
  • Loading branch information
Marius Isken committed Jan 5, 2024
1 parent c4b1200 commit 670a525
Show file tree
Hide file tree
Showing 8 changed files with 371 additions and 143 deletions.
52 changes: 51 additions & 1 deletion src/qseek/magnitudes/base.py
Original file line number Diff line number Diff line change
@@ -1,8 +1,14 @@
from __future__ import annotations

from typing import TYPE_CHECKING, Literal
from typing import TYPE_CHECKING, Literal, NamedTuple

import numpy as np
from pydantic import BaseModel, Field
from pyrocko.trace import Trace
from typing_extensions import Self

from qseek.models.detection import Receiver
from qseek.utils import NSL

if TYPE_CHECKING:
from pyrocko.squirrel import Squirrel
Expand Down Expand Up @@ -89,3 +95,47 @@ async def prepare(
NotImplementedError: This method must be implemented by subclasses.
"""
raise NotImplementedError


class StationAmplitudes(NamedTuple):
station_nsl: NSL
peak: float
noise: float
std_noise: float

distance_epi: float
distance_hypo: float

@property
def anr(self) -> float:
"""Amplitude to noise ratio."""
if self.noise == 0.0:
return 0.0
return self.peak / self.noise

@classmethod
def create(
cls,
receiver: Receiver,
traces: list[Trace],
event: EventDetection,
noise_padding: float = 3.0,
) -> Self:
time_arrival = min(receiver.get_arrivals_time_window()).timestamp()
noise_traces = [
tr.chop(tmin=tr.tmin, tmax=time_arrival - noise_padding, inplace=False)
for tr in traces
]

peak_amp = max(np.abs(tr.ydata).max() for tr in traces)
noise_amp = max(np.abs(tr.ydata).max() for tr in noise_traces)
std_noise = max(np.std(tr.ydata) for tr in noise_traces)

return cls(
station_nsl=receiver.nsl,
peak=peak_amp,
noise=noise_amp,
std_noise=std_noise,
distance_hypo=receiver.distance_to(event),
distance_epi=receiver.surface_distance_to(event),
)
52 changes: 4 additions & 48 deletions src/qseek/magnitudes/local_magnitude_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,11 @@
import numpy as np
from pyrocko.trace import PoleZeroResponse

from qseek.utils import ChannelSelector, ChannelSelectors, MeasurementUnit
from qseek.magnitudes.base import StationAmplitudes
from qseek.utils import NSL, ChannelSelector, ChannelSelectors, MeasurementUnit

if TYPE_CHECKING:
from pyrocko.trace import Trace
from typing_extensions import Self

from qseek.models.detection import EventDetection, Receiver

Expand Down Expand Up @@ -88,58 +88,14 @@ def inside(self, value: float) -> bool:


class StationLocalMagnitude(NamedTuple):
station_nsl: tuple[str, str, str]
station_nsl: NSL
magnitude: float
magnitude_error: float
peak_amp: float
distance_epi: float
distance_hypo: float


class StationAmplitudes(NamedTuple):
station_nsl: tuple[str, str, str]
peak: float
noise: float
std_noise: float

distance_epi: float
distance_hypo: float

@property
def anr(self) -> float:
"""Amplitude to noise ratio."""
if self.noise == 0.0:
return 0.0
return self.peak / self.noise

@classmethod
def create(
cls,
receiver: Receiver,
traces: list[Trace],
event: EventDetection,
noise_padding: float = 3.0,
) -> Self:
time_arrival = min(receiver.get_arrivals_time_window()).timestamp()
noise_traces = [
tr.chop(tmin=tr.tmin, tmax=time_arrival - noise_padding, inplace=False)
for tr in traces
]

peak_amp = max(np.abs(tr.ydata).max() for tr in traces)
noise_amp = max(np.abs(tr.ydata).max() for tr in noise_traces)
std_noise = max(np.std(tr.ydata) for tr in noise_traces)

return cls(
station_nsl=receiver.nsl,
peak=peak_amp,
noise=noise_amp,
std_noise=std_noise,
distance_hypo=receiver.distance_to(event),
distance_epi=receiver.surface_distance_to(event),
)


class LocalMagnitudeModel:
epicentral_range: ClassVar[Range | None] = None
hypocentral_range: ClassVar[Range | None] = None
Expand Down Expand Up @@ -203,7 +159,7 @@ def get_station_magnitude(
try:
traces = _COMPONENT_MAP[self.component](traces)
except KeyError:
logger.warning("Could not get channels for %s", receiver.nsl_pretty)
logger.warning("Could not get channels for %s", receiver.nsl.pretty)
return None
if not traces:
return None
Expand Down
Loading

0 comments on commit 670a525

Please sign in to comment.