You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
343 lines
12 KiB
343 lines
12 KiB
from __future__ import print_function
|
|
from builtins import str
|
|
import os
|
|
import copy
|
|
import numpy as num
|
|
from collections import defaultdict
|
|
from pyrocko.gui.snuffling import Snuffling, Param, PhaseMarker, Switch, Choice, \
|
|
EventMarker
|
|
from pyrocko import guts, orthodrome, trace
|
|
from pyrocko.gui.util import to01
|
|
from pyrocko.plot import graph_colors
|
|
|
|
__author__ = 'Catarina Matos (cpfcmatos@gmail.com)'
|
|
|
|
km = 1000.
|
|
wood_anderson_response = trace.PoleZeroResponse(
|
|
zeros=[0., 0.],
|
|
poles=[(-5.49779 - 5.60886j), (-5.49779 + 5.60886j)],
|
|
constant=1.
|
|
)
|
|
|
|
|
|
class LocalMagnitudeSnuffling(Snuffling):
|
|
|
|
'''
|
|
Local Magnitude Estimation
|
|
--------------------------
|
|
|
|
Main Control High- and Lowpass filters are applied to the data before
|
|
simulating the Wood-Anderson receiver.
|
|
|
|
The suggested default values for the geometrical spreading, anelastic
|
|
attenuation and static magnification are those recommended by IASPEI.
|
|
|
|
For correct estimates these values have to be calibrated for the region
|
|
under investigation.
|
|
|
|
For further information:
|
|
http://gfzpublic.gfz-potsdam.de/pubman/item/escidoc:816929:1/component/escidoc:816928/IS_3.3_rev1.pdf
|
|
|
|
Waveform data either need to have unit 'meters' or if instument
|
|
responses have to be removed select *Needs restitution* and use the file
|
|
browser to read responses.
|
|
|
|
References:
|
|
- Bormann, P. and Dewey J., 2012. The new IASPEI standards for determining
|
|
magnitudes from digital data and their relation to classical
|
|
magnitudes. doi:
|
|
- Hutton, K.L. and Boore D.M., 1987. The ML scale in southern California.
|
|
Bull. seism. Soc. Am., 77, 2074-2094
|
|
- Richter C.F., 1935. An instrumental earthquake magnitude scale, Bull.
|
|
seism. Soc. Am., 25, 1-32.
|
|
'''
|
|
def setup(self):
|
|
self._responses = None
|
|
self.add_parameter(Param(
|
|
'geom. spreading', 'const_a', 1.11, 1., 2.))
|
|
self.add_parameter(Param(
|
|
'anelastic attenuation', 'const_b', 0.00189, 0., 1.))
|
|
self.add_parameter(Param(
|
|
'static magnification', 'const_c', -2.09, -5., 5.))
|
|
self.add_parameter(Param(
|
|
'Duration for "fixed" time window',
|
|
'duration_fixed', 200., 1., 500.))
|
|
self.add_parameter(Choice(
|
|
'Time window', 'time_window', 'visible / selected',
|
|
['visible / selected', 'fixed', 'distance dependent']))
|
|
self.add_parameter(Choice(
|
|
'Apply to', 'apply_to', 'active event',
|
|
['active event', 'selected events', 'all events']))
|
|
self.add_parameter(Switch(
|
|
'Show restituted traces', 'show_restituded_traces', False))
|
|
self.add_parameter(Switch(
|
|
'Mark readings', 'show_markers', False))
|
|
self.add_parameter(Switch(
|
|
'Show plot', 'show_plot', False))
|
|
self.add_parameter(Switch(
|
|
'Show message', 'do_show_message', True))
|
|
self.add_parameter(Switch(
|
|
'Needs restitution', 'needs_restitution', False))
|
|
self.add_parameter(Switch(
|
|
'Set event magnitude', 'modify_inplace', False))
|
|
|
|
self.set_name('Local Magnitude')
|
|
self.set_live_update(False)
|
|
|
|
self.vmin = 1500.
|
|
self.vmax = 6000.
|
|
|
|
def read_responses(self, dirname):
|
|
responses = {}
|
|
entries = os.listdir(dirname)
|
|
for entry in entries:
|
|
if entry.endswith('.pf'):
|
|
key = tuple(entry[:-3].split('.'))
|
|
fn = os.path.join(dirname, entry)
|
|
resp = guts.load(filename=fn)
|
|
responses[key] = resp
|
|
|
|
return responses
|
|
|
|
def local_magnitude(self, distance, amplitude):
|
|
|
|
return num.log10(amplitude*1.0e9) + \
|
|
self.const_a*num.log10(distance/km) + \
|
|
self.const_b*distance/km + self.const_c
|
|
|
|
def get_response(self, nslc):
|
|
if self._responses is None:
|
|
self._responses = self.read_responses(self.input_directory())
|
|
|
|
n, s, l, c = nslc
|
|
for k in [(n, s, l, c), (n, s, c), (s, c), (s,)]:
|
|
if k in self._responses:
|
|
return self._responses[k]
|
|
|
|
self.fail(
|
|
'no response information available for trace %s.%s.%s.%s' % nslc)
|
|
|
|
def get_traces(self, event, stations, trace_selector, tpad):
|
|
p = self.get_pile()
|
|
trace_selector_viewer = self.get_viewer_trace_selector('visible')
|
|
if self.time_window == 'distance dependant':
|
|
for station in stations:
|
|
distance = orthodrome.distance_accurate50m(event, station)
|
|
tmin = distance / self.vmax
|
|
tmax = (distance + event.depth) / self.vmin
|
|
|
|
for trs in p.chopper(
|
|
tmin=event.time + tmin,
|
|
tmax=event.time + tmax,
|
|
tpad=tpad,
|
|
trace_selector=lambda tr: (
|
|
trace_selector(tr) and
|
|
trace_selector_viewer(tr) and
|
|
tr.nslc_id[:3] == station.nsl())):
|
|
|
|
for tr in trs:
|
|
yield tr
|
|
|
|
elif self.time_window == 'fixed':
|
|
tmin = 0.
|
|
tmax = self.duration_fixed
|
|
for trs in p.chopper(
|
|
tmin=event.time + tmin,
|
|
tmax=event.time + tmax,
|
|
tpad=tpad,
|
|
trace_selector=lambda tr: (
|
|
trace_selector(tr) and
|
|
trace_selector_viewer(tr))):
|
|
|
|
for tr in trs:
|
|
yield tr
|
|
|
|
else:
|
|
for trs in self.chopper_selected_traces(
|
|
fallback=True, tpad=tpad,
|
|
trace_selector=trace_selector, mode='inview'):
|
|
|
|
for tr in trs:
|
|
yield tr
|
|
|
|
def call(self):
|
|
self.cleanup()
|
|
|
|
if self.apply_to == 'active event':
|
|
event, _ = self.get_active_event_and_stations()
|
|
events = [event]
|
|
|
|
elif self.apply_to == 'selected events':
|
|
events = [m.get_event() for m in self.get_selected_event_markers()]
|
|
|
|
elif self.apply_to == 'all events':
|
|
events = []
|
|
for m in self.get_markers():
|
|
if isinstance(m, EventMarker):
|
|
events.append(m.get_event())
|
|
|
|
events.sort(key=lambda ev: ev.time)
|
|
|
|
if not events:
|
|
self.fail('no event selected')
|
|
|
|
if self.time_window == 'visible / selected' and len(events) != 1:
|
|
self.fail('cannot work with multiple events with "visible / '
|
|
'selected" time window setting')
|
|
|
|
stations_dict = dict((s.nsl(), s) for s in self.get_stations())
|
|
if len(stations_dict) == 0:
|
|
self.fail('no station information found')
|
|
|
|
markers = []
|
|
local_magnitudes = []
|
|
viewer = self.get_viewer()
|
|
fmin = viewer.highpass
|
|
fmax = viewer.lowpass
|
|
if not fmin or not fmax:
|
|
self.fail('Main Controls Highpass and Lowpass filters have to be set') # noqa
|
|
|
|
for event in events:
|
|
mags = defaultdict(list)
|
|
tpad = 2./fmin
|
|
|
|
def trace_selector(tr):
|
|
c = tr.channel.upper()
|
|
return c.endswith('E') or c.endswith('N') or \
|
|
tr.location.endswith('_rest')
|
|
|
|
distances = {}
|
|
rest_traces = []
|
|
|
|
event2 = copy.deepcopy(event)
|
|
|
|
for tr in self.get_traces(
|
|
event, stations_dict.values(), trace_selector, tpad):
|
|
|
|
nslc = tr.nslc_id
|
|
|
|
try:
|
|
tr.highpass(4, fmin, nyquist_exception=True)
|
|
tr.lowpass(4, fmax, nyquist_exception=True)
|
|
except trace.AboveNyquist as e:
|
|
self.fail(str(e))
|
|
|
|
try:
|
|
station = stations_dict[nslc[:3]]
|
|
except KeyError as e:
|
|
print(e)
|
|
continue
|
|
|
|
if self.needs_restitution:
|
|
resp = self.get_response(nslc)
|
|
|
|
try:
|
|
tr_vel = tr.transfer(
|
|
tfade=tpad,
|
|
freqlimits=(
|
|
fmin*0.5, fmin,
|
|
fmax, fmax*2.0),
|
|
transfer_function=resp,
|
|
invert=True)
|
|
except trace.TraceTooShort as e:
|
|
self.fail(str(e))
|
|
continue
|
|
|
|
else:
|
|
try:
|
|
tr_vel = tr.transfer(
|
|
tfade=tpad,
|
|
freqlimits=(
|
|
fmin*0.5, fmin,
|
|
fmax, fmax*2.0),
|
|
transfer_function=wood_anderson_response,
|
|
invert=False)
|
|
except trace.TraceTooShort as e:
|
|
self.fail(str(e))
|
|
continue
|
|
|
|
distance = orthodrome.distance_accurate50m(event, station)
|
|
|
|
tr_vel.set_codes(location=tr_vel.location+'_rest')
|
|
tr_vel.meta = dict(tabu=True)
|
|
t_of_max, amplitude = tr_vel.absmax()
|
|
|
|
if self.show_restituded_traces:
|
|
rest_traces.append(tr_vel)
|
|
m_nslc = tr_vel.nslc_id
|
|
else:
|
|
m_nslc = tr.nslc_id
|
|
|
|
mag = self.local_magnitude(distance, amplitude)
|
|
if self.show_markers:
|
|
markers.append(PhaseMarker(
|
|
[m_nslc],
|
|
t_of_max, t_of_max, 1, phasename='%3.1f' % mag,
|
|
event=event2))
|
|
|
|
mags[nslc[:2]].append(mag)
|
|
distances[nslc[:2]] = distance
|
|
|
|
if not mags:
|
|
continue
|
|
|
|
if rest_traces:
|
|
self.add_traces(rest_traces)
|
|
|
|
for k in mags:
|
|
mags[k] = max(mags[k])
|
|
|
|
local_magnitude = round(num.median(list(mags.values())), 1)
|
|
|
|
if self.show_plot:
|
|
data = []
|
|
for k in mags:
|
|
data.append((distances[k], mags[k]))
|
|
|
|
dists, mags_arr = num.array(data).T
|
|
|
|
dists /= km
|
|
fig = self.figure()
|
|
axes = fig.add_subplot(1, 1, 1)
|
|
axes.plot(dists, mags_arr, 'o', color=to01(graph_colors[0]))
|
|
for x, y, label in zip(dists, mags_arr, mags.keys()):
|
|
axes.text(x, y, '.'.join(label))
|
|
|
|
axes.axhline(local_magnitude, color=to01(graph_colors[0]))
|
|
mag_std = num.std(list(mags.values()))
|
|
|
|
msg = 'local magnitude: %s, std: %s' % \
|
|
(round(local_magnitude, 1),
|
|
round(mag_std, 1))
|
|
axes.text(max(dists), local_magnitude, msg,
|
|
verticalalignment='bottom',
|
|
horizontalalignment='right')
|
|
|
|
axes.axhspan(
|
|
local_magnitude-mag_std,
|
|
local_magnitude+mag_std,
|
|
alpha=0.1)
|
|
|
|
axes.set_xlabel('Distance [km]')
|
|
axes.set_ylabel('Local Magnitude')
|
|
fig.canvas.draw()
|
|
|
|
local_magnitudes.append(local_magnitude)
|
|
|
|
if self.modify_inplace:
|
|
event.magnitude = local_magnitude
|
|
|
|
if markers:
|
|
self.add_markers(markers)
|
|
|
|
if not local_magnitudes:
|
|
self.fail('no results')
|
|
|
|
if self.do_show_message:
|
|
self.show_message('Magnitude', ', '.join(
|
|
'%3.1f' % mag for mag in local_magnitudes))
|
|
|
|
|
|
def __snufflings__():
|
|
return [LocalMagnitudeSnuffling()]
|