1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5# python 2/3 

6from __future__ import absolute_import 

7 

8from math import pi as PI 

9import logging 

10import numpy as num 

11 

12from matplotlib.collections import PatchCollection 

13from matplotlib.patches import Polygon 

14from matplotlib.transforms import Transform 

15from matplotlib.colors import LinearSegmentedColormap 

16 

17from pyrocko import moment_tensor as mtm 

18from pyrocko.util import num_full 

19 

20logger = logging.getLogger('pyrocko.plot.beachball') 

21 

22NA = num.newaxis 

23 

24_view_south = num.array([[0, 0, -1], 

25 [0, 1, 0], 

26 [1, 0, 0]]) 

27 

28_view_north = _view_south.T 

29 

30_view_east = num.array([[1, 0, 0], 

31 [0, 0, -1], 

32 [0, 1, 0]]) 

33 

34_view_west = _view_east.T 

35 

36 

37class BeachballError(Exception): 

38 pass 

39 

40 

41class FixedPointOffsetTransform(Transform): 

42 def __init__(self, trans, dpi_scale_trans, fixed_point): 

43 Transform.__init__(self) 

44 self.input_dims = self.output_dims = 2 

45 self.has_inverse = False 

46 self.trans = trans 

47 self.dpi_scale_trans = dpi_scale_trans 

48 self.fixed_point = num.asarray(fixed_point, dtype=num.float64) 

49 

50 def transform_non_affine(self, values): 

51 fp = self.trans.transform(self.fixed_point) 

52 return fp + self.dpi_scale_trans.transform(values) 

53 

54 

55def vnorm(points): 

56 return num.sqrt(num.sum(points**2, axis=1)) 

57 

58 

59def clean_poly(points): 

60 if not num.all(points[0, :] == points[-1, :]): 

61 points = num.vstack((points, points[0:1, :])) 

62 

63 dupl = num.concatenate( 

64 (num.all(points[1:, :] == points[:-1, :], axis=1), [False])) 

65 points = points[num.logical_not(dupl)] 

66 return points 

67 

68 

69def close_poly(points): 

70 if not num.all(points[0, :] == points[-1, :]): 

71 points = num.vstack((points, points[0:1, :])) 

72 

73 return points 

74 

75 

76def circulation(points, axis): 

77 # assert num.all(points[:, axis] >= 0.0) or num.all(points[:, axis] <= 0.0) 

78 

79 points2 = points[:, ((axis+2) % 3, (axis+1) % 3)].copy() 

80 points2 *= 1.0 / num.sqrt(1.0 + num.abs(points[:, axis]))[:, num.newaxis] 

81 

82 result = -num.sum( 

83 (points2[1:, 0] - points2[:-1, 0]) * 

84 (points2[1:, 1] + points2[:-1, 1])) 

85 

86 result -= (points2[0, 0] - points2[-1, 0]) \ 

87 * (points2[0, 1] + points2[-1, 1]) 

88 return result 

89 

90 

91def spoly_cut(l_points, axis=0, nonsimple=True, arcres=181): 

92 dphi = 2.*PI / (2*arcres) 

93 

94 # cut sub-polygons and gather crossing point information 

95 crossings = [] 

96 snippets = {} 

97 for ipath, points in enumerate(l_points): 

98 if not num.all(points[0, :] == points[-1, :]): 

99 points = num.vstack((points, points[0:1, :])) 

100 

101 # get upward crossing points 

102 iup = num.where(num.logical_and(points[:-1, axis] <= 0., 

103 points[1:, axis] > 0.))[0] 

104 aup = - points[iup, axis] / (points[iup+1, axis] - points[iup, axis]) 

105 pup = points[iup, :] + aup[:, num.newaxis] * (points[iup+1, :] - 

106 points[iup, :]) 

107 phiup = num.arctan2(pup[:, (axis+2) % 3], pup[:, (axis+1) % 3]) 

108 

109 for i in range(len(iup)): 

110 crossings.append((phiup[i], ipath, iup[i], 1, pup[i], [1, -1])) 

111 

112 # get downward crossing points 

