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(

323 'No data. Time length too short?', action='status')

324

325 y = tr.get_ydata()

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

327

328 group.append(tr)

329

330 group = self.sorted_by_component_scheme(group)

331

332 tr_e = group[0]

333 tr_n = group[1]

334

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

336 cha_e = tr_e.channel

337 cha_n = tr_n.channel

338 if bazimuth != 0.0:

339 trs_rot = trace.rotate(

340 [tr_n, tr_e], bazimuth,

341 in_channels=[cha_n, cha_e],

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

343 else:

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

345 for tr in trs_rot:

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

347

348 if len(trs_rot) == 2:

349 group.extend(

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

351 for ch in [cha_e, cha_n])

352 else:

353 self.fail(

354 'Cannot rotate traces. '

355 'Are the samples properly aligned?')

356

357 d[k] = group

358

359 return d

360

361 def setup_figure_frame(self):

362 self.iframe += 1

363 self.fframe = self.figure_frame(

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

365

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

367

368 def get_figure(self):

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

370 self.setup_figure_frame()

371

372 return self.fframe.gcf()

373

374 def setup_figure(self, fig, nstations):

375 fig.clf()

376 new_axes = []

377

378 def iwrap(iy, ix):

379 return (ix + iy * 4) + 1

380

381 for istation in range(nstations):

382 axes_01 = fig.add_subplot(

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

384 axes_02 = fig.add_subplot(

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

386 axes_12 = fig.add_subplot(

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

388 axes_tr = fig.add_subplot(

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

390

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

392 axes.my_stuff = []

393

394 axes.my_line, = axes.plot(

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

396 axes.my_dot, = axes.plot(

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

398

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

400 axes.get_xaxis().set_tick_params(

401 labelbottom=False, bottom=False)

402 axes.get_yaxis().set_tick_params(

403 labelleft=False, left=False)

404

405 axes_tr.get_yaxis().set_tick_params(

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

407

408 # if istation != nstations - 1:

409 # axes_tr.get_xaxis().set_tick_params(

410 # bottom=False, labelbottom=False)

411

412 lines = []

413 dots = []

414 for _ in range(3):

415 lines.append(

416 axes_tr.plot(

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

418 dots.append(

419 axes_tr.plot(

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

421

422 axes_tr.my_lines = lines

423 axes_tr.my_dots = dots

424 axes_tr.my_stuff = []

425

426 new_axes.append(

427 (axes_01, axes_02, axes_12, axes_tr))

428

429 def resize_handler(*args):

430 self.layout(fig, new_axes)

431

432 if fig.my_disconnect:

433 fig.my_disconnect()

434

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

436 try:

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

438 except (AttributeError, ValueError):

439 pass # not available anymore in MPL >= 3.8

440

441 def disconnect():

442 fig.canvas.mpl_disconnect(cid_resize)

443 fig.callbacks.disconnect(cid_dpi)

444

445 fig.my_disconnect = disconnect

446

447 self.axes = new_axes

448

449 def get_trace(self, traces, pattern):

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

451 if len(trs) > 1:

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

453 elif len(trs) == 0:

454 return None

455 else:

456 return trs[0]

457

458 def get_vector_abs_max(self, traces):

459 tr_abs = None

460 for tr in traces:

461 if tr is not None:

462 tr = tr.copy()

463 tr.ydata **= 2

464 if tr_abs is None:

465 tr_abs = tr

466 else:

467 tr_abs.add(tr)

468

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

470 return num.max(tr_abs.ydata)

471

472 def set_labels(

473 self, istation, nstations, tref,

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

475

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

477 rcolor = self.colors['rot']

478 else:

479 rcolor = self.colors['fg']

480 chas_rot = chas

481

482 axes_01.set_ylabel(chas[1])

483 axes_02.set_ylabel(chas[2])

484 axes_12.set_ylabel(chas[2])

485

486 axes_01.set_xlabel(chas[0])

487 axes_02.set_xlabel(chas_rot[0])

488 axes_12.set_xlabel(chas_rot[1])

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

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

491 axes_tr.set_xlabel(

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

493

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

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

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

497 tlab.set_color(rcolor)

498

499 def pca(self, trs):

500

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

502 raise PCAError('Missing component')

503

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

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

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

507

508 ns = nss[0]

509

510 if ns < 3:

511 raise PCAError('Traces too short.')

512

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

514 for itr, tr in enumerate(trs):

515 data[itr, :] = tr.ydata

516

517 cov = num.cov(data)

518

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

520

521 pc = evecs[:, -1]

522

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

524

525 if len(trs) > 2:

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

527 else:

528 incidence = 90.

529

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

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

532

533 return cov, evals, evecs, azimuth, incidence

534

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

536 evals = num.sqrt(evals)

537 ell = patches.Ellipse(

538 xy=(0.0, 0.0),

539 width=evals[0] * 2.,

540 height=evals[1] * 2.,

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

542 zorder=-10,

543 fc=color,

544 ec=darken(color),

545 alpha=alpha)

546

547 return ell

548

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

550

551 def count(d, k):

552 if k not in d:

553 d[k] = 0

554

555 d[k] += 1

556

557 nsli_n = {}

558 for k in sorted(groups):

559 count(nsli_n, k[:4])

560

561 nsli_i = {}

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

563

564 tmin, tmax = nsl_to_tspan[k]

565 nsl = k[:3]

566 nsli = k[:4]

567 count(nsli_i, nsli)

568

569 bazimuth = self.azimuth

570

571 if self.rotate_to_rt:

572 if nsl not in nsl_to_bazi:

573 self.fail(

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

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

576 'coordinates must be available.')

577

578 bazimuth += nsl_to_bazi[nsl] + 90.

579

580 for axes in self.axes[igroup]:

581 while axes.my_stuff:

582 stuff = axes.my_stuff.pop()

583 stuff.remove()

584

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

586

587 axes_01.set_title(

588 '.'.join(nsli)

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

590

591 trs_all = groups[k]

592

593 if len(trs_all) == 5:

594 trs_orig = trs_all[:3]

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

596 elif len(trs_all) == 4:

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

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

599 else:

600 raise self.fail(

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

602 % len(trs_all))

603

604 trs_orig_chopped = [

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

606 for tr in trs_orig]

607

608 trs_rot_chopped = [

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

610 for tr in trs_rot]

611

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

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

614

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

616 amax = self.nsl_to_amax[nsl]

617 else:

618 amax = self.get_vector_abs_max(trs_orig_chopped)

619 self.nsl_to_amax[nsl] = amax

620

621 if nsl in nsl_to_bazi:

622 event_bazi = nsl_to_bazi[nsl]

623 l1, = axes_01.plot(

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

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

626 '--',

627 lw=3,

628 zorder=-20,

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

630

631 a1 = axes_01.annotate(

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

633 (0., 1.),

634 (self.font_size, -self.font_size),

635 xycoords='axes fraction',

636 textcoords='offset points',

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

638 ha='left',

639 va='top')

640

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

642

643 for ix, iy, axes, trs in [

644 (0, 1, axes_01, trs_orig_chopped),

645 (0, 2, axes_02, trs_rot_chopped),

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

647

648 if amax == 0.0:

649 amax = 1.0

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

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

652

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

654 continue

655

656 x = trs[ix].get_ydata()

657 y = trs[iy].get_ydata()

658

659 axes.my_line.set_data(x, y)

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

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

662

663 tref = tmin

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

665 trs_rot, trs_rot_chopped)):

666

667 if tr is None or tr_chopped is None:

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

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

670

671 else:

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

673 t = tr.get_xdata()

674 t = t - tref

675

676 ipos = int(round(

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

678

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

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

681

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

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

684

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

686 color = self.colors['rot']

687

688 xn = num.sin(bazimuth*d2r)

689 yn = num.cos(bazimuth*d2r)

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

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

692

693 l1, = axes_01.plot(

694 [0., amax*xn],

695 [0., amax*yn],

696 color=color)

697

698 font_size = self.font_size

699

700 a1 = axes_01.annotate(

701 chas_rot[1],

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

703 xycoords='data',

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

705 textcoords='offset points',

706 va='center',

707 ha='center',

708 color=color)

709

710 l2, = axes_01.plot(

711 [0., amax*xe],

712 [0., amax*ye],

713 color=color)

714

715 a2 = axes_01.annotate(

716 chas_rot[0],

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

718 xycoords='data',

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

720 textcoords='offset points',

721 va='center',

722 ha='center',

723 color=color)

724

725 aa = axes_01.annotate(

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

727 (1., 1.),

728 (-self.font_size, -self.font_size),

729 xycoords='axes fraction',

730 textcoords='offset points',

731 color=color,

732 ha='right',

733 va='top')

734

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

736

737 axes_tr.my_stuff.append(axes_tr.axvspan(

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

739

740 axes_tr.set_ylim(-1, 3)

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

742

743 self.set_labels(

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

745

746 if self.show_eigensystem_2d:

747

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

749 (0, 1, axes_01, trs_orig_chopped),

750 (0, 2, axes_02, trs_rot_chopped),

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

752

753 try:

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

755 [trs[ix], trs[iy]])

756

757 axes.my_stuff.append(

758 axes.annotate(

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

760 (0., 0.),

761 (self.font_size, self.font_size),

762 xycoords='axes fraction',

763 textcoords='offset points',

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

765 alpha=0.5))

766

767 ell = self.draw_cov_ellipse(

768 evals[:2], evecs[:2, :2],

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

770

771 axes.add_artist(ell)

772 axes.my_stuff.append(ell)

773

774 l1, = axes.plot(

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

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

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

778

779 l2, = axes.plot(

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

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

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

783

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

785

786 except PCAError as e:

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

788

789 if self.show_eigensystem_3d:

790 try:

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

792 trs_orig_chopped)

793 cosa = num.cos(bazimuth*d2r)

794 sina = num.sin(bazimuth*d2r)

795 rot = num.array(

796 [[cosa, -sina, 0.0],

797 [sina, cosa, 0.0],

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

799

800 evecs_rot = num.dot(rot, evecs)

801

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

803 (0, 1, axes_01, evecs),

804 (0, 2, axes_02, evecs_rot),

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

806

807 for (ie, alpha) in [

808 (-1, 0.5),

809 (-2, 0.5),

810 (-3, 0.5)]:

811

812 for vlen, lw in [

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

814 (amax, 1)]:

815

816 lv, = axes.plot(

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

818 vlen*evecs_[ix, ie]],

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

820 vlen*evecs_[iy, ie]],

821 lw=lw,

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

823

824 axes.my_stuff.extend([lv])

825

826 axes_01.my_stuff.append(

827 axes_01.annotate(

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

829 (1., 0.),

830 (-self.font_size, self.font_size),

831 xycoords='axes fraction',

832 textcoords='offset points',

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

834 alpha=0.8,

835 ha='right'))

836

837 axes_02.my_stuff.append(

838 axes_02.annotate(

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

840 (1., 0.),

841 (-self.font_size, self.font_size),

842 xycoords='axes fraction',

843 textcoords='offset points',

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

845 alpha=0.8,

846 ha='right'))

847

848 except PCAError as e:

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

850

851 def get_bazis(self):

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

853 if not event:

854 return {}

855

856 nsl_to_bazi = dict(

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

858 for station in self.get_stations())

859

860 return nsl_to_bazi

861

862 def layout(self, fig, axes):

863

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

865

866 def get_pixels_factor(fig):

867 try:

868 r = fig.canvas.get_renderer()

869 return 1.0 / r.points_to_pixels(1.0)

870 except AttributeError:

871 return 1.0

872

873 def rect_to_figure_coords(rect):

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

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

876

877 ny = len(axes)

878 if ny == 0:

879 raise LayoutError('No axes given.')

880

881 nx = len(axes[0])

882

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

884 pixels = get_pixels_factor(fig)

885

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

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

888

889 spacing_width = 3. * self.font_size / pixels

890 spacing_height = 5. * self.font_size / pixels

891

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

893 - margin_top - margin_bottom

894

895 a_min = 5.

896

897 if axes_height_avail < a_min * ny:

898 axes_height_avail = a_min * ny

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

900

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

902 - margin_left - margin_right

903

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

905 axes_width_avail = a_min * (nx + 2)

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

907

908 a_height = axes_height_avail / ny

909 a_width = axes_width_avail / (nx + 2)

910

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

912

913 pad_height = (a_height - a) * ny

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

915

916 for iy in range(ny):

917 y = height - 0.5 * pad_height - margin_top \

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

919 h = a

920 for ix in range(nx):

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

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

923 axes[iy][ix].set_position(

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

925

926 def call(self):

927

928 self.cleanup()

929

930 if self.rotate_to_rt != self.last_rotate_to_rt:

931 # reset delta azimuth to avoid confusion

932

933 self.set_parameter('azimuth', 0.0)

934

935 self.last_rotate_to_rt = self.rotate_to_rt

936

937 selected = self.get_selected()

938

939 nsl_to_bazi = self.get_bazis()

940

941 tpad = self.t_length

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

943

944 if not groups:

945 self.fail(

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

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

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

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

950 'artifacts.')

951

952 selection_markers = self.make_selection_markers(selected, groups)

953 self.add_markers(selection_markers)

954

955 fig = self.get_figure()

956

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

958

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

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

961 self.figure_key = figure_key

962

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

964

965 self.layout(fig, self.axes)

966

967 self.fframe.draw()

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

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

970 from pyrocko.gui.pile_viewer import PileViewer

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

972 tabs.setCurrentWidget(self.fframe)

973

974

975def __snufflings__():

976 return [Polarization()]