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

648 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 06:59 +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 

9import os.path as op 

10import re 

11import logging 

12import tempfile 

13import shutil 

14import inspect 

15 

16import numpy as num 

17from scipy.interpolate import RegularGridInterpolator as scrgi 

18 

19from matplotlib import pyplot as plt, patheffects 

20from matplotlib.ticker import FuncFormatter 

21 

22from pyrocko import orthodrome as pod 

23from pyrocko.guts import Object 

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

25 gmtpy, mpl_get_cmap 

26from pyrocko.plot.automap import Map, NoTopo 

27from pyrocko.gf import PseudoDynamicRupture 

28from pyrocko.gf.seismosizer import map_anchor 

29from pyrocko.dataset.topo.tile import Tile 

30 

31logger = logging.getLogger(__name__) 

32 

33km = 1e3 

34 

35d2m = 111180. 

36m2d = 1. / d2m 

37 

38cm2inch = gmtpy.cm/gmtpy.inch 

39 

40d2r = num.pi / 180. 

41r2d = 1. / d2r 

42 

43 

44def km_formatter(v, p): 

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

46 

47 

48def get_kwargs(cls): 

49 sig = inspect.signature(cls.__init__) 

50 kwargs = {} 

51 

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

53 for p in Map.T.properties: 

54 kwargs[p.name] = p._default 

55 

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

57 if par.default is inspect._empty: 

58 continue 

59 kwargs[par.name] = par.default 

60 

61 return kwargs 

62 

63 

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

65 ''' 

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

67 

68 :param lats: 

69 Grid latitude coordinates in [deg]. 

70 :type lats: 

71 iterable 

72 

73 :param lons: 

74 Grid longitude coordinates in [deg]. 

75 :type lons: 

76 iterable 

77 

78 :param data: 

79 Grid data of any kind. 

80 :type data: 

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

82 

83 :param filename: 

84 Filename of the written grid file. 

85 :type filename: 

86 str 

87 ''' 

88 

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

90 

91 

92def _mplcmap_to_gmtcpt_code(mplcmap, steps=256): 

93 ''' 

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

95 

96 :param mplcmap: 

97 Name of the demanded matplotlib colormap. 

98 :type mplcmap: 

99 str 

100 

101 :returns: 

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

103 :rtype: 

104 str 

105 ''' 

106 

107 cmap = mpl_get_cmap(mplcmap) 

108 

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

110 

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

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

113 

114 

115def make_colormap( 

116 gmt, 

117 vmin, 

118 vmax, 

119 C=None, 

120 cmap=None, 

121 space=False): 

122 ''' 

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

124 

125 :param vmin: 

126 Minimum value covered by the colormap. 

127 :type vmin: 

128 float 

129 

130 :param vmax: 

131 Maximum value covered by the colormap. 

132 :type vmax: 

133 float 

134 

135 :param C: 

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

137 :type C: 

138 str 

139 

140 :param cmap: 

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

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

143 extracted from this colormap. 

144 :type cmap: 

145 str 

146 

147 :param space: 

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

149 and above vmax. 

150 :type space: bool 

151 ''' 

152 

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

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

155 

156 incr = scale[2] 

157 margin = 0. 

158 

159 if vmin == vmax: 

160 space = True 

161 

162 if space: 

163 margin = incr 

164 

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

166 ' valid matplotlib colormap name.') 

167 

168 if C is None and cmap is None: 

169 raise ValueError(msg) 

170 

171 if C is None: 

172 try: 

173 C = _mplcmap_to_gmtcpt_code(cmap) 

174 except ValueError: 

175 raise ValueError(msg) 

176 

177 if cmap is None: 

