Skip to content

Commit

Permalink
Fixes/local magnitude (#8)
Browse files Browse the repository at this point in the history
* fixes local magnitudes

* adding docs
  • Loading branch information
miili authored Oct 31, 2023
1 parent a397482 commit 60916cd
Showing 1 changed file with 35 additions and 7 deletions.
42 changes: 35 additions & 7 deletions lassie/features/local_magnitude.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,15 @@
constant=2080.0,
)

WOOD_ANDERSON = trace.PoleZeroResponse(
poles=[
-6.283 - 4.7124j,
-6.283 + 4.7124j,
],
zeros=[0.0 + 0.0j, 0.0 + 0.0j],
constant=2080.0,
)

KM = 1e3
MM = 1e3

Expand Down Expand Up @@ -81,6 +90,16 @@ def calculate_magnitude(
receiver: Receiver,
traces: list[Trace],
) -> StationMagnitude | None:
"""Calculate the local magnitude for a given event and receiver.
Args:
event (EventDetection): The event to calculate the magnitude for.
receiver (Receiver): The seismic station to calculate the magnitude for.
traces (list[Trace]): The traces to calculate the magnitude for.
Returns:
StationMagnitude | None: The calculated magnitude or None if the magnitude.
"""
dist_hypo = event.distance_to(receiver)
dist_epi = event.surface_distance_to(receiver)
if not self._is_distance_valid(dist_hypo, dist_epi):
Expand Down Expand Up @@ -128,7 +147,7 @@ def calculate_magnitude(
logger.exception(exc)
return None

amp_max *= 1000000 # To nm
amp_max *= 1e6 # To nm
local_magnitude = (
np.log10(amp_max) + 1.11 * np.log10(dist_hypo) + 0.00189 * dist_hypo - 2.09
)
Expand Down Expand Up @@ -199,8 +218,9 @@ def get_amp_0(self, dist_hypo_km: float, dist_epi_km: float) -> float:
class LocalMagnitudeExtractor(FeatureExtractor):
feature: Literal["LocalMagnitude"] = "LocalMagnitude"

seconds_before: float = 5.0
seconds_after: float = 15.0
seconds_before: float = 3.0
seconds_after: float = 10.0
window_slack: float = 10.0
estimator: Union[
IASPEISouthernCalifornia,
SouthernCalifornia,
Expand All @@ -220,19 +240,27 @@ async def add_features(self, squirrel: Squirrel, event: EventDetection) -> None:
try:
traces = receiver.get_waveforms_restituted(
squirrel,
seconds_before=self.seconds_before,
seconds_after=self.seconds_after,
seconds_before=self.seconds_before - self.window_slack,
seconds_after=self.seconds_after + self.window_slack,
quantity="velocity",
)
except NotAvailable:
logger.error("cannot get responses for %s", receiver.pretty_nsl)
continue

restituded_traces = [
restituted_traces = [
tr.transfer(transfer_function=WOOD_ANDERSON) for tr in traces
]

for tr in restituted_traces:
tr.chop(
float(tr.tmin) + self.window_slack,
float(tr.tmax) - self.window_slack,
inplace=True,
)

magnitude = self.estimator.calculate_magnitude(
event, receiver, restituded_traces
event, receiver, restituted_traces
)
if magnitude is None:
continue
Expand Down

0 comments on commit 60916cd

Please sign in to comment.