Commit 1e0e750e authored by Sebastian Heimann's avatar Sebastian Heimann
Browse files

improved detections plot

parent f5a998ff
......@@ -335,6 +335,7 @@ def scan(
deltat_cf, imax, iframe, fsmooth_min, xpeak, ypeak,
zpeak,
tr_stackmax, tpeaks, apeaks, config.detector_threshold,
wmin, wmax,
pdata, trs, fmin, fmax, idetection,
grid_station_shift_max=shift_max,
movie=show_movie)
......
......@@ -84,6 +84,12 @@ class Carthesian3DGrid(Grid):
return distances
def distance_max(self):
return math.sqrt(
(self.xmax - self.xmin)**2 +
(self.ymax - self.ymin)**2 +
(self.zmax - self.zmin)**2)
def surface_points(self, system='latlon'):
x, y, z = self._get_coords()
xs = num.tile(x, y.size)
......@@ -108,7 +114,10 @@ class Carthesian3DGrid(Grid):
amax=None,
z_slice=None,
cmap=None,
system='latlon'):
system='latlon',
shading='gouraud',
units=1.,
artists=[]):
if system == 'latlon':
assert False, 'not implemented yet'
......@@ -127,7 +136,21 @@ class Carthesian3DGrid(Grid):
else:
a2d = num.max(a3d, axis=0)
axes.pcolormesh(y, x, a2d.T, vmin=amin, vmax=amax, cmap=cmap)
if artists:
if shading == 'gouraud':
artists[0].set_array(a2d.T.ravel())
elif shading == 'flat':
artists[0].set_array(a2d.T[:-1, :-1].ravel())
else:
assert False, 'unknown shading option'
return artists
else:
return [
axes.pcolormesh(
y/units, x/units, a2d.T,
vmin=amin, vmax=amax, cmap=cmap, shading=shading)]
def geometrical_normalization(grid, receivers):
......
import time
import numpy as num
from pyrocko import automap
from lassie import grid as gridmod, geo
km = 1000.
def map_gmt(
receivers, grid, center_lat, center_lon, radius, output_path,
......@@ -62,9 +63,10 @@ def map_geometry(config, output_path):
show_station_labels=False)
def plot_receivers(axes, receivers, system='latlon'):
def plot_receivers(axes, receivers, system='latlon', units=1.0, style=dict(
color='black', ms=10.)):
x, y = geo.points_coords(receivers, system=system)
axes.plot(y, x, '^', color='black', ms=10.)
return axes.plot(y/units, x/units, '^', **style)
def plot_geometry_carthesian(grid, receivers):
......@@ -92,15 +94,20 @@ def plot_geometry_carthesian(grid, receivers):
def plot_detection(
grid, receivers, frames, tmin_frames, deltat_cf, imax, iframe,
fsmooth_min, xpeak, ypeak, zpeak, tr_stackmax, tpeaks, apeaks,
detector_threshold, pdata, trs_raw, fmin, fmax, idetection,
grid_station_shift_max,
detector_threshold, wmin, wmax, pdata, trs_raw, fmin, fmax,
idetection, grid_station_shift_max,
movie=False):
from matplotlib import pyplot as plt
from matplotlib import cm
from matplotlib.animation import FuncAnimation
from pyrocko.cake_plot import mpl_init, labelspace, colors, \
str_to_mpl_color as scolor
nsls = set(tr.nslc_id[:3] for tr in trs_raw)
receivers_on = [r for r in receivers if r.codes in nsls]
receivers_off = [r for r in receivers if r.codes not in nsls]
distances = grid.distances(receivers)
mpl_init()
......@@ -116,8 +123,17 @@ def plot_detection(
axes3 = plt.subplot2grid((2, 3), (0, 1), rowspan=2)
axes4 = plt.subplot2grid((2, 3), (0, 0), rowspan=2)
axes.set_xlabel('X [m]')
axes.set_ylabel('Y [m]')
if grid.distance_max() > km:
dist_units = km
axes.set_xlabel('X [km]')
axes.set_ylabel('Y [km]')
axes4.set_ylabel('Distance [km]')
else:
dist_units = 1.0
axes.set_xlabel('X [m]')
axes.set_ylabel('Y [m]')
axes4.set_ylabel('Distance [m]')
axes.locator_params(nbins=6, tight=True)
axes2.set_xlabel('Time [s]')
......@@ -128,23 +144,27 @@ def plot_detection(
el.set_visible(False)
axes4.set_xlabel('Time [s]')
axes4.set_ylabel('Distance [m]')
tpeak_current = tmin_frames + deltat_cf * iframe
t0 = tpeak_current
tduration = 2.0*grid_station_shift_max + 1./fsmooth_min
for tpeak, apeak in zip(tpeaks, apeaks):
axes2.axvline(
tpeak-t0,
color=scolor('aluminium3'),
lw=2.)
axes2.axvspan(
tr_stackmax.tmin - t0, wmin - t0,
color=scolor('aluminium2'))
axes2.axvspan(
wmax - t0, tr_stackmax.tmax - t0,
color=scolor('aluminium2'))
axes2.axvspan(
tpeak_current-0.5*tduration - t0,
tpeak_current+0.5*tduration - t0,
color=scolor('aluminium2'),
lw=2.)
color=scolor('scarletred2'),
alpha=0.3,
lw=0.)
axes2.set_xlim(tr_stackmax.tmin - t0, tr_stackmax.tmax - t0)
axes2.axhline(
detector_threshold,
......@@ -153,7 +173,7 @@ def plot_detection(
t = tr_stackmax.get_xdata()
amp = tr_stackmax.get_ydata()
axes2.plot(t - t0, amp, color=scolor('scarletred2'), lw=2.)
axes2.plot(t - t0, amp, color=scolor('scarletred2'), lw=1.)
for tpeak, apeak in zip(tpeaks, apeaks):
axes2.plot(
......@@ -197,16 +217,21 @@ def plot_detection(
scalefactor = (dist_max - dist_min) / len(trs) * 0.5
axes3.set_xlim(-0.5*tduration + shift_min, 0.5*tduration + shift_max)
axes3.set_ylim(dist_min - scalefactor, dist_max + scalefactor)
axes3.set_ylim(
(dist_min - scalefactor)/dist_units,
(dist_max + scalefactor)/dist_units)
axes4.set_xlim(-0.5*tduration + shift_min, 0.5*tduration + shift_max)
axes4.set_ylim(dist_min - scalefactor, dist_max + scalefactor)
axes4.set_ylim(
(dist_min - scalefactor)/dist_units,
(dist_max + scalefactor)/dist_units)
axes3.axvline(
0.,
color=scolor('aluminium3'),
lw=2.)
nsl_have = set()
for ishifter, (trs, shift_table, shifter) in enumerate(pdata2):
color = colors[ishifter % len(colors)]
......@@ -215,18 +240,27 @@ def plot_detection(
tpeak_current - 0.5*tduration + shift_min,
tpeak_current + 0.5*tduration + shift_max, inplace=False)
istation = station_index[tr.nslc_id[:3]]
nsl = tr.nslc_id[:3]
istation = station_index[nsl]
shift = shift_table[imax, istation]
print ishifter, shift
axes3.plot(shift, dist, '|', mew=2, mec=color, ms=10)
axes3.plot(shift, dist/dist_units, '|', mew=2, mec=color, ms=10)
t = tr.get_xdata()
amp = tr.get_ydata() * shifter.weight
amp /= amp_max
axes3.plot(
t-t0,
dist + scalefactor*amp + ishifter*scalefactor*0.1,
(dist + scalefactor*amp + ishifter*scalefactor*0.1)/dist_units,
color=color)
axes3.text(t[0]-t0, dist, '.'.join(tr.nslc_id), verticalalignment='top')
if nsl not in nsl_have:
axes3.annotate(
'.'.join(nsl),
xy=(t[0]-t0, dist/dist_units),
xytext=(10., 0.),
textcoords='offset points',
verticalalignment='top')
nsl_have.add(nsl)
for tr in trs_raw:
istation = station_index[tr.nslc_id[:3]]
......@@ -245,8 +279,14 @@ def plot_detection(
t = tr.get_xdata()
amp = tr.get_ydata().astype(num.float)
amp = amp / num.max(num.abs(amp))
axes4.plot(t-t0, dist + scalefactor*amp, color='black', alpha=0.5)
axes4.plot(shift, dist, '|', mew=2, mec=color, ms=10)
axes4.plot(
t-t0, (dist + scalefactor*amp)/dist_units,
color='black', alpha=0.5)
axes4.plot(
shift, dist/dist_units, '|',
mew=2, mec=color, ms=10)
nframes = frames.shape[1]
......@@ -255,51 +295,75 @@ def plot_detection(
amax = frames[imax, iframe]
axes.set_xlim(grid.ymin, grid.ymax)
axes.set_ylim(grid.xmin, grid.xmax)
axes.set_xlim(grid.ymin/dist_units, grid.ymax/dist_units)
axes.set_ylim(grid.xmin/dist_units, grid.xmax/dist_units)
cmap = cm.YlOrBr
system = ('ne', grid.lat, grid.lon)
if movie:
plt.ion()
plt.show()
for iframe in xrange(iframe_min, iframe_max+1):
frame = frames[:, iframe]
axes.cla()
grid.plot(
axes, frame,
amin=0.0,
amax=amax,
cmap=cmap,
system=system)
plot_receivers(axes, receivers, system=system)
axes.plot(
ypeak, xpeak, '*',
ms=20.,
mfc='white',
mec='black')
plt.draw()
time.sleep(0.05)
plt.ioff()
static_artists = []
static_artists.extend(plot_receivers(
axes, receivers_on, system=system, units=dist_units, style=dict(
mfc='black',
ms=5.0)))
else:
frame = num.max(frames[:, iframe_min:iframe_max+1], axis=1)
grid.plot(axes, frame, amin=0.0, amax=amax, cmap=cmap, system=system)
static_artists.extend(plot_receivers(
axes, receivers_off, system=system, units=dist_units, style=dict(
mfc='none',
ms=5.0)))
plot_receivers(axes, receivers, system=system)
axes.plot(
ypeak, xpeak, '*',
ms=20.,
mec='black',
mfc='white')
static_artists.extend(axes.plot(
ypeak/dist_units, xpeak/dist_units, '*',
ms=20.,
mec='black',
mfc='white'))
fig.tight_layout()
frame_artists = []
progress_artists = []
plt.show()
def update(iframe):
if iframe is not None:
frame = frames[:, iframe]
if not progress_artists:
progress_artists[:] = [axes2.axvline(
tmin_frames - t0 + deltat_cf * iframe,
color=scolor('scarletred3'),
alpha=0.5,
lw=2.)]
else:
progress_artists[0].set_xdata(
tmin_frames - t0 + deltat_cf * iframe)
else:
frame = num.max(frames[:, iframe_min:iframe_max+1], axis=1)
frame_artists[:] = grid.plot(
axes, frame,
amin=0.0,
amax=amax,
cmap=cmap,
system=system,
artists=frame_artists,
units=dist_units,
shading='gouraud')
return frame_artists + progress_artists + static_artists
if movie:
ani = FuncAnimation(
fig, update,
frames=list(xrange(iframe_min, iframe_max+1))[::10] + [None],
interval=20.,
repeat=False,
blit=True)
else:
ani = None
update(None)
plt.show()
del ani
__all__ = [
......@@ -309,4 +373,3 @@ __all__ = [
'plot_geometry_carthesian',
'plot_receivers',
]
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment