Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/gui/snuffler/snufflings/polarization.py: 13%

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 def help(self):

44 return '''

45<html>

46<head>

47<style type="text/css">

48 body { margin-left:10px };

49</style>

50</head>

51<body>

53<h1>

54Investigate patterns of ground motion during the passage of seismic waves.

55</h1>

57

58This Snuffling can be used to analyze and visualize the polarization of seismic

59waves from 3-component seismic recordings or to check component orientations of

60a seismic sensors when used on signals with known directional properties. The

61spatial pattern of ground motion is shown in horizontal and vertical

62projections. Principal component analysis and rotation to radial/transverse

63components are available as tools. Time window and filter settings can be

64adjusted interactively.

65

67<h2>

68Usage

69</h2>

71

72Select one or more normal/phase markers as anchor points for the extraction and

73press Run. Multiple stations can be selected for direct

74comparison. Channels matching the pattern-triplet given in the

75Channels setting are selected for extraction of a time window

76around the anchor point. It is assumed that these three channels correspond to

77sensor components in the order east, north, vertical (upward), even if they are

78named differently.

79

81

82The time window can be adjusted with the Length and

83Offset parameters. Extracted waveforms are filtered according

84the Highpass and Lowpass parameters

85(Butterworth 4th order) and demeaned.

86

88

89To rotate the horizontal components around a vertical axis, use the

90\u0394 Azimuth setting. When station coordinates and an

91"active" event are available, the horizontal components can be rotated to

92radial (away from event) and transverse (leftwards) orientations using the

93computed event back-azimuth (dashed gray line).

94

96

97If Show 2D Eigensystems is selected, a principal component

98analysis is performed in each of the shown projects. Red lines are shown to

99depict the eigenvectors and the eigenvalues are visualized using ellipse

100symbols. If Show 3D Eigensystems is selected, a principal

101component analysis is performed using all three (non-rotated) components and

102the resulting eigenvectors are depicted with purple lines.

103

104

105

106By default the scaling is automatically adjusted to the maximum value (the

107maximum vector length of the three-component signal within the selected time

108window). The currently used scaling factor can be frozen by checking

109Fix Scale.

110

111</body>

112</html>

113'''

114

115 def setup(self):

116 self.set_name('Polarization')

117

118 self.add_parameter(

119 Param('Length [s]', 't_length', 1., 0.001, 1000.))

120

121 self.add_parameter(

122 Param('Offset (relative)', 'ft_offset', 0., -5., 5.))

123

124 self.add_parameter(

125 Param(u'\u0394 Azimuth', 'azimuth', 0., -180., 180.))

126

127 self.add_parameter(

128 Param(

129 'Highpass [Hz]', 'highpass', None, 0.001, 1000.,

130 low_is_none=True))

131

132 self.add_parameter(

133 Param(

134 'Lowpass [Hz]', 'lowpass', None, 0.001, 1000.,

135 high_is_none=True))

136

137 self.add_parameter(

138 Param('Dot Position', 'dot_position', 0., 0., 1.))

139

140 self.add_parameter(

141 Switch(

142 'Rotate to RT',

143 'rotate_to_rt',

144 False))

145

146 self.add_parameter(

147 Switch(

148 'Show 2D Eigensystems',

149 'show_eigensystem_2d',

150 False))

151

152 self.add_parameter(

153 Switch(

154 'Show 3D Eigensystem',

155 'show_eigensystem_3d',

156 False))

157

158 self.add_parameter(

159 Switch(

160 'Fix Scale',

161 'fix_scale',

162 False))

163

164 def new_figure():

165 self.setup_figure_frame()

166 self.call()

167

168 self.add_trigger('New Figure', new_figure)

169

170 self.fframe = None

171 self.iframe = 0

172 self.figure_key = None

173 self.nsl_to_amax = {}

174 self.last_rotate_to_rt = False

175

176 self.font_size = float(matplotlib.rcParams['font.size'])

177

