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

431 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 math import pi as PI 

7import logging 

8import numpy as num 

9 

10from matplotlib.collections import PatchCollection 

11from matplotlib.patches import Polygon 

12from matplotlib.transforms import Transform 

13from matplotlib.colors import LinearSegmentedColormap 

14 

15from pyrocko import moment_tensor as mtm 

16 

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

18 

19NA = num.newaxis 

20d2r = num.pi / 180. 

21 

22 

23def view_rotation(strike, dip): 

24 return mtm.euler_to_matrix( 

25 dip*d2r, strike*d2r, -90.*d2r) 

26 

27 

28_view_south = view_rotation(90., 90.) 

29_view_north = view_rotation(-90., 90.) 

30_view_east = view_rotation(0., 90.) 

31_view_west = view_rotation(180., 90.) 

32 

33 

34class BeachballError(Exception): 

35 pass 

36 

37 

38class _FixedPointOffsetTransform(Transform): 

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

40 Transform.__init__(self) 

41 self.input_dims = self.output_dims = 2 

42 self.has_inverse = False 

43 self.trans = trans 

44 self.dpi_scale_trans = dpi_scale_trans 

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

46 

47 def transform_non_affine(self, values): 

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

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

50 

51 

52def vnorm(points): 

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

54 

55 

56def clean_poly(points): 

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

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

59 

60 dupl = num.concatenate( 

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

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

63 return points 

64 

65 

66def close_poly(points): 

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

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

69 

70 return points 

71 

72 

73def circulation(points, axis): 

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

75 

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

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

78 

79 result = -num.sum( 

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

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

82 

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

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

85 return result 

86 

87 

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

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

90 

91 # cut sub-polygons and gather crossing point information 

92 crossings = [] 

93 snippets = {} 

94 for ipath, points in enumerate(l_points): 

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

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

97 

98 # get upward crossing points 

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

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

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

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

103 points[iup, :]) 

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

105 

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

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

108 

109 # get downward crossing points 

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

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

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

113 points[idown+1, axis]) 

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

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

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

117 

118 for i in range(idown.size): 

