User contributed plug-ins for Pyrocko's seismic waveform browser Snuffler.
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.
 
 
 
 

260 lines
8.8 KiB

from pyrocko.gui.snuffling import Snuffling, Param, Choice, Switch
from pyrocko.gui.util import EventMarker
from pyrocko import orthodrome
from pyrocko import moment_tensor
from pyrocko.orthodrome import distance_accurate50m as distance
from pyrocko import util, model
import matplotlib.dates as mdates
from matplotlib import cm
import numpy as num
from datetime import datetime
km = 1000.
cmap = cm.RdYlBu
cmaps = {'Red-Yellow-Blue': 'RdYlBu',
'Viridis': 'viridis',
'Magma': 'magma'}
save_cmaps = {}
for k, v in cmaps.items():
try:
x = getattr(cm, v)
save_cmaps[k] = x
except AttributeError:
continue
class TimeLine(Snuffling):
'''
<html>
<body>
<h1>Temporal Evolution of Seismicity</h1>
The considered region can be limited by defining one central coordinate
and a maximum epicentral distance.
</body>
</html>
'''
def setup(self):
'''Customization of the snuffling.'''
self.set_name('Time Line')
self.add_parameter(
Param('Latitude:', 'lat', 90, -90., 90., high_is_none=True))
self.add_parameter(
Param('Longitude:', 'lon', 180., -180, 180., high_is_none=True))
self.add_parameter(
Param('Maximum Distance [km]:', 'maxd', 20000., 0., 20000.,
high_is_none=True))
self.add_parameter(
Choice('Color by', 'color_by', 'time',
['time', 'longitude', 'latitude', 'magnitude', 'depth', 'kind']))
self.add_parameter(Choice('Colormap', 'cmap_selector',
'Red-Yellow-Blue', list(save_cmaps.keys())))
self.add_parameter(Choice('Coordinate system', 'coord_system',
'Lat/Lon', ['Lat/Lon', 'cartesian']))
self.add_parameter(Switch('Show stations', 'show_stations', False))
self.add_trigger('Save Figure', self.save_as)
self.set_live_update(False)
self.fig = None
self.cli_mode = False
def call(self):
'''Main work routine of the snuffling.'''
self.cleanup()
cmap = save_cmaps[self.cmap_selector]
viewer = self.get_viewer()
tmin, tmax = self.get_selected_time_range(fallback=True)
event_markers = filter(lambda x: x.tmin >= tmin and x.tmax <= tmax,
viewer.markers)
event_markers = filter(
lambda x: isinstance(x, EventMarker), event_markers)
if self.maxd:
event_markers = filter(
lambda x: distance(self, x._event) <= self.maxd*km,
event_markers)
if event_markers == []:
self.fail('No events in selected area found')
event_markers = list(event_markers)
stations = self.get_stations() if self.show_stations else None
self.make_time_line(event_markers, stations, cmap=cmap)
def make_time_line(self, markers, stations=None, cmap=cmap):
events = [m.get_event() for m in markers]
kinds = num.array([m.kind for m in markers])
if self.cli_mode:
self.fig = plt.figure()
else:
fframe = self.figure_frame()
self.fig = fframe.gcf()
ax = self.fig.add_subplot(311)
ax_cum = ax.twinx()
ax1 = self.fig.add_subplot(323)
ax2 = self.fig.add_subplot(325, sharex=ax1)
ax3 = self.fig.add_subplot(324, sharey=ax1)
num_events = len(events)
data = num.zeros((num_events, 6))
column_to_index = dict(zip(['magnitude', 'latitude', 'longitude', 'depth', 'time', 'kind'],
range(6)))
c2i = column_to_index
for i, e in enumerate(events):
if e.magnitude:
mag = e.magnitude
else:
mag = 0.
data[i, :] = mag, e.lat, e.lon, e.depth, e.time, kinds[i]
s_coords = num.array([])
s_labels = []
if stations is not None:
s_coords = num.array([(s.lon, s.lat, s.elevation-s.depth) for s in stations])
s_labels = ['.'.join(s.nsl()) for s in stations]
isorted = num.argsort(data[:, c2i['time']])
data = data[isorted]
def _D(key):
return data[:, c2i[key]]
tmin = _D('time').min()
tmax = _D('time').max()
lon_max = _D('longitude').max()
lon_min = _D('longitude').min()
lat_max = _D('latitude').max()
lat_min = _D('latitude').min()
depths_min = _D('depth').min()
depths_max = _D('depth').max()
mags_min = _D('magnitude').min()
mags_max = _D('magnitude').max()
moments = moment_tensor.magnitude_to_moment(_D('magnitude'))
dates = list(map(datetime.fromtimestamp, _D('time')))
fds = mdates.date2num(dates)
tday = 3600*24
tweek = tday*7
if tmax-tmin < 1*tday:
hfmt = mdates.DateFormatter('%Y-%m-%d %H:%M:%S')
elif tmax-tmin < tweek*52:
hfmt = mdates.DateFormatter('%Y-%m-%d')
else:
hfmt = mdates.DateFormatter('%Y/%m')
color_values = _D(self.color_by)
color_args = dict(c=color_values, vmin=color_values.min(),
vmax=color_values.max(), cmap=cmap)
ax.scatter(fds, _D('magnitude'), s=20, **color_args)
ax.xaxis.set_major_formatter(hfmt)
ax.spines['top'].set_color('none')
ax.spines['right'].set_color('none')
ax.set_ylim((mags_min, mags_max*1.10))
ax.set_xlim(map(datetime.fromtimestamp, (tmin, tmax)))
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
ax.set_ylabel('Magnitude')
init_pos = ax.get_position()
ax_cum.plot(fds, num.cumsum(moments), 'grey')
ax_cum.xaxis.set_major_formatter(hfmt)
ax_cum.spines['top'].set_color('none')
ax_cum.spines['right'].set_color('grey')
ax_cum.set_ylabel('Cumulative seismic moment')
lats_min = num.array([lat_min for x in range(num_events)])
lons_min = num.array([lon_min for x in range(num_events)])
if self.coord_system == 'cartesian':
lats, lons = orthodrome.latlon_to_ne_numpy(
lats_min, lons_min, _D('latitude'), _D('longitude'))
_x = num.empty((len(s_coords), 3))
for i, (slon, slat, sele) in enumerate(s_coords):
n, e = orthodrome.latlon_to_ne(lat_min, lon_min, slat, slon)
_x[i, :] = (e, n, sele)
s_coords = _x
else:
lats = _D('latitude')
lons = _D('longitude')
s_coords = s_coords.T
ax1.scatter(lons, lats, s=20, **color_args)
ax1.set_aspect('equal')
ax1.grid(True, which='both')
ax1.set_ylabel('Northing [m]')
ax1.get_yaxis().tick_left()
if len(s_coords):
ax1.scatter(s_coords[0], s_coords[1], marker='v', s=40, color='black')
for c, sl in zip(s_coords.T, s_labels):
ax1.text(c[0], c[1], sl, color='black')
# bottom left plot
ax2.scatter(lons, _D('depth'), s=20, **color_args)
ax2.grid(True)
ax2.set_xlabel('Easting [m]')
ax2.set_ylabel('Depth [m]')
ax2.get_yaxis().tick_left()
ax2.get_xaxis().tick_bottom()
ax2.invert_yaxis()
ax2.text(1.1, 0, 'Origin at:\nlat=%1.3f, lon=%1.3f' %
(lat_min, lon_min), transform=ax2.transAxes)
# top right plot
ax3.scatter(_D('depth'), lats, s=20, **color_args)
ax3.set_xlim((depths_min, depths_max))
ax3.grid(True)
ax3.set_xlabel('Depth [m]')
ax3.get_xaxis().tick_bottom()
ax3.get_yaxis().tick_right()
self.fig.subplots_adjust(
bottom=0.05, right=0.95, left=0.075, top=0.95, wspace=0.02, hspace=0.02)
init_pos.y0 += 0.05
ax.set_position(init_pos)
ax_cum.set_position(init_pos)
if self.cli_mode:
plt.show()
else:
self.fig.canvas.draw()
def save_as(self):
if self.fig:
fn = self.output_filename()
self.fig.savefig(fn,
pad_inches=0.05,
bbox_inches='tight')
def configure_cli_parser(self, parser):
parser.add_option(
'--events',
dest='events_filename',
default=None,
metavar='FILENAME',
help='Read events from FILENAME')
def __snufflings__():
'''Returns a list of snufflings to be exported by this module.'''
return [TimeLine()]
if __name__ == '__main__':
import matplotlib.pyplot as plt
util.setup_logging('time_line.py', 'info')
s = TimeLine()
options, args, parser = s.setup_cli()
s.cli_mode = True
if options.events_filename:
s.make_time_line(
list(model.Event.load_catalog(options.events_filename)))