178 self.colors = {

179 'fg': matplotlib.rcParams['axes.edgecolor'],

180 'lines': 'black',

181 'rot': plot.mpl_color('skyblue1'),

182 'pca_2d': plot.mpl_color('scarletred1'),

183 'pca_3d': plot.mpl_color('plum2'),

184 'backazimuth': plot.mpl_color('aluminium3'),

185 'time_window': plot.mpl_color('aluminium2')}

186

187 def panel_visibility_changed(self, visible):

188 viewer = self.get_viewer()

189 if visible:

190 viewer.pile_has_changed_signal.connect(self.adjust_controls)

191 self.adjust_controls()

192

193 else:

194 viewer.pile_has_changed_signal.disconnect(self.adjust_controls)

195

196 def adjust_controls(self):

197 viewer = self.get_viewer()

198 dtmin, dtmax = viewer.content_deltat_range()

199 maxfreq = 0.5/dtmin

200 minfreq = (0.5/dtmax)*0.001

201 self.set_parameter_range('lowpass', minfreq, maxfreq)

202 self.set_parameter_range('highpass', minfreq, maxfreq)

203

204 def get_selected(self):

205 markers = [

206 marker for marker in self.get_selected_markers()

207 if not isinstance(marker, EventMarker)]

208

209 if not markers:

210 self.fail(

211 'No selected markers.\n\nCreate and select markers at points '

212 'of interest. Normal and phase markers are accepted.')

213

214 d = {}

215 for marker in markers:

216 tspan = self.transform_time_span(marker.tmin, marker.tmax)

217 for nslc in marker.nslc_ids:

218 k = nslc[:3] + (nslc[3][:-1], marker.tmin)

219 d[k] = tspan

220

221 return d

222

223 def transform_time_span(self, tmin, tmax):

224 tmin = tmin + self.t_length * self.ft_offset

225 tmax = tmin + self.t_length

226 return (tmin, tmax)

227

228 def make_selection_markers(self, selected, groups):

229 selection_markers = []

230 for k in sorted(selected):

231 tmin, tmax = selected[k]

232 patterns = [tr.nslc_id for tr in groups[k]]

233 marker = Marker(

234 nslc_ids=patterns, tmin=tmin, tmax=tmax, kind=2)

235 selection_markers.append(marker)

236

237 return selection_markers

238

239 def sorted_by_component_scheme(self, group):

240 if len(group) not in (2, 3):

241 self.fail('Need 2 or 3 component recordings. Got %i.' % len(group))

242

243 comps = []

244 channels = []

245 for tr in group:

246 comps.append(tr.channel[-1])

247 channels.append(tr.channel)

248

249 # special case for channel naming convention of Cube dataloggers

250 for scheme in [['p2', 'p1', 'p0'], ['p2', 'p1']]:

251

252 if sorted(channels) == sorted(scheme):

253 return [

254 get_tr(group, lambda tr: tr.channel == ch)

255 for ch in scheme]

256

257 for scheme in [

258 ['E', 'N', 'Z'], ['1', '2', 'Z'], ['1', '2', '3'],

259 ['0', '1', '2'], ['R', 'T', 'Z'],

260 ['E', 'N'], ['1', '2'], ['0', '1']]:

261

262 if sorted(comps) == sorted(scheme):

263 return [

264 get_tr(group, lambda tr: tr.channel[-1] == comp)

265 for comp in scheme]

266

267 self.fail('Cannot handle component group: %s' % ', '.join(channels))

268

269 def get_traces(self, selected, nsl_to_bazi, tpad):

270 if self.highpass is not None:

271 tpad_filter = 3.0 / self.highpass

272 elif self.lowpass is not None:

273 tpad_filter = 3.0 / self.lowpass

274 else:

275 tpad_filter = 0.0

276

277 # prevent getting insanely long cutouts if e.g. if highpass is still

278 # very low, e.g. while moving the slider

279 tpad_filter = min(tpad_filter, 5.0 * self.t_length)

280

281 d = {}

282 for k in sorted(selected):

283 nsl = k[0:3]

284 tmin, tmax = selected[k]

285

