Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/plot/dynamic_rupture.py: 82%

633 statements  

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

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6from subprocess import check_call, CalledProcessError 

7 

8import os.path as op 

9import re 

10import logging 

11import tempfile 

12import shutil 

13import inspect 

14 

15import numpy as num 

16from scipy.interpolate import RegularGridInterpolator as scrgi 

17 

18from matplotlib import pyplot as plt, patheffects 

19from matplotlib.ticker import FuncFormatter 

20 

21from pyrocko import orthodrome as pod 

22from pyrocko.guts import Object 

23from pyrocko.plot import mpl_init, mpl_papersize, mpl_color, AutoScaler, \ 

24 gmtpy, mpl_get_cmap 

25from pyrocko.plot.automap import Map, NoTopo 

26from pyrocko.gf import PseudoDynamicRupture 

27from pyrocko.gf.seismosizer import map_anchor 

28from pyrocko.dataset.topo.tile import Tile 

29 

30logger = logging.getLogger(__name__) 

31 

32km = 1e3 

33 

34d2m = 111180. 

35m2d = 1. / d2m 

36 

37cm2inch = gmtpy.cm/gmtpy.inch 

38 

39d2r = num.pi / 180. 

40r2d = 1. / d2r 

41 

42 

43def km_formatter(v, p): 

44 return '%g' % (v / km) 

45 

46 

47def get_kwargs(cls): 

48 sig = inspect.signature(cls.__init__) 

49 kwargs = {} 

50 

51 if issubclass(cls.T._cls, RuptureMap): 

52 for p in Map.T.properties: 

53 kwargs[p.name] = p._default 

54 

55 for par in sig.parameters.values(): 

56 if par.default is inspect._empty: 

57 continue 

58 kwargs[par.name] = par.default 

59 

60 return kwargs 

61 

62 

63def _save_grid(lats, lons, data, filename): 

64 ''' 

65 Save lat-lon gridded data as gmt .grd file. 

66 

67 :param lats: 

68 Grid latitude coordinates in [deg]. 

69 :type lats: 

70 iterable 

71 

72 :param lons: 

73 Grid longitude coordinates in [deg]. 

74 :type lons: 

75 iterable 

76 

77 :param data: 

78 Grid data of any kind. 

79 :type data: 

80 :py:class:`~numpy.ndarray`, ``(n_lons, n_lats)`` 

81 

82 :param filename: 

83 Filename of the written grid file. 

84 :type filename: 

85 str 

86 ''' 

87 

88 gmtpy.savegrd(lons, lats, data, filename=filename, naming='lonlat') 

89 

90 

91def _mplcmap_to_gmtcpt_code(mplcmap, steps=256): 

92 ''' 

93 Get gmt readable R/G/A code from a given matplotlib colormap. 

94 

95 :param mplcmap: 

96 Name of the demanded matplotlib colormap. 

97 :type mplcmap: 

98 str 

99 

100 :returns: 

101 Series of comma seperate R/G/B values for gmtpy usage. 

102 :rtype: 

103 str 

104 ''' 

105 

106 cmap = mpl_get_cmap(mplcmap) 

107 

108 rgbas = [cmap(i) for i in num.linspace(0, 255, steps).astype(num.int64)] 

109 

110 return ','.join(['%d/%d/%d' % ( 

111 c[0] * 255, c[1] * 255, c[2] * 255) for c in rgbas]) 

112 

113 

114def make_cpt_path(gmt, base_cmap): 

115 return gmt.tempfilename('my_%s.cpt' % base_cmap) 

116 

117 

118def make_colormap( 

119 gmt, 

120 vmin, 

121 vmax, 

122 C=None, 

123 cmap=None, 

124 space=False): 

125 ''' 

126 Create gmt-readable colormap cpt file called my_<cmap>.cpt. 

127 

128 :param vmin: 

129 Minimum value covered by the colormap. 

130 :type vmin: 

131 float 

132 

133 :param vmax: 

134 Maximum value covered by the colormap. 

135 :type vmax: 

136 float 

137 

138 :param C: 

139 Comma seperated R/G/B values for cmap definition. 

140 :type C: 

141 str 

142 

143 :param cmap: 

144 Name of the colormap. Colormap is stored as "my_<cmap>.cpt". 

145 If name is equivalent to a matplotlib colormap, R/G/B strings are 

146 extracted from this colormap. 

147 :type cmap: 

148 str 

149 

150 :param space: 

151 If ``True``, the range of the colormap is broadened below vmin 

152 and above vmax. 

153 :type space: bool 

154 ''' 

155 

156 scaler = AutoScaler(mode='min-max') 

157 scale = scaler.make_scale((vmin, vmax)) 

158 

159 incr = scale[2] 

160 margin = 0. 

161 

162 if vmin == vmax: 

163 space = True 

164 

165 if space: 

166 margin = incr 

167 

168 msg = ('Please give either a valid color code or a' 

169 ' valid matplotlib colormap name.') 

170 

171 if C is None and cmap is None: 

172 raise ValueError(msg) 

173 

174 if C is None: 

175 try: 

176 C = _mplcmap_to_gmtcpt_code(cmap) 

177 except ValueError: 

178 raise ValueError(msg) 

179 

180 if cmap is None: 

181 logger.warning( 

182 'No colormap name given. Uses temporary filename instead') 

183 cmap = 'temp_cmap' 

184 

185 return gmt.makecpt( 

186 C=C, 

187 D='o', 

188 T='%g/%g/%g' % ( 

189 vmin - margin, vmax + margin, incr), 

190 Z=True, 

191 out_filename=make_cpt_path(gmt, cmap), 

192 suppress_defaults=True) 

193 

194 

195def xy_to_latlon(source, x, y): 

196 ''' 

197 Convert x and y relative coordinates on extended ruptures into latlon. 

198 

199 :param source: 

200 Extended source class, on which the given point is located. 

201 :type source: 

202 :py:class:`~pyrocko.gf.seismosizer.RectangularSource` or 

203 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture` 

204 

205 :param x: 

206 Relative point coordinate along strike (range: -1:1). 

207 :type x: 

208 :py:class:`float` or :py:class:`~numpy.ndarray` 

209 

210 :param y: 

211 Relative downdip point coordinate (range: -1:1). 

212 :type y: 

213 :py:class:`float` or :py:class:`~numpy.ndarray` 

214 

215 :returns: 

216 Latitude and Longitude of the given point in [deg]. 

217 :rtype: 

218 :py:class:`tuple` of :py:class:`float` 

219 ''' 

220 

221 s = source 

222 ax, ay = map_anchor[s.anchor] 

223 

224 length, width = (x - ax) / 2. * s.length, (y - ay) / 2. * s.width 

225 strike, dip = s.strike * d2r, s.dip * d2r 

226 

227 northing = num.cos(strike)*length - num.cos(dip) * num.sin(strike)*width 

228 easting = num.sin(strike)*length + num.cos(dip) * num.cos(strike)*width 

229 

230 northing += source.north_shift 

231 easting += source.east_shift 

232 

233 return pod.ne_to_latlon(s.lat, s.lon, northing, easting) 

234 

235 

236def xy_to_lw(source, x, y): 

237 ''' 

238 Convert relative coordinates on extended ruptures into length and width. 

239 

240 :param source: 

241 Extended source, on which the given points are located. 

242 :type source: 

243 :py:class:`~pyrocko.gf.seismosizer.RectangularSource` or 

244 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture` 

245 

246 :param x: 

247 Relative point coordinates along strike (range: -1:1). 

248 :type x: 

249 float or :py:class:`~numpy.ndarray` 

250 

251 :param y: 

252 Relative downdip point coordinates (range: -1:1). 

253 :type y: 

254 float or :py:class:`~numpy.ndarray` 

255 

256 :returns: 

257 Length and downdip width of the given points from the anchor in [m]. 

258 :rtype: 

259 :py:class:`tuple` of :py:class:`float` 

260 ''' 

261 

262 length, width = source.length, source.width 

263 

264 ax, ay = map_anchor[source.anchor] 

265 

266 lengths = (x - ax) / 2. * length 

267 widths = (y - ay) / 2. * width 

268 

269 return lengths, widths 

270 

271 

272cbar_anchor = { 

273 'center': 'MC', 

274 'center_left': 'ML', 

275 'center_right': 'MR', 

276 'top': 'TC', 

277 'top_left': 'TL', 

278 'top_right': 'TR', 

279 'bottom': 'BC', 

280 'bottom_left': 'BL', 

281 'bottom_right': 'BR'} 

282 

283 

284cbar_helper = { 

285 'traction': { 

286 'unit': 'MPa', 

287 'factor': 1e-6}, 

288 'tx': { 

289 'unit': 'MPa', 

290 'factor': 1e-6}, 

291 'ty': { 

292 'unit': 'MPa', 

293 'factor': 1e-6}, 

294 'tz': { 

295 'unit': 'MPa', 

296 'factor': 1e-6}, 

297 'time': { 

298 'unit': 's', 

299 'factor': 1.}, 

300 'strike': { 

301 'unit': '°', 

302 'factor': 1.}, 

303 'dip': { 

304 'unit': '°', 

305 'factor': 1.}, 

306 'vr': { 

307 'unit': 'km/s', 

308 'factor': 1e-3}, 

309 'length': { 

310 'unit': 'km', 

311 'factor': 1e-3}, 

312 'width': { 

313 'unit': 'km', 

314 'factor': 1e-3} 

315} 