178 logger.warning( 

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

180 cmap = 'temp_cmap' 

181 

182 return gmt.makecpt( 

183 C=C, 

184 D='o', 

185 T='%g/%g/%g' % ( 

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

187 Z=True, 

188 out_filename='my_%s.cpt' % cmap, 

189 suppress_defaults=True) 

190 

191 

192def clear_temp(gridfiles=[], cpts=[]): 

193 ''' 

194 Clear all temporary needed grid and colormap cpt files. 

195 

196 :param gridfiles: 

197 List of all "...grd" files, which shall be deleted. 

198 :type gridfiles: 

199 list 

200 

201 :param cpts: 

202 Cmaps, whose corresponding "my_<cmap>.cpt" file shall be deleted. 

203 :type cpts: 

204 list 

205 ''' 

206 

207 for fil in gridfiles: 

208 try: 

209 os.remove(fil) 

210 except OSError: 

211 continue 

212 for fil in cpts: 

213 try: 

214 os.remove('my_%s.cpt' % fil) 

215 except OSError: 

216 continue 

217 

218 

219def xy_to_latlon(source, x, y): 

220 ''' 

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

222 

223 :param source: 

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

225 :type source: 

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

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

228 

229 :param x: 

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

231 :type x: 

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

233 

234 :param y: 

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

236 :type y: 

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

238 

239 :returns: 

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

241 :rtype: 

242 :py:class:`tuple` of :py:class:`float` 

243 ''' 

244 

245 s = source 

246 ax, ay = map_anchor[s.anchor] 

247 

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

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

250 

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

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

253 

254 northing += source.north_shift 

255 easting += source.east_shift 

256 

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

258 

259 

260def xy_to_lw(source, x, y): 

261 ''' 

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

263 

264 :param source: 

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

266 :type source: 

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

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

269 

270 :param x: 

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

272 :type x: 

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

274 

275 :param y: 

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

277 :type y: 

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

279 

280 :returns: 

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

282 :rtype: 

283 :py:class:`tuple` of :py:class:`float` 

284 ''' 

285 

286 length, width = source.length, source.width 

287 

288 ax, ay = map_anchor[source.anchor] 

289 

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

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

292 

293 return lengths, widths 

294 

295 

296cbar_anchor = { 

297 'center': 'MC', 

298 'center_left': 'ML', 

299 'center_right': 'MR', 

300 'top': 'TC', 

301 'top_left': 'TL', 

302 'top_right': 'TR', 

303 'bottom': 'BC', 

304 'bottom_left': 'BL', 

305 'bottom_right': 'BR'} 

306 

307 

308cbar_helper = { 

309 'traction': { 

310 'unit': 'MPa', 

311 'factor': 1e-6}, 

312 'tx': { 

313 'unit': 'MPa', 

314 'factor': 1e-6}, 

315 'ty': { 

316 'unit': 'MPa', 

317 'factor': 1e-6}, 

318 'tz': { 

319 'unit': 'MPa', 

320 'factor': 1e-6}, 

321 'time': { 

322 'unit': 's', 

323 'factor': 1.}, 

324 'strike': { 

325 'unit': '°', 

326 'factor': 1.}, 

327 'dip': { 

328 'unit': '°', 

329 'factor': 1.}, 

330 'vr': { 

331 'unit': 'km/s', 

332 'factor': 1e-3}, 

333 'length': { 

334 'unit': 'km', 

335 'factor': 1e-3}, 

336 'width': { 

337 'unit': 'km', 

338 'factor': 1e-3} 

339} 

340 

341 

342fonttype = 'Helvetica' 

343 

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

345 

346 

347def _make_gmt_conf(fontcolor, size): 

348 ''' 

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

350 

351 :param fontcolor: 

352 GMT readable colorcode or colorstring for the font. 

353 :type fontcolor: 

354 str 

355 

356 :param size: 

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

358 :type size: 

359 :py:class:`tuple` of :py:class:`float` 

360 

361 :returns: 

362 Best fitting fontsize, GMT configuration parameter 

363 :rtype: 

364 float, dict 

365 ''' 

366 

367 color = fontcolor 

368 fontsize = num.max(size) 

369 

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

371 

372 pen_size = fontsize / 16. 

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

374 

375 return fontsize, dict( 

376 MAP_FRAME_PEN='%s' % color, 

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

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

379 MAP_TICK_LENGTH_PRIMARY='%gc' % tick_size, 

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

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

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

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

384 PS_CHAR_ENCODING='ISOLatin1+', 

385 MAP_FRAME_TYPE='fancy', 

386 FORMAT_GEO_MAP='D', 

387 PS_PAGE_ORIENTATION='portrait', 

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

389 MAP_ANNOT_OBLIQUE='6') 

390 

391 

392class SourceError(Exception): 

393 pass 

394 

395 

396class RuptureMap(Map): 

397 ''' 

398 Map plotting of attributes and results of the 

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

400 ''' 

401 

402 def __init__( 

403 self, 

404 source=None, 

405 fontcolor='darkslategrey', 

406 width=20., 

407 height=14., 

408 margins=None, 

409 color_wet=(216, 242, 254), 

410 color_dry=(238, 236, 230), 

411 topo_cpt_wet='light_sea_uniform', 

412 topo_cpt_dry='light_land_uniform', 

413 show_cities=False, 

414 *args, 

415 **kwargs): 

416 

417 size = (width, height) 

418 fontsize, gmt_config = _make_gmt_conf(fontcolor, size) 

419 

420 if margins is None: 

421 margins = [ 

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

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

424 

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

426 gmt_config=gmt_config, 

427 topo_cpt_dry=topo_cpt_dry, topo_cpt_wet=topo_cpt_wet, 

428 color_wet=color_wet, color_dry=color_dry, 

429 **kwargs) 

430 

431 if show_cities: 

432 self.draw_cities() 

433 

434 self._source = source 

435 self._fontcolor = fontcolor 

436 self._fontsize = fontsize 

437 self.axes_layout = 'WeSn' 

438 

439 @property 

440 def size(self): 

441 ''' 

442 Figure size in [cm]. 

443 ''' 

444 return (self.width, self.height) 

445 

446 @property 

447 def font(self): 

448 ''' 

449 Font style (size and type). 

450 ''' 

451 

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

453 

454 @property 

455 def source(self): 

456 ''' 

457 PseudoDynamicRupture whose attributes are plotted. 

458 

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

460 ''' 

461 

462 if self._source is None: 

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

464 

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

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

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

468 

469 if not self._source.patches: 

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

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

472 

473 return self._source 

474 

475 @source.setter 

476 def source(self, source): 

477 self._source = source 

478 

479 def _get_topotile(self): 

480 if self._dems is None: 

481 self._setup() 

482 

483 try: 

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

485 except NoTopo: 

486 wesn = self.wesn 

487 

488 nx = int(num.ceil( 

489 self.width * cm2inch * self.topo_resolution_max)) 

490 ny = int(num.ceil( 

491 self.height * cm2inch * self.topo_resolution_max)) 

492 

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

494 

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

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

497 data) 

498 

499 return t 

500 

501 def _patches_to_lw(self): 

502 ''' 

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

504 

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

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

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

508 

509 :returns: 

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

511 :rtype: 

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

513 ''' 

514 

515 source = self.source 

516 patches = source.patches 

517 

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

519 

520 patch_lengths = num.concatenate(( 

521 num.array([0.]), 

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

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

524 

525 patch_widths = num.concatenate(( 

526 num.array([0.]), 

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

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

529 

530 ax, ay = map_anchor[source.anchor] 

531 

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

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

534 

535 return patch_lengths, patch_widths 

536 

537 def _xy_to_lw(self, x, y): 

538 ''' 

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

540 

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

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

543 downdip. 

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

545 in [m]. 

546 

547 :returns: 

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

549 :rtype: 

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

551 ''' 

552 

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

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

555 

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

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

558 ' Please check the x coordinates.') 

559 

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

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

562 ' Please check the y coordinates.') 

563 

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

565 

566 def _tile_to_lw(self, ref_lat, ref_lon, 

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

568 

569 ''' 

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

571 

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

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

574 is used for interpolation of grid data. 

575 

576 :param ref_lat: 

577 Reference latitude, from which length-width relative coordinates 

578 grid are calculated. 

579 :type ref_lat: 

580 float 

581 

582 :param ref_lon: 

583 Reference longitude, from which length-width relative coordinates 

584 grid are calculated. 

585 :type ref_lon: 

586 float 

587 

588 :param north_shift: 

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

590 :type north_shift: 

591 float 

592 

593 :param east_shift: 

594 East shift of the reference point [m]. 

595 :type east_shift: 

596 float 

597 

598 :param strike: 

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

600 :type strike: 

601 float 

602 

603 :returns: 

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

605 :rtype: 

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

607 ''' 

608 

609 t = self._get_topotile() 

610 grid_lats = t.y() 

611 grid_lons = t.x() 

612 

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

614 

615 grid_northing, grid_easting = pod.latlon_to_ne_numpy( 

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

617 

618 grid_northing -= north_shift 

619 grid_easting -= east_shift 

620 

621 strike *= d2r 

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

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

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

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

626 

627 return points_out 

628 

629 def _prep_patch_grid_data(self, data): 

630 ''' 

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

632 

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

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

635 plane. 

636 

637 :param data: 

638 Patch wise data. 

639 :type data: 

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

641 

642 :returns: 

643 Extended data array. 

644 :rtype: 

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

646 ''' 

647 

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

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

650 

651 data_new = num.zeros(shape) 

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

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

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

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

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

657 

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

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

660 

661 return data_new 

662 

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

664 ''' 

665 Interpolate regular data onto topotile grid. 

666 

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

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

669 

670 :param lengths: 

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

672 :type lengths: 

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

674 

675 :param widths: 

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

677 :type widths: 

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

679 

680 :param data: 

681 Data grid array. 

682 :type data: 

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

684 

685 :param filename: 

686 Filename, where grid is stored. 

687 :type filename: 

688 str 

689 ''' 

690 

691 source = self.source 

692 

693 interpolator = scrgi( 

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

695 data.T, 

696 bounds_error=False, 

697 method='nearest') 

698 

699 points_out = self._tile_to_lw( 

700 ref_lat=source.lat, 

701 ref_lon=source.lon, 

702 north_shift=source.north_shift, 

703 east_shift=source.east_shift, 

704 strike=source.strike) 

705 

706 t = self._get_topotile() 

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

708 t.data[:] = num.nan 

709 

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

711 

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

713 

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

715 ''' 

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

717 

718 :param data: 

719 Patchwise grid data. 

720 :type data: 

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

722 ''' 

723 

724 lengths, widths = self._patches_to_lw() 

725 

726 data_new = self._prep_patch_grid_data(data) 

727 

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

729 

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

731 ''' 

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

733 

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

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

736 

737 :param x: 

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

739 :type x: 

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

741 

742 :param y: 

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

744 :type y: 

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

746 

747 :param data: 

748 Patchwise grid data. 

749 :type data: 

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

751 ''' 

752 

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

754 

755 self._regular_data_to_grid( 

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

757 *args, **kwargs) 

758 

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

760 ''' 

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

762 

763 :param gridfile: 

764 File of the grid which shall be plotted. 

765 :type gridfile: 

766 str 

767 

768 :param cmap: 

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

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

771 :type cmap: 

772 str 

773 

774 :param cbar: 

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

776 added. Keyword arguments are parsed to it. 

777 :type cbar: 

778 bool 

779 ''' 

780 

781 self.gmt.grdimage( 

782 gridfile, 

783 C='my_%s.cpt' % cmap, 

784 E='200', 

785 Q=True, 

786 n='+t0.0', 

787 *self.jxyr) 

788 

789 if cbar: 

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

791 

792 def draw_contour( 

793 self, 

794 gridfile, 

795 contour_int, 

796 anot_int, 

797 angle=None, 

798 unit='', 

799 color='', 

800 style='', 

801 **kwargs): 

802 

803 ''' 

804 Draw grid data as contour lines. 

805 

806 :param gridfile: 

807 File of the grid which shall be plotted. 

808 :type gridfile: 

809 str 

810 

811 :param contour_int: 

812 Interval of contour lines in units of the gridfile. 

813 :type contour_int: 

814 float 

815 

816 :param anot_int: 

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

818 be a integer multiple of contour_int. 

819 :type anot_int: 

820 float 

821 

822 :param angle: 

823 Rotation angle of the labels in [deg]. 

824 :type angle: 

825 float 

826 

827 :param unit: 

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

829 label on labelled contour lines. 

830 :type unit: 

831 str 

832 

833 :param color: 

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

835 :type color: 

836 str 

837 

838 :param style: 

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

840 plotted. 

841 :type style: 

842 str 

843 ''' 

844 

845 pen_size = self._fontsize / 40. 

846 

847 if not color: 

848 color = self._fontcolor 

849 

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

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

852 if angle: 

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

854 c_string = '%g' % contour_int 

855 

856 if kwargs: 

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

858 else: 

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

860 

861 if style: 

862 style = ',' + style 

863 

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

865 

866 self.gmt.grdcontour( 

867 gridfile, 

868 S='10', 

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

870 *self.jxyr + args, 

871 **kwargs) 

872 

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

874 ''' 

875 Draw a colorbar based on a existing colormap. 

876 

877 :param cmap: 

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

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

880 :type cmap: 

881 str 

882 

883 :param label: 

884 Title of the colorbar. 

885 :type label: 

886 str 

887 

888 :param anchor: 

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

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

891 :type anchor: 

892 str 

893 ''' 

894 

895 if not kwargs: 

896 kwargs = {} 

897 

898 if label: 

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

900 

901 kwargs['C'] = 'my_%s.cpt' % cmap 

902 a_str = cbar_anchor[anchor] 

903 

904 w = self.width / 3. 

905 h = w / 10. 

906 

907 lgap = rgap = w / 10. 

908 bgap, tgap = h, h / 10. 

909 

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

911 

912 if 'bottom' in anchor: 

913 dy += 4 * h 

914 

915 self.gmt.psscale( 

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

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

918 *self.jxyr, 

919 **kwargs) 

920 

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

922 ''' 

923 Draw vectors based on two grid files. 

924 

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

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

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

928 are applied if not defined differently. 

929 

930 :param x_gridfile: 

931 File of the grid defining vector lengths in x. 

932 :type x_gridfile: 

933 str 

934 

935 :param y_gridfile: 

936 File of the grid defining vector lengths in y. 

937 :type y_gridfile: 

938 str 

939 

940 :param vcolor: 

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

942 :type vcolor: 

943 str 

944 ''' 

945 

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

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

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

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

950 

951 self.gmt.grdvector( 

952 x_gridfile, y_gridfile, 

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

954 *self.jxyr, 

955 **kwargs) 

956 

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

958 ''' 

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

960 

961 :param data: 

962 Patchwise grid data. 

963 :type data: 

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

965 ''' 

966 

967 plot_data = data 

968 

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

970 

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

972 

973 cpt = [] 

974 if not op.exists('my_%s.cpt' % kwargs['cmap']): 

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

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

977 cpt = [kwargs['cmap']] 

978 

979 tmp_grd_file = 'tmpdata.grd' 

980 self.patch_data_to_grid(plot_data, tmp_grd_file) 

981 self.draw_image(tmp_grd_file, **kwargs) 

982 

983 clear_temp(gridfiles=[tmp_grd_file], cpts=cpt) 

984 

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

986 ''' 

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

988 

989 :param attribute: 

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

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

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

993 length or the single components of the traction vector. 

994 :type attribute: 

995 str 

996 ''' 

997 

998 a = attribute 

999 source = self.source 

1000 

1001 if a == 'traction': 

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

1003 elif a == 'tx': 

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

1005 elif a == 'ty': 

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

1007 elif a == 'tz': 

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

1009 else: 

1010 data = source.get_patch_attribute(attribute) 

1011 

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

1013 data *= factor 

1014 

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

1016 'label', 

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

1018 

1019 self.draw_dynamic_data(data, **kwargs) 

1020 

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

1022 ''' 

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

1024 

1025 :param store: 

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

1027 :type store: 

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

1029 :param clevel: 

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

1031 :type clevel: 

1032 :py:class:`list` of :py:class:`float` 

1033 ''' 

1034 

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

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

1037 

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

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

1040 

1041 if clevel: 

1042 if len(clevel) > 1: 

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

1044 else: 

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

1046 

1047 kwargs['contour_int'] = kwargs['anot_int'] 

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

1049 

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

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

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

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

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

1055 

1056 tmp_grd_file = 'tmpdata.grd' 

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

1058 times, tmp_grd_file) 

1059 self.draw_contour(tmp_grd_file, **kwargs) 

1060 

1061 clear_temp(gridfiles=[tmp_grd_file], cpts=[]) 

1062 

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

1064 ''' 

1065 Draw points at given locations. 

1066 

1067 :param lats: 

1068 Point latitude coordinates in [deg]. 

1069 :type lats: 

1070 iterable of :py:class:`float` 

1071 

1072 :param lons: 

1073 Point longitude coordinates in [deg]. 

1074 :type lons: 

1075 iterable of :py:class:`float` 

1076 

1077 :param symbol: 

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

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

1080 :type symbol: 

1081 str 

1082 

1083 :param size: 

1084 Size of the points in [points]. 

1085 :type size: 

1086 float 

1087 ''' 

1088 

1089 sym_to_gmt = dict( 

1090 star='a', 

1091 circle='c', 

1092 point='p', 

1093 triangle='t') 

1094 

1095 lats = num.atleast_1d(lats) 

1096 lons = num.atleast_1d(lons) 

1097 

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

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

1100 

1101 if size is None: 

1102 size = self._fontsize / 3. 

1103 

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

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

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

1107 

1108 self.gmt.psxy( 

1109 in_columns=[lons, lats], 

1110 *self.jxyr, 

1111 **kwargs) 

1112 

1113 def draw_nucleation_point(self, **kwargs): 

1114 ''' 

1115 Plot the nucleation point onto the map. 

1116 ''' 

1117 

1118 nlat, nlon = xy_to_latlon( 

1119 self.source, self.source.nucleation_x, self.source.nucleation_y) 

1120 

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

1122 

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

1124 ''' 

1125 Draw dislocation onto map at any time. 

1126 

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

1128 component the patchwise dislocation is plotted onto the map. 

1129 

1130 :param time: 

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

1132 ``tmax`` is taken. 

1133 :type time: 

1134 float 

1135 

1136 :param component: 

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

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

1139 length of the dislocation vector is plotted. 

1140 :type component: str 

1141 ''' 

1142 

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

1144 

1145 if component: 

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

1147 else: 

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

1149 

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

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

1152 

1153 self.draw_dynamic_data(data, **kwargs) 

1154 

1155 def draw_dislocation_contour( 

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

1157 ''' 

1158 Draw dislocation contour onto map at any time. 

1159 

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

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

1162 

1163 :param time: 

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

1165 ``tmax`` is taken. 

1166 :type time: 

1167 float 

1168 

1169 :param component: 

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

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

1172 dislocation vector is plotted. 

1173 :type component: str 

1174 

1175 :param clevel: 

1176 Times, for which contour lines are drawn. 

1177 :type clevel: 

1178 :py:class:`list` of :py:class:`float` 

1179 ''' 

1180 

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

1182 

1183 if component: 

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

1185 else: 

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

1187 

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

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

1190 

1191 if clevel: 

1192 if len(clevel) > 1: 

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

1194 else: 

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

1196 

1197 kwargs['contour_int'] = kwargs['anot_int'] 

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

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

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

1201 else: 

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

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

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

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

1206 

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

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

1209 

1210 tmp_grd_file = 'tmpdata.grd' 

1211 self.patch_data_to_grid(data, tmp_grd_file) 

1212 self.draw_contour(tmp_grd_file, **kwargs) 

1213 

1214 clear_temp(gridfiles=[tmp_grd_file], cpts=[]) 

1215 

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

1217 ''' 

1218 Draw vector arrows onto map indicating direction of dislocation. 

1219 

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

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

1222 

1223 :param time: 

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

1225 ``None``, ``tmax`` is used. 

1226 :type time: 

1227 float 

1228 ''' 

1229 

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

1231 

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

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

1234 

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

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

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

1238 

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

1240 

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

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

1243 

1244 self.draw_vector( 

1245 tmp_grd_files[1], tmp_grd_files[0], 

1246 **kwargs) 

1247 

1248 clear_temp(gridfiles=tmp_grd_files, cpts=[]) 

1249 

1250 def draw_top_edge(self, **kwargs): 

1251 ''' 

1252 Indicate rupture top edge on map. 

1253 ''' 

1254 

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

1256 top_edge = outline[:2, :] 

1257 

1258 kwargs = kwargs or {} 

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

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

1261 

1262 self.gmt.psxy( 

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

1264 *self.jxyr, 

1265 **kwargs) 

1266 

1267 

1268class RuptureView(Object): 

1269 ''' 

1270 Plot of attributes and results of the 

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

1272 ''' 

1273 

1274 _patches_to_lw = RuptureMap._patches_to_lw 

1275 

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

1277 self._source = source 

1278 self._axes = None 

1279 

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

1281 self.fontsize = fontsize or 10 

1282 mpl_init(fontsize=self.fontsize) 

1283 

1284 self._fig = None 

1285 self._axes = None 

1286 self._is_1d = False 

1287 

1288 @property 

1289 def source(self): 

1290 ''' 

1291 PseudoDynamicRupture whose attributes are plotted. 

1292 

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

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

1295 ''' 

1296 

1297 if self._source is None: 

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

1299 

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

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

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

1303 

1304 if not self._source.patches: 

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

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

1307 

1308 return self._source 

1309 

1310 @source.setter 

1311 def source(self, source): 

1312 self._source = source 

1313 

1314 def _setup( 

1315 self, 

1316 title='', 

1317 xlabel='', 

1318 ylabel='', 

1319 aspect=1., 

1320 spatial_plot=True, 

1321 **kwargs): 

1322 

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

1324 return 

1325 

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

1327 

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

1329 ax = self._axes 

1330 

1331 if ax is not None: 

1332 ax.set_title(title) 

1333 ax.grid(alpha=.3) 

1334 ax.set_xlabel(xlabel) 

1335 ax.set_ylabel(ylabel) 

1336 

1337 if spatial_plot: 

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

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

1340 

1341 def _clear_all(self): 

1342 plt.cla() 

1343 plt.clf() 

1344 plt.close() 

1345 

1346 self._fig, self._axes = None, None 

1347 

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

1349 default_kwargs = dict( 

1350 linewidth=0, 

1351 marker='o', 

1352 markerfacecolor=mpl_color('skyblue2'), 

1353 markersize=6., 

1354 markeredgecolor=mpl_color('skyblue3')) 

1355 default_kwargs.update(kwargs) 

1356 

1357 if self._axes is not None: 

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

1359 

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

1361 if self._axes is not None: 

1362 if 'extent' not in kwargs: 

1363 kwargs['extent'] = [ 

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

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

1366 

1367 im = self._axes.imshow( 

1368 data, 

1369 *args, 

1370 interpolation='none', 

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

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

1373 **kwargs) 

1374 

1375 del kwargs['extent'] 

1376 if 'aspect' in kwargs: 

1377 del kwargs['aspect'] 

1378 if 'clim' in kwargs: 

1379 del kwargs['clim'] 

1380 if 'cmap' in kwargs: 

1381 del kwargs['cmap'] 

1382 

1383 plt.colorbar( 

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

1385 

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

1387 setup_kwargs = dict( 

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

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

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

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

1392 

1393 self._setup(**setup_kwargs) 

1394 

1395 if self._axes is not None: 

1396 if clevel is None: 

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

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

1399 

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

1401 

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

1403 clevel = num.array([clevel]) 

1404 

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

1406 

1407 cont = self._axes.contour( 

1408 x, y, data, clevel, *args, 

1409 linewidths=1.5, **kwargs) 

1410 

1411 clabels = self._axes.clabel( 

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

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

1414 inline_spacing=15, rightside_up=True, use_clabeltext=True, 

1415 **kwargs) 

1416 

1417 plt.setp(clabels, path_effects=[ 

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

1419 patheffects.Normal()]) 

1420 

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

1422 ''' 

1423 Draw a point onto the figure. 

1424 

1425 Args and kwargs can be defined according to 

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

1427 

1428 :param length: 

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

1430 the anchor point in [m]. 

1431 :type length: 

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

1433 

1434 :param width: 

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

1436 the anchor point in [m]. 

1437 :type width: 

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

1439 ''' 

1440 

1441 if self._axes is not None: 

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

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

1444 

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

1446 

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

1448 ''' 

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

1450 

1451 :param data: 

1452 Patchwise grid data. 

1453 :type data: 

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

1455 ''' 

1456 

1457 plot_data = data 

1458 

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

1460 

1461 length, width = xy_to_lw( 

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

1463 

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

1465 

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

1467 

1468 setup_kwargs = dict( 

1469 xlabel='along strike [km]', 

1470 ylabel='down dip [km]', 

1471 title='', 

1472 aspect=1) 

1473 setup_kwargs.update(kwargs) 

1474 

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

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

1477 self._setup(**setup_kwargs) 

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

1479 

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

1481 ''' 

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

1483 

1484 :param attribute: 

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

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

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

1488 or the single components of the traction vector. 

1489 :type attribute: 

1490 str 

1491 ''' 

1492 

1493 a = attribute 

1494 source = self.source 

1495 

1496 if a == 'traction': 

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

1498 elif a == 'tx': 

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

1500 elif a == 'ty': 

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

1502 elif a == 'tz': 

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

1504 else: 

1505 data = source.get_patch_attribute(attribute) 

1506 

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

1508 data *= factor 

1509 

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

1511 'label', 

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

1513 

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

1515 

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

1517 ''' 

1518 Draw high resolution contours of the rupture front propgation time 

1519 

1520 :param store: 

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

1522 :type store: 

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

1524 

1525 :param clevel: 

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

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

1528 :type clevel: 

1529 list 

1530 ''' 

1531 

1532 source = self.source 

1533 default_kwargs = dict( 

1534 colors='#474747' 

1535 ) 

1536 default_kwargs.update(kwargs) 

1537 

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

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

1540 store) 

1541 

1542 times = time_interpolator.values 

1543 

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

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

1546 

1547 if len(clevel) == 0: 

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

1549 

1550 points_x = time_interpolator.grid[0] 

1551 points_y = time_interpolator.grid[1] 

1552 

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

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

1555 

1556 def draw_nucleation_point(self, **kwargs): 

1557 ''' 

1558 Draw the nucleation point onto the map. 

1559 ''' 

1560 

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

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

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

1564 

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

1566 ''' 

1567 Draw dislocation onto map at any time. 

1568 

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

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

1571 

1572 :param time: 

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

1574 If ``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``, the 

1581 length of the dislocation vector is plotted 

1582 :type component: str 

1583 ''' 

1584 

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

1586 

1587 if component: 

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

1589 else: 

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

1591 

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

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

1594 

1595 self.draw_dynamic_data(data, **kwargs) 

1596 

1597 def draw_dislocation_contour( 

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

1599 ''' 

1600 Draw dislocation contour onto map at any time. 

1601 

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

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

1604 

1605 :param time: 

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

1607 ``None``, ``tmax`` is taken. 

1608 :type time: 

1609 float 

1610 

1611 :param component: 

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

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

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

1615 :type component: 

1616 str 

1617 ''' 

1618 

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

1620 

1621 if component: 

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

1623 else: 

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

1625 

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

1627 

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

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

1630 

1631 if len(clevel) == 0: 

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

1633 

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

1635 

1636 length, width = self._patches_to_lw() 

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

1638 

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

1640 

1641 self._setup(**kwargs) 

1642 self._draw_contour( 

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

1644 

1645 def draw_source_dynamics( 

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

1647 ''' 

1648 Display dynamic source parameter. 

1649 

1650 Fast inspection possibility for the cumulative moment and the source 

1651 time function approximation (assuming equal paths between different 

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

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

1654 rate function. 

1655 

1656 :param variable: 

1657 Dynamic parameter, which shall be plotted. Choose 

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

1659 :type variable: 

1660 str 

1661 

1662 :param store: 

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

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

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

1666 :type store: 

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

1668 

1669 :param deltat: 

1670 Time increment between two parameter snapshots. If not 

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

1672 :type deltat: 

1673 float 

1674 ''' 

1675 

1676 v = variable 

1677 

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

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

1680 

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

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

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

1684 data = num.cumsum(data * deltat) 

1685 name, unit = 'M', 'Nm' 

1686 else: 

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

1688 

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

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

1691 aspect='auto', 

1692 spatial_plot=False) 

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

1694 self._is_1d = True 

1695 

1696 def draw_patch_dynamics( 

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

1698 ''' 

1699 Display dynamic boundary element / patch parameter. 

1700 

1701 Fast inspection possibility for different dynamic parameter for a 

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

1703 the chosen patch. 

1704 

1705 :param variable: 

1706 Dynamic parameter, which shall be plotted. Choose 

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

1708 :type variable: 

1709 str 

1710 

1711 :param nx: 

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

1713 :type nx: 

1714 int 

1715 

1716 :param nx: 

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

1718 :type nx: 

1719 int 

1720 

1721 :param store: 

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

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

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

1725 :type store: 

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

1727 

1728 :param deltat: 

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

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

1731 :type deltat: 

1732 float 

1733 ''' 

1734 

1735 v = variable 

1736 source = self.source 

1737 idx = nx * source.ny + nx 

1738 

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

1740 

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

1742 data, times = source.get_moment_rate_patches( 

1743 store=store, deltat=deltat) 

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

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

1746 

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

1748 data, times = source.get_moment_rate_patches( 

1749 store=store, deltat=deltat) 

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

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

1752 data, times = source.get_moment_rate_patches( 

1753 store=store, deltat=deltat) 

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

1755 name, unit = 'M', 'Nm' 

1756 elif v == 'slip_rate': 

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

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

1759 elif v == 'dislocation': 

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

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

1762 name, unit = 'du', 'm' 

1763 elif m: 

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

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

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

1767 else: 

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

1769 

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

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

1772 aspect='auto', 

1773 spatial_plot=False) 

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

1775 *args, **kwargs) 

1776 self._is_1d = True 

1777 

1778 def finalize(self): 

1779 if self._is_1d: 

1780 return 

1781 

1782 length, width = xy_to_lw( 

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

1784 

1785 self._axes.set_xlim(length) 

1786 self._axes.set_ylim(width) 

1787 

1788 def gcf(self): 

1789 self.finalize() 

1790 return self._fig 

1791 

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

1793 ''' 

1794 Save plot to file. 

1795 

1796 :param filename: 

1797 Filename and path, where the plot is stored. 

1798 :type filename: 

1799 str 

1800 

1801 :param dpi: 

1802 Resolution of the output plot in [dpi]. 

1803 :type dpi: 

1804 int 

1805 ''' 

1806 

1807 self.finalize() 

1808 try: 

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

1810 except TypeError: 

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

1812 

1813 self._clear_all() 

1814 

1815 def show_plot(self): 

1816 ''' 

1817 Show plot. 

1818 ''' 

1819 

1820 self.finalize() 

1821 plt.show() 

1822 self._clear_all() 

1823 

1824 

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

1826 ''' 

1827 Generate a mp4 movie based on given png files using 

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

1829 

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

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

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

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

1834 

1835 :param fn_path: 

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

1837 :type fn_path: 

1838 str 

1839 

1840 :param output_path: 

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

1842 :type output_path: 

1843 str 

1844 

1845 :param deltat: 

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

1847 :type deltat: 

1848 float 

1849 ''' 

1850 

1851 try: 

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

1853 except CalledProcessError: 

1854 pass 

1855 except (TypeError, IOError): 

1856 logger.warning( 

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

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

1859 return False 

1860 

1861 try: 

1862 check_call([ 

1863 'ffmpeg', '-loglevel', 'info', '-y', 

1864 '-framerate', '%g' % framerate, 

1865 '-i', fn_path, 

1866 '-vcodec', 'libx264', 

1867 '-preset', 'medium', 

1868 '-tune', 'animation', 

1869 '-pix_fmt', 'yuv420p', 

1870 '-movflags', '+faststart', 

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

1872 '-crf', '15', 

1873 output_path]) 

1874 

1875 return True 

1876 except CalledProcessError as e: 

1877 logger.warning(e) 

1878 return False 

1879 

1880 

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

1882 ''' 

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

1884 

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

1886 

1887 :param fn: 

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

1889 :type fn: 

1890 str 

1891 

1892 :param output_path: 

1893 Path and filename of the output animated gif file. 

1894 :type output_path: 

1895 str 

1896 :param loops: 

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

1898 infinite. 

1899 :type loops: 

1900 int 

1901 ''' 

1902 

1903 try: 

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

1905 except CalledProcessError: 

1906 pass 

1907 except (TypeError, IOError): 

1908 logger.warning( 

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

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

1911 return False 

1912 

1913 try: 

1914 check_call([ 

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

1916 fn, 

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

1918 '-loop', '%d' % loops, 

1919 output_path]) 

1920 

1921 return True 

1922 except CalledProcessError as e: 

1923 logger.warning(e) 

1924 return False 

1925 

1926 

1927def rupture_movie( 

1928 source, 

1929 store, 

1930 variable='dislocation', 

1931 draw_time_contours=False, 

1932 fn_path='.', 

1933 prefix='', 

1934 plot_type='map', 

1935 deltat=None, 

1936 framerate=None, 

1937 store_images=False, 

1938 render_as_gif=False, 

1939 gif_loops=-1, 

1940 **kwargs): 

1941 ''' 

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

1943 

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

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

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

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

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

1949 

1950 :param source: 

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

1952 :type source: 

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

1954 

1955 :param store: 

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

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

1958 :type store: 

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

1960 

1961 :param variable: 

1962 Dynamic parameter, which shall be plotted. Choose between 

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

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

1965 ``moment_rate``, default ``dislocation``. 

1966 :type variable: 

1967 str 

1968 

1969 :param draw_time_contours: 

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

1971 :type draw_time_contours: 

1972 bool 

1973 

1974 :param fn_path: 

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

1976 :type fn_path: 

1977 str 

1978 

1979 :param prefix: 

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

1981 :type prefix: 

1982 str 

1983 

1984 :param plot_type: 

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

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

1987 or plane view using 

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

1989 :type plot_type: 

1990 str 

1991 

1992 :param deltat: 

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

1994 used to define ``deltat``. 

1995 :type deltat: 

1996 float 

1997 

1998 :param store_images: 

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

2000 not. 

2001 :type store_images: 

2002 bool 

2003 

2004 :param render_as_gif: 

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

2006 is returned as mp4. 

2007 :type render_as_gif: 

2008 bool 

2009 

2010 :param gif_loops: 

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

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

2013 infinite. 

2014 :type gif_loops: 

2015 int 

2016 ''' 

2017 

2018 v = variable 

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

2020 

2021 if not source.patches: 

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

2023 

2024 if source.coef_mat is None: 

2025 source.calc_coef_mat() 

2026 

2027 prefix = prefix or v 

2028 deltat = deltat or store.config.deltat 

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

2030 

2031 if v == 'moment_rate': 

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

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

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

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

2036 else: 

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

2038 

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

2040 

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

2042 if m: 

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

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

2045 elif v == 'dislocation': 

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

2047 name, unit = 'du', 'm' 

2048 elif v == 'slip_rate': 

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

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

2051 

2052 if plot_type == 'map': 

2053 plt_base = RuptureMap 

2054 elif plot_type == 'view': 

2055 plt_base = RuptureView 

2056 else: 

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

2058 

2059 attrs_base = get_kwargs(plt_base) 

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

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

2062 

2063 if 'clim' in kwargs_plt: 

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

2065 else: 

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

2067 

2068 if 'label' not in kwargs_plt: 

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

2070 data /= vmax 

2071 

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

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

2074 

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

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

2077 for it, _ in enumerate(times)] 

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

2079 

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

2081 plot_data = data[:, it] 

2082 

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

2084 plt.draw_dynamic_data(plot_data, **kwargs_plt) 

2085 plt.draw_nucleation_point() 

2086 

2087 if draw_time_contours: 

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

2089 

2090 plt.save(ft) 

2091 

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

2093 return_code = render_movie( 

2094 fn_temp_path, 

2095 output_path=fn_mp4, 

2096 framerate=framerate) 

2097 

2098 if render_as_gif and return_code: 

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

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

2101 loops=gif_loops) 

2102 

2103 elif return_code: 

2104 shutil.move(fn_mp4, op.join( 

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

2106 

2107 else: 

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

2109 

2110 if store_images: 

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

2112 for t in times] 

2113 

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

2115 shutil.move(ft, f) 

2116 

2117 

2118__all__ = [ 

2119 'make_colormap', 

2120 'clear_temp', 

2121 'xy_to_latlon', 

2122 'xy_to_lw', 

2123 'SourceError', 

2124 'RuptureMap', 

2125 'RuptureView', 

2126 'rupture_movie', 

2127 'render_movie']