286 bazimuth = self.azimuth

287

288 if self.rotate_to_rt:

289 if nsl not in nsl_to_bazi:

290 self.fail(

291 'Cannot rotate to RT.\n\nStation coordinates must be '

292 'available and an event must be marked as the '

293 '"active" event (select event and press "e").')

294

295 bazimuth += nsl_to_bazi[nsl] + 90.

296

297 pattern = '.'.join(nsl + (k[3] + '?',))

298

299 group = []

300 for batch in self.get_pile().chopper(

301 tmin=tmin, tmax=tmax, tpad=tpad + tpad_filter,

302 want_incomplete=False,

303 style='batch',

304 trace_selector=lambda tr: util.match_nslc(

305 pattern, tr.nslc_id)):

306

307 for tr in batch.traces:

308 tr = tr.copy()

309 if self.lowpass is not None \

310 and self.lowpass < 0.5/tr.deltat:

311 tr.lowpass(4, self.lowpass)

312

313 if self.highpass is not None \

314 and self.highpass < 0.5/tr.deltat:

315

316 tr.highpass(4, self.highpass)

317

318 try:

319 tr.chop(tmin - tpad, tmax + tpad)

320 tr_chop = tr.chop(tmin, tmax, inplace=False)

321 except trace.NoData:

322 self.fail('No data. Time length too short?')

323

324 y = tr.get_ydata()

325 tr.set_ydata(y - num.mean(tr_chop.get_ydata()))

326

327 group.append(tr)

328

329 group = self.sorted_by_component_scheme(group)

330

331 tr_e = group[0]

332 tr_n = group[1]

333

334 if tr_e is not None and tr_n is not None:

335 cha_e = tr_e.channel

336 cha_n = tr_n.channel

337 if bazimuth != 0.0:

338 trs_rot = trace.rotate(

339 [tr_n, tr_e], bazimuth,

340 in_channels=[cha_n, cha_e],

341 out_channels=[ch+"'" for ch in [cha_n, cha_e]])

342 else:

343 trs_rot = [tr.copy() for tr in [tr_n, tr_e]]

344 for tr in trs_rot:

345 tr.set_codes(channel=tr.channel + "'")

346

347 if len(trs_rot) == 2:

348 group.extend(

349 get_tr(trs_rot, lambda tr: tr.channel == ch+"'")

350 for ch in [cha_e, cha_n])

351 else:

352 self.fail(

353 'Cannot rotate traces. '

354 'Are the samples properly aligned?')

355

356 d[k] = group

357

358 return d

359

360 def setup_figure_frame(self):

361 self.iframe += 1

362 self.fframe = self.figure_frame(

363 'Polarization (%i)' % self.iframe)

364

365 self.fframe.gcf().my_disconnect = None

366

367 def get_figure(self):

368 if not self.fframe or self.fframe.closed:

369 self.setup_figure_frame()

370

371 return self.fframe.gcf()

372

373 def setup_figure(self, fig, nstations):

374 fig.clf()

375 new_axes = []

376

377 def iwrap(iy, ix):

378 return (ix + iy * 4) + 1

379

380 for istation in range(nstations):

381 axes_01 = fig.add_subplot(

382 nstations, 4, iwrap(istation, 0), aspect=1.0)

383 axes_02 = fig.add_subplot(

384 nstations, 4, iwrap(istation, 1), aspect=1.0)

385 axes_12 = fig.add_subplot(

386 nstations, 4, iwrap(istation, 2), aspect=1.0)

387 axes_tr = fig.add_subplot(

388 nstations, 4, iwrap(istation, 3))

389

390 for axes in (axes_01, axes_02, axes_12, axes_tr):

391 axes.my_stuff = []

392

393 axes.my_line, = axes.plot(

394 [], [], color=self.colors['lines'], lw=1.0)

395 axes.my_dot, = axes.plot(

396 [], [], 'o', ms=4, color=self.colors['lines'])

397

398 for axes in (axes_01, axes_02, axes_12):

399 axes.get_xaxis().set_tick_params(

400 labelbottom=False, bottom=False)

