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