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

447 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2024-03-07 11:54 +0000

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 

7 

8logger = logging.getLogger('pyrocko.gui.snufflings.polarization') 

9 

10d2r = num.pi / 180. 

11r2d = 1.0 / d2r 

12 

13 

14def get_callbacks(obj): 

15 try: 

16 return obj.callbacks 

17 except AttributeError: 

18 return obj._callbacks 

19 

20 

21def get_tr(traces, f): 

22 for tr in traces: 

23 if f(tr): 

24 return tr 

25 

26 raise ValueError('Element not found.') 

27 

28 

29def darken(color, f=0.5): 

30 return tuple(c*f for c in color[:3]) + color[3:] 

31 

32 

33class PCAError(Exception): 

34 pass 

35 

36 

37class LayoutError(Exception): 

38 pass 

39 

40 

41class Polarization(Snuffling): 

42 

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> 

52 

53<h1> 

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

55</h1> 

56 

57<p> 

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</p> 

66 

67<h2> 

68Usage 

69</h2> 

70 

71<p> 

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

73press <strong>Run</strong>. Multiple stations can be selected for direct 

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

75<strong>Channels</strong> 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</p> 

80 

81<p> 

82The time window can be adjusted with the <strong>Length</strong> and 

83<strong>Offset</strong> parameters. Extracted waveforms are filtered according 

84the <strong>Highpass</strong> and <strong>Lowpass</strong> parameters 

85(Butterworth 4th order) and demeaned. 

86</p> 

87 

88<p> 

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

90<strong>\u0394 Azimuth</strong> 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</p> 

95 

96<p> 

97If <strong>Show 2D Eigensystems</strong> 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 <strong>Show 3D Eigensystems</strong> is selected, a principal 

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

102the resulting eigenvectors are depicted with purple lines. 

103</p> 

104 

105<p> 

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 

109<strong>Fix Scale</strong>. 

110</p> 

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()]