401 axes.get_yaxis().set_tick_params(

402 labelleft=False, left=False)

403

404 axes_tr.get_yaxis().set_tick_params(

405 left=False, labelleft=True, length=2.0)

406

407 # if istation != nstations - 1:

408 # axes_tr.get_xaxis().set_tick_params(

409 # bottom=False, labelbottom=False)

410

411 lines = []

412 dots = []

413 for _ in range(3):

414 lines.append(

415 axes_tr.plot(

416 [], [], color=self.colors['lines'], lw=1.0)[0])

417 dots.append(

418 axes_tr.plot(

419 [], [], 'o', ms=4, color=self.colors['lines'])[0])

420

421 axes_tr.my_lines = lines

422 axes_tr.my_dots = dots

423 axes_tr.my_stuff = []

424

425 new_axes.append(

426 (axes_01, axes_02, axes_12, axes_tr))

427

428 def resize_handler(*args):

429 self.layout(fig, new_axes)

430

431 if fig.my_disconnect:

432 fig.my_disconnect()

433

434 cid_resize = fig.canvas.mpl_connect('resize_event', resize_handler)

435 try:

436 cid_dpi = get_callbacks(fig).connect('dpi_changed', resize_handler)

437 except (AttributeError, ValueError):

438 pass # not available anymore in MPL >= 3.8

439

440 def disconnect():

441 fig.canvas.mpl_disconnect(cid_resize)

442 fig.callbacks.disconnect(cid_dpi)

443

444 fig.my_disconnect = disconnect

445

446 self.axes = new_axes

447

448 def get_trace(self, traces, pattern):

449 trs = [tr for tr in traces if util.match_nslc([pattern], tr.nslc_id)]

450 if len(trs) > 1:

451 self.fail('Multiple traces matching pattern %s' % pattern)

452 elif len(trs) == 0:

453 return None

454 else:

455 return trs[0]

456

457 def get_vector_abs_max(self, traces):

458 tr_abs = None

459 for tr in traces:

460 if tr is not None:

461 tr = tr.copy()

462 tr.ydata **= 2

463 if tr_abs is None:

464 tr_abs = tr

465 else:

466 tr_abs.add(tr)

467

468 tr_abs.set_ydata(num.sqrt(tr_abs.ydata))

469 return num.max(tr_abs.ydata)

470

471 def set_labels(

472 self, istation, nstations, tref,

473 axes_01, axes_02, axes_12, axes_tr, chas, chas_rot):

474

475 if self.azimuth != 0 or self.rotate_to_rt:

476 rcolor = self.colors['rot']

477 else:

478 rcolor = self.colors['fg']

479 chas_rot = chas

480

481 axes_01.set_ylabel(chas[1])

482 axes_02.set_ylabel(chas[2])

483 axes_12.set_ylabel(chas[2])

484

485 axes_01.set_xlabel(chas[0])

486 axes_02.set_xlabel(chas_rot[0])

487 axes_12.set_xlabel(chas_rot[1])

488 axes_02.get_xaxis().label.set_color(rcolor)

489 axes_12.get_xaxis().label.set_color(rcolor)

490 axes_tr.set_xlabel(

491 'Time [s] relative to %s' % util.time_to_str(tref))

492

493 axes_tr.set_yticks([0., 1., 2])

494 axes_tr.set_yticklabels([chas_rot[0], chas_rot[1], chas[2]])

495 for tlab in axes_tr.get_yticklabels()[:2]:

496 tlab.set_color(rcolor)

497

498 def pca(self, trs):

499

500 if any(tr is None for tr in trs):

501 raise PCAError('Missing component')

502

503 nss = [tr.data_len() for tr in trs]

504 if not all(ns == nss[0] for ns in nss):

505 raise PCAError('Traces have different lengths.')

506

507 ns = nss[0]

508

509 if ns < 3:

510 raise PCAError('Traces too short.')

511

512 data = num.zeros((len(trs), ns))

513 for itr, tr in enumerate(trs):