113 idown = num.where(num.logical_and(points[:-1, axis] > 0., 

114 points[1:, axis] <= 0.))[0] 

115 adown = - points[idown+1, axis] / (points[idown, axis] - 

116 points[idown+1, axis]) 

117 pdown = points[idown+1, :] + adown[:, num.newaxis] * ( 

118 points[idown, :] - points[idown+1, :]) 

119 phidown = num.arctan2(pdown[:, (axis+2) % 3], pdown[:, (axis+1) % 3]) 

120 

121 for i in range(idown.size): 

122 crossings.append( 

123 (phidown[i], ipath, idown[i], -1, pdown[i], [1, -1])) 

124 

125 icuts = num.sort(num.concatenate((iup, idown))) 

126 

127 for i in range(icuts.size-1): 

128 snippets[ipath, icuts[i]] = ( 

129 ipath, icuts[i+1], points[icuts[i]+1:icuts[i+1]+1]) 

130 

131 if icuts.size: 

132 points_last = num.concatenate(( 

133 points[icuts[-1]+1:], 

134 points[:icuts[0]+1])) 

135 

136 snippets[ipath, icuts[-1]] = (ipath, icuts[0], points_last) 

137 else: 

138 snippets[ipath, 0] = (ipath, 0, points) 

139 

140 crossings.sort() 

141 

142 # assemble new sub-polygons 

143 current = snippets.pop(list(snippets.keys())[0]) 

144 outs = [[]] 

145 while True: 

146 outs[-1].append(current[2]) 

147 for i, c1 in enumerate(crossings): 

148 if c1[1:3] == current[:2]: 

149 direction = -1 * c1[3] 

150 break 

151 else: 

152 if not snippets: 

153 break 

154 current = snippets.pop(list(snippets.keys())[0]) 

155 outs.append([]) 

156 continue 

157 

158 while True: 

159 i = (i + direction) % len(crossings) 

160 if crossings[i][3] == direction and direction in crossings[i][-1]: 

161 break 

162 

163 c2 = crossings[i] 

164 c2[-1].remove(direction) 

165 

166 phi1 = c1[0] 

167 phi2 = c2[0] 

168 if direction == 1: 

169 if phi1 > phi2: 

170 phi2 += PI * 2. 

171 

172 if direction == -1: 

173 if phi1 < phi2: 

174 phi2 -= PI * 2. 

175 

176 n = int(abs(phi2 - phi1) / dphi) + 2 

177 

178 phis = num.linspace(phi1, phi2, n) 

179 cpoints = num.zeros((n, 3)) 

180 cpoints[:, (axis+1) % 3] = num.cos(phis) 

181 cpoints[:, (axis+2) % 3] = num.sin(phis) 

182 cpoints[:, axis] = 0.0 

183 

184 outs[-1].append(cpoints) 

185 

186 try: 

187 current = snippets[c2[1:3]] 

188 del snippets[c2[1:3]] 

189 

190 except KeyError: 

191 if not snippets: 

192 break 

193 

194 current = snippets.pop(list(snippets.keys())[0]) 

195 outs.append([]) 

196 

197 # separate hemispheres, force polygons closed, remove duplicate points 

198 # remove polygons with less than 3 points (4, when counting repeated 

199 # endpoint) 

200 

201 outs_upper = [] 

202 outs_lower = [] 

203 for out in outs: 

204 if out: 

205 out = clean_poly(num.vstack(out)) 

206 if out.shape[0] >= 4: 

207 if num.sum(out[:, axis]) > 0.0: 

208 outs_upper.append(out) 

209 else: 

210 outs_lower.append(out) 

211 

212 if nonsimple and ( 

213 len(crossings) == 0 or 

214 len(outs_upper) == 0 or 

215 len(outs_lower) == 0): 

216 

217 # check if we are cutting between holes 

218 need_divider = False 

219 if outs_upper: 

220 candis = sorted( 

221 outs_upper, key=lambda out: num.min(out[:, axis])) 

222 

223 if circulation(candis[0], axis) > 0.0: 

224 need_divider = True 

225 

226 if outs_lower: 