119 crossings.append( 

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

121 

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

123 

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

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

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

127 

128 if icuts.size: 

129 points_last = num.concatenate(( 

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

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

132 

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

134 else: 

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

136 

137 crossings.sort() 

138 

139 # assemble new sub-polygons 

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

141 outs = [[]] 

142 while True: 

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

144 for i, c1 in enumerate(crossings): 

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

146 direction = -1 * c1[3] 

147 break 

148 else: 

149 if not snippets: 

150 break 

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

152 outs.append([]) 

153 continue 

154 

155 while True: 

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

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

158 break 

159 

160 c2 = crossings[i] 

161 c2[-1].remove(direction) 

162 

163 phi1 = c1[0] 

164 phi2 = c2[0] 

165 if direction == 1: 

166 if phi1 > phi2: 

167 phi2 += PI * 2. 

168 

169 if direction == -1: 

170 if phi1 < phi2: 

171 phi2 -= PI * 2. 

172 

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

174 

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

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

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

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

179 cpoints[:, axis] = 0.0 

180 

181 outs[-1].append(cpoints) 

182 

183 try: 

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

185 del snippets[c2[1:3]] 

186 

187 except KeyError: 

188 if not snippets: 

189 break 

190 

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

192 outs.append([]) 

193 

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

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

196 # endpoint) 

197 

198 outs_upper = [] 

199 outs_lower = [] 

200 for out in outs: 

201 if out: 

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

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

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

205 outs_upper.append(out) 

206 else: 

207 outs_lower.append(out) 

208 

209 if nonsimple and ( 

210 len(crossings) == 0 or 

211 len(outs_upper) == 0 or 

212 len(outs_lower) == 0): 

213 

214 # check if we are cutting between holes 

215 need_divider = False 

216 if outs_upper: 

217 candis = sorted( 

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

219 

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

221 need_divider = True 

222 

223 if outs_lower: 

224 candis = sorted( 

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

226 

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

228 need_divider = True 

229 

230 if need_divider: 

231 phi1 = 0. 

232 phi2 = PI*2. 

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

234 

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

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

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

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

239 cpoints[:, axis] = 0.0 

240 

241 outs_upper.append(cpoints) 

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

243 

244 return outs_lower, outs_upper 

245 

246 

247def numpy_rtp2xyz(rtp): 

248 r = rtp[:, 0] 

249 theta = rtp[:, 1] 

250 phi = rtp[:, 2] 

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

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

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

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

255 return vecs 

256 

257 

258def numpy_xyz2rtp(xyz): 

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

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

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

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

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

264 return vecs 

265 

266 

267def circle_points(aphi, sign=1.0): 

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

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

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

271 vecs[:, 2] = 0.0 

272 return vecs 

273 

274 

275def eig2gx(eig, arcres=181): 

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

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

278 

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

280 

281 groups = [] 

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

283 patches = [] 

284 patches_lower = [] 

285 patches_upper = [] 

286 lines = [] 

287 lines_lower = [] 

288 lines_upper = [] 

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

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

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

292 

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

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

295 from_e = to_e.T 

296 

297 poly_es = [] 

298 polys = [] 

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

300 xphi = perm_sign*pt_sign*sign*aphi 

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

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

303 continue 

304 

305 Y = -ea/denom 

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

307 continue 

308 

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

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

311 rtp[:, 0] = 1. 

312 if sign > 0: 

313 rtp[:, 1] = xtheta 

314 else: 

315 rtp[:, 1] = PI - xtheta 

316 

317 rtp[:, 2] = xphi 

318 poly_e = numpy_rtp2xyz(rtp) 

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

320 poly[:, 2] -= 0.001 

321 

322 poly_es.append(poly_e) 

323 polys.append(poly) 

324 

325 if polys: 

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

327 lines.extend(polys) 

328 lines_lower.extend(polys_lower) 

329 lines_upper.extend(polys_upper) 

330 

331 if poly_es: 

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

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

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

335 for poly_e in cc: 

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

337 poly[:, 2] -= 0.001 

338 polys_lower, polys_upper = spoly_cut( 

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

340 

341 patches.append(poly) 

342 patches_lower.extend(polys_lower) 

343 patches_upper.extend(polys_upper) 

344 

345 if not patches: 

346 if mt_sign * pt_sign == 1.: 

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

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

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

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

351 

352 groups.append(( 

353 pt_name, 

354 patches, patches_lower, patches_upper, 

355 lines, lines_lower, lines_upper)) 

356 

357 return groups 

358 

359 

360def extr(points): 

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

362 return points + pmean*0.05 

363 

364 

365def draw_eigenvectors_mpl(eig, axes): 

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

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

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

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

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

371 

372 

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

374 ''' 

375 Project 3D points to 2D. 

376 

377 :param projection: 

378 Projection to use. Choices: ``'lambert'``, ``'stereographic'``, 

379 ``'orthographic'``. 

380 :type projection: 

381 str 

382 ''' 

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

384 if projection == 'lambert': 

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

386 elif projection == 'stereographic': 

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

388 elif projection == 'orthographic': 

389 factor = None 

390 else: 

391 raise BeachballError( 

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

393 

394 if factor is not None: 

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

396 

397 return points_out 

398 

399 

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

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

402 

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

404 if projection == 'lambert': 

405 points_out[:, 2] = 1.0 - rsqr 

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

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

408 elif projection == 'stereographic': 

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

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

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

412 elif projection == 'orthographic': 

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

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

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

416 else: 

417 raise BeachballError( 

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

419 

420 return points_out 

421 

422 

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

424 mt = mtm.as_mt(mt) 

425 

426 if isinstance(view, str): 

427 if view == 'top': 

428 pass 

429 elif view == 'north': 

430 mt = mt.rotated(_view_north) 

431 elif view == 'south': 

432 mt = mt.rotated(_view_south) 

433 elif view == 'east': 

434 mt = mt.rotated(_view_east) 

435 elif view == 'west': 

436 mt = mt.rotated(_view_west) 

437 elif isinstance(view, tuple): 

438 mt = mt.rotated(view_rotation(*view)) 

439 else: 

440 raise BeachballError( 

441 'Invaild argument for `view`. Allowed values are "top", "north", ' 

442 '"south", "east", "west" or a tuple of angles `(strike, dip)` ' 

443 'orienting the view plane.') 

444 

445 if mt_type == 'full': 

446 return mt 

447 

448 res = mt.standard_decomposition() 

449 m = dict( 

450 dc=res[1][2], 

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

452 

453 return mtm.MomentTensor(m=m) 

454 

455 

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

457 

458 if size_units == 'points': 

459 transform = _FixedPointOffsetTransform( 

460 axes.transData, 

461 axes.figure.dpi_scale_trans, 

462 position) 

463 

464 if size is None: 

465 size = 12. 

466 

467 size = size * 0.5 / 72. 

468 position = (0., 0.) 

469 

470 elif size_units == 'data': 

471 transform = axes.transData 

472 

473 if size is None: 

474 size = 1.0 

475 

476 size = size * 0.5 

477 

478 elif size_units == 'axes': 

479 transform = axes.transAxes 

480 if size is None: 

481 size = 1. 

482 

483 size = size * .5 

484 

485 else: 

486 raise BeachballError( 

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

488 

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

490 

491 return transform, position, size 

492 

493 

494def mt2beachball( 

495 mt, 

496 beachball_type='deviatoric', 

497 position=(0., 0.), 

498 size=None, 

499 color_t='red', 

500 color_p='white', 

501 edgecolor='black', 

502 linewidth=2, 

503 projection='lambert', 

504 view='top'): 

505 

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

507 size = size or 1 

508 mt = deco_part(mt, beachball_type, view) 

509 

510 eig = mt.eigensystem() 

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

512 raise BeachballError('eigenvalues are zero') 

513 

514 data = [] 

515 for (group, patches, patches_lower, patches_upper, 

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

517 

518 if group == 'P': 

519 color = color_p 

520 else: 

521 color = color_t 

522 

523 for poly in patches_upper: 

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

525 position[NA, :] 

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

527 

528 for poly in lines_upper: 

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

530 position[NA, :] 

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

532 return data 

533 

534 

535def plot_beachball_mpl( 

536 mt, axes, 

537 beachball_type='deviatoric', 

538 position=(0., 0.), 

539 size=None, 

540 zorder=0, 

541 color_t='red', 

542 color_p='white', 

543 edgecolor='black', 

544 linewidth=2, 

545 alpha=1.0, 

546 arcres=181, 

547 decimation=1, 

548 projection='lambert', 

549 size_units='points', 

550 view='top'): 

551 

552 ''' 

553 Plot beachball diagram to a Matplotlib plot 

554 

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

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

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

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

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

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

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

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

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

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

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

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

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

568 ``'orthographic'`` 

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

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

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

572 beachball will be shown distorted when using this). 

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

574 ``'east'`` or ``'west'``, or project onto plane given by 

575 ``(strike, dip)``. Useful to show beachballs in cross-sections. 

576 Default is ``'top'``. 

577 ''' 

578 

579 transform, position, size = choose_transform( 

580 axes, size_units, position, size) 

581 

582 mt = deco_part(mt, beachball_type, view) 

583 

584 eig = mt.eigensystem() 

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

586 raise BeachballError('eigenvalues are zero') 

587 

588 data = [] 

589 for (group, patches, patches_lower, patches_upper, 

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

591 

592 if group == 'P': 

593 color = color_p 

594 else: 

595 color = color_t 

596 

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

598 # is NED 

599 

600 for poly in patches_upper: 

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

602 if alpha == 1.0: 

603 data.append( 

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

605 else: 

606 data.append( 

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

608 

609 for poly in lines_upper: 

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

611 data.append( 

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

613 

614 patches = [] 

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

616 patches.append(Polygon( 

617 xy=path, facecolor=facecolor, 

618 edgecolor=edgecolor, 

619 linewidth=linewidth, 

620 alpha=alpha)) 

621 

622 collection = PatchCollection( 

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

624 

625 axes.add_artist(collection) 

626 return collection 

627 

628 

629def amplitudes_ned(mt, vecs): 

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

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

632 vecs_e = num.dot(to_e, vecs.T).T 

633 rtp = numpy_xyz2rtp(vecs_e) 

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

635 return ep * num.cos(atheta)**2 + ( 

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

637 

638 

639def amplitudes(mt, azimuths, takeoff_angles): 

640 ''' 

641 Get beachball amplitude values for selected azimuths and takeoff angles. 

642 

643 :param azimuths: 

644 Azimuths, measured clockwise from north [deg]. 

645 :type azimuths: 

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

647 

648 :param takeoff_angles: 

649 Takeoff angles, measured from downward vertical [deg]. 

650 :type takeoff_angles: 

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

652 ''' 

653 azimuths = num.asarray(azimuths, dtype=float) 

654 takeoff_angles = num.asarray(takeoff_angles, dtype=float) 

655 assert azimuths.size == takeoff_angles.size 

656 rtps = num.vstack( 

657 (num.ones(azimuths.size), takeoff_angles*d2r, azimuths*d2r)).T 

658 vecs = numpy_rtp2xyz(rtps) 

659 return amplitudes_ned(mt, vecs) 

660 

661 

662def mts2amps( 

663 mts, 

664 projection, 

665 beachball_type, 

666 grid_resolution=200, 

667 mask=True, 

668 view='top'): 

669 

670 n_balls = len(mts) 

671 nx = grid_resolution 

672 ny = grid_resolution 

673 

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

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

676 

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

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

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

680 

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

682 amps = num.full(nx * ny, num.nan, dtype=num.float64) 

683 

684 amps[ii_ok] = 0. 

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

686 

687 for mt in mts: 

688 amps_ok = amplitudes_ned(deco_part(mt, beachball_type, view), vecs3_ok) 

689 if mask: 

690 amps_ok[amps_ok > 0] = 1. 

691 amps_ok[amps_ok < 0] = 0. 

692 

693 amps[ii_ok] += amps_ok 

694 

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

696 

697 

698def plot_fuzzy_beachball_mpl_pixmap( 

699 mts, axes, 

700 best_mt=None, 

701 beachball_type='deviatoric', 

702 position=(0., 0.), 

703 size=None, 

704 zorder=0, 

705 color_t='red', 

706 color_p='white', 

707 edgecolor='black', 

708 best_color='red', 

709 linewidth=2, 

710 alpha=1.0, 

711 projection='lambert', 

712 size_units='data', 

713 grid_resolution=200, 

714 method='imshow', 

715 view='top'): 

716 ''' 

717 Plot fuzzy beachball from a list of given MomentTensors 

718 

719 :param mts: list of 

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

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

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

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

724 of most likely or minimum misfit solution to extra highlight 

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

726 polygons are not plotted 

727 

728 See plot_beachball_mpl for other arguments 

729 ''' 

730 if size_units == 'points': 

731 raise BeachballError( 

732 'size_units="points" not supported in ' 

733 'plot_fuzzy_beachball_mpl_pixmap') 

734 

735 transform, position, size = choose_transform( 

736 axes, size_units, position, size) 

737 

738 amps, x, y = mts2amps( 

739 mts, 

740 grid_resolution=grid_resolution, 

741 projection=projection, 

742 beachball_type=beachball_type, 

743 mask=True, 

744 view=view) 

745 

746 ncolors = 256 

747 cmap = LinearSegmentedColormap.from_list( 

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

749 

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

751 if method == 'contourf': 

752 axes.contourf( 

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

754 levels=levels, 

755 cmap=cmap, 

756 transform=transform, 

757 zorder=zorder, 

758 alpha=alpha) 

759 

760 elif method == 'imshow': 

761 axes.imshow( 

762 amps.T, 

763 extent=( 

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

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

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

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

768 cmap=cmap, 

769 transform=transform, 

770 zorder=zorder-0.1, 

771 alpha=alpha) 

772 else: 

773 assert False, 'invalid `method` argument' 

774 

775 # draw optimum edges 

776 if best_mt is not None: 

777 best_amps, bx, by = mts2amps( 

778 [best_mt], 

779 grid_resolution=grid_resolution, 

780 projection=projection, 

781 beachball_type=beachball_type, 

782 mask=False) 

783 

784 axes.contour( 

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

786 levels=[0.], 

787 colors=[best_color], 

788 linewidths=linewidth, 

789 transform=transform, 

790 zorder=zorder, 

791 alpha=alpha) 

792 

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

794 x = num.cos(phi) 

795 y = num.sin(phi) 

796 axes.plot( 

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

798 linewidth=linewidth, 

799 color=edgecolor, 

800 transform=transform, 

801 zorder=zorder, 

802 alpha=alpha) 

803 

804 

805def plot_beachball_mpl_construction( 

806 mt, axes, 

807 show='patches', 

808 beachball_type='deviatoric', 

809 view='top'): 

810 

811 mt = deco_part(mt, beachball_type, view) 

812 eig = mt.eigensystem() 

813 

814 for (group, patches, patches_lower, patches_upper, 

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

816 

817 if group == 'P': 

818 color = 'blue' 

819 lw = 1 

820 else: 

821 color = 'red' 

822 lw = 1 

823 

824 if show == 'patches': 

825 for poly in patches_upper: 

826 px, py, pz = poly.T 

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

828 

829 if show == 'lines': 

830 for poly in lines_upper: 

831 px, py, pz = poly.T 

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

833 

834 

835def plot_beachball_mpl_pixmap( 

836 mt, axes, 

837 beachball_type='deviatoric', 

838 position=(0., 0.), 

839 size=None, 

840 zorder=0, 

841 color_t='red', 

842 color_p='white', 

843 edgecolor='black', 

844 linewidth=2, 

845 alpha=1.0, 

846 projection='lambert', 

847 size_units='data', 

848 view='top'): 

849 

850 if size_units == 'points': 

851 raise BeachballError( 

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

853 

854 transform, position, size = choose_transform( 

855 axes, size_units, position, size) 

856 

857 mt = deco_part(mt, beachball_type, view) 

858 

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

860 

861 amps, x, y = mts2amps( 

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

863 

864 axes.contourf( 

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

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

867 colors=[color_p, color_t], 

868 transform=transform, 

869 zorder=zorder, 

870 alpha=alpha) 

871 

872 axes.contour( 

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

874 levels=[0.], 

875 colors=[edgecolor], 

876 linewidths=linewidth, 

877 transform=transform, 

878 zorder=zorder, 

879 alpha=alpha) 

880 

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

882 x = num.cos(phi) 

883 y = num.sin(phi) 

884 axes.plot( 

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

886 linewidth=linewidth, 

887 color=edgecolor, 

888 transform=transform, 

889 zorder=zorder, 

890 alpha=alpha) 

891 

892 

893if __name__ == '__main__': 

894 import sys 

895 import os 

896 import matplotlib.pyplot as plt 

897 from pyrocko import model 

898 

899 args = sys.argv[1:] 

900 

901 data = [] 

902 for iarg, arg in enumerate(args): 

903 

904 if os.path.exists(arg): 

905 events = model.load_events(arg) 

906 for ev in events: 

907 if not ev.moment_tensor: 

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

909 continue 

910 

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

912 else: 

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

914 mt = mtm.as_mt(vals) 

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

916 

917 n = len(data) 

918 

919 ncols = 1 

920 while ncols**2 < n: 

921 ncols += 1 

922 

923 nrows = ncols 

924 

925 fig = plt.figure() 

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

927 axes.axison = False 

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

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

930 

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

932 irow = ibeach // ncols 

933 icol = ibeach % ncols 

934 plot_beachball_mpl( 

935 mt, axes, 

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

937 size_units='data') 

938 

939 axes.annotate( 

940 name, 

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

942 xycoords='data', 

943 xytext=(0, 0), 

944 textcoords='offset points', 

945 verticalalignment='center', 

946 horizontalalignment='center', 

947 rotation=0.) 

948 

949 plt.show()