514 data[itr, :] = tr.ydata

515

516 cov = num.cov(data)

517

518 evals, evecs = num.linalg.eigh(cov)

519

520 pc = evecs[:, -1]

521

522 eh = num.sqrt(pc[1]**2 + pc[0]**2)

523

524 if len(trs) > 2:

525 incidence = r2d * num.arctan2(eh, abs(pc[2]))

526 else:

527 incidence = 90.

528

529 azimuth = r2d*num.arctan2(pc[1], pc[0])

530 azimuth = (-azimuth) % 180. - 90.

531

532 return cov, evals, evecs, azimuth, incidence

533

534 def draw_cov_ellipse(self, evals, evecs, color, alpha=1.0):

535 evals = num.sqrt(evals)

536 ell = patches.Ellipse(

537 xy=(0.0, 0.0),

538 width=evals[0] * 2.,

539 height=evals[1] * 2.,

540 angle=r2d*num.arctan2(evecs[1, 0], evecs[0, 0]),

541 zorder=-10,

542 fc=color,

543 ec=darken(color),

544 alpha=alpha)

545

546 return ell

547

548 def draw(self, groups, nsl_to_tspan, tpad, nsl_to_bazi):

549

550 def count(d, k):

551 if k not in d:

552 d[k] = 0

553

554 d[k] += 1

555

556 nsli_n = {}

557 for k in sorted(groups):

558 count(nsli_n, k[:4])

559

560 nsli_i = {}

561 for igroup, k in enumerate(sorted(groups)):

562

563 tmin, tmax = nsl_to_tspan[k]

564 nsl = k[:3]

565 nsli = k[:4]

566 count(nsli_i, nsli)

567

568 bazimuth = self.azimuth

569

570 if self.rotate_to_rt:

571 if nsl not in nsl_to_bazi:

572 self.fail(

573 'Cannot rotate to RT.\n\nActive event must '

574 'available (select event and press "e"). Station '

575 'coordinates must be available.')

576

577 bazimuth += nsl_to_bazi[nsl] + 90.

578

579 for axes in self.axes[igroup]:

580 while axes.my_stuff:

581 stuff = axes.my_stuff.pop()

582 stuff.remove()

583

584 axes_01, axes_02, axes_12, axes_tr = self.axes[igroup]

585

586 axes_01.set_title(

587 '.'.join(nsli)

588 + ('' if nsli_n[nsli] == 1 else ' (%i)' % nsli_i[nsli]))

589

590 trs_all = groups[k]

591

592 if len(trs_all) == 5:

593 trs_orig = trs_all[:3]

594 trs_rot = trs_all[3:] + trs_all[2:3]

595 elif len(trs_all) == 4:

596 trs_orig = trs_all[:2] + [None]

597 trs_rot = trs_all[2:] + [None]

598 else:

599 raise self.fail(

600 'Unexpectedly got other that 4 or 6 traces: %i'

601 % len(trs_all))

602

603 trs_orig_chopped = [

604 (tr.chop(tmin, tmax, inplace=False) if tr else None)

605 for tr in trs_orig]

606

607 trs_rot_chopped = [

608 (tr.chop(tmin, tmax, inplace=False) if tr else None)

609 for tr in trs_rot]

610

611 chas = [tr.channel[-1:] for tr in trs_orig_chopped]

612 chas_rot = [tr.channel[-2:] for tr in trs_rot_chopped]

613

614 if self.fix_scale and nsl in self.nsl_to_amax:

615 amax = self.nsl_to_amax[nsl]

616 else:

617 amax = self.get_vector_abs_max(trs_orig_chopped)

618 self.nsl_to_amax[nsl] = amax

619

620 if nsl in nsl_to_bazi:

621 event_bazi = nsl_to_bazi[nsl]

622 l1, = axes_01.plot(

623 [0., amax*num.cos((90. - event_bazi)*d2r)],

624 [0., amax*num.sin((90. - event_bazi)*d2r)],

625 '--',

626 lw=3,

627 zorder=-20,

628 color=self.colors['backazimuth'])

