1import logging
2import numpy as num
3import matplotlib
4from matplotlib import patches
5from pyrocko import util, trace, plot
6from ..snuffling import Snuffling, Param, Marker, Switch, EventMarker
8logger = logging.getLogger('pyrocko.gui.snufflings.polarization')
10d2r = num.pi / 180.
11r2d = 1.0 / d2r
14def get_tr(traces, f):
15 for tr in traces:
16 if f(tr):
17 return tr
19 raise ValueError('Element not found.')
22def darken(color, f=0.5):
23 return tuple(c*f for c in color[:3]) + color[3:]
26class PCAError(Exception):
27 pass
30class LayoutError(Exception):
31 pass
34class Polarization(Snuffling):
36 u'''
37Polarization
38============
40Investigate patterns of ground motion during the passage of seismic waves.
42This Snuffling can be used to analyze and visualize the polarization of seismic
43waves from 3-component seismic recordings or to check component orientations of
44a seismic sensors when used on signals with known directional properties. The
45spatial pattern of ground motion is shown in horizontal and vertical
46projections. Principal component analysis and rotation to radial/transverse
47components are available as tools. Time window and filter settings can be
48adjusted interactively.
50Usage
51-----
53Select one or more normal/phase markers as anchor points for the extraction and
54press *Run*. Multiple stations can be selected for direct comparison. Channels
55matching the pattern-triplet given in the *Channels* setting are selected for
56extraction of a time window around the anchor point. It is assumed that these
57three channels correspond to sensor components in the order east, north,
58vertical (upward), even if they are named differently.
60The time window can be adjusted with the *Length* and *Offset* parameters.
61Extracted waveforms are filtered according the *Highpass* and *Lowpass*
62parameters (Butterworth 4th order) and demeaned.
64To rotate the horizontal components around a vertical axis, use the *\u0394
65Azimuth* setting. When station coordinates and an "active" event are available,
66the horizontal components can be rotated to radial (away from event) and
67transverse (leftwards) orientations using the computed event back-azimuth
68(dashed gray line).
70If *Show 2D Eigensystems* is selected, a principal component analysis is
71performed in each of the shown projects. Red lines are shown to depict the
72eigenvectors and the eigenvalues are visualized using ellipse symbols. If *Show
733D Eigensystems* is selected, a principal component analysis is performed using
74all three (non-rotated) components and the resulting eigenvectors are depicted
75with purple lines.
77By default the scaling is automatically adjusted to the maximum value (the
78maximum vector length of the three-component signal within the selected time
79window). The currently used scaling factor can be frozen by checking *Fix
80Scale*.
81'''
83 def setup(self):
84 self.set_name('Polarization')
86 self.add_parameter(
87 Param('Length [s]', 't_length', 1., 0.001, 1000.))
89 self.add_parameter(
90 Param('Offset (relative)', 'ft_offset', 0., -5., 5.))
92 self.add_parameter(
93 Param(u'\u0394 Azimuth', 'azimuth', 0., -180., 180.))
95 self.add_parameter(
96 Param(
97 'Highpass [Hz]', 'highpass', None, 0.001, 1000.,
98 low_is_none=True))
100 self.add_parameter(
101 Param(
102 'Lowpass [Hz]', 'lowpass', None, 0.001, 1000.,
103 high_is_none=True))
105 self.add_parameter(
106 Param('Dot Position', 'dot_position', 0., 0., 1.))
108 self.add_parameter(
109 Switch(
110 'Rotate to RT',
111 'rotate_to_rt',
112 False))
114 self.add_parameter(
115 Switch(
116 'Show 2D Eigensystems',
117 'show_eigensystem_2d',
118 False))
120 self.add_parameter(
121 Switch(
122 'Show 3D Eigensystem',
123 'show_eigensystem_3d',
124 False))
126 self.add_parameter(
127 Switch(
128 'Fix Scale',
129 'fix_scale',
130 False))
132 def new_figure():
133 self.setup_figure_frame()
134 self.call()
136 self.add_trigger('New Figure', new_figure)
138 self.fframe = None
139 self.iframe = 0
140 self.figure_key = None
141 self.nsl_to_amax = {}
142 self.last_rotate_to_rt = False
144 self.font_size = float(matplotlib.rcParams['font.size'])
146 self.colors = {
147 'fg': matplotlib.rcParams['axes.edgecolor'],
148 'lines': 'black',
149 'rot': plot.mpl_color('skyblue1'),
150 'pca_2d': plot.mpl_color('scarletred1'),
151 'pca_3d': plot.mpl_color('plum2'),
152 'backazimuth': plot.mpl_color('aluminium3'),
153 'time_window': plot.mpl_color('aluminium2')}
155 def panel_visibility_changed(self, visible):
156 viewer = self.get_viewer()
157 if visible:
158 viewer.pile_has_changed_signal.connect(self.adjust_controls)
159 self.adjust_controls()
161 else:
162 viewer.pile_has_changed_signal.disconnect(self.adjust_controls)
164 def adjust_controls(self):
165 viewer = self.get_viewer()
166 dtmin, dtmax = viewer.content_deltat_range()
167 maxfreq = 0.5/dtmin
168 minfreq = (0.5/dtmax)*0.001
169 self.set_parameter_range('lowpass', minfreq, maxfreq)
170 self.set_parameter_range('highpass', minfreq, maxfreq)
172 def get_selected(self):
173 markers = [
174 marker for marker in self.get_selected_markers()
175 if not isinstance(marker, EventMarker)]
177 if not markers:
178 self.fail(
179 'No selected markers.\n\nCreate and select markers at points '
180 'of interest. Normal and phase markers are accepted.')
182 d = {}
183 for marker in markers:
184 tspan = self.transform_time_span(marker.tmin, marker.tmax)
185 for nslc in marker.nslc_ids:
186 k = nslc[:3] + (nslc[3][:-1], marker.tmin)
187 d[k] = tspan
189 return d
191 def transform_time_span(self, tmin, tmax):
192 tmin = tmin + self.t_length * self.ft_offset
193 tmax = tmin + self.t_length
194 return (tmin, tmax)
196 def make_selection_markers(self, selected, groups):
197 selection_markers = []
198 for k in sorted(selected):
199 tmin, tmax = selected[k]
200 patterns = [tr.nslc_id for tr in groups[k]]
201 marker = Marker(
202 nslc_ids=patterns, tmin=tmin, tmax=tmax, kind=2)
203 selection_markers.append(marker)
205 return selection_markers
207 def sorted_by_component_scheme(self, group):
208 if len(group) not in (2, 3):
209 self.fail('Need 2 or 3 component recordings. Got %i.' % len(group))
211 comps = []
212 channels = []
213 for tr in group:
214 comps.append(tr.channel[-1])
215 channels.append(tr.channel)
217 # special case for channel naming convention of Cube dataloggers
218 for scheme in [['p2', 'p1', 'p0'], ['p2', 'p1']]:
220 if sorted(channels) == sorted(scheme):
221 return [
222 get_tr(group, lambda tr: tr.channel == ch)
223 for ch in scheme]
225 for scheme in [
226 ['E', 'N', 'Z'], ['1', '2', 'Z'], ['1', '2', '3'],
227 ['0', '1', '2'], ['R', 'T', 'Z'],
228 ['E', 'N'], ['1', '2'], ['0', '1']]:
230 if sorted(comps) == sorted(scheme):
231 return [
232 get_tr(group, lambda tr: tr.channel[-1] == comp)
233 for comp in scheme]
235 self.fail('Cannot handle component group: %s' % ', '.join(channels))
237 def get_traces(self, selected, nsl_to_bazi, tpad):
238 if self.highpass is not None:
239 tpad_filter = 3.0 / self.highpass
240 elif self.lowpass is not None:
241 tpad_filter = 3.0 / self.lowpass
242 else:
243 tpad_filter = 0.0
245 # prevent getting insanely long cutouts if e.g. if highpass is still
246 # very low, e.g. while moving the slider
247 tpad_filter = min(tpad_filter, 5.0 * self.t_length)
249 d = {}
250 for k in sorted(selected):
251 nsl = k[0:3]
252 tmin, tmax = selected[k]
254 bazimuth = self.azimuth
256 if self.rotate_to_rt:
257 if nsl not in nsl_to_bazi:
258 self.fail(
259 'Cannot rotate to RT.\n\nStation coordinates must be '
260 'available and an event must be marked as the '
261 '"active" event (select event and press "e").')
263 bazimuth += nsl_to_bazi[nsl] + 90.
265 pattern = '.'.join(nsl + (k[3] + '?',))
267 group = []
268 for batch in self.get_pile().chopper(
269 tmin=tmin, tmax=tmax, tpad=tpad + tpad_filter,
270 want_incomplete=False,
271 style='batch',
272 trace_selector=lambda tr: util.match_nslc(
273 pattern, tr.nslc_id)):
275 for tr in batch.traces:
276 tr = tr.copy()
277 if self.lowpass is not None \
278 and self.lowpass < 0.5/tr.deltat:
279 tr.lowpass(4, self.lowpass)
281 if self.highpass is not None \
282 and self.highpass < 0.5/tr.deltat:
284 tr.highpass(4, self.highpass)
286 try:
287 tr.chop(tmin - tpad, tmax + tpad)
288 tr_chop = tr.chop(tmin, tmax, inplace=False)
289 except trace.NoData:
290 self.fail('No data. Time length too short?')
292 y = tr.get_ydata()
293 tr.set_ydata(y - num.mean(tr_chop.get_ydata()))
295 group.append(tr)
297 group = self.sorted_by_component_scheme(group)
299 tr_e = group[0]
300 tr_n = group[1]
302 if tr_e is not None and tr_n is not None:
303 cha_e = tr_e.channel
304 cha_n = tr_n.channel
305 if bazimuth != 0.0:
306 trs_rot = trace.rotate(
307 [tr_n, tr_e], bazimuth,
308 in_channels=[cha_n, cha_e],
309 out_channels=[ch+"'" for ch in [cha_n, cha_e]])
310 else:
311 trs_rot = [tr.copy() for tr in [tr_n, tr_e]]
312 for tr in trs_rot:
313 tr.set_codes(channel=tr.channel + "'")
315 if len(trs_rot) == 2:
316 group.extend(
317 get_tr(trs_rot, lambda tr: tr.channel == ch+"'")
318 for ch in [cha_e, cha_n])
319 else:
320 self.fail(
321 'Cannot rotate traces. '
322 'Are the samples properly aligned?')
324 d[k] = group
326 return d
328 def setup_figure_frame(self):
329 self.iframe += 1
330 self.fframe = self.figure_frame(
331 'Polarization (%i)' % self.iframe)
333 self.fframe.gcf().my_disconnect = None
335 def get_figure(self):
336 if not self.fframe or self.fframe.closed:
337 self.setup_figure_frame()
339 return self.fframe.gcf()
341 def setup_figure(self, fig, nstations):
342 fig.clf()
343 new_axes = []
345 def iwrap(iy, ix):
346 return (ix + iy * 4) + 1
348 for istation in range(nstations):
349 axes_01 = fig.add_subplot(
350 nstations, 4, iwrap(istation, 0), aspect=1.0)
351 axes_02 = fig.add_subplot(
352 nstations, 4, iwrap(istation, 1), aspect=1.0)
353 axes_12 = fig.add_subplot(
354 nstations, 4, iwrap(istation, 2), aspect=1.0)
355 axes_tr = fig.add_subplot(
356 nstations, 4, iwrap(istation, 3))
358 for axes in (axes_01, axes_02, axes_12, axes_tr):
359 axes.my_stuff = []
361 axes.my_line, = axes.plot(
362 [], [], color=self.colors['lines'], lw=1.0)
363 axes.my_dot, = axes.plot(
364 [], [], 'o', ms=4, color=self.colors['lines'])
366 for axes in (axes_01, axes_02, axes_12):
367 axes.get_xaxis().set_tick_params(
368 labelbottom=False, bottom=False)
369 axes.get_yaxis().set_tick_params(
370 labelleft=False, left=False)
372 axes_tr.get_yaxis().set_tick_params(
373 left=False, labelleft=True, length=2.0)
375 # if istation != nstations - 1:
376 # axes_tr.get_xaxis().set_tick_params(
377 # bottom=False, labelbottom=False)
379 lines = []
380 dots = []
381 for _ in range(3):
382 lines.append(
383 axes_tr.plot(
384 [], [], color=self.colors['lines'], lw=1.0)[0])
385 dots.append(
386 axes_tr.plot(
387 [], [], 'o', ms=4, color=self.colors['lines'])[0])
389 axes_tr.my_lines = lines
390 axes_tr.my_dots = dots
391 axes_tr.my_stuff = []
393 new_axes.append(
394 (axes_01, axes_02, axes_12, axes_tr))
396 def resize_handler(*args):
397 self.layout(fig, new_axes)
399 if fig.my_disconnect:
400 fig.my_disconnect()
402 cid_resize = fig.canvas.mpl_connect('resize_event', resize_handler)
403 cid_dpi = fig.callbacks.connect('dpi_changed', resize_handler)
405 def disconnect():
406 fig.canvas.mpl_disconnect(cid_resize)
407 fig.callbacks.disconnect(cid_dpi)
409 fig.my_disconnect = disconnect
411 self.axes = new_axes
413 def get_trace(self, traces, pattern):
414 trs = [tr for tr in traces if util.match_nslc([pattern], tr.nslc_id)]
415 if len(trs) > 1:
416 self.fail('Multiple traces matching pattern %s' % pattern)
417 elif len(trs) == 0:
418 return None
419 else:
420 return trs[0]
422 def get_vector_abs_max(self, traces):
423 tr_abs = None
424 for tr in traces:
425 if tr is not None:
426 tr = tr.copy()
427 tr.ydata **= 2
428 if tr_abs is None:
429 tr_abs = tr
430 else:
431 tr_abs.add(tr)
433 tr_abs.set_ydata(num.sqrt(tr_abs.ydata))
434 return num.max(tr_abs.ydata)
436 def set_labels(
437 self, istation, nstations, tref,
438 axes_01, axes_02, axes_12, axes_tr, chas, chas_rot):
440 if self.azimuth != 0 or self.rotate_to_rt:
441 rcolor = self.colors['rot']
442 else:
443 rcolor = self.colors['fg']
444 chas_rot = chas
446 axes_01.set_ylabel(chas[1])
447 axes_02.set_ylabel(chas[2])
448 axes_12.set_ylabel(chas[2])
450 axes_01.set_xlabel(chas[0])
451 axes_02.set_xlabel(chas_rot[0])
452 axes_12.set_xlabel(chas_rot[1])
453 axes_02.get_xaxis().label.set_color(rcolor)
454 axes_12.get_xaxis().label.set_color(rcolor)
455 axes_tr.set_xlabel(
456 'Time [s] relative to %s' % util.time_to_str(tref))
458 axes_tr.set_yticks([0., 1., 2])
459 axes_tr.set_yticklabels([chas_rot[0], chas_rot[1], chas[2]])
460 for tlab in axes_tr.get_yticklabels()[:2]:
461 tlab.set_color(rcolor)
463 def pca(self, trs):
465 if any(tr is None for tr in trs):
466 raise PCAError('Missing component')
468 nss = [tr.data_len() for tr in trs]
469 if not all(ns == nss[0] for ns in nss):
470 raise PCAError('Traces have different lengths.')
472 ns = nss[0]
474 if ns < 3:
475 raise PCAError('Traces too short.')
477 data = num.zeros((len(trs), ns))
478 for itr, tr in enumerate(trs):
479 data[itr, :] = tr.ydata
481 cov = num.cov(data)
483 evals, evecs = num.linalg.eigh(cov)
485 pc = evecs[:, -1]
487 eh = num.sqrt(pc[1]**2 + pc[0]**2)
489 if len(trs) > 2:
490 incidence = r2d * num.arctan2(eh, abs(pc[2]))
491 else:
492 incidence = 90.
494 azimuth = r2d*num.arctan2(pc[1], pc[0])
495 azimuth = (-azimuth) % 180. - 90.
497 return cov, evals, evecs, azimuth, incidence
499 def draw_cov_ellipse(self, evals, evecs, color, alpha=1.0):
500 evals = num.sqrt(evals)
501 ell = patches.Ellipse(
502 xy=(0.0, 0.0),
503 width=evals[0] * 2.,
504 height=evals[1] * 2.,
505 angle=r2d*num.arctan2(evecs[1, 0], evecs[0, 0]),
506 zorder=-10,
507 fc=color,
508 ec=darken(color),
509 alpha=alpha)
511 return ell
513 def draw(self, groups, nsl_to_tspan, tpad, nsl_to_bazi):
515 def count(d, k):
516 if k not in d:
517 d[k] = 0
519 d[k] += 1
521 nsli_n = {}
522 for k in sorted(groups):
523 count(nsli_n, k[:4])
525 nsli_i = {}
526 for igroup, k in enumerate(sorted(groups)):
528 tmin, tmax = nsl_to_tspan[k]
529 nsl = k[:3]
530 nsli = k[:4]
531 count(nsli_i, nsli)
533 bazimuth = self.azimuth
535 if self.rotate_to_rt:
536 if nsl not in nsl_to_bazi:
537 self.fail(
538 'Cannot rotate to RT.\n\nActive event must '
539 'available (select event and press "e"). Station '
540 'coordinates must be available.')
542 bazimuth += nsl_to_bazi[nsl] + 90.
544 for axes in self.axes[igroup]:
545 while axes.my_stuff:
546 stuff = axes.my_stuff.pop()
547 stuff.remove()
549 axes_01, axes_02, axes_12, axes_tr = self.axes[igroup]
551 axes_01.set_title(
552 '.'.join(nsli)
553 + ('' if nsli_n[nsli] == 1 else ' (%i)' % nsli_i[nsli]))
555 trs_all = groups[k]
557 if len(trs_all) == 5:
558 trs_orig = trs_all[:3]
559 trs_rot = trs_all[3:] + trs_all[2:3]
560 elif len(trs_all) == 4:
561 trs_orig = trs_all[:2] + [None]
562 trs_rot = trs_all[2:] + [None]
563 else:
564 raise self.fail(
565 'Unexpectedly got other that 4 or 6 traces: %i'
566 % len(trs_all))
568 trs_orig_chopped = [
569 (tr.chop(tmin, tmax, inplace=False) if tr else None)
570 for tr in trs_orig]
572 trs_rot_chopped = [
573 (tr.chop(tmin, tmax, inplace=False) if tr else None)
574 for tr in trs_rot]
576 chas = [tr.channel[-1:] for tr in trs_orig_chopped]
577 chas_rot = [tr.channel[-2:] for tr in trs_rot_chopped]
579 if self.fix_scale and nsl in self.nsl_to_amax:
580 amax = self.nsl_to_amax[nsl]
581 else:
582 amax = self.get_vector_abs_max(trs_orig_chopped)
583 self.nsl_to_amax[nsl] = amax
585 if nsl in nsl_to_bazi:
586 event_bazi = nsl_to_bazi[nsl]
587 l1, = axes_01.plot(
588 [0., amax*num.cos((90. - event_bazi)*d2r)],
589 [0., amax*num.sin((90. - event_bazi)*d2r)],
590 '--',
591 lw=3,
592 zorder=-20,
593 color=self.colors['backazimuth'])
595 a1 = axes_01.annotate(
596 '%.0f\u00b0' % event_bazi,
597 (0., 1.),
598 (self.font_size, -self.font_size),
599 xycoords='axes fraction',
600 textcoords='offset points',
601 color=self.colors['backazimuth'],
602 ha='left',
603 va='top')
605 axes_01.my_stuff.extend([l1, a1])
607 for ix, iy, axes, trs in [
608 (0, 1, axes_01, trs_orig_chopped),
609 (0, 2, axes_02, trs_rot_chopped),
610 (1, 2, axes_12, trs_rot_chopped)]:
612 axes.set_xlim(-amax*1.05, amax*1.05)
613 axes.set_ylim(-amax*1.05, amax*1.05)
615 if not (trs[ix] and trs[iy]):
616 continue
618 x = trs[ix].get_ydata()
619 y = trs[iy].get_ydata()
621 axes.my_line.set_data(x, y)
622 ipos = int(round(self.dot_position * (x.size-1)))
623 axes.my_dot.set_data(x[ipos], y[ipos])
625 tref = tmin
626 for itr, (tr, tr_chopped) in enumerate(zip(
627 trs_rot, trs_rot_chopped)):
629 if tr is None or tr_chopped is None:
630 axes_tr.my_lines[itr].set_data([], [])
631 axes_tr.my_dots[itr].set_data([], [])
633 else:
634 y = tr.get_ydata() / (2.*amax) + itr
635 t = tr.get_xdata()
636 t = t - tref
638 ipos = int(round(
639 self.dot_position * (tr_chopped.data_len()-1)))
641 yp = tr_chopped.ydata[ipos] / (2.*amax) + itr
642 tp = tr_chopped.tmin - tref + tr_chopped.deltat*ipos
644 axes_tr.my_lines[itr].set_data(t, y)
645 axes_tr.my_dots[itr].set_data(tp, yp)
647 if self.azimuth != 0.0 or self.rotate_to_rt:
648 color = self.colors['rot']
650 xn = num.sin(bazimuth*d2r)
651 yn = num.cos(bazimuth*d2r)
652 xe = num.sin(bazimuth*d2r + 0.5*num.pi)
653 ye = num.cos(bazimuth*d2r + 0.5*num.pi)
655 l1, = axes_01.plot(
656 [0., amax*xn],
657 [0., amax*yn],
658 color=color)
660 font_size = self.font_size
662 a1 = axes_01.annotate(
663 chas_rot[1],
664 xy=(amax*xn, amax*yn),
665 xycoords='data',
666 xytext=(-font_size*(xe+.5*xn), -font_size*(ye+.5*yn)),
667 textcoords='offset points',
668 va='center',
669 ha='center',
670 color=color)
672 l2, = axes_01.plot(
673 [0., amax*xe],
674 [0., amax*ye],
675 color=color)
677 a2 = axes_01.annotate(
678 chas_rot[0],
679 xy=(amax*xe, amax*ye),
680 xycoords='data',
681 xytext=(-font_size*(xn+.5*xe), -font_size*(yn+.5*ye)),
682 textcoords='offset points',
683 va='center',
684 ha='center',
685 color=color)
687 aa = axes_01.annotate(
688 '%.0f\u00b0' % bazimuth,
689 (1., 1.),
690 (-self.font_size, -self.font_size),
691 xycoords='axes fraction',
692 textcoords='offset points',
693 color=color,
694 ha='right',
695 va='top')
697 axes_01.my_stuff.extend([l1, a1, l2, a2, aa])
699 axes_tr.my_stuff.append(axes_tr.axvspan(
700 tmin - tref, tmax - tref, color=self.colors['time_window']))
702 axes_tr.set_ylim(-1, 3)
703 axes_tr.set_xlim(tmin - tref - tpad, tmax - tref + tpad)
705 self.set_labels(
706 igroup, len(groups), tref, *self.axes[igroup], chas, chas_rot)
708 if self.show_eigensystem_2d:
710 for (ix, iy, axes, trs) in [
711 (0, 1, axes_01, trs_orig_chopped),
712 (0, 2, axes_02, trs_rot_chopped),
713 (1, 2, axes_12, trs_rot_chopped)]:
715 try:
716 cov, evals, evecs, azimuth, _ = self.pca(
717 [trs[ix], trs[iy]])
719 axes.my_stuff.append(
720 axes.annotate(
721 '%.0f\u00b0' % azimuth,
722 (0., 0.),
723 (self.font_size, self.font_size),
724 xycoords='axes fraction',
725 textcoords='offset points',
726 color=self.colors['pca_2d'],
727 alpha=0.5))
729 ell = self.draw_cov_ellipse(
730 evals[:2], evecs[:2, :2],
731 color=self.colors['pca_2d'], alpha=0.5)
733 axes.add_artist(ell)
734 axes.my_stuff.append(ell)
736 l1, = axes.plot(
737 [-amax*evecs[0, -1], amax*evecs[0, -1]],
738 [-amax*evecs[1, -1], amax*evecs[1, -1]],
739 color=self.colors['pca_2d'], alpha=0.8)
741 l2, = axes.plot(
742 [-amax*evecs[0, -2], amax*evecs[0, -2]],
743 [-amax*evecs[1, -2], amax*evecs[1, -2]],
744 color=self.colors['pca_2d'], alpha=0.2)
746 axes.my_stuff.extend([l1, l2])
748 except PCAError as e:
749 logger.warning('PCA failed: %s' % e)
751 if self.show_eigensystem_3d:
752 try:
753 cov, evals, evecs, azimuth, incidence = self.pca(
754 trs_orig_chopped)
755 cosa = num.cos(bazimuth*d2r)
756 sina = num.sin(bazimuth*d2r)
757 rot = num.array(
758 [[cosa, -sina, 0.0],
759 [sina, cosa, 0.0],
760 [0.0, 0.0, 1.0]], dtype=float)
762 evecs_rot = num.dot(rot, evecs)
764 for (ix, iy, axes, evecs_) in [
765 (0, 1, axes_01, evecs),
766 (0, 2, axes_02, evecs_rot),
767 (1, 2, axes_12, evecs_rot)]:
769 for (ie, alpha) in [
770 (-1, 0.5),
771 (-2, 0.5),
772 (-3, 0.5)]:
774 for vlen, lw in [
775 (num.sqrt(evals[ie]), 3),
776 (amax, 1)]:
778 lv, = axes.plot(
779 [-vlen*evecs_[ix, ie],
780 vlen*evecs_[ix, ie]],
781 [-vlen*evecs_[iy, ie],
782 vlen*evecs_[iy, ie]],
783 lw=lw,
784 color=self.colors['pca_3d'], alpha=alpha)
786 axes.my_stuff.extend([lv])
788 axes_01.my_stuff.append(
789 axes_01.annotate(
790 '%.0f\u00b0' % azimuth,
791 (1., 0.),
792 (-self.font_size, self.font_size),
793 xycoords='axes fraction',
794 textcoords='offset points',
795 color=self.colors['pca_3d'],
796 alpha=0.8,
797 ha='right'))
799 axes_02.my_stuff.append(
800 axes_02.annotate(
801 '%.0f\u00b0' % incidence,
802 (1., 0.),
803 (-self.font_size, self.font_size),
804 xycoords='axes fraction',
805 textcoords='offset points',
806 color=self.colors['pca_3d'],
807 alpha=0.8,
808 ha='right'))
810 except PCAError as e:
811 logger.warning('PCA failed: %s' % e)
813 def get_bazis(self):
814 event = self.get_viewer().get_active_event()
815 if not event:
816 return {}
818 nsl_to_bazi = dict(
819 (station.nsl(), event.azibazi_to(station)[1])
820 for station in self.get_stations())
822 return nsl_to_bazi
824 def layout(self, fig, axes):
826 # Do not access self in here. Called from resize in finished plots.
828 def get_pixels_factor(fig):
829 try:
830 r = fig.canvas.get_renderer()
831 return 1.0 / r.points_to_pixels(1.0)
832 except AttributeError:
833 return 1.0
835 def rect_to_figure_coords(rect):
836 l, b, w, h = rect # noqa
837 return (l / width, b / height, w / width, h / height) # noqa
839 ny = len(axes)
840 if ny == 0:
841 raise LayoutError('No axes given.')
843 nx = len(axes[0])
845 width, height = fig.canvas.get_width_height()
846 pixels = get_pixels_factor(fig)
848 margin_left = margin_right = 4. * self.font_size / pixels
849 margin_top = margin_bottom = 5. * self.font_size / pixels
851 spacing_width = 3. * self.font_size / pixels
852 spacing_height = 5. * self.font_size / pixels
854 axes_height_avail = height - (ny - 1) * spacing_height \
855 - margin_top - margin_bottom
857 a_min = 5.
859 if axes_height_avail < a_min * ny:
860 axes_height_avail = a_min * ny
861 logger.warning('Not enough space vertically.')
863 axes_width_avail = width - (nx - 1) * spacing_width \
864 - margin_left - margin_right
866 if axes_width_avail < a_min * (nx + 2):
867 axes_width_avail = a_min * (nx + 2)
868 logger.warning('Not enough space horizontally.')
870 a_height = axes_height_avail / ny
871 a_width = axes_width_avail / (nx + 2)
873 a = max(min(a_height, a_width), a_min)
875 pad_height = (a_height - a) * ny
876 pad_width = (a_width - a) * (nx + 2)
878 for iy in range(ny):
879 y = height - 0.5 * pad_height - margin_top \
880 - (iy + 1) * a - iy * spacing_height
881 h = a
882 for ix in range(nx):
883 x = margin_right + 0.5 * pad_width + ix * (a + spacing_width)
884 w = a if ix != (nx - 1) else a * 3.0
885 axes[iy][ix].set_position(
886 rect_to_figure_coords((x, y, w, h)), which='both')
888 def call(self):
890 self.cleanup()
892 if self.rotate_to_rt != self.last_rotate_to_rt:
893 # reset delta azimuth to avoid confusion
895 self.set_parameter('azimuth', 0.0)
897 self.last_rotate_to_rt = self.rotate_to_rt
899 selected = self.get_selected()
901 nsl_to_bazi = self.get_bazis()
903 tpad = self.t_length
904 groups = self.get_traces(selected, nsl_to_bazi, tpad)
906 if not groups:
907 self.fail(
908 'No matching traces. Check time and channel settings. Traces '
909 'may not contain gaps within the extracted time window and in '
910 'the padding areas left and right. Traces are extracted with '
911 'additional padding of 3 x filter corner period to eliminate '
912 'artifacts.')
914 selection_markers = self.make_selection_markers(selected, groups)
915 self.add_markers(selection_markers)
917 fig = self.get_figure()
919 figure_key = (len(groups), self.iframe)
921 if not self.figure_key or self.figure_key != figure_key:
922 self.setup_figure(fig, len(groups))
923 self.figure_key = figure_key
925 self.draw(groups, selected, tpad, nsl_to_bazi)
927 self.layout(fig, self.axes)
929 self.fframe.draw()
930 tabs = self.fframe.parent().parent()
931 # bring plot to front if we are not looking at the markers
932 from pyrocko.gui.pile_viewer import PileViewer
933 if not isinstance(tabs.currentWidget(), PileViewer):
934 tabs.setCurrentWidget(self.fframe)
937def __snufflings__():
938 return [Polarization()]