316 

317 

318fonttype = 'Helvetica' 

319 

320c2disl = dict([('x', 0), ('y', 1), ('z', 2)]) 

321 

322 

323def _make_gmt_conf(fontcolor, size): 

324 ''' 

325 Update different gmt parameters depending on figure size and fontcolor. 

326 

327 :param fontcolor: 

328 GMT readable colorcode or colorstring for the font. 

329 :type fontcolor: 

330 str 

331 

332 :param size: 

333 Tuple of the figure size (width, height) in [cm]. 

334 :type size: 

335 :py:class:`tuple` of :py:class:`float` 

336 

337 :returns: 

338 Best fitting fontsize, GMT configuration parameter 

339 :rtype: 

340 float, dict 

341 ''' 

342 

343 color = fontcolor 

344 fontsize = num.max(size) 

345 

346 font = '%gp,%s' % (fontsize, fonttype) 

347 

348 pen_size = fontsize / 16. 

349 tick_size = num.min(size) / 200. 

350 

351 return fontsize, dict( 

352 MAP_FRAME_PEN='%s' % color, 

353 MAP_TICK_PEN_PRIMARY='%gp,%s' % (pen_size, color), 

354 MAP_TICK_PEN_SECONDARY='%gp,%s' % (pen_size, color), 

355 MAP_TICK_LENGTH_PRIMARY='%gc' % tick_size, 

356 MAP_TICK_LENGTH_SECONDARY='%gc' % (tick_size * 3), 

357 FONT_ANNOT_PRIMARY='%s-Bold,%s' % (font, color), 

358 FONT_LABEL='%s-Bold,%s' % (font, color), 

359 FONT_TITLE='%s-Bold,%s' % (font, color), 

360 PS_CHAR_ENCODING='ISOLatin1+', 

361 MAP_FRAME_TYPE='fancy', 

362 FORMAT_GEO_MAP='D', 

363 PS_PAGE_ORIENTATION='portrait', 

364 MAP_GRID_PEN_PRIMARY='thinnest,%s' % color, 

365 MAP_ANNOT_OBLIQUE='6') 

366 

367 

368class SourceError(Exception): 

369 pass 

370 

371 

372class RuptureMap(Map): 

373 ''' 

374 Map plotting of attributes and results of the 

375 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`. 

376 ''' 

377 

378 def __init__( 

379 self, 

380 source=None, 

381 fontcolor='darkslategrey', 

382 width=20., 

383 height=14., 

384 margins=None, 

385 color_wet=(216, 242, 254), 

386 color_dry=(238, 236, 230), 

387 topo_cpt_wet='light_sea_uniform', 

388 topo_cpt_dry='light_land_uniform', 

389 show_cities=False, 

390 *args, 

391 **kwargs): 

392 

393 size = (width, height) 

394 fontsize, gmt_config = _make_gmt_conf(fontcolor, size) 

395 

396 if margins is None: 

397 margins = [ 

398 fontsize * 0.15, num.min(size) / 200., 

399 num.min(size) / 200., fontsize * 0.05] 

400 

401 Map.__init__(self, *args, margins=margins, width=width, height=height, 

402 gmt_config=gmt_config, 

403 topo_cpt_dry=topo_cpt_dry, topo_cpt_wet=topo_cpt_wet, 

404 color_wet=color_wet, color_dry=color_dry, 

405 **kwargs) 

406 

407 if show_cities: 

408 self.draw_cities() 

409 

410 self._source = source 

411 self._fontcolor = fontcolor 

412 self._fontsize = fontsize 

413 self.axes_layout = 'WeSn' 

414 

415 @property 

416 def size(self): 

417 ''' 

418 Figure size in [cm]. 

419 ''' 

420 return (self.width, self.height) 

421 

422 @property 

423 def font(self): 

424 ''' 

425 Font style (size and type). 

426 ''' 

427 

428 return '%sp,%s' % (self._fontsize, fonttype) 

429 

430 @property 

431 def source(self): 

432 ''' 

433 PseudoDynamicRupture whose attributes are plotted. 

434 

435 Note, that source.patches attribute needs to be calculated in advance. 

436 ''' 

437 

438 if self._source is None: 

439 raise SourceError('No source given. Please define it!') 

440 

441 if not isinstance(self._source, PseudoDynamicRupture): 

442 raise SourceError('This function is only capable for a source of' 

443 ' type: %s' % type(PseudoDynamicRupture())) 

444 

445 if not self._source.patches: 

446 raise TypeError('No source patches are defined. Please run' 

447 '"source.discretize_patches()" on your source') 

448 

449 return self._source 

450 

451 @source.setter 

452 def source(self, source): 

453 self._source = source 

454 

455 def _get_topotile(self): 

456 if self._dems is None: 

457 self._setup() 

458 

459 try: 

460 t, _ = self._get_topo_tile('land') 

461 except NoTopo: 

462 wesn = self.wesn 

463 

464 nx = int(num.ceil( 

465 self.width * cm2inch * self.topo_resolution_max)) 

466 ny = int(num.ceil( 

467 self.height * cm2inch * self.topo_resolution_max)) 

468 

469 data = num.zeros((nx, ny)) 

470 

471 t = Tile(wesn[0], wesn[2], 

472 (wesn[1] - wesn[0]) / nx, (wesn[3] - wesn[2]) / ny, 

473 data) 

474 

475 return t 

476 

477 def _patches_to_lw(self): 

478 ''' 

479 Generate regular rect. length-width grid based on the patch distance. 

480 

481 Prerequesite is a regular grid of patches (constant lengths and widths) 

482 Both coordinates are given relative to the source anchor point in [m] 

483 The grid is extended from the patch centres to the edges of the source. 

484 

485 :returns: 

486 Lengths along strike, widths downdip in [m]. 

487 :rtype: 

488 :py:class:`~numpy.ndarray`, :py:class:`~numpy.ndarray` 

489 ''' 

490 

491 source = self.source 

492 patches = source.patches 

493 

494 patch_l, patch_w = patches[0].length, patches[0].width 

495 

496 patch_lengths = num.concatenate(( 

497 num.array([0.]), 

498 num.array([il*patch_l+patch_l/2. for il in range(source.nx)]), 

499 num.array([patch_l * source.nx]))) 

500 

501 patch_widths = num.concatenate(( 

502 num.array([0.]), 

503 num.array([iw*patch_w+patch_w/2. for iw in range(source.ny)]), 

504 num.array([patch_w * source.ny]))) 

505 

506 ax, ay = map_anchor[source.anchor] 

507 

508 patch_lengths -= source.length * (ax + 1.) / 2. 

509 patch_widths -= source.width * (ay + 1.) / 2. 

510 

511 return patch_lengths, patch_widths 

512 

513 def _xy_to_lw(self, x, y): 

514 ''' 

515 Generate regular rect. length-width grid based on the xy coordinates. 

516 

517 Prerequesite is a regular grid with constant dx and dy. x and y are 

518 relative coordinates on the rupture plane (range -1:1) along strike and 

519 downdip. 

520 Length and width are obtained relative to the source anchor point 

521 in [m]. 

522 

523 :returns: 

524 Lengths along strike, widths downdip in [m]. 

525 :rtype: 

526 :py:class:`~numpy.ndarray`, :py:class:`~numpy.ndarray` 

527 ''' 

528 

529 x, y = num.unique(x), num.unique(y) 

530 dx, dy = x[1] - x[0], y[1] - y[0] 

531 

532 if any(num.abs(num.diff(x) - dx) >= 1e-6): 

533 raise ValueError('Regular grid with constant spacing needed.' 

534 ' Please check the x coordinates.') 

535 

536 if any(num.abs(num.diff(y) - dy) >= 1e-6): 

537 raise ValueError('Regular grid with constant spacing needed.' 

538 ' Please check the y coordinates.') 

539 

540 return xy_to_lw(self.source, x, y) 

541 

542 def _tile_to_lw(self, ref_lat, ref_lon, 

543 north_shift=0., east_shift=0., strike=0.): 

544 

545 ''' 

546 Coordinate transformation from the topo tile grid into length-width. 

547 

548 The topotile latlon grid is rotated into the length width grid. The 

549 width is defined here as its horizontal component. The resulting grid 

550 is used for interpolation of grid data. 

551 

552 :param ref_lat: 

553 Reference latitude, from which length-width relative coordinates 

554 grid are calculated. 

555 :type ref_lat: 

556 float 

557 

558 :param ref_lon: 

559 Reference longitude, from which length-width relative coordinates 

560 grid are calculated. 

561 :type ref_lon: 

562 float 

563 

564 :param north_shift: 

565 North shift of the reference point in [m]. 

566 :type north_shift: 

567 float 

568 

569 :param east_shift: 

570 East shift of the reference point [m]. 

571 :type east_shift: 

572 float 

573 

574 :param strike: 

575 striking of the length axis compared to the North axis in [deg]. 

576 :type strike: 

577 float 

578 

579 :returns: 

580 Topotile grid nodes as array of length-width coordinates. 

581 :rtype: 

582 :py:class:`~numpy.ndarray`, ``(n_nodes, 2)`` 

583 ''' 