629

630 a1 = axes_01.annotate(

631 '%.0f\u00b0' % event_bazi,

632 (0., 1.),

633 (self.font_size, -self.font_size),

634 xycoords='axes fraction',

635 textcoords='offset points',

636 color=self.colors['backazimuth'],

637 ha='left',

638 va='top')

639

640 axes_01.my_stuff.extend([l1, a1])

641

642 for ix, iy, axes, trs in [

643 (0, 1, axes_01, trs_orig_chopped),

644 (0, 2, axes_02, trs_rot_chopped),

645 (1, 2, axes_12, trs_rot_chopped)]:

646

647 axes.set_xlim(-amax*1.05, amax*1.05)

648 axes.set_ylim(-amax*1.05, amax*1.05)

649

650 if not (trs[ix] and trs[iy]):

651 continue

652

653 x = trs[ix].get_ydata()

654 y = trs[iy].get_ydata()

655

656 axes.my_line.set_data(x, y)

657 ipos = int(round(self.dot_position * (x.size-1)))

658 axes.my_dot.set_data(x[ipos], y[ipos])

659

660 tref = tmin

661 for itr, (tr, tr_chopped) in enumerate(zip(

662 trs_rot, trs_rot_chopped)):

663

664 if tr is None or tr_chopped is None:

665 axes_tr.my_lines[itr].set_data([], [])

666 axes_tr.my_dots[itr].set_data([], [])

667

668 else:

669 y = tr.get_ydata() / (2.*amax) + itr

670 t = tr.get_xdata()

671 t = t - tref

672

673 ipos = int(round(

674 self.dot_position * (tr_chopped.data_len()-1)))

675

676 yp = tr_chopped.ydata[ipos] / (2.*amax) + itr

677 tp = tr_chopped.tmin - tref + tr_chopped.deltat*ipos

678

679 axes_tr.my_lines[itr].set_data(t, y)

680 axes_tr.my_dots[itr].set_data(tp, yp)

681

682 if self.azimuth != 0.0 or self.rotate_to_rt:

683 color = self.colors['rot']

684

685 xn = num.sin(bazimuth*d2r)

686 yn = num.cos(bazimuth*d2r)

687 xe = num.sin(bazimuth*d2r + 0.5*num.pi)

688 ye = num.cos(bazimuth*d2r + 0.5*num.pi)

689

690 l1, = axes_01.plot(

691 [0., amax*xn],

692 [0., amax*yn],

693 color=color)

694

695 font_size = self.font_size

696

697 a1 = axes_01.annotate(

698 chas_rot[1],

699 xy=(amax*xn, amax*yn),

700 xycoords='data',

701 xytext=(-font_size*(xe+.5*xn), -font_size*(ye+.5*yn)),

702 textcoords='offset points',

703 va='center',

704 ha='center',

705 color=color)

706

707 l2, = axes_01.plot(

708 [0., amax*xe],

709 [0., amax*ye],

710 color=color)

711

712 a2 = axes_01.annotate(

713 chas_rot[0],

714 xy=(amax*xe, amax*ye),

715 xycoords='data',

716 xytext=(-font_size*(xn+.5*xe), -font_size*(yn+.5*ye)),

717 textcoords='offset points',

718 va='center',

719 ha='center',

720 color=color)

721

722 aa = axes_01.annotate(

723 '%.0f\u00b0' % bazimuth,

724 (1., 1.),

725 (-self.font_size, -self.font_size),

726 xycoords='axes fraction',

727 textcoords='offset points',

728 color=color,

729 ha='right',

730 va='top')

731

732 axes_01.my_stuff.extend([l1, a1, l2, a2, aa])

733

734 axes_tr.my_stuff.append(axes_tr.axvspan(

735 tmin - tref, tmax - tref, color=self.colors['time_window']))

736

737 axes_tr.set_ylim(-1, 3)

738 axes_tr.set_xlim(tmin - tref - tpad, tmax - tref + tpad)

739

