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

445 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-04 09:52 +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('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()]