584 

585 t = self._get_topotile() 

586 grid_lats = t.y() 

587 grid_lons = t.x() 

588 

589 meshgrid_lons, meshgrid_lats = num.meshgrid(grid_lons, grid_lats) 

590 

591 grid_northing, grid_easting = pod.latlon_to_ne_numpy( 

592 ref_lat, ref_lon, meshgrid_lats.flatten(), meshgrid_lons.flatten()) 

593 

594 grid_northing -= north_shift 

595 grid_easting -= east_shift 

596 

597 strike *= d2r 

598 sin, cos = num.sin(strike), num.cos(strike) 

599 points_out = num.zeros((grid_northing.shape[0], 2)) 

600 points_out[:, 0] = -sin * grid_northing + cos * grid_easting 

601 points_out[:, 1] = cos * grid_northing + sin * grid_easting 

602 

603 return points_out 

604 

605 def _prep_patch_grid_data(self, data): 

606 ''' 

607 Extend patch data from patch centres to the outer source edges. 

608 

609 Patch data is always defined in the centre of the patches. For 

610 interpolation the data is extended here to the edges of the rupture 

611 plane. 

612 

613 :param data: 

614 Patch wise data. 

615 :type data: 

616 :py:class:`~numpy.ndarray` 

617 

618 :returns: 

619 Extended data array. 

620 :rtype: 

621 :py:class:`~numpy.ndarray` 

622 ''' 

623 

624 shape = (self.source.nx + 2, self.source.ny + 2) 

625 data = data.reshape(self.source.nx, self.source.ny).copy() 

626 

627 data_new = num.zeros(shape) 

628 data_new[1:-1, 1:-1] = data 

629 data_new[1:-1, 0] = data[:, 0] 

630 data_new[1:-1, -1] = data[:, -1] 

631 data_new[0, 1:-1] = data[0, :] 

632 data_new[-1, 1:-1] = data[-1, :] 

633 

634 for i, j in zip([-1, -1, 0, 0], [-1, 0, -1, 0]): 

635 data_new[i, j] = data[i, j] 

636 

637 return data_new 

638 

639 def _regular_data_to_grid(self, lengths, widths, data, filename): 

640 ''' 

641 Interpolate regular data onto topotile grid. 

642 

643 Regular gridded data is interpolated onto the latlon grid of the 

644 topotile. It is then stored as a gmt-readable .grd-file. 

645 

646 :param lengths: 

647 Grid coordinates along strike relative to anchor in [m]. 

648 :type lengths: 

649 :py:class:`~numpy.ndarray` 

650 

651 :param widths: 

652 Grid coordinates downdip relative to anchor in [m]. 

653 :type widths: 

654 :py:class:`~numpy.ndarray` 

655 

656 :param data: 

657 Data grid array. 

658 :type data: 

659 :py:class:`~numpy.ndarray` 

660 

661 :param filename: 

662 Filename, where grid is stored. 

663 :type filename: 

664 str 

665 ''' 

666 

667 source = self.source 

668 

669 interpolator = scrgi( 

670 (widths * num.cos(d2r * source.dip), lengths), 

671 data.T, 

672 bounds_error=False, 

673 method='nearest') 

674 

675 points_out = self._tile_to_lw( 

676 ref_lat=source.lat, 

677 ref_lon=source.lon, 

678 north_shift=source.north_shift, 

679 east_shift=source.east_shift, 

680 strike=source.strike) 

681 

682 t = self._get_topotile() 

683 t.data = num.zeros_like(t.data, dtype=float) 

684 t.data[:] = num.nan 

685 

686 t.data = interpolator(points_out).reshape(t.data.shape) 

687 

688 _save_grid(t.y(), t.x(), t.data, filename=filename) 

689 

690 def patch_data_to_grid(self, data, *args, **kwargs): 

691 ''' 

692 Generate a grid file based on regular patch wise data. 

693 

694 :param data: 

695 Patchwise grid data. 

696 :type data: 

697 :py:class:`~numpy.ndarray` 

698 ''' 

699 

700 lengths, widths = self._patches_to_lw() 

701 

702 data_new = self._prep_patch_grid_data(data) 

703 

704 self._regular_data_to_grid(lengths, widths, data_new, *args, **kwargs) 

705 

706 def xy_data_to_grid(self, x, y, data, *args, **kwargs): 

707 ''' 

708 Generate a grid file based on gridded data using xy coordinates. 

709 

710 Convert a grid based on relative fault coordinates (range -1:1) along 

711 strike (x) and downdip (y) into a .grd file. 

712 

713 :param x: 

714 Relative point coordinate along strike (range: -1:1). 

715 :type x: 

716 float or :py:class:`~numpy.ndarray` 

717 

718 :param y: 

719 Relative downdip point coordinate (range: -1:1). 

720 :type y: 

721 float or :py:class:`~numpy.ndarray` 

722 

723 :param data: 

724 Patchwise grid data. 

725 :type data: 

726 :py:class:`~numpy.ndarray` 

727 ''' 

728 

729 lengths, widths = self._xy_to_lw(x, y) 

730 

731 self._regular_data_to_grid( 

732 lengths, widths, data.reshape((lengths.shape[0], widths.shape[0])), 

733 *args, **kwargs) 

734 

735 def draw_image(self, gridfile, cmap, cbar=True, **kwargs): 

736 ''' 

737 Draw grid data as image and include, if whished, a colorbar. 

738 

739 :param gridfile: 

740 File of the grid which shall be plotted. 

741 :type gridfile: 

742 str 

743 

744 :param cmap: 

745 Name of the colormap, which shall be used. A .cpt-file 

746 "my_<cmap>.cpt" must exist. 

747 :type cmap: 

748 str 

749 

750 :param cbar: 

751 If ``True``, a colorbar corresponding to the grid data is 

752 added. Keyword arguments are parsed to it. 

753 :type cbar: 

754 bool 

755 ''' 

756 

757 self.gmt.grdimage( 

758 gridfile, 

759 C=make_cpt_path(self.gmt, cmap), 

760 E='200', 

761 Q=True, 

762 n='+t0.0', 

763 *self.jxyr) 

764 

765 if cbar: 

766 self.draw_colorbar(cmap=cmap, **kwargs) 

767 

768 def draw_contour( 

769 self, 

770 gridfile, 

771 contour_int, 

772 anot_int, 

773 angle=None, 

774 unit='', 

775 color='', 

776 style='', 

777 **kwargs): 

778 

779 ''' 

780 Draw grid data as contour lines. 

781 

782 :param gridfile: 

783 File of the grid which shall be plotted. 

784 :type gridfile: 

785 str 

786 

787 :param contour_int: 

788 Interval of contour lines in units of the gridfile. 

789 :type contour_int: 

790 float 

791 

792 :param anot_int: 

793 Interval of labelled contour lines in units of the gridfile. Must 

794 be a integer multiple of contour_int. 

795 :type anot_int: 

796 float 

797 

798 :param angle: 

799 Rotation angle of the labels in [deg]. 

800 :type angle: 

801 float 

802 

803 :param unit: 

804 Name of the unit in the grid file. It will be displayed behind the 

805 label on labelled contour lines. 

806 :type unit: 

807 str 

808 

809 :param color: 

810 GMT readable color code or string of the contour lines. 

811 :type color: 

812 str 

813 

814 :param style: 

815 Line style of the contour lines. If not given, solid lines are 

816 plotted. 

817 :type style: 

818 str 

819 ''' 

820 

821 pen_size = self._fontsize / 40. 

822 

823 if not color: 

824 color = self._fontcolor 

825 

826 a_string = '%g+f%s,%s+r%gc+u%s' % ( 

827 anot_int, self.font, color, pen_size*4, unit) 

828 if angle: 

829 a_string += '+a%g' % angle 

830 c_string = '%g' % contour_int 

831 

832 if kwargs: 

833 kwargs['A'], kwargs['C'] = a_string, c_string 

834 else: 

835 kwargs = dict(A=a_string, C=c_string) 

836 

837 if style: 

838 style = ',' + style 

839 

840 args = ['-Wc%gp,%s%s+s' % (pen_size, color, style)] 

841 

842 self.gmt.grdcontour( 

843 gridfile, 

844 S='10', 

845 W='a%gp,%s%s+s' % (pen_size*4, color, style), 

846 *self.jxyr + args, 

847 **kwargs) 

848 

849 def draw_colorbar(self, cmap, label='', anchor='top_right', **kwargs): 

850 ''' 

851 Draw a colorbar based on a existing colormap. 

852 

853 :param cmap: 

854 Name of the colormap, which shall be used. A .cpt-file 

855 "my_<cmap>.cpt" must exist. 

856 :type cmap: 

857 str 

858 

859 :param label: 

860 Title of the colorbar. 

861 :type label: 

862 str 

863 

864 :param anchor: 

865 Placement of the colorbar. Combine 'top', 'center' and 'bottom' 

866 with 'left', None for middle and 'right'. 

867 :type anchor: 

868 str 

869 ''' 

870 

871 if not kwargs: 

872 kwargs = {} 

873 

874 if label: 

875 kwargs['B'] = 'af+l%s' % label 

876 

877 kwargs['C'] = make_cpt_path(self.gmt, cmap) 