227 candis = sorted( 

228 outs_lower, key=lambda out: num.max(out[:, axis])) 

229 

230 if circulation(candis[0], axis) < 0.0: 

231 need_divider = True 

232 

233 if need_divider: 

234 phi1 = 0. 

235 phi2 = PI*2. 

236 n = int(abs(phi2 - phi1) / dphi) + 2 

237 

238 phis = num.linspace(phi1, phi2, n) 

239 cpoints = num.zeros((n, 3)) 

240 cpoints[:, (axis+1) % 3] = num.cos(phis) 

241 cpoints[:, (axis+2) % 3] = num.sin(phis) 

242 cpoints[:, axis] = 0.0 

243 

244 outs_upper.append(cpoints) 

245 outs_lower.append(cpoints[::-1, :]) 

246 

247 return outs_lower, outs_upper 

248 

249 

250def numpy_rtp2xyz(rtp): 

251 r = rtp[:, 0] 

252 theta = rtp[:, 1] 

253 phi = rtp[:, 2] 

254 vecs = num.empty(rtp.shape, dtype=num.float64) 

255 vecs[:, 0] = r*num.sin(theta)*num.cos(phi) 

256 vecs[:, 1] = r*num.sin(theta)*num.sin(phi) 

257 vecs[:, 2] = r*num.cos(theta) 

258 return vecs 

259 

260 

261def numpy_xyz2rtp(xyz): 

262 x, y, z = xyz[:, 0], xyz[:, 1], xyz[:, 2] 

263 vecs = num.empty(xyz.shape, dtype=num.float64) 

264 vecs[:, 0] = num.sqrt(x**2+y**2+z**2) 

265 vecs[:, 1] = num.arctan2(num.sqrt(x**2+y**2), z) 

266 vecs[:, 2] = num.arctan2(y, x) 

267 return vecs 

268 

269 

270def circle_points(aphi, sign=1.0): 

271 vecs = num.empty((aphi.size, 3), dtype=num.float64) 

272 vecs[:, 0] = num.cos(sign*aphi) 

273 vecs[:, 1] = num.sin(sign*aphi) 

274 vecs[:, 2] = 0.0 

275 return vecs 

276 

277 

278def eig2gx(eig, arcres=181): 

279 aphi = num.linspace(0., 2.*PI, arcres) 

280 ep, en, et, vp, vn, vt = eig 

281 

282 mt_sign = num.sign(ep + en + et) 

283 

284 groups = [] 

285 for (pt_name, pt_sign) in [('P', -1.), ('T', 1.)]: 

286 patches = [] 

287 patches_lower = [] 

288 patches_upper = [] 

289 lines = [] 

290 lines_lower = [] 

291 lines_upper = [] 

292 for iperm, (va, vb, vc, ea, eb, ec) in enumerate([ 

293 (vp, vn, vt, ep, en, et), 

294 (vt, vp, vn, et, ep, en)]): # (vn, vt, vp, en, et, ep)]): 

295 

296 perm_sign = [-1.0, 1.0][iperm] 

297 to_e = num.vstack((vb, vc, va)) 

298 from_e = to_e.T 

299 

300 poly_es = [] 

301 polys = [] 

302 for sign in (-1., 1.): 

303 xphi = perm_sign*pt_sign*sign*aphi 

304 denom = eb*num.cos(xphi)**2 + ec*num.sin(xphi)**2 

305 if num.any(denom == 0.): 

306 continue 

307 

308 Y = -ea/denom 

309 if num.any(Y < 0.): 

310 continue 

311 

312 xtheta = num.arctan(num.sqrt(Y)) 

313 rtp = num.empty(xphi.shape+(3,), dtype=num.float64) 

314 rtp[:, 0] = 1. 

315 if sign > 0: 

316 rtp[:, 1] = xtheta 

317 else: 

318 rtp[:, 1] = PI - xtheta 

319 

320 rtp[:, 2] = xphi 

321 poly_e = numpy_rtp2xyz(rtp) 

322 poly = num.dot(from_e, poly_e.T).T 

323 poly[:, 2] -= 0.001 

