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

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