878 a_str = cbar_anchor[anchor] 

879 

880 w = self.width / 3. 

881 h = w / 10. 

882 

883 lgap = rgap = w / 10. 

884 bgap, tgap = h, h / 10. 

885 

886 dx, dy = 2.5 * lgap, 2. * tgap 

887 

888 if 'bottom' in anchor: 

889 dy += 4 * h 

890 

891 self.gmt.psscale( 

892 D='j%s+w%gc/%gc+h+o%gc/%gc' % (a_str, w, h, dx, dy), 

893 F='+g238/236/230+c%g/%g/%g/%g' % (lgap, rgap, bgap, tgap), 

894 *self.jxyr, 

895 **kwargs) 

896 

897 def draw_vector(self, x_gridfile, y_gridfile, vcolor='', **kwargs): 

898 ''' 

899 Draw vectors based on two grid files. 

900 

901 Two grid files for vector lengths in x and y need to be given. The 

902 function calls gmt.grdvector. All arguments defined for this function 

903 in gmt can be passed as keyword arguments. Different standard settings 

904 are applied if not defined differently. 

905 

906 :param x_gridfile: 

907 File of the grid defining vector lengths in x. 

908 :type x_gridfile: 

909 str 

910 

911 :param y_gridfile: 

912 File of the grid defining vector lengths in y. 

913 :type y_gridfile: 

914 str 

915 

916 :param vcolor: 

917 Vector face color as defined in "G" option. 

918 :type vcolor: 

919 str 

920 ''' 

921 

922 kwargs['S'] = kwargs.get('S', 'il1.') 

923 kwargs['I'] = kwargs.get('I', 'x20') 

924 kwargs['W'] = kwargs.get('W', '0.3p,%s' % 'black') 

925 kwargs['Q'] = kwargs.get('Q', '4c+e+n5km+h1') 

926 

927 self.gmt.grdvector( 

928 x_gridfile, y_gridfile, 

929 G='%s' % 'lightgrey' if not vcolor else vcolor, 

930 *self.jxyr, 

931 **kwargs) 

932 

933 def draw_dynamic_data(self, data, **kwargs): 

934 ''' 

935 Draw an image of any data gridded on the patches e.g dislocation. 

936 

937 :param data: 

938 Patchwise grid data. 

939 :type data: 

940 :py:class:`~numpy.ndarray` 

941 ''' 

942 

943 plot_data = data 

944 

945 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r') 

946 

947 clim = kwargs.pop('clim', (plot_data.min(), plot_data.max())) 

948 

949 cpt_path = make_cpt_path(self.gmt, kwargs['cmap']) 

950 if not op.exists(cpt_path): 

951 make_colormap(self.gmt, clim[0], clim[1], 

952 cmap=kwargs['cmap'], space=False) 

953 

954 tmp_grd_file = self.gmt.tempfilename('tmpdata.grd') 

955 self.patch_data_to_grid(plot_data, tmp_grd_file) 

956 self.draw_image(tmp_grd_file, **kwargs) 

957 

958 def draw_patch_parameter(self, attribute, **kwargs): 

959 ''' 

960 Draw an image of a chosen patch attribute e.g traction. 

961 

962 :param attribute: 

963 Patch attribute, which is plotted. All patch attributes can be 

964 taken (see doc of :py:class:`~pyrocko.modelling.okada.OkadaSource`) 

965 and also ``traction``, ``tx``, ``ty`` or ``tz`` to display the 

966 length or the single components of the traction vector. 

967 :type attribute: 

968 str 

969 ''' 

970 

971 a = attribute 

972 source = self.source 

973 

974 if a == 'traction': 

975 data = num.linalg.norm(source.get_tractions(), axis=1) 

976 elif a == 'tx': 

977 data = source.get_tractions()[:, 0] 

978 elif a == 'ty': 

979 data = source.get_tractions()[:, 1] 

980 elif a == 'tz': 

981 data = source.get_tractions()[:, 2] 

982 else: 

983 data = source.get_patch_attribute(attribute) 

984 

985 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor'] 

986 data *= factor 

987 

988 kwargs['label'] = kwargs.get( 

989 'label', 

990 '%s [%s]' % (a, cbar_helper[a]['unit'])) 

991 

992 self.draw_dynamic_data(data, **kwargs) 

993 

994 def draw_time_contour(self, store, clevel=[], **kwargs): 

995 ''' 

996 Draw high contour lines of the rupture front propgation time. 

997 

998 :param store: 

999 Greens function store, which is used for time calculation. 

1000 :type store: 

1001 :py:class:`~pyrocko.gf.store.Store` 

1002 :param clevel: 

1003 List of times, for which contour lines are drawn. 

1004 :type clevel: 

1005 :py:class:`list` of :py:class:`float` 

1006 ''' 

1007 

1008 _, _, _, _, points_xy = self.source._discretize_points(store, cs='xyz') 

1009 _, _, times, _, _, _ = self.source.get_vr_time_interpolators(store) 

1010 

1011 scaler = AutoScaler(mode='0-max', approx_ticks=8) 

1012 scale = scaler.make_scale([num.min(times), num.max(times)]) 

1013 

1014 if clevel: 

1015 if len(clevel) > 1: 

1016 kwargs['anot_int'] = num.min(num.diff(clevel)) 

1017 else: 

1018 kwargs['anot_int'] = clevel[0] 

1019 

1020 kwargs['contour_int'] = kwargs['anot_int'] 

1021 kwargs['L'] = '0/%g' % num.max(clevel) 

1022 

1023 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.) 

1024 kwargs['contour_int'] = kwargs.get('contour_int', scale[2]) 

1025 kwargs['unit'] = kwargs.get('unit', cbar_helper['time']['unit']) 

1026 kwargs['L'] = kwargs.get('L', '0/%g' % (num.max(times) + 1.)) 

1027 kwargs['G'] = kwargs.get('G', 'n2/3c') 

1028 

1029 tmp_grd_file = self.gmt.tempfilename('tmpdata.grd') 

1030 self.xy_data_to_grid(points_xy[:, 0], points_xy[:, 1], 

1031 times, tmp_grd_file) 

1032 self.draw_contour(tmp_grd_file, **kwargs) 

1033 

1034 def draw_points(self, lats, lons, symbol='point', size=None, **kwargs): 

1035 ''' 

1036 Draw points at given locations. 

1037 

1038 :param lats: 

1039 Point latitude coordinates in [deg]. 

1040 :type lats: 

1041 iterable of :py:class:`float` 

1042 

1043 :param lons: 

1044 Point longitude coordinates in [deg]. 

1045 :type lons: 

1046 iterable of :py:class:`float` 

1047 

1048 :param symbol: 

1049 Define symbol of the points (``'star', 'circle', 'point', 

1050 'triangle'``) - default is ``point``. 

1051 :type symbol: 

1052 str 

1053 

1054 :param size: 

1055 Size of the points in [points]. 

1056 :type size: 

1057 float 

1058 ''' 

1059 

1060 sym_to_gmt = dict( 

1061 star='a', 

1062 circle='c', 

1063 point='p', 

1064 triangle='t') 

1065 

1066 lats = num.atleast_1d(lats) 

1067 lons = num.atleast_1d(lons) 

1068 

1069 if lats.shape[0] != lons.shape[0]: 

1070 raise IndexError('lats and lons do not have the same shape!') 

1071 

1072 if size is None: 

1073 size = self._fontsize / 3. 

1074 

1075 kwargs['S'] = kwargs.get('S', sym_to_gmt[symbol] + '%gp' % size) 

1076 kwargs['G'] = kwargs.get('G', gmtpy.color('scarletred2')) 

1077 kwargs['W'] = kwargs.get('W', '2p,%s' % self._fontcolor) 

1078 

1079 self.gmt.psxy( 

1080 in_columns=[lons, lats], 

1081 *self.jxyr, 

1082 **kwargs) 

1083 

1084 def draw_nucleation_point(self, **kwargs): 

1085 ''' 

1086 Plot the nucleation point onto the map. 

1087 ''' 

1088 

1089 nlat, nlon = xy_to_latlon( 

1090 self.source, self.source.nucleation_x, self.source.nucleation_y) 

1091 

1092 self.draw_points(nlat, nlon, **kwargs) 

1093 

1094 def draw_dislocation(self, time=None, component='', **kwargs): 

1095 ''' 

1096 Draw dislocation onto map at any time. 

1097 

1098 For a given time (if ``time`` is ``None``, ``tmax`` is used) and given 

1099 component the patchwise dislocation is plotted onto the map. 

1100 

1101 :param time: 

1102 Time after origin, for which dislocation is computed. If ``None``, 

1103 ``tmax`` is taken. 

1104 :type time: 

1105 float 

1106 

1107 :param component: 

1108 Dislocation component, which shall be plotted: ``x`` along strike, 

1109 ``y`` along updip, ``z`` normal. If ``None``, the 

1110 length of the dislocation vector is plotted. 

1111 :type component: str 

1112 ''' 

1113 

1114 disl = self.source.get_slip(time=time) 

1115 

1116 if component: 

1117 data = disl[:, c2disl[component]] 

1118 else: 

1119 data = num.linalg.norm(disl, axis=1) 

1120 

1121 kwargs['label'] = kwargs.get( 

1122 'label', 'u%s [m]' % (component)) 