324 

325 poly_es.append(poly_e) 

326 polys.append(poly) 

327 

328 if polys: 

329 polys_lower, polys_upper = spoly_cut(polys, 2, arcres=arcres) 

330 lines.extend(polys) 

331 lines_lower.extend(polys_lower) 

332 lines_upper.extend(polys_upper) 

333 

334 if poly_es: 

335 for aa in spoly_cut(poly_es, 0, arcres=arcres): 

336 for bb in spoly_cut(aa, 1, arcres=arcres): 

337 for cc in spoly_cut(bb, 2, arcres=arcres): 

338 for poly_e in cc: 

339 poly = num.dot(from_e, poly_e.T).T 

340 poly[:, 2] -= 0.001 

341 polys_lower, polys_upper = spoly_cut( 

342 [poly], 2, nonsimple=False, arcres=arcres) 

343 

344 patches.append(poly) 

345 patches_lower.extend(polys_lower) 

346 patches_upper.extend(polys_upper) 

347 

348 if not patches: 

349 if mt_sign * pt_sign == 1.: 

350 patches_lower.append(circle_points(aphi, -1.0)) 

351 patches_upper.append(circle_points(aphi, 1.0)) 

352 lines_lower.append(circle_points(aphi, -1.0)) 

353 lines_upper.append(circle_points(aphi, 1.0)) 

354 

355 groups.append(( 

356 pt_name, 

357 patches, patches_lower, patches_upper, 

358 lines, lines_lower, lines_upper)) 

359 

360 return groups 

361 

362 

363def extr(points): 

364 pmean = num.mean(points, axis=0) 

365 return points + pmean*0.05 

366 

367 

368def draw_eigenvectors_mpl(eig, axes): 

369 vp, vn, vt = eig[3:] 

370 for lab, v in [('P', vp), ('N', vn), ('T', vt)]: 

371 sign = num.sign(v[2]) + (v[2] == 0.0) 

372 axes.plot(sign*v[1], sign*v[0], 'o', color='black') 

373 axes.text(sign*v[1], sign*v[0], ' '+lab) 

374 

375 

376def project(points, projection='lambert'): 

377 points_out = points[:, :2].copy() 

378 if projection == 'lambert': 

379 factor = 1.0 / num.sqrt(1.0 + points[:, 2]) 

380 elif projection == 'stereographic': 

381 factor = 1.0 / (1.0 + points[:, 2]) 

382 elif projection == 'orthographic': 

383 factor = None 

384 else: 

385 raise BeachballError( 

386 'invalid argument for projection: %s' % projection) 

387 

388 if factor is not None: 

389 points_out *= factor[:, num.newaxis] 

390 

391 return points_out 

392 

393 

394def inverse_project(points, projection='lambert'): 

395 points_out = num.zeros((points.shape[0], 3)) 

396 

397 rsqr = points[:, 0]**2 + points[:, 1]**2 

398 if projection == 'lambert': 

399 points_out[:, 2] = 1.0 - rsqr 

400 points_out[:, 1] = num.sqrt(2.0 - rsqr) * points[:, 1] 

401 points_out[:, 0] = num.sqrt(2.0 - rsqr) * points[:, 0] 

402 elif projection == 'stereographic': 

403 points_out[:, 2] = - (rsqr - 1.0) / (rsqr + 1.0) 

404 points_out[:, 1] = 2.0 * points[:, 1] / (rsqr + 1.0) 

405 points_out[:, 0] = 2.0 * points[:, 0] / (rsqr + 1.0) 

406 elif projection == 'orthographic': 

407 points_out[:, 2] = num.sqrt(num.maximum(1.0 - rsqr, 0.0)) 

408 points_out[:, 1] = points[:, 1] 

409 points_out[:, 0] = points[:, 0] 

410 else: 

411 raise BeachballError( 

412 'invalid argument for projection: %s' % projection) 

413 

414 return points_out 

415 

416 

417def deco_part(mt, mt_type='full', view='top'): 

418 assert view in ('top', 'north', 'south', 'east', 'west'),\ 

