1import logging
2import numpy as num
3import matplotlib
4from matplotlib import patches
5from pyrocko import util, trace, plot
6from pyrocko.gui.snuffling import Snuffling, Param, Marker, Switch, \
7 EventMarker
9logger = logging.getLogger('pyrocko.gui.snufflings.polarization')
11d2r = num.pi / 180.
12r2d = 1.0 / d2r
15def get_tr(traces, f):
16 for tr in traces:
17 if f(tr):
18 return tr
20 raise ValueError('Element not found.')
23def darken(color, f=0.5):
24 return tuple(c*f for c in color[:3]) + color[3:]
27class PCAError(Exception):
28 pass
31class LayoutError(Exception):
32 pass
35class Polarization(Snuffling):
37 u'''
38Polarization
39============
41Investigate patterns of ground motion during the passage of seismic waves.
43This Snuffling can be used to analyze and visualize the polarization of seismic
44waves from 3-component seismic recordings or to check component orientations of
45a seismic sensors when used on signals with known directional properties. The
46spatial pattern of ground motion is shown in horizontal and vertical
47projections. Principal component analysis and rotation to radial/transverse
48components are available as tools. Time window and filter settings can be
49adjusted interactively.
51Usage
52-----
54Select one or more normal/phase markers as anchor points for the extraction and
55press *Run*. Multiple stations can be selected for direct comparison. Channels
56matching the pattern-triplet given in the *Channels* setting are selected for
57extraction of a time window around the anchor point. It is assumed that these
58three channels correspond to sensor components in the order east, north,
59vertical (upward), even if they are named differently.
61The time window can be adjusted with the *Length* and *Offset* parameters.
62Extracted waveforms are filtered according the *Highpass* and *Lowpass*
63parameters (Butterworth 4th order) and demeaned.
65To rotate the horizontal components around a vertical axis, use the *\u0394
66Azimuth* setting. When station coordinates and an "active" event are available,
67the horizontal components can be rotated to radial (away from event) and
68transverse (leftwards) orientations using the computed event back-azimuth
69(dashed gray line).
71If *Show 2D Eigensystems* is selected, a principal component analysis is
72performed in each of the shown projects. Red lines are shown to depict the
73eigenvectors and the eigenvalues are visualized using ellipse symbols. If *Show
743D Eigensystems* is selected, a principal component analysis is performed using
75all three (non-rotated) components and the resulting eigenvectors are depicted
76with purple lines.
78By default the scaling is automatically adjusted to the maximum value (the
79maximum vector length of the three-component signal within the selected time
80window). The currently used scaling factor can be frozen by checking *Fix
81Scale*.
82'''
84 def setup(self):
85 self.set_name('Polarization')
87 self.add_parameter(
88 Param('Length [s]', 't_length', 1., 0.001, 1000.))
90 self.add_parameter(
91 Param('Offset (relative)', 'ft_offset', 0., -5., 5.))
93 self.add_parameter(
94 Param(u'\u0394 Azimuth', 'azimuth', 0., -180., 180.))
96 self.add_parameter(
97 Param(
98 'Highpass [Hz]', 'highpass', None, 0.001, 1000.,
99 low_is_none=True))
101 self.add_parameter(
102 Param(
103 'Lowpass [Hz]', 'lowpass', None, 0.001, 1000.,
104 high_is_none=True))
106 self.add_parameter(
107 Param('Dot Position', 'dot_position', 0., 0., 1.))
109 self.add_parameter(
110 Switch(
111 'Rotate to RT',
112 'rotate_to_rt',
113 False))
115 self.add_parameter(
116 Switch(
117 'Show 2D Eigensystems',
118 'show_eigensystem_2d',
119 False))
121 self.add_parameter(
122 Switch(
123 'Show 3D Eigensystem',
124 'show_eigensystem_3d',
125 False))
127 self.add_parameter(
128 Switch(
129 'Fix Scale',
130 'fix_scale',
131 False))
133 def new_figure():
134 self.setup_figure_frame()
135 self.call()
137 self.add_trigger('New Figure', new_figure)
139 self.fframe = None
140 self.iframe = 0
141 self.figure_key = None
142 self.nsl_to_amax = {}
143 self.last_rotate_to_rt = False
145 self.font_size = float(matplotlib.rcParams['font.size'])
147 self.colors = {
148 'fg': matplotlib.rcParams['axes.edgecolor'],
149 'lines': 'black',
150 'rot': plot.mpl_color('skyblue1'),
151 'pca_2d': plot.mpl_color('scarletred1'),
152 'pca_3d': plot.mpl_color('plum2'),
153 'backazimuth': plot.mpl_color('aluminium3'),
154 'time_window': plot.mpl_color('aluminium2')}
156 def panel_visibility_changed(self, visible):
157 viewer = self.get_viewer()
158 if visible:
159 viewer.pile_has_changed_signal.connect(self.adjust_controls)
160 self.adjust_controls()
162 else:
163 viewer.pile_has_changed_signal.disconnect(self.adjust_controls)
165 def adjust_controls(self):
166 viewer = self.get_viewer()
167 dtmin, dtmax = viewer.content_deltat_range()
168 maxfreq = 0.5/dtmin
169 minfreq = (0.5/dtmax)*0.001
170 self.set_parameter_range('lowpass', minfreq, maxfreq)
171 self.set_parameter_range('highpass', minfreq, maxfreq)
173 def get_selected(self):
174 markers = [
175 marker for marker in self.get_selected_markers()
176 if not isinstance(marker, EventMarker)]
178 if not markers:
179 self.fail(
180 'No selected markers.\n\nCreate and select markers at points '
181 'of interest. Normal and phase markers are accepted.')
183 d = {}
184 for marker in markers:
185 tspan = self.transform_time_span(marker.tmin, marker.tmax)
186 for nslc in marker.nslc_ids:
187 k = nslc[:3] + (nslc[3][:-1], marker.tmin)
188 d[k] = tspan
190 return d
192 def transform_time_span(self, tmin, tmax):
193 tmin = tmin + self.t_length * self.ft_offset
194 tmax = tmin + self.t_length
195 return (tmin, tmax)
197 def make_selection_markers(self, selected, groups):
198 selection_markers = []
199 for k in sorted(selected):
200 tmin, tmax = selected[k]
201 patterns = [tr.nslc_id for tr in groups[k]]
202 marker = Marker(
203 nslc_ids=patterns, tmin=tmin, tmax=tmax, kind=2)
204 selection_markers.append(marker)
206 return selection_markers
208 def sorted_by_component_scheme(self, group):
209 if len(group) not in (2, 3):
210 self.fail('Need 2 or 3 component recordings. Got %i.' % len(group))
212 comps = []
213 channels = []
214 for tr in group:
215 comps.append(tr.channel[-1])
216 channels.append(tr.channel)
218 # special case for channel naming convention of Cube dataloggers
219 for scheme in [['p2', 'p1', 'p0'], ['p2', 'p1']]:
221 if sorted(channels) == sorted(scheme):
222 return [
223 get_tr(group, lambda tr: tr.channel == ch)
224 for ch in scheme]
226 for scheme in [
227 ['E', 'N', 'Z'], ['1', '2', 'Z'], ['1', '2', '3'],
228 ['0', '1', '2'], ['R', 'T', 'Z'],
229 ['E', 'N'], ['1', '2'], ['0', '1']]:
231 if sorted(comps) == sorted(scheme):
232 return [
233 get_tr(group, lambda tr: tr.channel[-1] == comp)
234 for comp in scheme]
236 self.fail('Cannot handle component group: %s' % ', '.join(channels))
238 def get_traces(self, selected, nsl_to_bazi, tpad):
239 if self.highpass is not None:
240 tpad_filter = 3.0 / self.highpass
241 elif self.lowpass is not None:
242 tpad_filter = 3.0 / self.lowpass
243 else:
244 tpad_filter = 0.0
246 # prevent getting insanely long cutouts if e.g. if highpass is still
247 # very low, e.g. while moving the slider
248 tpad_filter = min(tpad_filter, 5.0 * self.t_length)
250 d = {}
251 for k in sorted(selected):
252 nsl = k[0:3]
253 tmin, tmax = selected[k]
255 bazimuth = self.azimuth
257 if self.rotate_to_rt:
258 if nsl not in nsl_to_bazi:
259 self.fail(
260 'Cannot rotate to RT.\n\nStation coordinates must be '
261 'available and an event must be marked as the '
262 '"active" event (select event and press "e").')
264 bazimuth += nsl_to_bazi[nsl] + 90.
266 pattern = '.'.join(nsl + (k[3] + '?',))
268 group = []
269 for batch in self.get_pile().chopper(
270 tmin=tmin, tmax=tmax, tpad=tpad + tpad_filter,
271 want_incomplete=False,
272 style='batch',
273 trace_selector=lambda tr: util.match_nslc(
274 pattern, tr.nslc_id)):
276 for tr in batch.traces:
277 tr = tr.copy()
278 if self.lowpass is not None \
279 and self.lowpass < 0.5/tr.deltat:
280 tr.lowpass(4, self.lowpass)
282 if self.highpass is not None \
283 and self.highpass < 0.5/tr.deltat:
285 tr.highpass(4, self.highpass)
287 try:
288 tr.chop(tmin - tpad, tmax + tpad)
289 tr_chop = tr.chop(tmin, tmax, inplace=False)
290 except trace.NoData:
291 self.fail('No data. Time length too short?')
293 y = tr.get_ydata()
294 tr.set_ydata(y - num.mean(tr_chop.get_ydata()))
296 group.append(tr)
298 group = self.sorted_by_component_scheme(group)
300 tr_e = group[0]
301 tr_n = group[1]
303 if tr_e is not None and tr_n is not None:
304 cha_e = tr_e.channel
305 cha_n = tr_n.channel
306 if bazimuth != 0.0:
307 trs_rot = trace.rotate(
308 [tr_n, tr_e], bazimuth,
309 in_channels=[cha_n, cha_e],
310 out_channels=[ch+"'" for ch in [cha_n, cha_e]])
311 else:
312 trs_rot = [tr.copy() for tr in [tr_n, tr_e]]
313 for tr in trs_rot:
314 tr.set_codes(channel=tr.channel + "'")
316 if len(trs_rot) == 2:
317 group.extend(
318 get_tr(trs_rot, lambda tr: tr.channel == ch+"'")
319 for ch in [cha_e, cha_n])
320 else:
321 self.fail(
322 'Cannot rotate traces. '
323 'Are the samples properly aligned?')
325 d[k] = group
327 return d
329 def setup_figure_frame(self):
330 self.iframe += 1
331 self.fframe = self.figure_frame(
332 'Polarization (%i)' % self.iframe)
334 self.fframe.gcf().my_disconnect = None
336 def get_figure(self):
337 if not self.fframe or self.fframe.closed:
338 self.setup_figure_frame()
340 return self.fframe.gcf()
342 def setup_figure(self, fig, nstations):
343 fig.clf()
344 new_axes = []
346 def iwrap(iy, ix):
347 return (ix + iy * 4) + 1
349 for istation in range(nstations):
350 axes_01 = fig.add_subplot(
351 nstations, 4, iwrap(istation, 0), aspect=1.0)
352 axes_02 = fig.add_subplot(
353 nstations, 4, iwrap(istation, 1), aspect=1.0)
354 axes_12 = fig.add_subplot(
355 nstations, 4, iwrap(istation, 2), aspect=1.0)
356 axes_tr = fig.add_subplot(
357 nstations, 4, iwrap(istation, 3))
359 for axes in (axes_01, axes_02, axes_12, axes_tr):
360 axes.my_stuff = []
362 axes.my_line, = axes.plot(
363 [], [], color=self.colors['lines'], lw=1.0)
364 axes.my_dot, = axes.plot(
365 [], [], 'o', ms=4, color=self.colors['lines'])
367 for axes in (axes_01, axes_02, axes_12):
368 axes.get_xaxis().set_tick_params(
369 labelbottom=False, bottom=False)
370 axes.get_yaxis().set_tick_params(
371 labelleft=False, left=False)
373 axes_tr.get_yaxis().set_tick_params(
374 left=False, labelleft=True, length=2.0)
376 # if istation != nstations - 1:
377 # axes_tr.get_xaxis().set_tick_params(
378 # bottom=False, labelbottom=False)
380 lines = []
381 dots = []
382 for _ in range(3):
383 lines.append(
384 axes_tr.plot(
385 [], [], color=self.colors['lines'], lw=1.0)[0])
386 dots.append(
387 axes_tr.plot(
388 [], [], 'o', ms=4, color=self.colors['lines'])[0])
390 axes_tr.my_lines = lines
391 axes_tr.my_dots = dots
392 axes_tr.my_stuff = []
394 new_axes.append(
395 (axes_01, axes_02, axes_12, axes_tr))
397 def resize_handler(*args):
398 self.layout(fig, new_axes)
400 if fig.my_disconnect:
401 fig.my_disconnect()
403 cid_resize = fig.canvas.mpl_connect('resize_event', resize_handler)
404 cid_dpi = fig.callbacks.connect('dpi_changed', resize_handler)
406 def disconnect():
407 fig.canvas.mpl_disconnect(cid_resize)
408 fig.callbacks.disconnect(cid_dpi)
410 fig.my_disconnect = disconnect
412 self.axes = new_axes
414 def get_trace(self, traces, pattern):
415 trs = [tr for tr in traces if util.match_nslc([pattern], tr.nslc_id)]
416 if len(trs) > 1:
417 self.fail('Multiple traces matching pattern %s' % pattern)
418 elif len(trs) == 0:
419 return None
420 else:
421 return trs[0]
423 def get_vector_abs_max(self, traces):
424 tr_abs = None
425 for tr in traces:
426 if tr is not None:
427 tr = tr.copy()
428 tr.ydata **= 2
429 if tr_abs is None:
430 tr_abs = tr
431 else:
432 tr_abs.add(tr)
434 tr_abs.set_ydata(num.sqrt(tr_abs.ydata))
435 return num.max(tr_abs.ydata)
437 def set_labels(
438 self, istation, nstations, tref,
439 axes_01, axes_02, axes_12, axes_tr, chas, chas_rot):
441 if self.azimuth != 0 or self.rotate_to_rt:
442 rcolor = self.colors['rot']
443 else:
444 rcolor = self.colors['fg']
445 chas_rot = chas
447 axes_01.set_ylabel(chas[1])
448 axes_02.set_ylabel(chas[2])
449 axes_12.set_ylabel(chas[2])
451 axes_01.set_xlabel(chas[0])
452 axes_02.set_xlabel(chas_rot[0])
453 axes_12.set_xlabel(chas_rot[1])
454 axes_02.get_xaxis().label.set_color(rcolor)
455 axes_12.get_xaxis().label.set_color(rcolor)
456 axes_tr.set_xlabel(
457 'Time [s] relative to %s' % util.time_to_str(tref))
459 axes_tr.set_yticks([0., 1., 2])
460 axes_tr.set_yticklabels([chas_rot[0], chas_rot[1], chas[2]])
461 for tlab in axes_tr.get_yticklabels()[:2]:
462 tlab.set_color(rcolor)
464 def pca(self, trs):
466 if any(tr is None for tr in trs):
467 raise PCAError('Missing component')
469 nss = [tr.data_len() for tr in trs]
470 if not all(ns == nss[0] for ns in nss):
471 raise PCAError('Traces have different lengths.')
473 ns = nss[0]
475 if ns < 3:
476 raise PCAError('Traces too short.')
478 data = num.zeros((len(trs), ns))
479 for itr, tr in enumerate(trs):
480 data[itr, :] = tr.ydata
482 cov = num.cov(data)
484 evals, evecs = num.linalg.eigh(cov)
486 pc = evecs[:, -1]
488 eh = num.sqrt(pc[1]**2 + pc[0]**2)
490 if len(trs) > 2:
491 incidence = r2d * num.arctan2(eh, abs(pc[2]))
492 else:
493 incidence = 90.
495 azimuth = r2d*num.arctan2(pc[1], pc[0])
496 azimuth = (-azimuth) % 180. - 90.
498 return cov, evals, evecs, azimuth, incidence
500 def draw_cov_ellipse(self, evals, evecs, color, alpha=1.0):
501 evals = num.sqrt(evals)
502 ell = patches.Ellipse(
503 xy=(0.0, 0.0),
504 width=evals[0] * 2.,
505 height=evals[1] * 2.,
506 angle=r2d*num.arctan2(evecs[1, 0], evecs[0, 0]),
507 zorder=-10,
508 fc=color,
509 ec=darken(color),
510 alpha=alpha)
512 return ell
514 def draw(self, groups, nsl_to_tspan, tpad, nsl_to_bazi):
516 def count(d, k):
517 if k not in d:
518 d[k] = 0
520 d[k] += 1
522 nsli_n = {}
523 for k in sorted(groups):
524 count(nsli_n, k[:4])
526 nsli_i = {}
527 for igroup, k in enumerate(sorted(groups)):
529 tmin, tmax = nsl_to_tspan[k]
530 nsl = k[:3]
531 nsli = k[:4]
532 count(nsli_i, nsli)
534 bazimuth = self.azimuth
536 if self.rotate_to_rt:
537 if nsl not in nsl_to_bazi:
538 self.fail(
539 'Cannot rotate to RT.\n\nActive event must '
540 'available (select event and press "e"). Station '
541 'coordinates must be available.')
543 bazimuth += nsl_to_bazi[nsl] + 90.
545 for axes in self.axes[igroup]:
546 while axes.my_stuff:
547 stuff = axes.my_stuff.pop()
548 stuff.remove()
550 axes_01, axes_02, axes_12, axes_tr = self.axes[igroup]
552 axes_01.set_title(
553 '.'.join(nsli)
554 + ('' if nsli_n[nsli] == 1 else ' (%i)' % nsli_i[nsli]))
556 trs_all = groups[k]
558 if len(trs_all) == 5:
559 trs_orig = trs_all[:3]
560 trs_rot = trs_all[3:] + trs_all[2:3]
561 elif len(trs_all) == 4:
562 trs_orig = trs_all[:2] + [None]
563 trs_rot = trs_all[2:] + [None]
564 else:
565 raise self.fail(
566 'Unexpectedly got other that 4 or 6 traces: %i'
567 % len(trs_all))
569 trs_orig_chopped = [
570 (tr.chop(tmin, tmax, inplace=False) if tr else None)
571 for tr in trs_orig]
573 trs_rot_chopped = [
574 (tr.chop(tmin, tmax, inplace=False) if tr else None)
575 for tr in trs_rot]
577 chas = [tr.channel[-1:] for tr in trs_orig_chopped]
578 chas_rot = [tr.channel[-2:] for tr in trs_rot_chopped]
580 if self.fix_scale and nsl in self.nsl_to_amax:
581 amax = self.nsl_to_amax[nsl]
582 else:
583 amax = self.get_vector_abs_max(trs_orig_chopped)
584 self.nsl_to_amax[nsl] = amax
586 if nsl in nsl_to_bazi:
587 event_bazi = nsl_to_bazi[nsl]
588 l1, = axes_01.plot(
589 [0., amax*num.cos((90. - event_bazi)*d2r)],
590 [0., amax*num.sin((90. - event_bazi)*d2r)],
591 '--',
592 lw=3,
593 zorder=-20,
594 color=self.colors['backazimuth'])
596 a1 = axes_01.annotate(
597 '%.0f\u00b0' % event_bazi,
598 (0., 1.),
599 (self.font_size, -self.font_size),
600 xycoords='axes fraction',
601 textcoords='offset points',
602 color=self.colors['backazimuth'],
603 ha='left',
604 va='top')
606 axes_01.my_stuff.extend([l1, a1])
608 for ix, iy, axes, trs in [
609 (0, 1, axes_01, trs_orig_chopped),
610 (0, 2, axes_02, trs_rot_chopped),
611 (1, 2, axes_12, trs_rot_chopped)]:
613 axes.set_xlim(-amax*1.05, amax*1.05)
614 axes.set_ylim(-amax*1.05, amax*1.05)
616 if not (trs[ix] and trs[iy]):
617 continue
619 x = trs[ix].get_ydata()
620 y = trs[iy].get_ydata()
622 axes.my_line.set_data(x, y)
623 ipos = int(round(self.dot_position * (x.size-1)))
624 axes.my_dot.set_data(x[ipos], y[ipos])
626 tref = tmin
627 for itr, (tr, tr_chopped) in enumerate(zip(
628 trs_rot, trs_rot_chopped)):
630 if tr is None or tr_chopped is None:
631 axes_tr.my_lines[itr].set_data([], [])
632 axes_tr.my_dots[itr].set_data([], [])
634 else:
635 y = tr.get_ydata() / (2.*amax) + itr
636 t = tr.get_xdata()
637 t = t - tref
639 ipos = int(round(
640 self.dot_position * (tr_chopped.data_len()-1)))
642 yp = tr_chopped.ydata[ipos] / (2.*amax) + itr
643 tp = tr_chopped.tmin - tref + tr_chopped.deltat*ipos
645 axes_tr.my_lines[itr].set_data(t, y)
646 axes_tr.my_dots[itr].set_data(tp, yp)
648 if self.azimuth != 0.0 or self.rotate_to_rt:
649 color = self.colors['rot']
651 xn = num.sin(bazimuth*d2r)
652 yn = num.cos(bazimuth*d2r)
653 xe = num.sin(bazimuth*d2r + 0.5*num.pi)
654 ye = num.cos(bazimuth*d2r + 0.5*num.pi)
656 l1, = axes_01.plot(
657 [0., amax*xn],
658 [0., amax*yn],
659 color=color)
661 font_size = self.font_size
663 a1 = axes_01.annotate(
664 chas_rot[1],
665 xy=(amax*xn, amax*yn),
666 xycoords='data',
667 xytext=(-font_size*(xe+.5*xn), -font_size*(ye+.5*yn)),
668 textcoords='offset points',
669 va='center',
670 ha='center',
671 color=color)
673 l2, = axes_01.plot(
674 [0., amax*xe],
675 [0., amax*ye],
676 color=color)
678 a2 = axes_01.annotate(
679 chas_rot[0],
680 xy=(amax*xe, amax*ye),
681 xycoords='data',
682 xytext=(-font_size*(xn+.5*xe), -font_size*(yn+.5*ye)),
683 textcoords='offset points',
684 va='center',
685 ha='center',
686 color=color)
688 aa = axes_01.annotate(
689 '%.0f\u00b0' % bazimuth,
690 (1., 1.),
691 (-self.font_size, -self.font_size),
692 xycoords='axes fraction',
693 textcoords='offset points',
694 color=color,
695 ha='right',
696 va='top')
698 axes_01.my_stuff.extend([l1, a1, l2, a2, aa])
700 axes_tr.my_stuff.append(axes_tr.axvspan(
701 tmin - tref, tmax - tref, color=self.colors['time_window']))
703 axes_tr.set_ylim(-1, 3)
704 axes_tr.set_xlim(tmin - tref - tpad, tmax - tref + tpad)
706 self.set_labels(
707 igroup, len(groups), tref, *self.axes[igroup], chas, chas_rot)
709 if self.show_eigensystem_2d:
711 for (ix, iy, axes, trs) in [
712 (0, 1, axes_01, trs_orig_chopped),
713 (0, 2, axes_02, trs_rot_chopped),
714 (1, 2, axes_12, trs_rot_chopped)]:
716 try:
717 cov, evals, evecs, azimuth, _ = self.pca(
718 [trs[ix], trs[iy]])
720 axes.my_stuff.append(
721 axes.annotate(
722 '%.0f\u00b0' % azimuth,
723 (0., 0.),
724 (self.font_size, self.font_size),
725 xycoords='axes fraction',
726 textcoords='offset points',
727 color=self.colors['pca_2d'],
728 alpha=0.5))
730 ell = self.draw_cov_ellipse(
731 evals[:2], evecs[:2, :2],
732 color=self.colors['pca_2d'], alpha=0.5)
734 axes.add_artist(ell)
735 axes.my_stuff.append(ell)
737 l1, = axes.plot(
738 [-amax*evecs[0, -1], amax*evecs[0, -1]],
739 [-amax*evecs[1, -1], amax*evecs[1, -1]],
740 color=self.colors['pca_2d'], alpha=0.8)
742 l2, = axes.plot(
743 [-amax*evecs[0, -2], amax*evecs[0, -2]],
744 [-amax*evecs[1, -2], amax*evecs[1, -2]],
745 color=self.colors['pca_2d'], alpha=0.2)
747 axes.my_stuff.extend([l1, l2])
749 except PCAError as e:
750 logger.warn('PCA failed: %s' % e)
752 if self.show_eigensystem_3d:
753 try:
754 cov, evals, evecs, azimuth, incidence = self.pca(
755 trs_orig_chopped)
756 cosa = num.cos(bazimuth*d2r)
757 sina = num.sin(bazimuth*d2r)
758 rot = num.array(
759 [[cosa, -sina, 0.0],
760 [sina, cosa, 0.0],
761 [0.0, 0.0, 1.0]], dtype=num.float)
763 evecs_rot = num.dot(rot, evecs)
765 for (ix, iy, axes, evecs_) in [
766 (0, 1, axes_01, evecs),
767 (0, 2, axes_02, evecs_rot),
768 (1, 2, axes_12, evecs_rot)]:
770 for (ie, alpha) in [
771 (-1, 0.5),
772 (-2, 0.5),
773 (-3, 0.5)]:
775 for vlen, lw in [
776 (num.sqrt(evals[ie]), 3),
777 (amax, 1)]:
779 lv, = axes.plot(
780 [-vlen*evecs_[ix, ie],
781 vlen*evecs_[ix, ie]],
782 [-vlen*evecs_[iy, ie],
783 vlen*evecs_[iy, ie]],
784 lw=lw,
785 color=self.colors['pca_3d'], alpha=alpha)
787 axes.my_stuff.extend([lv])
789 axes_01.my_stuff.append(
790 axes_01.annotate(
791 '%.0f\u00b0' % azimuth,
792 (1., 0.),
793 (-self.font_size, self.font_size),
794 xycoords='axes fraction',
795 textcoords='offset points',
796 color=self.colors['pca_3d'],
797 alpha=0.8,
798 ha='right'))
800 axes_02.my_stuff.append(
801 axes_02.annotate(
802 '%.0f\u00b0' % incidence,
803 (1., 0.),
804 (-self.font_size, self.font_size),
805 xycoords='axes fraction',
806 textcoords='offset points',
807 color=self.colors['pca_3d'],
808 alpha=0.8,
809 ha='right'))
811 except PCAError as e:
812 logger.warn('PCA failed: %s' % e)
814 def get_bazis(self):
815 event = self.get_viewer().get_active_event()
816 if not event:
817 return {}
819 nsl_to_bazi = dict(
820 (station.nsl(), event.azibazi_to(station)[1])
821 for station in self.get_stations())
823 return nsl_to_bazi
825 def layout(self, fig, axes):
827 # Do not access self in here. Called from resize in finished plots.
829 def get_pixels_factor(fig):
830 try:
831 r = fig.canvas.get_renderer()
832 return 1.0 / r.points_to_pixels(1.0)
833 except AttributeError:
834 return 1.0
836 def rect_to_figure_coords(rect):
837 l, b, w, h = rect
838 return (l / width, b / height, w / width, h / height)
840 ny = len(axes)
841 if ny == 0:
842 raise LayoutError('No axes given.')
844 nx = len(axes[0])
846 width, height = fig.canvas.get_width_height()
847 pixels = get_pixels_factor(fig)
849 margin_left = margin_right = 4. * self.font_size / pixels
850 margin_top = margin_bottom = 5. * self.font_size / pixels
852 spacing_width = 3. * self.font_size / pixels
853 spacing_height = 5. * self.font_size / pixels
855 axes_height_avail = height - (ny - 1) * spacing_height \
856 - margin_top - margin_bottom
858 a_min = 5.
860 if axes_height_avail < a_min * ny:
861 axes_height_avail = a_min * ny
862 logger.warn('Not enough space vertically.')
864 axes_width_avail = width - (nx - 1) * spacing_width \
865 - margin_left - margin_right
867 if axes_width_avail < a_min * (nx + 2):
868 axes_width_avail = a_min * (nx + 2)
869 logger.warn('Not enough space horizontally.')
871 a_height = axes_height_avail / ny
872 a_width = axes_width_avail / (nx + 2)
874 a = max(min(a_height, a_width), a_min)
876 pad_height = (a_height - a) * ny
877 pad_width = (a_width - a) * (nx + 2)
879 for iy in range(ny):
880 y = height - 0.5 * pad_height - margin_top \
881 - (iy + 1) * a - iy * spacing_height
882 h = a
883 for ix in range(nx):
884 x = margin_right + 0.5 * pad_width + ix * (a + spacing_width)
885 w = a if ix != (nx - 1) else a * 3.0
886 axes[iy][ix].set_position(
887 rect_to_figure_coords((x, y, w, h)), which='both')
889 def call(self):
891 self.cleanup()
893 if self.rotate_to_rt != self.last_rotate_to_rt:
894 # reset delta azimuth to avoid confusion
896 self.set_parameter('azimuth', 0.0)
898 self.last_rotate_to_rt = self.rotate_to_rt
900 selected = self.get_selected()
902 nsl_to_bazi = self.get_bazis()
904 tpad = self.t_length
905 groups = self.get_traces(selected, nsl_to_bazi, tpad)
907 if not groups:
908 self.fail(
909 'No matching traces. Check time and channel settings. Traces '
910 'may not contain gaps within the extracted time window and in '
911 'the padding areas left and right. Traces are extracted with '
912 'additional padding of 3 x filter corner period to eliminate '
913 'artifacts.')
915 selection_markers = self.make_selection_markers(selected, groups)
916 self.add_markers(selection_markers)
918 fig = self.get_figure()
920 figure_key = (len(groups), self.iframe)
922 if not self.figure_key or self.figure_key != figure_key:
923 self.setup_figure(fig, len(groups))
924 self.figure_key = figure_key
926 self.draw(groups, selected, tpad, nsl_to_bazi)
928 self.layout(fig, self.axes)
930 self.fframe.draw()
931 tabs = self.fframe.parent().parent()
932 # bring plot to front if we are not looking at the markers
933 from pyrocko.gui.pile_viewer import PileViewer
934 if not isinstance(tabs.currentWidget(), PileViewer):
935 tabs.setCurrentWidget(self.fframe)
938def __snufflings__():
939 return [Polarization()]