1123 

1124 self.draw_dynamic_data(data, **kwargs) 

1125 

1126 def draw_dislocation_contour( 

1127 self, time=None, component=None, clevel=[], **kwargs): 

1128 ''' 

1129 Draw dislocation contour onto map at any time. 

1130 

1131 For a given time (if ``time`` is ``None``, ``tmax`` is used) and given 

1132 component the patchwise dislocation is plotted as contour onto the map. 

1133 

1134 :param time: 

1135 Time after origin, for which dislocation is computed. If ``None``, 

1136 ``tmax`` is taken. 

1137 :type time: 

1138 float 

1139 

1140 :param component: 

1141 Dislocation component, which shall be plotted: ``x`` along strike, 

1142 ``y`` along updip, ``z`` normal``. If ``None``, the length of the 

1143 dislocation vector is plotted. 

1144 :type component: str 

1145 

1146 :param clevel: 

1147 Times, for which contour lines are drawn. 

1148 :type clevel: 

1149 :py:class:`list` of :py:class:`float` 

1150 ''' 

1151 

1152 disl = self.source.get_slip(time=time) 

1153 

1154 if component: 

1155 data = disl[:, c2disl[component]] 

1156 else: 

1157 data = num.linalg.norm(disl, axis=1) 

1158 

1159 scaler = AutoScaler(mode='min-max', approx_ticks=7) 

1160 scale = scaler.make_scale([num.min(data), num.max(data)]) 

1161 

1162 if clevel: 

1163 if len(clevel) > 1: 

1164 kwargs['anot_int'] = num.min(num.diff(clevel)) 

1165 else: 

1166 kwargs['anot_int'] = clevel[0] 

1167 

1168 kwargs['contour_int'] = kwargs['anot_int'] 

1169 kwargs['L'] = '%g/%g' % ( 

1170 num.min(clevel) - kwargs['contour_int'], 

1171 num.max(clevel) + kwargs['contour_int']) 

1172 else: 

1173 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.) 

1174 kwargs['contour_int'] = kwargs.get('contour_int', scale[2]) 

1175 kwargs['L'] = kwargs.get('L', '%g/%g' % ( 

1176 num.min(data) - 1., num.max(data) + 1.)) 

1177 

1178 kwargs['unit'] = kwargs.get('unit', ' m') 

1179 kwargs['G'] = kwargs.get('G', 'n2/3c') 

1180 

1181 tmp_grd_file = self.gmt.tempfilename('tmpdata.grd') 

1182 self.patch_data_to_grid(data, tmp_grd_file) 

1183 self.draw_contour(tmp_grd_file, **kwargs) 

1184 

1185 def draw_dislocation_vector(self, time=None, **kwargs): 

1186 ''' 

1187 Draw vector arrows onto map indicating direction of dislocation. 

1188 

1189 For a given time (if ``time`` is ``None``, ``tmax`` is used) 

1190 and given component the dislocation is plotted as vectors onto the map. 

1191 

1192 :param time: 

1193 Time after origin [s], for which dislocation is computed. If 

1194 ``None``, ``tmax`` is used. 

1195 :type time: 

1196 float 

1197 ''' 

1198 

1199 disl = self.source.get_slip(time=time) 

1200 

1201 p_strike = self.source.get_patch_attribute('strike') * d2r 

1202 p_dip = self.source.get_patch_attribute('dip') * d2r 

1203 

1204 disl_dh = num.cos(p_dip) * disl[:, 1] 

1205 disl_n = num.cos(p_strike) * disl[:, 0] + num.sin(p_strike) * disl_dh 

1206 disl_e = num.sin(p_strike) * disl[:, 0] - num.cos(p_strike) * disl_dh 

1207 

1208 tmp_grd_files = ['tmpdata_%s.grd' % c for c in ('n', 'e')] 

1209 

1210 self.patch_data_to_grid(disl_n, tmp_grd_files[0]) 

1211 self.patch_data_to_grid(disl_e, tmp_grd_files[1]) 

1212 

1213 self.draw_vector( 

1214 tmp_grd_files[1], tmp_grd_files[0], 

1215 **kwargs) 

1216 

1217 def draw_top_edge(self, **kwargs): 

1218 ''' 

1219 Indicate rupture top edge on map. 

1220 ''' 

1221 

1222 outline = self.source.outline(cs='latlondepth') 

1223 top_edge = outline[:2, :] 

1224 

1225 kwargs = kwargs or {} 

1226 kwargs['W'] = kwargs.get( 

1227 'W', '%gp,%s' % (self._fontsize / 10., gmtpy.color('orange3'))) 

1228 

1229 self.gmt.psxy( 

1230 in_columns=[top_edge[:, 1], top_edge[:, 0]], 

1231 *self.jxyr, 

1232 **kwargs) 

1233 

1234 

1235class RuptureView(Object): 

1236 ''' 

1237 Plot of attributes and results of the 

1238 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`. 

1239 ''' 

1240 

1241 _patches_to_lw = RuptureMap._patches_to_lw 

1242 

1243 def __init__(self, source=None, figsize=None, fontsize=None): 

1244 self._source = source 

1245 self._axes = None 

1246 

1247 self.figsize = figsize or mpl_papersize('halfletter', 'landscape') 

1248 self.fontsize = fontsize or 10 

1249 mpl_init(fontsize=self.fontsize) 

1250 

1251 self._fig = None 

1252 self._axes = None 

1253 self._is_1d = False 

1254 

1255 @property 

1256 def source(self): 

1257 ''' 

1258 PseudoDynamicRupture whose attributes are plotted. 

1259 

1260 Note, that source.patches attribute needs to be calculated for 

1261 :type source: :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`. 

1262 ''' 

1263 

1264 if self._source is None: 

1265 raise SourceError('No source given. Please define it!') 

1266 

1267 if not isinstance(self._source, PseudoDynamicRupture): 

1268 raise SourceError('This function is only capable for a source of' 

1269 ' type: %s' % type(PseudoDynamicRupture())) 

1270 

1271 if not self._source.patches: 

1272 raise TypeError('No source patches are defined. Please run' 

1273 '\"discretize_patches\" on your source') 

1274 

1275 return self._source 

1276 

1277 @source.setter 

1278 def source(self, source): 

1279 self._source = source 

1280 

1281 def _setup( 

1282 self, 

1283 title='', 

1284 xlabel='', 

1285 ylabel='', 

1286 aspect=1., 

1287 spatial_plot=True, 

1288 **kwargs): 

1289 

1290 if self._fig is not None and self._axes is not None: 

1291 return 

1292 

1293 self._fig = plt.figure(figsize=self.figsize) 

1294 

1295 self._axes = self._fig.add_subplot(1, 1, 1, aspect=aspect) 

1296 ax = self._axes 

1297 

1298 if ax is not None: 

1299 ax.set_title(title) 

1300 ax.grid(alpha=.3) 

1301 ax.set_xlabel(xlabel) 

1302 ax.set_ylabel(ylabel) 

1303 

1304 if spatial_plot: 

1305 ax.xaxis.set_major_formatter(FuncFormatter(km_formatter)) 

1306 ax.yaxis.set_major_formatter(FuncFormatter(km_formatter)) 

1307 

1308 def _clear_all(self): 

1309 plt.cla() 

1310 plt.clf() 

1311 plt.close() 

1312 

1313 self._fig, self._axes = None, None 

1314 

1315 def _draw_scatter(self, x, y, *args, **kwargs): 

1316 default_kwargs = dict( 

1317 linewidth=0, 

1318 marker='o', 

1319 markerfacecolor=mpl_color('skyblue2'), 

1320 markersize=6., 

1321 markeredgecolor=mpl_color('skyblue3')) 

1322 default_kwargs.update(kwargs) 

1323 

1324 if self._axes is not None: 

1325 self._axes.plot(x, y, *args, **default_kwargs) 

1326 

1327 def _draw_image(self, length, width, data, *args, **kwargs): 

1328 if self._axes is not None: 

1329 if 'extent' not in kwargs: 

1330 kwargs['extent'] = [ 

1331 num.min(length), num.max(length), 

1332 num.max(width), num.min(width)] 

1333 

1334 im = self._axes.imshow( 

1335 data, 

1336 *args, 

1337 interpolation='none', 

1338 vmin=kwargs.get('clim', [None])[0], 

1339 vmax=kwargs.get('clim', [None, None])[1], 

1340 **kwargs) 

1341 

1342 del kwargs['extent'] 

1343 if 'aspect' in kwargs: 

1344 del kwargs['aspect'] 

1345 if 'clim' in kwargs: 

1346 del kwargs['clim'] 

1347 if 'cmap' in kwargs: 

1348 del kwargs['cmap'] 

1349 

1350 plt.colorbar( 

1351 im, *args, shrink=0.9, pad=0.03, aspect=15., **kwargs) 

1352 

1353 def _draw_contour(self, x, y, data, clevel=None, unit='', *args, **kwargs): 

1354 setup_kwargs = dict( 

1355 xlabel=kwargs.pop('xlabel', 'along strike [km]'), 

1356 ylabel=kwargs.pop('ylabel', 'down dip [km]'), 

1357 title=kwargs.pop('title', ''), 

1358 aspect=kwargs.pop('aspect', 1)) 

1359 