419 'Allowed views are top, north, south, east and west' 

420 mt = mtm.as_mt(mt) 

421 

422 if view == 'top': 

423 pass 

424 elif view == 'north': 

425 mt = mt.rotated(_view_north) 

426 elif view == 'south': 

427 mt = mt.rotated(_view_south) 

428 elif view == 'east': 

429 mt = mt.rotated(_view_east) 

430 elif view == 'west': 

431 mt = mt.rotated(_view_west) 

432 

433 if mt_type == 'full': 

434 return mt 

435 

436 res = mt.standard_decomposition() 

437 m = dict( 

438 dc=res[1][2], 

439 deviatoric=res[3][2])[mt_type] 

440 

441 return mtm.MomentTensor(m=m) 

442 

443 

444def choose_transform(axes, size_units, position, size): 

445 

446 if size_units == 'points': 

447 transform = FixedPointOffsetTransform( 

448 axes.transData, 

449 axes.figure.dpi_scale_trans, 

450 position) 

451 

452 if size is None: 

453 size = 12. 

454 

455 size = size * 0.5 / 72. 

456 position = (0., 0.) 

457 

458 elif size_units == 'data': 

459 transform = axes.transData 

460 

461 if size is None: 

462 size = 1.0 

463 

464 size = size * 0.5 

465 

466 elif size_units == 'axes': 

467 transform = axes.transAxes 

468 if size is None: 

469 size = 1. 

470 

471 size = size * .5 

472 

473 else: 

474 raise BeachballError( 

475 'invalid argument for size_units: %s' % size_units) 

476 

477 position = num.asarray(position, dtype=num.float64) 

478 

479 return transform, position, size 

480 

481 

482def mt2beachball( 

483 mt, 

484 beachball_type='deviatoric', 

485 position=(0., 0.), 

486 size=None, 

487 color_t='red', 

488 color_p='white', 

489 edgecolor='black', 

490 linewidth=2, 

491 projection='lambert', 

492 view='top'): 

493 

494 position = num.asarray(position, dtype=num.float64) 

495 size = size or 1 

496 mt = deco_part(mt, beachball_type, view) 

497 

498 eig = mt.eigensystem() 

499 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0: 

500 raise BeachballError('eigenvalues are zero') 

501 

502 data = [] 

503 for (group, patches, patches_lower, patches_upper, 

504 lines, lines_lower, lines_upper) in eig2gx(eig): 

505 

506 if group == 'P': 

507 color = color_p 

508 else: 

509 color = color_t 

510 

511 for poly in patches_upper: 

512 verts = project(poly, projection)[:, ::-1] * size + \ 

513 position[NA, :] 

514 data.append((verts, color, color, 1.0)) 

515 

516 for poly in lines_upper: 

517 verts = project(poly, projection)[:, ::-1] * size + \ 

518 position[NA, :] 

519 data.append((verts, 'none', edgecolor, linewidth)) 

520 return data 

521 

522 

523def plot_beachball_mpl( 

524 mt, axes, 

525 beachball_type='deviatoric', 

526 position=(0., 0.), 

527 size=None, 

528 zorder=0, 

529 color_t='red', 

530 color_p='white', 

531 edgecolor='black', 

532 linewidth=2, 

533 alpha=1.0, 

534 arcres=181, 

535 decimation=1, 

536 projection='lambert', 

537 size_units='points', 

538 view='top'): 

539 

540 ''' 

541 Plot beachball diagram to a Matplotlib plot 

542 

543 :param mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or an 

544 array or sequence which can be converted into an MT object 

545 :param beachball_type: ``'deviatoric'`` (default), ``'full'``, or ``'dc'`` 

546 :param position: position of the beachball in data coordinates 

547 :param size: diameter of the beachball either in points or in data 

548 coordinates, depending on the ``size_units`` setting 

549 :param zorder: (passed through to matplotlib drawing functions) 

550 :param color_t: color for compressional quadrants (default: ``'red'``) 

551 :param color_p: color for extensive quadrants (default: ``'white'``) 

552 :param edgecolor: color for lines (default: ``'black'``) 

553 :param linewidth: linewidth in points (default: ``2``) 

554 :param alpha: (passed through to matplotlib drawing functions) 

555 :param projection: ``'lambert'`` (default), ``'stereographic'``, or 

556 ``'orthographic'`` 

557 :param size_units: ``'points'`` (default) or ``'data'``, where the 

558 latter causes the beachball to be projected in the plots data 

559 coordinates (axes must have an aspect ratio of 1.0 or the 

560 beachball will be shown distorted when using this). 

561 :param view: View the beachball from ``top``, ``north``, ``south``, 

562 ``east`` or ``west``. Useful for to show beachballs in cross-sections. 

563 Default is ``top``. 

564 ''' 

