1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6from math import pi as PI
7import logging
8import numpy as num
10from matplotlib.collections import PatchCollection
11from matplotlib.patches import Polygon
12from matplotlib.transforms import Transform
13from matplotlib.colors import LinearSegmentedColormap
15from pyrocko import moment_tensor as mtm
16from pyrocko.util import num_full
18logger = logging.getLogger('pyrocko.plot.beachball')
20NA = num.newaxis
22_view_south = num.array([[0, 0, -1],
23 [0, 1, 0],
24 [1, 0, 0]])
26_view_north = _view_south.T
28_view_east = num.array([[1, 0, 0],
29 [0, 0, -1],
30 [0, 1, 0]])
32_view_west = _view_east.T
35class BeachballError(Exception):
36 pass
39class _FixedPointOffsetTransform(Transform):
40 def __init__(self, trans, dpi_scale_trans, fixed_point):
41 Transform.__init__(self)
42 self.input_dims = self.output_dims = 2
43 self.has_inverse = False
44 self.trans = trans
45 self.dpi_scale_trans = dpi_scale_trans
46 self.fixed_point = num.asarray(fixed_point, dtype=num.float64)
48 def transform_non_affine(self, values):
49 fp = self.trans.transform(self.fixed_point)
50 return fp + self.dpi_scale_trans.transform(values)
53def vnorm(points):
54 return num.sqrt(num.sum(points**2, axis=1))
57def clean_poly(points):
58 if not num.all(points[0, :] == points[-1, :]):
59 points = num.vstack((points, points[0:1, :]))
61 dupl = num.concatenate(
62 (num.all(points[1:, :] == points[:-1, :], axis=1), [False]))
63 points = points[num.logical_not(dupl)]
64 return points
67def close_poly(points):
68 if not num.all(points[0, :] == points[-1, :]):
69 points = num.vstack((points, points[0:1, :]))
71 return points
74def circulation(points, axis):
75 # assert num.all(points[:, axis] >= 0.0) or num.all(points[:, axis] <= 0.0)
77 points2 = points[:, ((axis+2) % 3, (axis+1) % 3)].copy()
78 points2 *= 1.0 / num.sqrt(1.0 + num.abs(points[:, axis]))[:, num.newaxis]
80 result = -num.sum(
81 (points2[1:, 0] - points2[:-1, 0]) *
82 (points2[1:, 1] + points2[:-1, 1]))
84 result -= (points2[0, 0] - points2[-1, 0]) \
85 * (points2[0, 1] + points2[-1, 1])
86 return result
89def spoly_cut(l_points, axis=0, nonsimple=True, arcres=181):
90 dphi = 2.*PI / (2*arcres)
92 # cut sub-polygons and gather crossing point information
93 crossings = []
94 snippets = {}
95 for ipath, points in enumerate(l_points):
96 if not num.all(points[0, :] == points[-1, :]):
97 points = num.vstack((points, points[0:1, :]))
99 # get upward crossing points
100 iup = num.where(num.logical_and(points[:-1, axis] <= 0.,
101 points[1:, axis] > 0.))[0]
102 aup = - points[iup, axis] / (points[iup+1, axis] - points[iup, axis])
103 pup = points[iup, :] + aup[:, num.newaxis] * (points[iup+1, :] -
104 points[iup, :])
105 phiup = num.arctan2(pup[:, (axis+2) % 3], pup[:, (axis+1) % 3])
107 for i in range(len(iup)):
108 crossings.append((phiup[i], ipath, iup[i], 1, pup[i], [1, -1]))
110 # get downward crossing points
111 idown = num.where(num.logical_and(points[:-1, axis] > 0.,
112 points[1:, axis] <= 0.))[0]
113 adown = - points[idown+1, axis] / (points[idown, axis] -
114 points[idown+1, axis])
115 pdown = points[idown+1, :] + adown[:, num.newaxis] * (
116 points[idown, :] - points[idown+1, :])
117 phidown = num.arctan2(pdown[:, (axis+2) % 3], pdown[:, (axis+1) % 3])
119 for i in range(idown.size):
120 crossings.append(
121 (phidown[i], ipath, idown[i], -1, pdown[i], [1, -1]))
123 icuts = num.sort(num.concatenate((iup, idown)))
125 for i in range(icuts.size-1):
126 snippets[ipath, icuts[i]] = (
127 ipath, icuts[i+1], points[icuts[i]+1:icuts[i+1]+1])
129 if icuts.size:
130 points_last = num.concatenate((
131 points[icuts[-1]+1:],
132 points[:icuts[0]+1]))
134 snippets[ipath, icuts[-1]] = (ipath, icuts[0], points_last)
135 else:
136 snippets[ipath, 0] = (ipath, 0, points)
138 crossings.sort()
140 # assemble new sub-polygons
141 current = snippets.pop(list(snippets.keys())[0])
142 outs = [[]]
143 while True:
144 outs[-1].append(current[2])
145 for i, c1 in enumerate(crossings):
146 if c1[1:3] == current[:2]:
147 direction = -1 * c1[3]
148 break
149 else:
150 if not snippets:
151 break
152 current = snippets.pop(list(snippets.keys())[0])
153 outs.append([])
154 continue
156 while True:
157 i = (i + direction) % len(crossings)
158 if crossings[i][3] == direction and direction in crossings[i][-1]:
159 break
161 c2 = crossings[i]
162 c2[-1].remove(direction)
164 phi1 = c1[0]
165 phi2 = c2[0]
166 if direction == 1:
167 if phi1 > phi2:
168 phi2 += PI * 2.
170 if direction == -1:
171 if phi1 < phi2:
172 phi2 -= PI * 2.
174 n = int(abs(phi2 - phi1) / dphi) + 2
176 phis = num.linspace(phi1, phi2, n)
177 cpoints = num.zeros((n, 3))
178 cpoints[:, (axis+1) % 3] = num.cos(phis)
179 cpoints[:, (axis+2) % 3] = num.sin(phis)
180 cpoints[:, axis] = 0.0
182 outs[-1].append(cpoints)
184 try:
185 current = snippets[c2[1:3]]
186 del snippets[c2[1:3]]
188 except KeyError:
189 if not snippets:
190 break
192 current = snippets.pop(list(snippets.keys())[0])
193 outs.append([])
195 # separate hemispheres, force polygons closed, remove duplicate points
196 # remove polygons with less than 3 points (4, when counting repeated
197 # endpoint)
199 outs_upper = []
200 outs_lower = []
201 for out in outs:
202 if out:
203 out = clean_poly(num.vstack(out))
204 if out.shape[0] >= 4:
205 if num.sum(out[:, axis]) > 0.0:
206 outs_upper.append(out)
207 else:
208 outs_lower.append(out)
210 if nonsimple and (
211 len(crossings) == 0 or
212 len(outs_upper) == 0 or
213 len(outs_lower) == 0):
215 # check if we are cutting between holes
216 need_divider = False
217 if outs_upper:
218 candis = sorted(
219 outs_upper, key=lambda out: num.min(out[:, axis]))
221 if circulation(candis[0], axis) > 0.0:
222 need_divider = True
224 if outs_lower:
225 candis = sorted(
226 outs_lower, key=lambda out: num.max(out[:, axis]))
228 if circulation(candis[0], axis) < 0.0:
229 need_divider = True
231 if need_divider:
232 phi1 = 0.
233 phi2 = PI*2.
234 n = int(abs(phi2 - phi1) / dphi) + 2
236 phis = num.linspace(phi1, phi2, n)
237 cpoints = num.zeros((n, 3))
238 cpoints[:, (axis+1) % 3] = num.cos(phis)
239 cpoints[:, (axis+2) % 3] = num.sin(phis)
240 cpoints[:, axis] = 0.0
242 outs_upper.append(cpoints)
243 outs_lower.append(cpoints[::-1, :])
245 return outs_lower, outs_upper
248def numpy_rtp2xyz(rtp):
249 r = rtp[:, 0]
250 theta = rtp[:, 1]
251 phi = rtp[:, 2]
252 vecs = num.empty(rtp.shape, dtype=num.float64)
253 vecs[:, 0] = r*num.sin(theta)*num.cos(phi)
254 vecs[:, 1] = r*num.sin(theta)*num.sin(phi)
255 vecs[:, 2] = r*num.cos(theta)
256 return vecs
259def numpy_xyz2rtp(xyz):
260 x, y, z = xyz[:, 0], xyz[:, 1], xyz[:, 2]
261 vecs = num.empty(xyz.shape, dtype=num.float64)
262 vecs[:, 0] = num.sqrt(x**2+y**2+z**2)
263 vecs[:, 1] = num.arctan2(num.sqrt(x**2+y**2), z)
264 vecs[:, 2] = num.arctan2(y, x)
265 return vecs
268def circle_points(aphi, sign=1.0):
269 vecs = num.empty((aphi.size, 3), dtype=num.float64)
270 vecs[:, 0] = num.cos(sign*aphi)
271 vecs[:, 1] = num.sin(sign*aphi)
272 vecs[:, 2] = 0.0
273 return vecs
276def eig2gx(eig, arcres=181):
277 aphi = num.linspace(0., 2.*PI, arcres)
278 ep, en, et, vp, vn, vt = eig
280 mt_sign = num.sign(ep + en + et)
282 groups = []
283 for (pt_name, pt_sign) in [('P', -1.), ('T', 1.)]:
284 patches = []
285 patches_lower = []
286 patches_upper = []
287 lines = []
288 lines_lower = []
289 lines_upper = []
290 for iperm, (va, vb, vc, ea, eb, ec) in enumerate([
291 (vp, vn, vt, ep, en, et),
292 (vt, vp, vn, et, ep, en)]): # (vn, vt, vp, en, et, ep)]):
294 perm_sign = [-1.0, 1.0][iperm]
295 to_e = num.vstack((vb, vc, va))
296 from_e = to_e.T
298 poly_es = []
299 polys = []
300 for sign in (-1., 1.):
301 xphi = perm_sign*pt_sign*sign*aphi
302 denom = eb*num.cos(xphi)**2 + ec*num.sin(xphi)**2
303 if num.any(denom == 0.):
304 continue
306 Y = -ea/denom
307 if num.any(Y < 0.):
308 continue
310 xtheta = num.arctan(num.sqrt(Y))
311 rtp = num.empty(xphi.shape+(3,), dtype=num.float64)
312 rtp[:, 0] = 1.
313 if sign > 0:
314 rtp[:, 1] = xtheta
315 else:
316 rtp[:, 1] = PI - xtheta
318 rtp[:, 2] = xphi
319 poly_e = numpy_rtp2xyz(rtp)
320 poly = num.dot(from_e, poly_e.T).T
321 poly[:, 2] -= 0.001
323 poly_es.append(poly_e)
324 polys.append(poly)
326 if polys:
327 polys_lower, polys_upper = spoly_cut(polys, 2, arcres=arcres)
328 lines.extend(polys)
329 lines_lower.extend(polys_lower)
330 lines_upper.extend(polys_upper)
332 if poly_es:
333 for aa in spoly_cut(poly_es, 0, arcres=arcres):
334 for bb in spoly_cut(aa, 1, arcres=arcres):
335 for cc in spoly_cut(bb, 2, arcres=arcres):
336 for poly_e in cc:
337 poly = num.dot(from_e, poly_e.T).T
338 poly[:, 2] -= 0.001
339 polys_lower, polys_upper = spoly_cut(
340 [poly], 2, nonsimple=False, arcres=arcres)
342 patches.append(poly)
343 patches_lower.extend(polys_lower)
344 patches_upper.extend(polys_upper)
346 if not patches:
347 if mt_sign * pt_sign == 1.:
348 patches_lower.append(circle_points(aphi, -1.0))
349 patches_upper.append(circle_points(aphi, 1.0))
350 lines_lower.append(circle_points(aphi, -1.0))
351 lines_upper.append(circle_points(aphi, 1.0))
353 groups.append((
354 pt_name,
355 patches, patches_lower, patches_upper,
356 lines, lines_lower, lines_upper))
358 return groups
361def extr(points):
362 pmean = num.mean(points, axis=0)
363 return points + pmean*0.05
366def draw_eigenvectors_mpl(eig, axes):
367 vp, vn, vt = eig[3:]
368 for lab, v in [('P', vp), ('N', vn), ('T', vt)]:
369 sign = num.sign(v[2]) + (v[2] == 0.0)
370 axes.plot(sign*v[1], sign*v[0], 'o', color='black')
371 axes.text(sign*v[1], sign*v[0], ' '+lab)
374def project(points, projection='lambert'):
375 points_out = points[:, :2].copy()
376 if projection == 'lambert':
377 factor = 1.0 / num.sqrt(1.0 + points[:, 2])
378 elif projection == 'stereographic':
379 factor = 1.0 / (1.0 + points[:, 2])
380 elif projection == 'orthographic':
381 factor = None
382 else:
383 raise BeachballError(
384 'invalid argument for projection: %s' % projection)
386 if factor is not None:
387 points_out *= factor[:, num.newaxis]
389 return points_out
392def inverse_project(points, projection='lambert'):
393 points_out = num.zeros((points.shape[0], 3))
395 rsqr = points[:, 0]**2 + points[:, 1]**2
396 if projection == 'lambert':
397 points_out[:, 2] = 1.0 - rsqr
398 points_out[:, 1] = num.sqrt(2.0 - rsqr) * points[:, 1]
399 points_out[:, 0] = num.sqrt(2.0 - rsqr) * points[:, 0]
400 elif projection == 'stereographic':
401 points_out[:, 2] = - (rsqr - 1.0) / (rsqr + 1.0)
402 points_out[:, 1] = 2.0 * points[:, 1] / (rsqr + 1.0)
403 points_out[:, 0] = 2.0 * points[:, 0] / (rsqr + 1.0)
404 elif projection == 'orthographic':
405 points_out[:, 2] = num.sqrt(num.maximum(1.0 - rsqr, 0.0))
406 points_out[:, 1] = points[:, 1]
407 points_out[:, 0] = points[:, 0]
408 else:
409 raise BeachballError(
410 'invalid argument for projection: %s' % projection)
412 return points_out
415def deco_part(mt, mt_type='full', view='top'):
416 assert view in ('top', 'north', 'south', 'east', 'west'),\
417 'Allowed views are top, north, south, east and west'
418 mt = mtm.as_mt(mt)
420 if view == 'top':
421 pass
422 elif view == 'north':
423 mt = mt.rotated(_view_north)
424 elif view == 'south':
425 mt = mt.rotated(_view_south)
426 elif view == 'east':
427 mt = mt.rotated(_view_east)
428 elif view == 'west':
429 mt = mt.rotated(_view_west)
431 if mt_type == 'full':
432 return mt
434 res = mt.standard_decomposition()
435 m = dict(
436 dc=res[1][2],
437 deviatoric=res[3][2])[mt_type]
439 return mtm.MomentTensor(m=m)
442def choose_transform(axes, size_units, position, size):
444 if size_units == 'points':
445 transform = _FixedPointOffsetTransform(
446 axes.transData,
447 axes.figure.dpi_scale_trans,
448 position)
450 if size is None:
451 size = 12.
453 size = size * 0.5 / 72.
454 position = (0., 0.)
456 elif size_units == 'data':
457 transform = axes.transData
459 if size is None:
460 size = 1.0
462 size = size * 0.5
464 elif size_units == 'axes':
465 transform = axes.transAxes
466 if size is None:
467 size = 1.
469 size = size * .5
471 else:
472 raise BeachballError(
473 'invalid argument for size_units: %s' % size_units)
475 position = num.asarray(position, dtype=num.float64)
477 return transform, position, size
480def mt2beachball(
481 mt,
482 beachball_type='deviatoric',
483 position=(0., 0.),
484 size=None,
485 color_t='red',
486 color_p='white',
487 edgecolor='black',
488 linewidth=2,
489 projection='lambert',
490 view='top'):
492 position = num.asarray(position, dtype=num.float64)
493 size = size or 1
494 mt = deco_part(mt, beachball_type, view)
496 eig = mt.eigensystem()
497 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0:
498 raise BeachballError('eigenvalues are zero')
500 data = []
501 for (group, patches, patches_lower, patches_upper,
502 lines, lines_lower, lines_upper) in eig2gx(eig):
504 if group == 'P':
505 color = color_p
506 else:
507 color = color_t
509 for poly in patches_upper:
510 verts = project(poly, projection)[:, ::-1] * size + \
511 position[NA, :]
512 data.append((verts, color, color, 1.0))
514 for poly in lines_upper:
515 verts = project(poly, projection)[:, ::-1] * size + \
516 position[NA, :]
517 data.append((verts, 'none', edgecolor, linewidth))
518 return data
521def plot_beachball_mpl(
522 mt, axes,
523 beachball_type='deviatoric',
524 position=(0., 0.),
525 size=None,
526 zorder=0,
527 color_t='red',
528 color_p='white',
529 edgecolor='black',
530 linewidth=2,
531 alpha=1.0,
532 arcres=181,
533 decimation=1,
534 projection='lambert',
535 size_units='points',
536 view='top'):
538 '''
539 Plot beachball diagram to a Matplotlib plot
541 :param mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
542 array or sequence which can be converted into an MT object
543 :param beachball_type: ``'deviatoric'`` (default), ``'full'``, or ``'dc'``
544 :param position: position of the beachball in data coordinates
545 :param size: diameter of the beachball either in points or in data
546 coordinates, depending on the ``size_units`` setting
547 :param zorder: (passed through to matplotlib drawing functions)
548 :param color_t: color for compressional quadrants (default: ``'red'``)
549 :param color_p: color for extensive quadrants (default: ``'white'``)
550 :param edgecolor: color for lines (default: ``'black'``)
551 :param linewidth: linewidth in points (default: ``2``)
552 :param alpha: (passed through to matplotlib drawing functions)
553 :param projection: ``'lambert'`` (default), ``'stereographic'``, or
554 ``'orthographic'``
555 :param size_units: ``'points'`` (default) or ``'data'``, where the
556 latter causes the beachball to be projected in the plots data
557 coordinates (axes must have an aspect ratio of 1.0 or the
558 beachball will be shown distorted when using this).
559 :param view: View the beachball from ``top``, ``north``, ``south``,
560 ``east`` or ``west``. Useful for to show beachballs in cross-sections.
561 Default is ``top``.
562 '''
564 transform, position, size = choose_transform(
565 axes, size_units, position, size)
567 mt = deco_part(mt, beachball_type, view)
569 eig = mt.eigensystem()
570 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0:
571 raise BeachballError('eigenvalues are zero')
573 data = []
574 for (group, patches, patches_lower, patches_upper,
575 lines, lines_lower, lines_upper) in eig2gx(eig, arcres):
577 if group == 'P':
578 color = color_p
579 else:
580 color = color_t
582 # plot "upper" features for lower hemisphere, because coordinate system
583 # is NED
585 for poly in patches_upper:
586 verts = project(poly, projection)[:, ::-1] * size + position[NA, :]
587 if alpha == 1.0:
588 data.append(
589 (verts[::decimation], color, color, linewidth))
590 else:
591 data.append(
592 (verts[::decimation], color, 'none', 0.0))
594 for poly in lines_upper:
595 verts = project(poly, projection)[:, ::-1] * size + position[NA, :]
596 data.append(
597 (verts[::decimation], 'none', edgecolor, linewidth))
599 patches = []
600 for (path, facecolor, edgecolor, linewidth) in data:
601 patches.append(Polygon(
602 xy=path, facecolor=facecolor,
603 edgecolor=edgecolor,
604 linewidth=linewidth,
605 alpha=alpha))
607 collection = PatchCollection(
608 patches, zorder=zorder, transform=transform, match_original=True)
610 axes.add_artist(collection)
611 return collection
614def amplitudes_ned(mt, vecs):
615 ep, en, et, vp, vn, vt = mt.eigensystem()
616 to_e = num.vstack((vn, vt, vp))
617 vecs_e = num.dot(to_e, vecs.T).T
618 rtp = numpy_xyz2rtp(vecs_e)
619 atheta, aphi = rtp[:, 1], rtp[:, 2]
620 return ep * num.cos(atheta)**2 + (
621 en * num.cos(aphi)**2 + et * num.sin(aphi)**2) * num.sin(atheta)**2
624def amplitudes(mt, azimuths, takeoff_angles):
625 azimuths = num.asarray(azimuths, dtype=float)
626 takeoff_angles = num.asarray(takeoff_angles, dtype=float)
627 assert azimuths.size == takeoff_angles.size
628 rtps = num.vstack((num.ones(azimuths.size), takeoff_angles, azimuths)).T
629 vecs = numpy_rtp2xyz(rtps)
630 return amplitudes_ned(mt, vecs)
633def mts2amps(
634 mts,
635 projection,
636 beachball_type,
637 grid_resolution=200,
638 mask=True,
639 view='top'):
641 n_balls = len(mts)
642 nx = grid_resolution
643 ny = grid_resolution
645 x = num.linspace(-1., 1., nx)
646 y = num.linspace(-1., 1., ny)
648 vecs2 = num.zeros((nx * ny, 2), dtype=num.float64)
649 vecs2[:, 0] = num.tile(x, ny)
650 vecs2[:, 1] = num.repeat(y, nx)
652 ii_ok = vecs2[:, 0]**2 + vecs2[:, 1]**2 <= 1.0
653 amps = num_full(nx * ny, num.nan, dtype=num.float64)
655 amps[ii_ok] = 0.
656 vecs3_ok = inverse_project(vecs2[ii_ok, :], projection)
658 for mt in mts:
659 amps_ok = amplitudes_ned(deco_part(mt, beachball_type, view), vecs3_ok)
660 if mask:
661 amps_ok[amps_ok > 0] = 1.
662 amps_ok[amps_ok < 0] = 0.
664 amps[ii_ok] += amps_ok
666 return num.reshape(amps, (ny, nx)) / n_balls, x, y
669def plot_fuzzy_beachball_mpl_pixmap(
670 mts, axes,
671 best_mt=None,
672 beachball_type='deviatoric',
673 position=(0., 0.),
674 size=None,
675 zorder=0,
676 color_t='red',
677 color_p='white',
678 edgecolor='black',
679 best_color='red',
680 linewidth=2,
681 alpha=1.0,
682 projection='lambert',
683 size_units='data',
684 grid_resolution=200,
685 method='imshow',
686 view='top'):
687 '''
688 Plot fuzzy beachball from a list of given MomentTensors
690 :param mts: list of
691 :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
692 array or sequence which can be converted into an MT object
693 :param best_mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or
694 an array or sequence which can be converted into an MT object
695 of most likely or minimum misfit solution to extra highlight
696 :param best_color: mpl color for best MomentTensor edges,
697 polygons are not plotted
699 See plot_beachball_mpl for other arguments
700 '''
701 if size_units == 'points':
702 raise BeachballError(
703 'size_units="points" not supported in '
704 'plot_fuzzy_beachball_mpl_pixmap')
706 transform, position, size = choose_transform(
707 axes, size_units, position, size)
709 amps, x, y = mts2amps(
710 mts,
711 grid_resolution=grid_resolution,
712 projection=projection,
713 beachball_type=beachball_type,
714 mask=True,
715 view=view)
717 ncolors = 256
718 cmap = LinearSegmentedColormap.from_list(
719 'dummy', [color_p, color_t], N=ncolors)
721 levels = num.linspace(0, 1., ncolors)
722 if method == 'contourf':
723 axes.contourf(
724 position[0] + y * size, position[1] + x * size, amps.T,
725 levels=levels,
726 cmap=cmap,
727 transform=transform,
728 zorder=zorder,
729 alpha=alpha)
731 elif method == 'imshow':
732 axes.imshow(
733 amps.T,
734 extent=(
735 position[0] + y[0] * size,
736 position[0] + y[-1] * size,
737 position[1] - x[0] * size,
738 position[1] - x[-1] * size),
739 cmap=cmap,
740 transform=transform,
741 zorder=zorder-0.1,
742 alpha=alpha)
743 else:
744 assert False, 'invalid `method` argument'
746 # draw optimum edges
747 if best_mt is not None:
748 best_amps, bx, by = mts2amps(
749 [best_mt],
750 grid_resolution=grid_resolution,
751 projection=projection,
752 beachball_type=beachball_type,
753 mask=False)
755 axes.contour(
756 position[0] + by * size, position[1] + bx * size, best_amps.T,
757 levels=[0.],
758 colors=[best_color],
759 linewidths=linewidth,
760 transform=transform,
761 zorder=zorder,
762 alpha=alpha)
764 phi = num.linspace(0., 2 * PI, 361)
765 x = num.cos(phi)
766 y = num.sin(phi)
767 axes.plot(
768 position[0] + x * size, position[1] + y * size,
769 linewidth=linewidth,
770 color=edgecolor,
771 transform=transform,
772 zorder=zorder,
773 alpha=alpha)
776def plot_beachball_mpl_construction(
777 mt, axes,
778 show='patches',
779 beachball_type='deviatoric',
780 view='top'):
782 mt = deco_part(mt, beachball_type, view)
783 eig = mt.eigensystem()
785 for (group, patches, patches_lower, patches_upper,
786 lines, lines_lower, lines_upper) in eig2gx(eig):
788 if group == 'P':
789 color = 'blue'
790 lw = 1
791 else:
792 color = 'red'
793 lw = 1
795 if show == 'patches':
796 for poly in patches_upper:
797 px, py, pz = poly.T
798 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
800 if show == 'lines':
801 for poly in lines_upper:
802 px, py, pz = poly.T
803 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
806def plot_beachball_mpl_pixmap(
807 mt, axes,
808 beachball_type='deviatoric',
809 position=(0., 0.),
810 size=None,
811 zorder=0,
812 color_t='red',
813 color_p='white',
814 edgecolor='black',
815 linewidth=2,
816 alpha=1.0,
817 projection='lambert',
818 size_units='data',
819 view='top'):
821 if size_units == 'points':
822 raise BeachballError(
823 'size_units="points" not supported in plot_beachball_mpl_pixmap')
825 transform, position, size = choose_transform(
826 axes, size_units, position, size)
828 mt = deco_part(mt, beachball_type, view)
830 ep, en, et, vp, vn, vt = mt.eigensystem()
832 amps, x, y = mts2amps(
833 [mt], projection, beachball_type, grid_resolution=200, mask=False)
835 axes.contourf(
836 position[0] + y * size, position[1] + x * size, amps.T,
837 levels=[-num.inf, 0., num.inf],
838 colors=[color_p, color_t],
839 transform=transform,
840 zorder=zorder,
841 alpha=alpha)
843 axes.contour(
844 position[0] + y * size, position[1] + x * size, amps.T,
845 levels=[0.],
846 colors=[edgecolor],
847 linewidths=linewidth,
848 transform=transform,
849 zorder=zorder,
850 alpha=alpha)
852 phi = num.linspace(0., 2 * PI, 361)
853 x = num.cos(phi)
854 y = num.sin(phi)
855 axes.plot(
856 position[0] + x * size, position[1] + y * size,
857 linewidth=linewidth,
858 color=edgecolor,
859 transform=transform,
860 zorder=zorder,
861 alpha=alpha)
864if __name__ == '__main__':
865 import sys
866 import os
867 import matplotlib.pyplot as plt
868 from pyrocko import model
870 args = sys.argv[1:]
872 data = []
873 for iarg, arg in enumerate(args):
875 if os.path.exists(arg):
876 events = model.load_events(arg)
877 for ev in events:
878 if not ev.moment_tensor:
879 logger.warning('no moment tensor given for event')
880 continue
882 data.append((ev.name, ev.moment_tensor))
883 else:
884 vals = list(map(float, arg.split(',')))
885 mt = mtm.as_mt(vals)
886 data.append(('%i' % (iarg+1), mt))
888 n = len(data)
890 ncols = 1
891 while ncols**2 < n:
892 ncols += 1
894 nrows = ncols
896 fig = plt.figure()
897 axes = fig.add_subplot(1, 1, 1, aspect=1.)
898 axes.axison = False
899 axes.set_xlim(-0.05 - ncols, ncols + 0.05)
900 axes.set_ylim(-0.05 - nrows, nrows + 0.05)
902 for ibeach, (name, mt) in enumerate(data):
903 irow = ibeach // ncols
904 icol = ibeach % ncols
905 plot_beachball_mpl(
906 mt, axes,
907 position=(icol*2-ncols+1, -irow*2+nrows-1),
908 size_units='data')
910 axes.annotate(
911 name,
912 xy=(icol*2-ncols+1, -irow*2+nrows-2),
913 xycoords='data',
914 xytext=(0, 0),
915 textcoords='offset points',
916 verticalalignment='center',
917 horizontalalignment='center',
918 rotation=0.)
920 plt.show()