1360 self._setup(**setup_kwargs) 

1361 

1362 if self._axes is not None: 

1363 if clevel is None: 

1364 scaler = AutoScaler(mode='min-max') 

1365 scale = scaler.make_scale([data.min(), data.max()]) 

1366 

1367 clevel = num.arange(scale[0], scale[1] + scale[2], scale[2]) 

1368 

1369 if not isinstance(clevel, num.ndarray): 

1370 clevel = num.array([clevel]) 

1371 

1372 clevel = clevel[clevel < data.max()] 

1373 

1374 cont = self._axes.contour( 

1375 x, y, data, clevel, *args, 

1376 linewidths=1.5, **kwargs) 

1377 

1378 clabels = self._axes.clabel( 

1379 cont, clevel[::2], *args, 

1380 inline=1, fmt='%g' + '%s' % unit, 

1381 inline_spacing=15, rightside_up=True, use_clabeltext=True, 

1382 **kwargs) 

1383 

1384 plt.setp(clabels, path_effects=[ 

1385 patheffects.withStroke(linewidth=1.25, foreground='beige'), 

1386 patheffects.Normal()]) 

1387 

1388 def draw_points(self, length, width, *args, **kwargs): 

1389 ''' 

1390 Draw a point onto the figure. 

1391 

1392 Args and kwargs can be defined according to 

1393 :py:func:`matplotlib.pyplot.scatter`. 

1394 

1395 :param length: 

1396 Point(s) coordinate on the rupture plane along strike relative to 

1397 the anchor point in [m]. 

1398 :type length: 

1399 float, :py:class:`~numpy.ndarray` 

1400 

1401 :param width: 

1402 Point(s) coordinate on the rupture plane along downdip relative to 

1403 the anchor point in [m]. 

1404 :type width: 

1405 float, :py:class:`~numpy.ndarray` 

1406 ''' 

1407 

1408 if self._axes is not None: 

1409 kwargs['color'] = kwargs.get('color', mpl_color('scarletred2')) 

1410 kwargs['s'] = kwargs.get('s', 100.) 

1411 

1412 self._axes.scatter(length, width, *args, **kwargs) 

1413 

1414 def draw_dynamic_data(self, data, **kwargs): 

1415 ''' 

1416 Draw an image of any data gridded on the patches e.g dislocation. 

1417 

1418 :param data: 

1419 Patchwise grid data. 

1420 :type data: 

1421 :py:class:`~numpy.ndarray` 

1422 ''' 

1423 

1424 plot_data = data 

1425 

1426 anchor_x, anchor_y = map_anchor[self.source.anchor] 

1427 

1428 length, width = xy_to_lw( 

1429 self.source, num.array([-1., 1.]), num.array([-1., 1.])) 

1430 

1431 data = plot_data.reshape(self.source.nx, self.source.ny).T 

1432 

1433 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r') 

1434 

1435 setup_kwargs = dict( 

1436 xlabel='along strike [km]', 

1437 ylabel='down dip [km]', 

1438 title='', 

1439 aspect=1) 

1440 setup_kwargs.update(kwargs) 

1441 

1442 kwargs = {k: v for k, v in kwargs.items() if k not in 

1443 ('xlabel', 'ylabel', 'title')} 

1444 self._setup(**setup_kwargs) 

1445 self._draw_image(length=length, width=width, data=data, **kwargs) 

1446 

1447 def draw_patch_parameter(self, attribute, **kwargs): 

1448 ''' 

1449 Draw an image of a chosen patch attribute e.g traction. 

1450 

1451 :param attribute: 

1452 Patch attribute, which is plotted. All patch attributes can be 

1453 taken (see doc of :py:class:`~pyrocko.modelling.okada.OkadaSource`) 

1454 and also ``'traction', 'tx', 'ty', 'tz'`` to display the length 

1455 or the single components of the traction vector. 

1456 :type attribute: 

1457 str 

1458 ''' 

1459 

1460 a = attribute 

1461 source = self.source 

1462 

1463 if a == 'traction': 

1464 data = num.linalg.norm(source.get_tractions(), axis=1) 

1465 elif a == 'tx': 

1466 data = source.get_tractions()[:, 0] 

1467 elif a == 'ty': 

1468 data = source.get_tractions()[:, 1] 

1469 elif a == 'tz': 

1470 data = source.get_tractions()[:, 2] 

1471 else: 

1472 data = source.get_patch_attribute(attribute) 

1473 

1474 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor'] 

1475 data *= factor 

1476 

1477 kwargs['label'] = kwargs.get( 

1478 'label', 

1479 '%s [%s]' % (a, cbar_helper[a]['unit'])) 

1480 

1481 return self.draw_dynamic_data(data=data, **kwargs) 

1482 

1483 def draw_time_contour(self, store, clevel=[], **kwargs): 

1484 ''' 

1485 Draw high resolution contours of the rupture front propgation time 

1486 

1487 :param store: 

1488 Greens function store, which is used for time calculation. 

1489 :type store: 

1490 :py:class:`~pyrocko.gf.store.Store` 

1491 

1492 :param clevel: 

1493 Levels of the contour lines. If no levels are given, they are 

1494 automatically computed based on ``tmin`` and ``tmax``. 

1495 :type clevel: 

1496 list 

1497 ''' 

1498 

1499 source = self.source 

1500 default_kwargs = dict( 

1501 colors='#474747' 

1502 ) 

1503 default_kwargs.update(kwargs) 

1504 

1505 _, _, _, _, points_xy = source._discretize_points(store, cs='xyz') 

1506 _, _, _, _, time_interpolator, _ = source.get_vr_time_interpolators( 

1507 store) 

1508 

1509 times = time_interpolator.values 

1510 

1511 scaler = AutoScaler(mode='min-max', approx_ticks=8) 

1512 scale = scaler.make_scale([times.min(), times.max()]) 

1513 

1514 if len(clevel) == 0: 

1515 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2]) 

1516 

1517 points_x = time_interpolator.grid[0] 

1518 points_y = time_interpolator.grid[1] 

1519 

1520 self._draw_contour(points_x, points_y, data=times.T, 

1521 clevel=clevel, unit='s', **default_kwargs) 

1522 

1523 def draw_nucleation_point(self, **kwargs): 

1524 ''' 

1525 Draw the nucleation point onto the map. 

1526 ''' 

1527 

1528 nuc_x, nuc_y = self.source.nucleation_x, self.source.nucleation_y 

1529 length, width = xy_to_lw(self.source, nuc_x, nuc_y) 

1530 self.draw_points(length, width, marker='o', **kwargs) 

1531 

1532 def draw_dislocation(self, time=None, component='', **kwargs): 

1533 ''' 

1534 Draw dislocation onto map at any time. 

1535 

1536 For a given time (if ``time`` is ``None``, ``tmax`` is used) 

1537 and given component the patchwise dislocation is plotted onto the map. 

1538 

1539 :param time: 

1540 Time after origin [s], for which dislocation is computed. 

1541 If ``None``, ``tmax`` is taken. 

1542 :type time: 

1543 float 

1544 

1545 :param component: 

1546 Dislocation component, which shall be plotted: ``x`` 

1547 along strike, ``y`` along updip, ``z`` normal. If ``None``, the 

1548 length of the dislocation vector is plotted 

1549 :type component: str 

1550 ''' 

1551 

1552 disl = self.source.get_slip(time=time) 

1553 

1554 if component: 

1555 data = disl[:, c2disl[component]] 

1556 else: 

1557 data = num.linalg.norm(disl, axis=1) 

1558 

1559 kwargs['label'] = kwargs.get( 

1560 'label', 'u%s [m]' % (component)) 

1561 

1562 self.draw_dynamic_data(data, **kwargs) 

1563 

1564 def draw_dislocation_contour( 

1565 self, time=None, component=None, clevel=[], **kwargs): 

1566 ''' 

1567 Draw dislocation contour onto map at any time. 

1568 

1569 For a given time (if time is ``None``, ``tmax`` is used) and given 

1570 component the patchwise dislocation is plotted as contour onto the map. 

1571 

1572 :param time: 

1573 Time after origin, for which dislocation is computed. If 

1574 ``None``, ``tmax`` is taken. 

1575 :type time: 

1576 float 

1577 

1578 :param component: 

1579 Dislocation component, which shall be plotted. ``x`` 

1580 along strike, ``y`` along updip, ``z`` - normal. If ``None`` 

1581 is given, the length of the dislocation vector is plotted. 

1582 :type component: 

1583 str 

1584 ''' 

1585 

1586 disl = self.source.get_slip(time=time) 

1587 

1588 if component: 

1589 data = disl[:, c2disl[component]] 

1590 else: 

1591 data = num.linalg.norm(disl, axis=1) 

1592 

1593 data = data.reshape(self.source.ny, self.source.nx, order='F') 

1594 

1595 scaler = AutoScaler(mode='min-max', approx_ticks=7) 

1596 scale = scaler.make_scale([data.min(), data.max()]) 

1597 

1598 if len(clevel) == 0: 

1599 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2]) 

1600 

1601 anchor_x, anchor_y = map_anchor[self.source.anchor] 

1602 

1603 length, width = self._patches_to_lw() 

1604 length, width = length[1:-1], width[1:-1] 

1605 