740 self.set_labels(

741 igroup, len(groups), tref, *self.axes[igroup], chas, chas_rot)

742

743 if self.show_eigensystem_2d:

744

745 for (ix, iy, axes, trs) in [

746 (0, 1, axes_01, trs_orig_chopped),

747 (0, 2, axes_02, trs_rot_chopped),

748 (1, 2, axes_12, trs_rot_chopped)]:

749

750 try:

751 cov, evals, evecs, azimuth, _ = self.pca(

752 [trs[ix], trs[iy]])

753

754 axes.my_stuff.append(

755 axes.annotate(

756 '%.0f\u00b0' % azimuth,

757 (0., 0.),

758 (self.font_size, self.font_size),

759 xycoords='axes fraction',

760 textcoords='offset points',

761 color=self.colors['pca_2d'],

762 alpha=0.5))

763

764 ell = self.draw_cov_ellipse(

765 evals[:2], evecs[:2, :2],

766 color=self.colors['pca_2d'], alpha=0.5)

767

768 axes.add_artist(ell)

769 axes.my_stuff.append(ell)

770

771 l1, = axes.plot(

772 [-amax*evecs[0, -1], amax*evecs[0, -1]],

773 [-amax*evecs[1, -1], amax*evecs[1, -1]],

774 color=self.colors['pca_2d'], alpha=0.8)

775

776 l2, = axes.plot(

777 [-amax*evecs[0, -2], amax*evecs[0, -2]],

778 [-amax*evecs[1, -2], amax*evecs[1, -2]],

779 color=self.colors['pca_2d'], alpha=0.2)

780

781 axes.my_stuff.extend([l1, l2])

782

783 except PCAError as e:

784 logger.warning('PCA failed: %s' % e)

785

786 if self.show_eigensystem_3d:

787 try:

788 cov, evals, evecs, azimuth, incidence = self.pca(

789 trs_orig_chopped)

790 cosa = num.cos(bazimuth*d2r)

791 sina = num.sin(bazimuth*d2r)

792 rot = num.array(

793 [[cosa, -sina, 0.0],

794 [sina, cosa, 0.0],

795 [0.0, 0.0, 1.0]], dtype=float)

796

797 evecs_rot = num.dot(rot, evecs)

798

799 for (ix, iy, axes, evecs_) in [

800 (0, 1, axes_01, evecs),

801 (0, 2, axes_02, evecs_rot),

802 (1, 2, axes_12, evecs_rot)]:

803

804 for (ie, alpha) in [

805 (-1, 0.5),

806 (-2, 0.5),

807 (-3, 0.5)]:

808

809 for vlen, lw in [

810 (num.sqrt(evals[ie]), 3),

811 (amax, 1)]:

812

813 lv, = axes.plot(

814 [-vlen*evecs_[ix, ie],

815 vlen*evecs_[ix, ie]],

816 [-vlen*evecs_[iy, ie],

817 vlen*evecs_[iy, ie]],

818 lw=lw,

819 color=self.colors['pca_3d'], alpha=alpha)

820

821 axes.my_stuff.extend([lv])

822

823 axes_01.my_stuff.append(

824 axes_01.annotate(

825 '%.0f\u00b0' % azimuth,

826 (1., 0.),

827 (-self.font_size, self.font_size),

828 xycoords='axes fraction',

829 textcoords='offset points',

830 color=self.colors['pca_3d'],

831 alpha=0.8,

832 ha='right'))

833

834 axes_02.my_stuff.append(

835 axes_02.annotate(

836 '%.0f\u00b0' % incidence,

837 (1., 0.),

838 (-self.font_size, self.font_size),

839 xycoords='axes fraction',

840 textcoords='offset points',

841 color=self.colors['pca_3d'],

842 alpha=0.8,

843 ha='right'))

844

845 except PCAError as e:

846 logger.warning('PCA failed: %s' % e)

847

848 def get_bazis(self):

849 event = self.get_viewer().get_active_event()

850 if not event:

851 return {}

852