565 

566 transform, position, size = choose_transform( 

567 axes, size_units, position, size) 

568 

569 mt = deco_part(mt, beachball_type, view) 

570 

571 eig = mt.eigensystem() 

572 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0: 

573 raise BeachballError('eigenvalues are zero') 

574 

575 data = [] 

576 for (group, patches, patches_lower, patches_upper, 

577 lines, lines_lower, lines_upper) in eig2gx(eig, arcres): 

578 

579 if group == 'P': 

580 color = color_p 

581 else: 

582 color = color_t 

583 

584 # plot "upper" features for lower hemisphere, because coordinate system 

585 # is NED 

586 

587 for poly in patches_upper: 

588 verts = project(poly, projection)[:, ::-1] * size + position[NA, :] 

589 if alpha == 1.0: 

590 data.append( 

591 (verts[::decimation], color, color, linewidth)) 

592 else: 

593 data.append( 

594 (verts[::decimation], color, 'none', 0.0)) 

595 

596 for poly in lines_upper: 

597 verts = project(poly, projection)[:, ::-1] * size + position[NA, :] 

598 data.append( 

599 (verts[::decimation], 'none', edgecolor, linewidth)) 

600 

601 patches = [] 

602 for (path, facecolor, edgecolor, linewidth) in data: 

603 patches.append(Polygon( 

604 xy=path, facecolor=facecolor, 

605 edgecolor=edgecolor, 

606 linewidth=linewidth, 

607 alpha=alpha)) 

608 

609 collection = PatchCollection( 

610 patches, zorder=zorder, transform=transform, match_original=True) 

611 

612 axes.add_artist(collection) 

613 return collection 

614 

615 

616def mts2amps(mts, projection, beachball_type, grid_resolution=200, mask=True, 

617 view='top'): 

618 

619 n_balls = len(mts) 

620 nx = grid_resolution 

621 ny = grid_resolution 

622 

623 x = num.linspace(-1., 1., nx) 

624 y = num.linspace(-1., 1., ny) 

625 

626 vecs2 = num.zeros((nx * ny, 2), dtype=num.float64) 

627 vecs2[:, 0] = num.tile(x, ny) 

628 vecs2[:, 1] = num.repeat(y, nx) 

629 

630 ii_ok = vecs2[:, 0]**2 + vecs2[:, 1]**2 <= 1.0 

631 amps = num_full(nx * ny, num.nan, dtype=num.float64) 

632 

633 amps[ii_ok] = 0. 

634 for mt in mts: 

635 mt = deco_part(mt, beachball_type, view) 

636 

637 ep, en, et, vp, vn, vt = mt.eigensystem() 

638 

639 vecs3_ok = inverse_project(vecs2[ii_ok, :], projection) 

640 

641 to_e = num.vstack((vn, vt, vp)) 

642 

643 vecs_e = num.dot(to_e, vecs3_ok.T).T 

644 rtp = numpy_xyz2rtp(vecs_e) 

645 

646 atheta, aphi = rtp[:, 1], rtp[:, 2] 

647 amps_ok = ep * num.cos(atheta)**2 + ( 

648 en * num.cos(aphi)**2 + et * num.sin(aphi)**2) * num.sin(atheta)**2 

649 

650 if mask: 

651 amps_ok[amps_ok > 0] = 1. 

652 amps_ok[amps_ok < 0] = 0. 

653 