1606 kwargs['colors'] = kwargs.get('colors', '#474747') 

1607 

1608 self._setup(**kwargs) 

1609 self._draw_contour( 

1610 length, width, data=data, clevel=clevel, unit='m', **kwargs) 

1611 

1612 def draw_source_dynamics( 

1613 self, variable, store, deltat=None, *args, **kwargs): 

1614 ''' 

1615 Display dynamic source parameter. 

1616 

1617 Fast inspection possibility for the cumulative moment and the source 

1618 time function approximation (assuming equal paths between different 

1619 patches and observation point - valid for an observation point in the 

1620 far field perpendicular to the source strike), so the cumulative moment 

1621 rate function. 

1622 

1623 :param variable: 

1624 Dynamic parameter, which shall be plotted. Choose 

1625 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment') 

1626 :type variable: 

1627 str 

1628 

1629 :param store: 

1630 Greens function store, whose store.config.deltat defines 

1631 the time increment between two parameter snapshots. If ``store`` is 

1632 not given, the time increment is defined is taken from ``deltat``. 

1633 :type store: 

1634 :py:class:`~pyrocko.gf.store.Store` 

1635 

1636 :param deltat: 

1637 Time increment between two parameter snapshots. If not 

1638 given, store.config.deltat is used to define ``deltat``. 

1639 :type deltat: 

1640 float 

1641 ''' 

1642 

1643 v = variable 

1644 

1645 data, times = self.source.get_moment_rate(store=store, deltat=deltat) 

1646 deltat = num.concatenate([(num.diff(times)[0], ), num.diff(times)]) 

1647 

1648 if v in ('moment_rate', 'stf'): 

1649 name, unit = 'dM/dt', 'Nm/s' 

1650 elif v in ('cumulative_moment', 'moment'): 

1651 data = num.cumsum(data * deltat) 

1652 name, unit = 'M', 'Nm' 

1653 else: 

1654 raise ValueError('No dynamic data for given variable %s found' % v) 

1655 

1656 self._setup(xlabel='time [s]', 

1657 ylabel='%s / %.2g %s' % (name, data.max(), unit), 

1658 aspect='auto', 

1659 spatial_plot=False) 

1660 self._draw_scatter(times, data/num.max(data), *args, **kwargs) 

1661 self._is_1d = True 

1662 

1663 def draw_patch_dynamics( 

1664 self, variable, nx, ny, store=None, deltat=None, *args, **kwargs): 

1665 ''' 

1666 Display dynamic boundary element / patch parameter. 

1667 

1668 Fast inspection possibility for different dynamic parameter for a 

1669 single patch / boundary element. The chosen parameter is plotted for 

1670 the chosen patch. 

1671 

1672 :param variable: 

1673 Dynamic parameter, which shall be plotted. Choose 

1674 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment'). 

1675 :type variable: 

1676 str 

1677 

1678 :param nx: 

1679 Patch index along strike (range: 0:source.nx - 1). 

1680 :type nx: 

1681 int 

1682 

1683 :param nx: 

1684 Patch index downdip (range: 0:source.ny - 1). 

1685 :type nx: 

1686 int 

1687 

1688 :param store: 

1689 Greens function store, whose store.config.deltat defines 

1690 the time increment between two parameter snapshots. If ``store`` is 

1691 not given, the time increment is defined is taken from ``deltat``. 

1692 :type store: 

1693 :py:class:`~pyrocko.gf.store.Store` 

1694 

1695 :param deltat: 

1696 Time increment between two parameter snapshots. If not given, 

1697 store.config.deltat is used to define ``deltat``. 

1698 :type deltat: 

1699 float 

1700 ''' 

1701 

1702 v = variable 

1703 source = self.source 

1704 idx = nx * source.ny + nx 

1705 

1706 m = re.match(r'dislocation_([xyz])', v) 

1707 

1708 if v in ('moment_rate', 'cumulative_moment', 'moment', 'stf'): 

1709 data, times = source.get_moment_rate_patches( 

1710 store=store, deltat=deltat) 

1711 elif 'dislocation' in v or 'slip_rate' == v: 

1712 ddisloc, times = source.get_delta_slip(store=store, deltat=deltat) 

1713 

1714 if v in ('moment_rate', 'stf'): 

1715 data, times = source.get_moment_rate_patches( 

1716 store=store, deltat=deltat) 

1717 name, unit = 'dM/dt', 'Nm/s' 

1718 elif v in ('cumulative_moment', 'moment'): 

1719 data, times = source.get_moment_rate_patches( 

1720 store=store, deltat=deltat) 

1721 data = num.cumsum(data, axis=1) 

1722 name, unit = 'M', 'Nm' 

1723 elif v == 'slip_rate': 

1724 data, times = source.get_slip_rate(store=store, deltat=deltat) 

1725 name, unit = 'du/dt', 'm/s' 

1726 elif v == 'dislocation': 

1727 data, times = source.get_delta_slip(store=store, deltat=deltat) 

1728 data = num.linalg.norm(num.cumsum(data, axis=2), axis=1) 

1729 name, unit = 'du', 'm' 

1730 elif m: 

1731 data, times = source.get_delta_slip(store=store, deltat=deltat) 

1732 data = num.cumsum(data, axis=2)[:, c2disl[m.group(1)], :] 

1733 name, unit = 'du%s' % m.group(1), 'm' 

1734 else: 

1735 raise ValueError('No dynamic data for given variable %s found' % v) 

1736 

1737 self._setup(xlabel='time [s]', 

1738 ylabel='%s / %.2g %s' % (name, num.max(data), unit), 

1739 aspect='auto', 

1740 spatial_plot=False) 

1741 self._draw_scatter(times, data[idx, :]/num.max(data), 

1742 *args, **kwargs) 

1743 self._is_1d = True 

1744 

1745 def finalize(self): 

1746 if self._is_1d: 

1747 return 

1748 

1749 length, width = xy_to_lw( 

1750 self.source, num.array([-1., 1.]), num.array([1., -1.])) 

1751 

1752 self._axes.set_xlim(length) 

1753 self._axes.set_ylim(width) 

1754 

1755 def gcf(self): 

1756 self.finalize() 

1757 return self._fig 

1758 

1759 def save(self, filename, dpi=None): 

1760 ''' 

1761 Save plot to file. 

1762 

1763 :param filename: 

1764 Filename and path, where the plot is stored. 

1765 :type filename: 

1766 str 

1767 

1768 :param dpi: 

1769 Resolution of the output plot in [dpi]. 

1770 :type dpi: 

1771 int 

1772 ''' 

1773 

1774 self.finalize() 

1775 try: 

1776 self._fig.savefig(fname=filename, dpi=dpi, bbox_inches='tight') 

1777 except TypeError: 

1778 self._fig.savefig(filename=filename, dpi=dpi, bbox_inches='tight') 

1779 

1780 self._clear_all() 

1781 

1782 def show_plot(self): 

1783 ''' 

1784 Show plot. 

1785 ''' 

1786 

1787 self.finalize() 

1788 plt.show() 

1789 self._clear_all() 

1790 

1791 

1792def render_movie(fn_path, output_path, framerate=20): 

1793 ''' 

1794 Generate a mp4 movie based on given png files using 

1795 `ffmpeg <https://ffmpeg.org>`_. 

1796 

1797 Render a movie based on a set of given .png files in fn_path. All files 

1798 must have a filename specified by ``fn_path`` (e.g. giving ``fn_path`` with 

1799 ``/temp/f%04.png`` a valid png filename would be ``/temp/f0001.png``). The 

1800 files must have a numbering, indicating their order in the movie. 

1801 

1802 :param fn_path: 

1803 Path and fileformat specification of the input .png files. 

1804 :type fn_path: 

1805 str 

1806 

1807 :param output_path: 

1808 Path and filename of the output ``.mp4`` movie file. 

1809 :type output_path: 

1810 str 

1811 

1812 :param deltat: 

1813 Time between individual frames (``1 / framerate``) in [s]. 

1814 :type deltat: 

1815 float 

1816 ''' 

1817 

1818 try: 

1819 check_call(['ffmpeg', '-loglevel', 'panic']) 

1820 except CalledProcessError: 

1821 pass 

1822 except (TypeError, IOError): 

1823 logger.warning( 

1824 'Package ffmpeg needed for movie rendering. Please install it ' 

1825 '(e.g. on linux distr. via sudo apt-get ffmpeg) and retry.') 

1826 return False 

1827 

1828 try: 

1829 check_call([ 

1830 'ffmpeg', '-loglevel', 'info', '-y', 

1831 '-framerate', '%g' % framerate, 

1832 '-i', fn_path, 

1833 '-vcodec', 'libx264', 

1834 '-preset', 'medium', 

1835 '-tune', 'animation', 

1836 '-pix_fmt', 'yuv420p', 

1837 '-movflags', '+faststart', 

1838 '-filter:v', "crop='floor(in_w/2)*2:floor(in_h/2)*2'", 

1839 '-crf', '15', 

1840 output_path]) 

1841 

1842 return True 

1843 except CalledProcessError as e: 

1844 logger.warning(e) 

1845 return False 

1846 

1847 

1848def render_gif(fn, output_path, loops=-1): 

