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
21d2r = num.pi / 180.
24def view_rotation(strike, dip):
25 return mtm.euler_to_matrix(
26 dip*d2r, strike*d2r, -90.*d2r)
29_view_south = view_rotation(90., 90.)
30_view_north = view_rotation(-90., 90.)
31_view_east = view_rotation(0., 90.)
32_view_west = view_rotation(180., 90.)
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 mt = mtm.as_mt(mt)
418 if isinstance(view, str):
419 if view == 'top':
420 pass
421 elif view == 'north':
422 mt = mt.rotated(_view_north)
423 elif view == 'south':
424 mt = mt.rotated(_view_south)
425 elif view == 'east':
426 mt = mt.rotated(_view_east)
427 elif view == 'west':
428 mt = mt.rotated(_view_west)
429 elif isinstance(view, tuple):
430 mt = mt.rotated(view_rotation(*view))
431 else:
432 raise BeachballError(
433 'Invaild argument for `view`. Allowed values are "top", "north", '
434 '"south", "east", "west" or a tuple of angles `(strike, dip)` '
435 'orienting the view plane.')
437 if mt_type == 'full':
438 return mt
440 res = mt.standard_decomposition()
441 m = dict(
442 dc=res[1][2],
443 deviatoric=res[3][2])[mt_type]
445 return mtm.MomentTensor(m=m)
448def choose_transform(axes, size_units, position, size):
450 if size_units == 'points':
451 transform = _FixedPointOffsetTransform(
452 axes.transData,
453 axes.figure.dpi_scale_trans,
454 position)
456 if size is None:
457 size = 12.
459 size = size * 0.5 / 72.
460 position = (0., 0.)
462 elif size_units == 'data':
463 transform = axes.transData
465 if size is None:
466 size = 1.0
468 size = size * 0.5
470 elif size_units == 'axes':
471 transform = axes.transAxes
472 if size is None:
473 size = 1.
475 size = size * .5
477 else:
478 raise BeachballError(
479 'invalid argument for size_units: %s' % size_units)
481 position = num.asarray(position, dtype=num.float64)
483 return transform, position, size
486def mt2beachball(
487 mt,
488 beachball_type='deviatoric',
489 position=(0., 0.),
490 size=None,
491 color_t='red',
492 color_p='white',
493 edgecolor='black',
494 linewidth=2,
495 projection='lambert',
496 view='top'):
498 position = num.asarray(position, dtype=num.float64)
499 size = size or 1
500 mt = deco_part(mt, beachball_type, view)
502 eig = mt.eigensystem()
503 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0:
504 raise BeachballError('eigenvalues are zero')
506 data = []
507 for (group, patches, patches_lower, patches_upper,
508 lines, lines_lower, lines_upper) in eig2gx(eig):
510 if group == 'P':
511 color = color_p
512 else:
513 color = color_t
515 for poly in patches_upper:
516 verts = project(poly, projection)[:, ::-1] * size + \
517 position[NA, :]
518 data.append((verts, color, color, 1.0))
520 for poly in lines_upper:
521 verts = project(poly, projection)[:, ::-1] * size + \
522 position[NA, :]
523 data.append((verts, 'none', edgecolor, linewidth))
524 return data
527def plot_beachball_mpl(
528 mt, axes,
529 beachball_type='deviatoric',
530 position=(0., 0.),
531 size=None,
532 zorder=0,
533 color_t='red',
534 color_p='white',
535 edgecolor='black',
536 linewidth=2,
537 alpha=1.0,
538 arcres=181,
539 decimation=1,
540 projection='lambert',
541 size_units='points',
542 view='top'):
544 '''
545 Plot beachball diagram to a Matplotlib plot
547 :param mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
548 array or sequence which can be converted into an MT object
549 :param beachball_type: ``'deviatoric'`` (default), ``'full'``, or ``'dc'``
550 :param position: position of the beachball in data coordinates
551 :param size: diameter of the beachball either in points or in data
552 coordinates, depending on the ``size_units`` setting
553 :param zorder: (passed through to matplotlib drawing functions)
554 :param color_t: color for compressional quadrants (default: ``'red'``)
555 :param color_p: color for extensive quadrants (default: ``'white'``)
556 :param edgecolor: color for lines (default: ``'black'``)
557 :param linewidth: linewidth in points (default: ``2``)
558 :param alpha: (passed through to matplotlib drawing functions)
559 :param projection: ``'lambert'`` (default), ``'stereographic'``, or
560 ``'orthographic'``
561 :param size_units: ``'points'`` (default) or ``'data'``, where the
562 latter causes the beachball to be projected in the plots data
563 coordinates (axes must have an aspect ratio of 1.0 or the
564 beachball will be shown distorted when using this).
565 :param view: View the beachball from ``'top'``, ``'north'``, ``'south'``,
566 ``'east'`` or ``'west'``, or project onto plane given by
567 ``(strike, dip)``. Useful to show beachballs in cross-sections.
568 Default is ``'top'``.
569 '''
571 transform, position, size = choose_transform(
572 axes, size_units, position, size)
574 mt = deco_part(mt, beachball_type, view)
576 eig = mt.eigensystem()
577 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0:
578 raise BeachballError('eigenvalues are zero')
580 data = []
581 for (group, patches, patches_lower, patches_upper,
582 lines, lines_lower, lines_upper) in eig2gx(eig, arcres):
584 if group == 'P':
585 color = color_p
586 else:
587 color = color_t
589 # plot "upper" features for lower hemisphere, because coordinate system
590 # is NED
592 for poly in patches_upper:
593 verts = project(poly, projection)[:, ::-1] * size + position[NA, :]
594 if alpha == 1.0:
595 data.append(
596 (verts[::decimation], color, color, linewidth))
597 else:
598 data.append(
599 (verts[::decimation], color, 'none', 0.0))
601 for poly in lines_upper:
602 verts = project(poly, projection)[:, ::-1] * size + position[NA, :]
603 data.append(
604 (verts[::decimation], 'none', edgecolor, linewidth))
606 patches = []
607 for (path, facecolor, edgecolor, linewidth) in data:
608 patches.append(Polygon(
609 xy=path, facecolor=facecolor,
610 edgecolor=edgecolor,
611 linewidth=linewidth,
612 alpha=alpha))
614 collection = PatchCollection(
615 patches, zorder=zorder, transform=transform, match_original=True)
617 axes.add_artist(collection)
618 return collection
621def amplitudes_ned(mt, vecs):
622 ep, en, et, vp, vn, vt = mt.eigensystem()
623 to_e = num.vstack((vn, vt, vp))
624 vecs_e = num.dot(to_e, vecs.T).T
625 rtp = numpy_xyz2rtp(vecs_e)
626 atheta, aphi = rtp[:, 1], rtp[:, 2]
627 return ep * num.cos(atheta)**2 + (
628 en * num.cos(aphi)**2 + et * num.sin(aphi)**2) * num.sin(atheta)**2
631def amplitudes(mt, azimuths, takeoff_angles):
632 azimuths = num.asarray(azimuths, dtype=float)
633 takeoff_angles = num.asarray(takeoff_angles, dtype=float)
634 assert azimuths.size == takeoff_angles.size
635 rtps = num.vstack((num.ones(azimuths.size), takeoff_angles, azimuths)).T
636 vecs = numpy_rtp2xyz(rtps)
637 return amplitudes_ned(mt, vecs)
640def mts2amps(
641 mts,
642 projection,
643 beachball_type,
644 grid_resolution=200,
645 mask=True,
646 view='top'):
648 n_balls = len(mts)
649 nx = grid_resolution
650 ny = grid_resolution
652 x = num.linspace(-1., 1., nx)
653 y = num.linspace(-1., 1., ny)
655 vecs2 = num.zeros((nx * ny, 2), dtype=num.float64)
656 vecs2[:, 0] = num.tile(x, ny)
657 vecs2[:, 1] = num.repeat(y, nx)
659 ii_ok = vecs2[:, 0]**2 + vecs2[:, 1]**2 <= 1.0
660 amps = num_full(nx * ny, num.nan, dtype=num.float64)
662 amps[ii_ok] = 0.
663 vecs3_ok = inverse_project(vecs2[ii_ok, :], projection)
665 for mt in mts:
666 amps_ok = amplitudes_ned(deco_part(mt, beachball_type, view), vecs3_ok)
667 if mask:
668 amps_ok[amps_ok > 0] = 1.
669 amps_ok[amps_ok < 0] = 0.
671 amps[ii_ok] += amps_ok
673 return num.reshape(amps, (ny, nx)) / n_balls, x, y
676def plot_fuzzy_beachball_mpl_pixmap(
677 mts, axes,
678 best_mt=None,
679 beachball_type='deviatoric',
680 position=(0., 0.),
681 size=None,
682 zorder=0,
683 color_t='red',
684 color_p='white',
685 edgecolor='black',
686 best_color='red',
687 linewidth=2,
688 alpha=1.0,
689 projection='lambert',
690 size_units='data',
691 grid_resolution=200,
692 method='imshow',
693 view='top'):
694 '''
695 Plot fuzzy beachball from a list of given MomentTensors
697 :param mts: list of
698 :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
699 array or sequence which can be converted into an MT object
700 :param best_mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or
701 an array or sequence which can be converted into an MT object
702 of most likely or minimum misfit solution to extra highlight
703 :param best_color: mpl color for best MomentTensor edges,
704 polygons are not plotted
706 See plot_beachball_mpl for other arguments
707 '''
708 if size_units == 'points':
709 raise BeachballError(
710 'size_units="points" not supported in '
711 'plot_fuzzy_beachball_mpl_pixmap')
713 transform, position, size = choose_transform(
714 axes, size_units, position, size)
716 amps, x, y = mts2amps(
717 mts,
718 grid_resolution=grid_resolution,
719 projection=projection,
720 beachball_type=beachball_type,
721 mask=True,
722 view=view)
724 ncolors = 256
725 cmap = LinearSegmentedColormap.from_list(
726 'dummy', [color_p, color_t], N=ncolors)
728 levels = num.linspace(0, 1., ncolors)
729 if method == 'contourf':
730 axes.contourf(
731 position[0] + y * size, position[1] + x * size, amps.T,
732 levels=levels,
733 cmap=cmap,
734 transform=transform,
735 zorder=zorder,
736 alpha=alpha)
738 elif method == 'imshow':
739 axes.imshow(
740 amps.T,
741 extent=(
742 position[0] + y[0] * size,
743 position[0] + y[-1] * size,
744 position[1] - x[0] * size,
745 position[1] - x[-1] * size),
746 cmap=cmap,
747 transform=transform,
748 zorder=zorder-0.1,
749 alpha=alpha)
750 else:
751 assert False, 'invalid `method` argument'
753 # draw optimum edges
754 if best_mt is not None:
755 best_amps, bx, by = mts2amps(
756 [best_mt],
757 grid_resolution=grid_resolution,
758 projection=projection,
759 beachball_type=beachball_type,
760 mask=False)
762 axes.contour(
763 position[0] + by * size, position[1] + bx * size, best_amps.T,
764 levels=[0.],
765 colors=[best_color],
766 linewidths=linewidth,
767 transform=transform,
768 zorder=zorder,
769 alpha=alpha)
771 phi = num.linspace(0., 2 * PI, 361)
772 x = num.cos(phi)
773 y = num.sin(phi)
774 axes.plot(
775 position[0] + x * size, position[1] + y * size,
776 linewidth=linewidth,
777 color=edgecolor,
778 transform=transform,
779 zorder=zorder,
780 alpha=alpha)
783def plot_beachball_mpl_construction(
784 mt, axes,
785 show='patches',
786 beachball_type='deviatoric',
787 view='top'):
789 mt = deco_part(mt, beachball_type, view)
790 eig = mt.eigensystem()
792 for (group, patches, patches_lower, patches_upper,
793 lines, lines_lower, lines_upper) in eig2gx(eig):
795 if group == 'P':
796 color = 'blue'
797 lw = 1
798 else:
799 color = 'red'
800 lw = 1
802 if show == 'patches':
803 for poly in patches_upper:
804 px, py, pz = poly.T
805 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
807 if show == 'lines':
808 for poly in lines_upper:
809 px, py, pz = poly.T
810 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
813def plot_beachball_mpl_pixmap(
814 mt, axes,
815 beachball_type='deviatoric',
816 position=(0., 0.),
817 size=None,
818 zorder=0,
819 color_t='red',
820 color_p='white',
821 edgecolor='black',
822 linewidth=2,
823 alpha=1.0,
824 projection='lambert',
825 size_units='data',
826 view='top'):
828 if size_units == 'points':
829 raise BeachballError(
830 'size_units="points" not supported in plot_beachball_mpl_pixmap')
832 transform, position, size = choose_transform(
833 axes, size_units, position, size)
835 mt = deco_part(mt, beachball_type, view)
837 ep, en, et, vp, vn, vt = mt.eigensystem()
839 amps, x, y = mts2amps(
840 [mt], projection, beachball_type, grid_resolution=200, mask=False)
842 axes.contourf(
843 position[0] + y * size, position[1] + x * size, amps.T,
844 levels=[-num.inf, 0., num.inf],
845 colors=[color_p, color_t],
846 transform=transform,
847 zorder=zorder,
848 alpha=alpha)
850 axes.contour(
851 position[0] + y * size, position[1] + x * size, amps.T,
852 levels=[0.],
853 colors=[edgecolor],
854 linewidths=linewidth,
855 transform=transform,
856 zorder=zorder,
857 alpha=alpha)
859 phi = num.linspace(0., 2 * PI, 361)
860 x = num.cos(phi)
861 y = num.sin(phi)
862 axes.plot(
863 position[0] + x * size, position[1] + y * size,
864 linewidth=linewidth,
865 color=edgecolor,
866 transform=transform,
867 zorder=zorder,
868 alpha=alpha)
871if __name__ == '__main__':
872 import sys
873 import os
874 import matplotlib.pyplot as plt
875 from pyrocko import model
877 args = sys.argv[1:]
879 data = []
880 for iarg, arg in enumerate(args):
882 if os.path.exists(arg):
883 events = model.load_events(arg)
884 for ev in events:
885 if not ev.moment_tensor:
886 logger.warning('no moment tensor given for event')
887 continue
889 data.append((ev.name, ev.moment_tensor))
890 else:
891 vals = list(map(float, arg.split(',')))
892 mt = mtm.as_mt(vals)
893 data.append(('%i' % (iarg+1), mt))
895 n = len(data)
897 ncols = 1
898 while ncols**2 < n:
899 ncols += 1
901 nrows = ncols
903 fig = plt.figure()
904 axes = fig.add_subplot(1, 1, 1, aspect=1.)
905 axes.axison = False
906 axes.set_xlim(-0.05 - ncols, ncols + 0.05)
907 axes.set_ylim(-0.05 - nrows, nrows + 0.05)
909 for ibeach, (name, mt) in enumerate(data):
910 irow = ibeach // ncols
911 icol = ibeach % ncols
912 plot_beachball_mpl(
913 mt, axes,
914 position=(icol*2-ncols+1, -irow*2+nrows-1),
915 size_units='data')
917 axes.annotate(
918 name,
919 xy=(icol*2-ncols+1, -irow*2+nrows-2),
920 xycoords='data',
921 xytext=(0, 0),
922 textcoords='offset points',
923 verticalalignment='center',
924 horizontalalignment='center',
925 rotation=0.)
927 plt.show()