654 amps[ii_ok] += amps_ok 

655 

656 return num.reshape(amps, (ny, nx)) / n_balls, x, y 

657 

658 

659def plot_fuzzy_beachball_mpl_pixmap( 

660 mts, axes, 

661 best_mt=None, 

662 beachball_type='deviatoric', 

663 position=(0., 0.), 

664 size=None, 

665 zorder=0, 

666 color_t='red', 

667 color_p='white', 

668 edgecolor='black', 

669 best_color='red', 

670 linewidth=2, 

671 alpha=1.0, 

672 projection='lambert', 

673 size_units='data', 

674 grid_resolution=200, 

675 method='imshow', 

676 view='top'): 

677 ''' 

678 Plot fuzzy beachball from a list of given MomentTensors 

679 

680 :param mts: list of 

681 :py:class:`pyrocko.moment_tensor.MomentTensor` object or an 

682 array or sequence which can be converted into an MT object 

683 :param best_mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or 

684 an array or sequence which can be converted into an MT object 

685 of most likely or minimum misfit solution to extra highlight 

686 :param best_color: mpl color for best MomentTensor edges, 

687 polygons are not plotted 

688 

689 See plot_beachball_mpl for other arguments 

690 ''' 

691 if size_units == 'points': 

692 raise BeachballError( 

693 'size_units="points" not supported in ' 

694 'plot_fuzzy_beachball_mpl_pixmap') 

695 

696 transform, position, size = choose_transform( 

697 axes, size_units, position, size) 

698 

699 amps, x, y = mts2amps( 

700 mts, 

701 grid_resolution=grid_resolution, 

702 projection=projection, 

703 beachball_type=beachball_type, 

704 mask=True, 

705 view=view) 

706 

707 ncolors = 256 

708 cmap = LinearSegmentedColormap.from_list( 

709 'dummy', [color_p, color_t], N=ncolors) 

710 

711 levels = num.linspace(0, 1., ncolors) 

712 if method == 'contourf': 

713 axes.contourf( 

714 position[0] + y * size, position[1] + x * size, amps.T, 

715 levels=levels, 

716 cmap=cmap, 

717 transform=transform, 

718 zorder=zorder, 

719 alpha=alpha) 

720 

721 elif method == 'imshow': 

722 axes.imshow( 

723 amps.T, 

724 extent=( 

725 position[0] + y[0] * size, 

726 position[0] + y[-1] * size, 

727 position[1] - x[0] * size, 

728 position[1] - x[-1] * size), 

729 cmap=cmap, 

730 transform=transform, 

731 zorder=zorder-0.1, 

732 alpha=alpha) 

733 else: 

734 assert False, 'invalid `method` argument' 

735 

736 # draw optimum edges 

737 if best_mt is not None: 

738 best_amps, bx, by = mts2amps( 

739 [best_mt], 

740 grid_resolution=grid_resolution, 

741 projection=projection, 

742 beachball_type=beachball_type, 

743 mask=False) 

744 

745 axes.contour( 

746 position[0] + by * size, position[1] + bx * size, best_amps.T, 

747 levels=[0.], 

748 colors=[best_color], 

749 linewidths=linewidth, 

750 transform=transform, 

751 zorder=zorder, 

752 alpha=alpha) 

753 

754 phi = num.linspace(0., 2 * PI, 361) 

755 x = num.cos(phi) 

756 y = num.sin(phi) 

757 axes.plot( 

758 position[0] + x * size, position[1] + y * size, 

759 linewidth=linewidth, 

760 color=edgecolor, 

761 transform=transform, 

762 zorder=zorder, 

763 alpha=alpha) 

764 

765 

766def plot_beachball_mpl_construction( 

767 mt, axes, 

768 show='patches', 

769 beachball_type='deviatoric', 

770 view='top'): 

771 

772 mt = deco_part(mt, beachball_type, view) 

773 eig = mt.eigensystem() 

774 

775 for (group, patches, patches_lower, patches_upper, 

776 lines, lines_lower, lines_upper) in eig2gx(eig): 

777 

778 if group == 'P': 

779 color = 'blue' 