853 nsl_to_bazi = dict(

854 (station.nsl(), event.azibazi_to(station)[1])

855 for station in self.get_stations())

856

857 return nsl_to_bazi

858

859 def layout(self, fig, axes):

860

861 # Do not access self in here. Called from resize in finished plots.

862

863 def get_pixels_factor(fig):

864 try:

865 r = fig.canvas.get_renderer()

866 return 1.0 / r.points_to_pixels(1.0)

867 except AttributeError:

868 return 1.0

869

870 def rect_to_figure_coords(rect):

871 l, b, w, h = rect # noqa

872 return (l / width, b / height, w / width, h / height) # noqa

873

874 ny = len(axes)

875 if ny == 0:

876 raise LayoutError('No axes given.')

877

878 nx = len(axes[0])

879

880 width, height = fig.canvas.get_width_height()

881 pixels = get_pixels_factor(fig)

882

883 margin_left = margin_right = 4. * self.font_size / pixels

884 margin_top = margin_bottom = 5. * self.font_size / pixels

885

886 spacing_width = 3. * self.font_size / pixels

887 spacing_height = 5. * self.font_size / pixels

888

889 axes_height_avail = height - (ny - 1) * spacing_height \

890 - margin_top - margin_bottom

891

892 a_min = 5.

893

894 if axes_height_avail < a_min * ny:

895 axes_height_avail = a_min * ny

896 logger.warning('Not enough space vertically.')

897

898 axes_width_avail = width - (nx - 1) * spacing_width \

899 - margin_left - margin_right

900

901 if axes_width_avail < a_min * (nx + 2):

902 axes_width_avail = a_min * (nx + 2)

903 logger.warning('Not enough space horizontally.')

904

905 a_height = axes_height_avail / ny

906 a_width = axes_width_avail / (nx + 2)

907

908 a = max(min(a_height, a_width), a_min)

909

910 pad_height = (a_height - a) * ny

911 pad_width = (a_width - a) * (nx + 2)

912

913 for iy in range(ny):

914 y = height - 0.5 * pad_height - margin_top \

915 - (iy + 1) * a - iy * spacing_height

916 h = a

917 for ix in range(nx):

918 x = margin_right + 0.5 * pad_width + ix * (a + spacing_width)

919 w = a if ix != (nx - 1) else a * 3.0

920 axes[iy][ix].set_position(

921 rect_to_figure_coords((x, y, w, h)), which='both')

922

923 def call(self):

924

925 self.cleanup()

926

927 if self.rotate_to_rt != self.last_rotate_to_rt:

928 # reset delta azimuth to avoid confusion

929

930 self.set_parameter('azimuth', 0.0)

931

932 self.last_rotate_to_rt = self.rotate_to_rt

933

934 selected = self.get_selected()

935

936 nsl_to_bazi = self.get_bazis()

937

938 tpad = self.t_length

939 groups = self.get_traces(selected, nsl_to_bazi, tpad)

940

941 if not groups:

942 self.fail(

943 'No matching traces. Check time and channel settings. Traces '

944 'may not contain gaps within the extracted time window and in '

945 'the padding areas left and right. Traces are extracted with '

946 'additional padding of 3 x filter corner period to eliminate '

947 'artifacts.')

948

949 selection_markers = self.make_selection_markers(selected, groups)

950 self.add_markers(selection_markers)

951

952 fig = self.get_figure()

953

954 figure_key = (len(groups), self.iframe)

955

956 if not self.figure_key or self.figure_key != figure_key:

957 self.setup_figure(fig, len(groups))

958 self.figure_key = figure_key

959

960 self.draw(groups, selected, tpad, nsl_to_bazi)

961

962 self.layout(fig, self.axes)

963

964 self.fframe.draw()

965 tabs = self.fframe.parent().parent()

966 # bring plot to front if we are not looking at the markers

967 from pyrocko.gui.pile_viewer import PileViewer

968 if not isinstance(tabs.currentWidget(), PileViewer):

969 tabs.setCurrentWidget(self.fframe)

970

971

972def __snufflings__():

973 return [Polarization()]