1849 ''' 

1850 Generate a gif based on a given mp4 using ffmpeg. 

1851 

1852 Render a gif based on a given .mp4 movie file in ``fn`` path. 

1853 

1854 :param fn: 

1855 Path and file name of the input .mp4 file. 

1856 :type fn: 

1857 str 

1858 

1859 :param output_path: 

1860 Path and filename of the output animated gif file. 

1861 :type output_path: 

1862 str 

1863 :param loops: 

1864 Number of gif repeats (loops). ``-1`` is not repetition, ``0`` 

1865 infinite. 

1866 :type loops: 

1867 int 

1868 ''' 

1869 

1870 try: 

1871 check_call(['ffmpeg', '-loglevel', 'panic']) 

1872 except CalledProcessError: 

1873 pass 

1874 except (TypeError, IOError): 

1875 logger.warning( 

1876 'Package ffmpeg needed for movie rendering. Please install it ' 

1877 '(e.g. on linux distr. via sudo apt-get ffmpeg.) and retry.') 

1878 return False 

1879 

1880 try: 

1881 check_call([ 

1882 'ffmpeg', '-hide_banner', '-loglevel', 'panic', '-y', '-i', 

1883 fn, 

1884 '-filter_complex', 'palettegen[v1];[0:v][v1]paletteuse', 

1885 '-loop', '%d' % loops, 

1886 output_path]) 

1887 

1888 return True 

1889 except CalledProcessError as e: 

1890 logger.warning(e) 

1891 return False 

1892 

1893 

1894def rupture_movie( 

1895 source, 

1896 store, 

1897 variable='dislocation', 

1898 draw_time_contours=False, 

1899 fn_path='.', 

1900 prefix='', 

1901 plot_type='map', 

1902 deltat=None, 

1903 framerate=None, 

1904 store_images=False, 

1905 render_as_gif=False, 

1906 gif_loops=-1, 

1907 **kwargs): 

1908 ''' 

1909 Generate a movie based on a given source for dynamic parameter. 

1910 

1911 Create a MPEG-4 movie or gif of one of the following dynamic parameters 

1912 (``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y`` 

1913 (along updip), ``dislocation_z`` (normal), ``slip_rate``, ``moment_rate``). 

1914 If desired, the single snap shots can be stored as images as well. 

1915 ``kwargs`` have to be given according to the chosen ``plot_type``. 

1916 

1917 :param source: 

1918 Pseudo dynamic rupture, for which the movie is produced. 

1919 :type source: 

1920 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture` 

1921 

1922 :param store: 

1923 Greens function store, which is used for time calculation. If 

1924 ``deltat`` is not given, it is taken from the store.config.deltat 

1925 :type store: 

1926 :py:class:`~pyrocko.gf.store.Store` 

1927 

1928 :param variable: 

1929 Dynamic parameter, which shall be plotted. Choose between 

1930 ``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y`` 

1931 (along updip), ``dislocation_z`` (normal), ``slip_rate`` and 

1932 ``moment_rate``, default ``dislocation``. 

1933 :type variable: 

1934 str 

1935 

1936 :param draw_time_contours: 

1937 If ``True``, corresponding isochrones are drawn on the each plots. 

1938 :type draw_time_contours: 

1939 bool 

1940 

1941 :param fn_path: 

1942 Absolut or relative path, where movie (and optional images) are stored. 

1943 :type fn_path: 

1944 str 

1945 

1946 :param prefix: 

1947 File prefix used for the movie (and optional image) files. 

1948 :type prefix: 

1949 str 

1950 

1951 :param plot_type: 

1952 Choice of plot type: ``map``, ``view`` (map plot using 

1953 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureMap` 

1954 or plane view using 

1955 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureView`). 

1956 :type plot_type: 

1957 str 

1958 

1959 :param deltat: 

1960 Time between parameter snapshots. If not given, store.config.deltat is 

1961 used to define ``deltat``. 

1962 :type deltat: 

1963 float 

1964 

1965 :param store_images: 

1966 Choice to store the single .png parameter snapshots in ``fn_path`` or 

1967 not. 

1968 :type store_images: 

1969 bool 

1970 

1971 :param render_as_gif: 

1972 If ``True``, the movie is converted into a gif. If ``False``, the movie 

1973 is returned as mp4. 

1974 :type render_as_gif: 

1975 bool 

1976 

1977 :param gif_loops: 

1978 If ``render_as_gif`` is ``True``, a gif with ``gif_loops`` number of 

1979 loops (repetitions) is returned. ``-1`` is no repetition, ``0`` 

1980 infinite. 

1981 :type gif_loops: 

1982 int 

1983 ''' 

1984 

1985 v = variable 

1986 assert plot_type in ('map', 'view') 

1987 

1988 if not source.patches: 

1989 source.discretize_patches(store, interpolation='multilinear') 

1990 

1991 if source.coef_mat is None: 

1992 source.calc_coef_mat() 

1993 

1994 prefix = prefix or v 

1995 deltat = deltat or store.config.deltat 

1996 framerate = max(framerate or int(1./deltat), 1) 

1997 

1998 if v == 'moment_rate': 

1999 data, times = source.get_moment_rate_patches(deltat=deltat) 

2000 name, unit = 'dM/dt', 'Nm/s' 

2001 elif 'dislocation' in v or 'slip_rate' == v: 

2002 ddisloc, times = source.get_delta_slip(deltat=deltat) 

2003 else: 

2004 raise ValueError('No dynamic data for given variable %s found' % v) 

2005 

2006 deltat = times[1] - times[0] 

2007 

2008 m = re.match(r'dislocation_([xyz])', v) 

2009 if m: 

2010 data = num.cumsum(ddisloc, axis=1)[:, :, c2disl[m.group(1)]] 

2011 name, unit = 'du%s' % m.group(1), 'm' 

2012 elif v == 'dislocation': 

2013 data = num.linalg.norm(num.cumsum(ddisloc, axis=1), axis=2) 

2014 name, unit = 'du', 'm' 

2015 elif v == 'slip_rate': 

2016 data = num.linalg.norm(ddisloc, axis=2) / deltat 

2017 name, unit = 'du/dt', 'm/s' 

2018 

2019 if plot_type == 'map': 

2020 plt_base = RuptureMap 

2021 elif plot_type == 'view': 

2022 plt_base = RuptureView 

2023 else: 

2024 raise AttributeError('invalid type: %s' % plot_type) 

2025 

2026 attrs_base = get_kwargs(plt_base) 

2027 kwargs_base = dict([k for k in kwargs.items() if k[0] in attrs_base]) 

2028 kwargs_plt = dict([k for k in kwargs.items() if k[0] not in kwargs_base]) 

2029 

2030 if 'clim' in kwargs_plt: 

2031 data = num.clip(data, kwargs_plt['clim'][0], kwargs_plt['clim'][1]) 

2032 else: 

2033 kwargs_plt['clim'] = [num.min(data), num.max(data)] 

2034 

2035 if 'label' not in kwargs_plt: 

2036 vmax = num.max(num.abs(kwargs_plt['clim'])) 

2037 data /= vmax 

2038 

2039 kwargs_plt['label'] = '%s / %.2g %s' % (name, vmax, unit) 

2040 kwargs_plt['clim'] = [i / vmax for i in kwargs_plt['clim']] 

2041 

2042 with tempfile.TemporaryDirectory(suffix='pyrocko') as temp_path: 

2043 fns_temp = [op.join(temp_path, 'f%09d.png' % (it + 1)) 

2044 for it, _ in enumerate(times)] 

2045 fn_temp_path = op.join(temp_path, 'f%09d.png') 

2046 

2047 for it, (t, ft) in enumerate(zip(times, fns_temp)): 

2048 plot_data = data[:, it] 

2049 

2050 plt = plt_base(source=source, **kwargs_base) 

2051 plt.draw_dynamic_data(plot_data, **kwargs_plt) 

2052 plt.draw_nucleation_point() 

2053 

2054 if draw_time_contours: 

2055 plt.draw_time_contour(store, clevel=[t]) 

2056 

2057 plt.save(ft) 

2058 

2059 fn_mp4 = op.join(temp_path, 'movie.mp4') 

2060 return_code = render_movie( 

2061 fn_temp_path, 

2062 output_path=fn_mp4, 

2063 framerate=framerate) 

2064 

2065 if render_as_gif and return_code: 

2066 render_gif(fn=fn_mp4, output_path=op.join( 

2067 fn_path, '%s_%s_gif.gif' % (prefix, plot_type)), 

2068 loops=gif_loops) 

2069 

2070 elif return_code: 

2071 shutil.move(fn_mp4, op.join( 

2072 fn_path, '%s_%s_movie.mp4' % (prefix, plot_type))) 

2073 

2074 else: 

2075 logger.error('ffmpeg failed. Exit') 

2076 

2077 if store_images: 

2078 fns = [op.join(fn_path, '%s_%s_%g.png' % (prefix, plot_type, t)) 

2079 for t in times] 

2080 

2081 for ft, f in zip(fns_temp, fns): 

2082 shutil.move(ft, f) 

2083 

2084 

2085__all__ = [ 

2086 'make_colormap', 

2087 'xy_to_latlon', 

2088 'xy_to_lw', 

2089 'SourceError', 

2090 'RuptureMap', 

2091 'RuptureView', 

2092 'rupture_movie', 

2093 'render_movie']