780 lw = 1 

781 else: 

782 color = 'red' 

783 lw = 1 

784 

785 if show == 'patches': 

786 for poly in patches_upper: 

787 px, py, pz = poly.T 

788 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5) 

789 

790 if show == 'lines': 

791 for poly in lines_upper: 

792 px, py, pz = poly.T 

793 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5) 

794 

795 

796def plot_beachball_mpl_pixmap( 

797 mt, axes, 

798 beachball_type='deviatoric', 

799 position=(0., 0.), 

800 size=None, 

801 zorder=0, 

802 color_t='red', 

803 color_p='white', 

804 edgecolor='black', 

805 linewidth=2, 

806 alpha=1.0, 

807 projection='lambert', 

808 size_units='data', 

809 view='top'): 

810 

811 if size_units == 'points': 

812 raise BeachballError( 

813 'size_units="points" not supported in plot_beachball_mpl_pixmap') 

814 

815 transform, position, size = choose_transform( 

816 axes, size_units, position, size) 

817 

818 mt = deco_part(mt, beachball_type, view) 

819 

820 ep, en, et, vp, vn, vt = mt.eigensystem() 

821 

822 amps, x, y = mts2amps( 

823 [mt], projection, beachball_type, grid_resolution=200, mask=False) 

824 

825 axes.contourf( 

826 position[0] + y * size, position[1] + x * size, amps.T, 

827 levels=[-num.inf, 0., num.inf], 

828 colors=[color_p, color_t], 

829 transform=transform, 

830 zorder=zorder, 

831 alpha=alpha) 

832 

833 axes.contour( 

834 position[0] + y * size, position[1] + x * size, amps.T, 

835 levels=[0.], 

836 colors=[edgecolor], 

837 linewidths=linewidth, 

838 transform=transform, 

839 zorder=zorder, 

840 alpha=alpha) 

841 

842 phi = num.linspace(0., 2 * PI, 361) 

843 x = num.cos(phi) 

844 y = num.sin(phi) 

845 axes.plot( 

846 position[0] + x * size, position[1] + y * size, 

847 linewidth=linewidth, 

848 color=edgecolor, 

849 transform=transform, 

850 zorder=zorder, 

851 alpha=alpha) 

852 

853 

854if __name__ == '__main__': 

855 import sys 

856 import os 

857 import matplotlib.pyplot as plt 

858 from pyrocko import model 

859 

860 args = sys.argv[1:] 

861 

862 data = [] 

863 for iarg, arg in enumerate(args): 

864 

865 if os.path.exists(arg): 

866 events = model.load_events(arg) 

867 for ev in events: 

868 if not ev.moment_tensor: 

869 logger.warning('no moment tensor given for event') 

870 continue 

871 

872 data.append((ev.name, ev.moment_tensor)) 

873 else: 

874 vals = list(map(float, arg.split(','))) 

875 mt = mtm.as_mt(vals) 

876 data.append(('%i' % (iarg+1), mt)) 

877 

878 n = len(data) 

879 

880 ncols = 1 

881 while ncols**2 < n: 

882 ncols += 1 

883 

884 nrows = ncols 

885 

886 fig = plt.figure() 

887 axes = fig.add_subplot(1, 1, 1, aspect=1.) 

888 axes.axison = False 

889 axes.set_xlim(-0.05 - ncols, ncols + 0.05) 

890 axes.set_ylim(-0.05 - nrows, nrows + 0.05) 

891 

892 for ibeach, (name, mt) in enumerate(data): 

893 irow = ibeach // ncols 

894 icol = ibeach % ncols 

895 plot_beachball_mpl( 

896 mt, axes, 

897 position=(icol*2-ncols+1, -irow*2+nrows-1), 

898 size_units='data') 

899 

900 axes.annotate( 

901 name, 

902 xy=(icol*2-ncols+1, -irow*2+nrows-2), 

903 xycoords='data', 

904 xytext=(0, 0), 

905 textcoords='offset points', 

906 verticalalignment='center', 

907 horizontalalignment='center', 

908 rotation=0.) 

909 

910 plt.show()