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_callbacks(obj):
15 try:
16 return obj.callbacks
17 except AttributeError:
18 return obj._callbacks
21def get_tr(traces, f):
22 for tr in traces:
23 if f(tr):
24 return tr
26 raise ValueError('Element not found.')
29def darken(color, f=0.5):
30 return tuple(c*f for c in color[:3]) + color[3:]
33class PCAError(Exception):
34 pass
37class LayoutError(Exception):
38 pass
41class Polarization(Snuffling):
43 u'''
44Polarization
45============
47Investigate patterns of ground motion during the passage of seismic waves.
49This Snuffling can be used to analyze and visualize the polarization of seismic
50waves from 3-component seismic recordings or to check component orientations of
51a seismic sensors when used on signals with known directional properties. The
52spatial pattern of ground motion is shown in horizontal and vertical
53projections. Principal component analysis and rotation to radial/transverse
54components are available as tools. Time window and filter settings can be
55adjusted interactively.
57Usage
58-----
60Select one or more normal/phase markers as anchor points for the extraction and
61press *Run*. Multiple stations can be selected for direct comparison. Channels
62matching the pattern-triplet given in the *Channels* setting are selected for
63extraction of a time window around the anchor point. It is assumed that these
64three channels correspond to sensor components in the order east, north,
65vertical (upward), even if they are named differently.
67The time window can be adjusted with the *Length* and *Offset* parameters.
68Extracted waveforms are filtered according the *Highpass* and *Lowpass*
69parameters (Butterworth 4th order) and demeaned.
71To rotate the horizontal components around a vertical axis, use the *\u0394
72Azimuth* setting. When station coordinates and an "active" event are available,
73the horizontal components can be rotated to radial (away from event) and
74transverse (leftwards) orientations using the computed event back-azimuth
75(dashed gray line).
77If *Show 2D Eigensystems* is selected, a principal component analysis is
78performed in each of the shown projects. Red lines are shown to depict the
79eigenvectors and the eigenvalues are visualized using ellipse symbols. If *Show
803D Eigensystems* is selected, a principal component analysis is performed using
81all three (non-rotated) components and the resulting eigenvectors are depicted
82with purple lines.
84By default the scaling is automatically adjusted to the maximum value (the
85maximum vector length of the three-component signal within the selected time
86window). The currently used scaling factor can be frozen by checking *Fix
87Scale*.
88'''
90 def setup(self):
91 self.set_name('Polarization')
93 self.add_parameter(
94 Param('Length [s]', 't_length', 1., 0.001, 1000.))
96 self.add_parameter(
97 Param('Offset (relative)', 'ft_offset', 0., -5., 5.))
99 self.add_parameter(
100 Param(u'\u0394 Azimuth', 'azimuth', 0., -180., 180.))
102 self.add_parameter(
103 Param(
104 'Highpass [Hz]', 'highpass', None, 0.001, 1000.,
105 low_is_none=True))
107 self.add_parameter(
108 Param(
109 'Lowpass [Hz]', 'lowpass', None, 0.001, 1000.,
110 high_is_none=True))
112 self.add_parameter(
113 Param('Dot Position', 'dot_position', 0., 0., 1.))
115 self.add_parameter(
116 Switch(
117 'Rotate to RT',
118 'rotate_to_rt',
119 False))
121 self.add_parameter(
122 Switch(
123 'Show 2D Eigensystems',
124 'show_eigensystem_2d',
125 False))
127 self.add_parameter(
128 Switch(
129 'Show 3D Eigensystem',
130 'show_eigensystem_3d',
131 False))
133 self.add_parameter(
134 Switch(
135 'Fix Scale',
136 'fix_scale',
137 False))
139 def new_figure():
140 self.setup_figure_frame()
141 self.call()
143 self.add_trigger('New Figure', new_figure)
145 self.fframe = None
146 self.iframe = 0
147 self.figure_key = None
148 self.nsl_to_amax = {}
149 self.last_rotate_to_rt = False
151 self.font_size = float(matplotlib.rcParams['font.size'])
153 self.colors = {
154 'fg': matplotlib.rcParams['axes.edgecolor'],
155 'lines': 'black',
156 'rot': plot.mpl_color('skyblue1'),
157 'pca_2d': plot.mpl_color('scarletred1'),
158 'pca_3d': plot.mpl_color('plum2'),
159 'backazimuth': plot.mpl_color('aluminium3'),
160 'time_window': plot.mpl_color('aluminium2')}
162 def panel_visibility_changed(self, visible):
163 viewer = self.get_viewer()
164 if visible:
165 viewer.pile_has_changed_signal.connect(self.adjust_controls)
166 self.adjust_controls()
168 else:
169 viewer.pile_has_changed_signal.disconnect(self.adjust_controls)
171 def adjust_controls(self):
172 viewer = self.get_viewer()
173 dtmin, dtmax = viewer.content_deltat_range()
174 maxfreq = 0.5/dtmin
175 minfreq = (0.5/dtmax)*0.001
176 self.set_parameter_range('lowpass', minfreq, maxfreq)
177 self.set_parameter_range('highpass', minfreq, maxfreq)
179 def get_selected(self):
180 markers = [
181 marker for marker in self.get_selected_markers()
182 if not isinstance(marker, EventMarker)]
184 if not markers:
185 self.fail(
186 'No selected markers.\n\nCreate and select markers at points '
187 'of interest. Normal and phase markers are accepted.')
189 d = {}
190 for marker in markers:
191 tspan = self.transform_time_span(marker.tmin, marker.tmax)
192 for nslc in marker.nslc_ids:
193 k = nslc[:3] + (nslc[3][:-1], marker.tmin)
194 d[k] = tspan
196 return d
198 def transform_time_span(self, tmin, tmax):
199 tmin = tmin + self.t_length * self.ft_offset
200 tmax = tmin + self.t_length
201 return (tmin, tmax)
203 def make_selection_markers(self, selected, groups):
204 selection_markers = []
205 for k in sorted(selected):
206 tmin, tmax = selected[k]
207 patterns = [tr.nslc_id for tr in groups[k]]
208 marker = Marker(
209 nslc_ids=patterns, tmin=tmin, tmax=tmax, kind=2)
210 selection_markers.append(marker)
212 return selection_markers
214 def sorted_by_component_scheme(self, group):
215 if len(group) not in (2, 3):
216 self.fail('Need 2 or 3 component recordings. Got %i.' % len(group))
218 comps = []
219 channels = []
220 for tr in group:
221 comps.append(tr.channel[-1])
222 channels.append(tr.channel)
224 # special case for channel naming convention of Cube dataloggers
225 for scheme in [['p2', 'p1', 'p0'], ['p2', 'p1']]:
227 if sorted(channels) == sorted(scheme):
228 return [
229 get_tr(group, lambda tr: tr.channel == ch)
230 for ch in scheme]
232 for scheme in [
233 ['E', 'N', 'Z'], ['1', '2', 'Z'], ['1', '2', '3'],
234 ['0', '1', '2'], ['R', 'T', 'Z'],
235 ['E', 'N'], ['1', '2'], ['0', '1']]:
237 if sorted(comps) == sorted(scheme):
238 return [
239 get_tr(group, lambda tr: tr.channel[-1] == comp)
240 for comp in scheme]
242 self.fail('Cannot handle component group: %s' % ', '.join(channels))
244 def get_traces(self, selected, nsl_to_bazi, tpad):
245 if self.highpass is not None:
246 tpad_filter = 3.0 / self.highpass
247 elif self.lowpass is not None:
248 tpad_filter = 3.0 / self.lowpass
249 else:
250 tpad_filter = 0.0
252 # prevent getting insanely long cutouts if e.g. if highpass is still
253 # very low, e.g. while moving the slider
254 tpad_filter = min(tpad_filter, 5.0 * self.t_length)
256 d = {}
257 for k in sorted(selected):
258 nsl = k[0:3]
259 tmin, tmax = selected[k]
261 bazimuth = self.azimuth
263 if self.rotate_to_rt:
264 if nsl not in nsl_to_bazi:
265 self.fail(
266 'Cannot rotate to RT.\n\nStation coordinates must be '
267 'available and an event must be marked as the '
268 '"active" event (select event and press "e").')
270 bazimuth += nsl_to_bazi[nsl] + 90.
272 pattern = '.'.join(nsl + (k[3] + '?',))
274 group = []
275 for batch in self.get_pile().chopper(
276 tmin=tmin, tmax=tmax, tpad=tpad + tpad_filter,
277 want_incomplete=False,
278 style='batch',
279 trace_selector=lambda tr: util.match_nslc(
280 pattern, tr.nslc_id)):
282 for tr in batch.traces:
283 tr = tr.copy()
284 if self.lowpass is not None \
285 and self.lowpass < 0.5/tr.deltat:
286 tr.lowpass(4, self.lowpass)
288 if self.highpass is not None \
289 and self.highpass < 0.5/tr.deltat:
291 tr.highpass(4, self.highpass)
293 try:
294 tr.chop(tmin - tpad, tmax + tpad)
295 tr_chop = tr.chop(tmin, tmax, inplace=False)
296 except trace.NoData:
297 self.fail('No data. Time length too short?')
299 y = tr.get_ydata()
300 tr.set_ydata(y - num.mean(tr_chop.get_ydata()))
302 group.append(tr)
304 group = self.sorted_by_component_scheme(group)
306 tr_e = group[0]
307 tr_n = group[1]
309 if tr_e is not None and tr_n is not None:
310 cha_e = tr_e.channel
311 cha_n = tr_n.channel
312 if bazimuth != 0.0:
313 trs_rot = trace.rotate(
314 [tr_n, tr_e], bazimuth,
315 in_channels=[cha_n, cha_e],
316 out_channels=[ch+"'" for ch in [cha_n, cha_e]])
317 else:
318 trs_rot = [tr.copy() for tr in [tr_n, tr_e]]
319 for tr in trs_rot:
320 tr.set_codes(channel=tr.channel + "'")
322 if len(trs_rot) == 2:
323 group.extend(
324 get_tr(trs_rot, lambda tr: tr.channel == ch+"'")
325 for ch in [cha_e, cha_n])
326 else:
327 self.fail(
328 'Cannot rotate traces. '
329 'Are the samples properly aligned?')
331 d[k] = group
333 return d
335 def setup_figure_frame(self):
336 self.iframe += 1
337 self.fframe = self.figure_frame(
338 'Polarization (%i)' % self.iframe)
340 self.fframe.gcf().my_disconnect = None
342 def get_figure(self):
343 if not self.fframe or self.fframe.closed:
344 self.setup_figure_frame()
346 return self.fframe.gcf()
348 def setup_figure(self, fig, nstations):
349 fig.clf()
350 new_axes = []
352 def iwrap(iy, ix):
353 return (ix + iy * 4) + 1
355 for istation in range(nstations):
356 axes_01 = fig.add_subplot(
357 nstations, 4, iwrap(istation, 0), aspect=1.0)
358 axes_02 = fig.add_subplot(
359 nstations, 4, iwrap(istation, 1), aspect=1.0)
360 axes_12 = fig.add_subplot(
361 nstations, 4, iwrap(istation, 2), aspect=1.0)
362 axes_tr = fig.add_subplot(
363 nstations, 4, iwrap(istation, 3))
365 for axes in (axes_01, axes_02, axes_12, axes_tr):
366 axes.my_stuff = []
368 axes.my_line, = axes.plot(
369 [], [], color=self.colors['lines'], lw=1.0)
370 axes.my_dot, = axes.plot(
371 [], [], 'o', ms=4, color=self.colors['lines'])
373 for axes in (axes_01, axes_02, axes_12):
374 axes.get_xaxis().set_tick_params(
375 labelbottom=False, bottom=False)
376 axes.get_yaxis().set_tick_params(
377 labelleft=False, left=False)
379 axes_tr.get_yaxis().set_tick_params(
380 left=False, labelleft=True, length=2.0)
382 # if istation != nstations - 1:
383 # axes_tr.get_xaxis().set_tick_params(
384 # bottom=False, labelbottom=False)
386 lines = []
387 dots = []
388 for _ in range(3):
389 lines.append(
390 axes_tr.plot(
391 [], [], color=self.colors['lines'], lw=1.0)[0])
392 dots.append(
393 axes_tr.plot(
394 [], [], 'o', ms=4, color=self.colors['lines'])[0])
396 axes_tr.my_lines = lines
397 axes_tr.my_dots = dots
398 axes_tr.my_stuff = []
400 new_axes.append(
401 (axes_01, axes_02, axes_12, axes_tr))
403 def resize_handler(*args):
404 self.layout(fig, new_axes)
406 if fig.my_disconnect:
407 fig.my_disconnect()
409 cid_resize = fig.canvas.mpl_connect('resize_event', resize_handler)
410 try:
411 cid_dpi = get_callbacks(fig).connect('dpi_changed', resize_handler)
412 except (AttributeError, ValueError):
413 pass # not available anymore in MPL >= 3.8
415 def disconnect():
416 fig.canvas.mpl_disconnect(cid_resize)
417 fig.callbacks.disconnect(cid_dpi)
419 fig.my_disconnect = disconnect
421 self.axes = new_axes
423 def get_trace(self, traces, pattern):
424 trs = [tr for tr in traces if util.match_nslc([pattern], tr.nslc_id)]
425 if len(trs) > 1:
426 self.fail('Multiple traces matching pattern %s' % pattern)
427 elif len(trs) == 0:
428 return None
429 else:
430 return trs[0]
432 def get_vector_abs_max(self, traces):
433 tr_abs = None
434 for tr in traces:
435 if tr is not None:
436 tr = tr.copy()
437 tr.ydata **= 2
438 if tr_abs is None:
439 tr_abs = tr
440 else:
441 tr_abs.add(tr)
443 tr_abs.set_ydata(num.sqrt(tr_abs.ydata))
444 return num.max(tr_abs.ydata)
446 def set_labels(
447 self, istation, nstations, tref,
448 axes_01, axes_02, axes_12, axes_tr, chas, chas_rot):
450 if self.azimuth != 0 or self.rotate_to_rt:
451 rcolor = self.colors['rot']
452 else:
453 rcolor = self.colors['fg']
454 chas_rot = chas
456 axes_01.set_ylabel(chas[1])
457 axes_02.set_ylabel(chas[2])
458 axes_12.set_ylabel(chas[2])
460 axes_01.set_xlabel(chas[0])
461 axes_02.set_xlabel(chas_rot[0])
462 axes_12.set_xlabel(chas_rot[1])
463 axes_02.get_xaxis().label.set_color(rcolor)
464 axes_12.get_xaxis().label.set_color(rcolor)
465 axes_tr.set_xlabel(
466 'Time [s] relative to %s' % util.time_to_str(tref))
468 axes_tr.set_yticks([0., 1., 2])
469 axes_tr.set_yticklabels([chas_rot[0], chas_rot[1], chas[2]])
470 for tlab in axes_tr.get_yticklabels()[:2]:
471 tlab.set_color(rcolor)
473 def pca(self, trs):
475 if any(tr is None for tr in trs):
476 raise PCAError('Missing component')
478 nss = [tr.data_len() for tr in trs]
479 if not all(ns == nss[0] for ns in nss):
480 raise PCAError('Traces have different lengths.')
482 ns = nss[0]
484 if ns < 3:
485 raise PCAError('Traces too short.')
487 data = num.zeros((len(trs), ns))
488 for itr, tr in enumerate(trs):
489 data[itr, :] = tr.ydata
491 cov = num.cov(data)
493 evals, evecs = num.linalg.eigh(cov)
495 pc = evecs[:, -1]
497 eh = num.sqrt(pc[1]**2 + pc[0]**2)
499 if len(trs) > 2:
500 incidence = r2d * num.arctan2(eh, abs(pc[2]))
501 else:
502 incidence = 90.
504 azimuth = r2d*num.arctan2(pc[1], pc[0])
505 azimuth = (-azimuth) % 180. - 90.
507 return cov, evals, evecs, azimuth, incidence
509 def draw_cov_ellipse(self, evals, evecs, color, alpha=1.0):
510 evals = num.sqrt(evals)
511 ell = patches.Ellipse(
512 xy=(0.0, 0.0),
513 width=evals[0] * 2.,
514 height=evals[1] * 2.,
515 angle=r2d*num.arctan2(evecs[1, 0], evecs[0, 0]),
516 zorder=-10,
517 fc=color,
518 ec=darken(color),
519 alpha=alpha)
521 return ell
523 def draw(self, groups, nsl_to_tspan, tpad, nsl_to_bazi):
525 def count(d, k):
526 if k not in d:
527 d[k] = 0
529 d[k] += 1
531 nsli_n = {}
532 for k in sorted(groups):
533 count(nsli_n, k[:4])
535 nsli_i = {}
536 for igroup, k in enumerate(sorted(groups)):
538 tmin, tmax = nsl_to_tspan[k]
539 nsl = k[:3]
540 nsli = k[:4]
541 count(nsli_i, nsli)
543 bazimuth = self.azimuth
545 if self.rotate_to_rt:
546 if nsl not in nsl_to_bazi:
547 self.fail(
548 'Cannot rotate to RT.\n\nActive event must '
549 'available (select event and press "e"). Station '
550 'coordinates must be available.')
552 bazimuth += nsl_to_bazi[nsl] + 90.
554 for axes in self.axes[igroup]:
555 while axes.my_stuff:
556 stuff = axes.my_stuff.pop()
557 stuff.remove()
559 axes_01, axes_02, axes_12, axes_tr = self.axes[igroup]
561 axes_01.set_title(
562 '.'.join(nsli)
563 + ('' if nsli_n[nsli] == 1 else ' (%i)' % nsli_i[nsli]))
565 trs_all = groups[k]
567 if len(trs_all) == 5:
568 trs_orig = trs_all[:3]
569 trs_rot = trs_all[3:] + trs_all[2:3]
570 elif len(trs_all) == 4:
571 trs_orig = trs_all[:2] + [None]
572 trs_rot = trs_all[2:] + [None]
573 else:
574 raise self.fail(
575 'Unexpectedly got other that 4 or 6 traces: %i'
576 % len(trs_all))
578 trs_orig_chopped = [
579 (tr.chop(tmin, tmax, inplace=False) if tr else None)
580 for tr in trs_orig]
582 trs_rot_chopped = [
583 (tr.chop(tmin, tmax, inplace=False) if tr else None)
584 for tr in trs_rot]
586 chas = [tr.channel[-1:] for tr in trs_orig_chopped]
587 chas_rot = [tr.channel[-2:] for tr in trs_rot_chopped]
589 if self.fix_scale and nsl in self.nsl_to_amax:
590 amax = self.nsl_to_amax[nsl]
591 else:
592 amax = self.get_vector_abs_max(trs_orig_chopped)
593 self.nsl_to_amax[nsl] = amax
595 if nsl in nsl_to_bazi:
596 event_bazi = nsl_to_bazi[nsl]
597 l1, = axes_01.plot(
598 [0., amax*num.cos((90. - event_bazi)*d2r)],
599 [0., amax*num.sin((90. - event_bazi)*d2r)],
600 '--',
601 lw=3,
602 zorder=-20,
603 color=self.colors['backazimuth'])
605 a1 = axes_01.annotate(
606 '%.0f\u00b0' % event_bazi,
607 (0., 1.),
608 (self.font_size, -self.font_size),
609 xycoords='axes fraction',
610 textcoords='offset points',
611 color=self.colors['backazimuth'],
612 ha='left',
613 va='top')
615 axes_01.my_stuff.extend([l1, a1])
617 for ix, iy, axes, trs in [
618 (0, 1, axes_01, trs_orig_chopped),
619 (0, 2, axes_02, trs_rot_chopped),
620 (1, 2, axes_12, trs_rot_chopped)]:
622 axes.set_xlim(-amax*1.05, amax*1.05)
623 axes.set_ylim(-amax*1.05, amax*1.05)
625 if not (trs[ix] and trs[iy]):
626 continue
628 x = trs[ix].get_ydata()
629 y = trs[iy].get_ydata()
631 axes.my_line.set_data(x, y)
632 ipos = int(round(self.dot_position * (x.size-1)))
633 axes.my_dot.set_data(x[ipos], y[ipos])
635 tref = tmin
636 for itr, (tr, tr_chopped) in enumerate(zip(
637 trs_rot, trs_rot_chopped)):
639 if tr is None or tr_chopped is None:
640 axes_tr.my_lines[itr].set_data([], [])
641 axes_tr.my_dots[itr].set_data([], [])
643 else:
644 y = tr.get_ydata() / (2.*amax) + itr
645 t = tr.get_xdata()
646 t = t - tref
648 ipos = int(round(
649 self.dot_position * (tr_chopped.data_len()-1)))
651 yp = tr_chopped.ydata[ipos] / (2.*amax) + itr
652 tp = tr_chopped.tmin - tref + tr_chopped.deltat*ipos
654 axes_tr.my_lines[itr].set_data(t, y)
655 axes_tr.my_dots[itr].set_data(tp, yp)
657 if self.azimuth != 0.0 or self.rotate_to_rt:
658 color = self.colors['rot']
660 xn = num.sin(bazimuth*d2r)
661 yn = num.cos(bazimuth*d2r)
662 xe = num.sin(bazimuth*d2r + 0.5*num.pi)
663 ye = num.cos(bazimuth*d2r + 0.5*num.pi)
665 l1, = axes_01.plot(
666 [0., amax*xn],
667 [0., amax*yn],
668 color=color)
670 font_size = self.font_size
672 a1 = axes_01.annotate(
673 chas_rot[1],
674 xy=(amax*xn, amax*yn),
675 xycoords='data',
676 xytext=(-font_size*(xe+.5*xn), -font_size*(ye+.5*yn)),
677 textcoords='offset points',
678 va='center',
679 ha='center',
680 color=color)
682 l2, = axes_01.plot(
683 [0., amax*xe],
684 [0., amax*ye],
685 color=color)
687 a2 = axes_01.annotate(
688 chas_rot[0],
689 xy=(amax*xe, amax*ye),
690 xycoords='data',
691 xytext=(-font_size*(xn+.5*xe), -font_size*(yn+.5*ye)),
692 textcoords='offset points',
693 va='center',
694 ha='center',
695 color=color)
697 aa = axes_01.annotate(
698 '%.0f\u00b0' % bazimuth,
699 (1., 1.),
700 (-self.font_size, -self.font_size),
701 xycoords='axes fraction',
702 textcoords='offset points',
703 color=color,
704 ha='right',
705 va='top')
707 axes_01.my_stuff.extend([l1, a1, l2, a2, aa])
709 axes_tr.my_stuff.append(axes_tr.axvspan(
710 tmin - tref, tmax - tref, color=self.colors['time_window']))
712 axes_tr.set_ylim(-1, 3)
713 axes_tr.set_xlim(tmin - tref - tpad, tmax - tref + tpad)
715 self.set_labels(
716 igroup, len(groups), tref, *self.axes[igroup], chas, chas_rot)
718 if self.show_eigensystem_2d:
720 for (ix, iy, axes, trs) in [
721 (0, 1, axes_01, trs_orig_chopped),
722 (0, 2, axes_02, trs_rot_chopped),
723 (1, 2, axes_12, trs_rot_chopped)]:
725 try:
726 cov, evals, evecs, azimuth, _ = self.pca(
727 [trs[ix], trs[iy]])
729 axes.my_stuff.append(
730 axes.annotate(
731 '%.0f\u00b0' % azimuth,
732 (0., 0.),
733 (self.font_size, self.font_size),
734 xycoords='axes fraction',
735 textcoords='offset points',
736 color=self.colors['pca_2d'],
737 alpha=0.5))
739 ell = self.draw_cov_ellipse(
740 evals[:2], evecs[:2, :2],
741 color=self.colors['pca_2d'], alpha=0.5)
743 axes.add_artist(ell)
744 axes.my_stuff.append(ell)
746 l1, = axes.plot(
747 [-amax*evecs[0, -1], amax*evecs[0, -1]],
748 [-amax*evecs[1, -1], amax*evecs[1, -1]],
749 color=self.colors['pca_2d'], alpha=0.8)
751 l2, = axes.plot(
752 [-amax*evecs[0, -2], amax*evecs[0, -2]],
753 [-amax*evecs[1, -2], amax*evecs[1, -2]],
754 color=self.colors['pca_2d'], alpha=0.2)
756 axes.my_stuff.extend([l1, l2])
758 except PCAError as e:
759 logger.warning('PCA failed: %s' % e)
761 if self.show_eigensystem_3d:
762 try:
763 cov, evals, evecs, azimuth, incidence = self.pca(
764 trs_orig_chopped)
765 cosa = num.cos(bazimuth*d2r)
766 sina = num.sin(bazimuth*d2r)
767 rot = num.array(
768 [[cosa, -sina, 0.0],
769 [sina, cosa, 0.0],
770 [0.0, 0.0, 1.0]], dtype=float)
772 evecs_rot = num.dot(rot, evecs)
774 for (ix, iy, axes, evecs_) in [
775 (0, 1, axes_01, evecs),
776 (0, 2, axes_02, evecs_rot),
777 (1, 2, axes_12, evecs_rot)]:
779 for (ie, alpha) in [
780 (-1, 0.5),
781 (-2, 0.5),
782 (-3, 0.5)]:
784 for vlen, lw in [
785 (num.sqrt(evals[ie]), 3),
786 (amax, 1)]:
788 lv, = axes.plot(
789 [-vlen*evecs_[ix, ie],
790 vlen*evecs_[ix, ie]],
791 [-vlen*evecs_[iy, ie],
792 vlen*evecs_[iy, ie]],
793 lw=lw,
794 color=self.colors['pca_3d'], alpha=alpha)
796 axes.my_stuff.extend([lv])
798 axes_01.my_stuff.append(
799 axes_01.annotate(
800 '%.0f\u00b0' % azimuth,
801 (1., 0.),
802 (-self.font_size, self.font_size),
803 xycoords='axes fraction',
804 textcoords='offset points',
805 color=self.colors['pca_3d'],
806 alpha=0.8,
807 ha='right'))
809 axes_02.my_stuff.append(
810 axes_02.annotate(
811 '%.0f\u00b0' % incidence,
812 (1., 0.),
813 (-self.font_size, self.font_size),
814 xycoords='axes fraction',
815 textcoords='offset points',
816 color=self.colors['pca_3d'],
817 alpha=0.8,
818 ha='right'))
820 except PCAError as e:
821 logger.warning('PCA failed: %s' % e)
823 def get_bazis(self):
824 event = self.get_viewer().get_active_event()
825 if not event:
826 return {}
828 nsl_to_bazi = dict(
829 (station.nsl(), event.azibazi_to(station)[1])
830 for station in self.get_stations())
832 return nsl_to_bazi
834 def layout(self, fig, axes):
836 # Do not access self in here. Called from resize in finished plots.
838 def get_pixels_factor(fig):
839 try:
840 r = fig.canvas.get_renderer()
841 return 1.0 / r.points_to_pixels(1.0)
842 except AttributeError:
843 return 1.0
845 def rect_to_figure_coords(rect):
846 l, b, w, h = rect # noqa
847 return (l / width, b / height, w / width, h / height) # noqa
849 ny = len(axes)
850 if ny == 0:
851 raise LayoutError('No axes given.')
853 nx = len(axes[0])
855 width, height = fig.canvas.get_width_height()
856 pixels = get_pixels_factor(fig)
858 margin_left = margin_right = 4. * self.font_size / pixels
859 margin_top = margin_bottom = 5. * self.font_size / pixels
861 spacing_width = 3. * self.font_size / pixels
862 spacing_height = 5. * self.font_size / pixels
864 axes_height_avail = height - (ny - 1) * spacing_height \
865 - margin_top - margin_bottom
867 a_min = 5.
869 if axes_height_avail < a_min * ny:
870 axes_height_avail = a_min * ny
871 logger.warning('Not enough space vertically.')
873 axes_width_avail = width - (nx - 1) * spacing_width \
874 - margin_left - margin_right
876 if axes_width_avail < a_min * (nx + 2):
877 axes_width_avail = a_min * (nx + 2)
878 logger.warning('Not enough space horizontally.')
880 a_height = axes_height_avail / ny
881 a_width = axes_width_avail / (nx + 2)
883 a = max(min(a_height, a_width), a_min)
885 pad_height = (a_height - a) * ny
886 pad_width = (a_width - a) * (nx + 2)
888 for iy in range(ny):
889 y = height - 0.5 * pad_height - margin_top \
890 - (iy + 1) * a - iy * spacing_height
891 h = a
892 for ix in range(nx):
893 x = margin_right + 0.5 * pad_width + ix * (a + spacing_width)
894 w = a if ix != (nx - 1) else a * 3.0
895 axes[iy][ix].set_position(
896 rect_to_figure_coords((x, y, w, h)), which='both')
898 def call(self):
900 self.cleanup()
902 if self.rotate_to_rt != self.last_rotate_to_rt:
903 # reset delta azimuth to avoid confusion
905 self.set_parameter('azimuth', 0.0)
907 self.last_rotate_to_rt = self.rotate_to_rt
909 selected = self.get_selected()
911 nsl_to_bazi = self.get_bazis()
913 tpad = self.t_length
914 groups = self.get_traces(selected, nsl_to_bazi, tpad)
916 if not groups:
917 self.fail(
918 'No matching traces. Check time and channel settings. Traces '
919 'may not contain gaps within the extracted time window and in '
920 'the padding areas left and right. Traces are extracted with '
921 'additional padding of 3 x filter corner period to eliminate '
922 'artifacts.')
924 selection_markers = self.make_selection_markers(selected, groups)
925 self.add_markers(selection_markers)
927 fig = self.get_figure()
929 figure_key = (len(groups), self.iframe)
931 if not self.figure_key or self.figure_key != figure_key:
932 self.setup_figure(fig, len(groups))
933 self.figure_key = figure_key
935 self.draw(groups, selected, tpad, nsl_to_bazi)
937 self.layout(fig, self.axes)
939 self.fframe.draw()
940 tabs = self.fframe.parent().parent()
941 # bring plot to front if we are not looking at the markers
942 from pyrocko.gui.pile_viewer import PileViewer
943 if not isinstance(tabs.currentWidget(), PileViewer):
944 tabs.setCurrentWidget(self.fframe)
947def __snufflings__():
948 return [Polarization()]