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 mts2amps(mts, projection, beachball_type, grid_resolution=200, mask=True,
622 view='top'):
624 n_balls = len(mts)
625 nx = grid_resolution
626 ny = grid_resolution
628 x = num.linspace(-1., 1., nx)
629 y = num.linspace(-1., 1., ny)
631 vecs2 = num.zeros((nx * ny, 2), dtype=num.float64)
632 vecs2[:, 0] = num.tile(x, ny)
633 vecs2[:, 1] = num.repeat(y, nx)
635 ii_ok = vecs2[:, 0]**2 + vecs2[:, 1]**2 <= 1.0
636 amps = num_full(nx * ny, num.nan, dtype=num.float64)
638 amps[ii_ok] = 0.
639 for mt in mts:
640 mt = deco_part(mt, beachball_type, view)
642 ep, en, et, vp, vn, vt = mt.eigensystem()
644 vecs3_ok = inverse_project(vecs2[ii_ok, :], projection)
646 to_e = num.vstack((vn, vt, vp))
648 vecs_e = num.dot(to_e, vecs3_ok.T).T
649 rtp = numpy_xyz2rtp(vecs_e)
651 atheta, aphi = rtp[:, 1], rtp[:, 2]
652 amps_ok = ep * num.cos(atheta)**2 + (
653 en * num.cos(aphi)**2 + et * num.sin(aphi)**2) * num.sin(atheta)**2
655 if mask:
656 amps_ok[amps_ok > 0] = 1.
657 amps_ok[amps_ok < 0] = 0.
659 amps[ii_ok] += amps_ok
661 return num.reshape(amps, (ny, nx)) / n_balls, x, y
664def plot_fuzzy_beachball_mpl_pixmap(
665 mts, axes,
666 best_mt=None,
667 beachball_type='deviatoric',
668 position=(0., 0.),
669 size=None,
670 zorder=0,
671 color_t='red',
672 color_p='white',
673 edgecolor='black',
674 best_color='red',
675 linewidth=2,
676 alpha=1.0,
677 projection='lambert',
678 size_units='data',
679 grid_resolution=200,
680 method='imshow',
681 view='top'):
682 '''
683 Plot fuzzy beachball from a list of given MomentTensors
685 :param mts: list of
686 :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
687 array or sequence which can be converted into an MT object
688 :param best_mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or
689 an array or sequence which can be converted into an MT object
690 of most likely or minimum misfit solution to extra highlight
691 :param best_color: mpl color for best MomentTensor edges,
692 polygons are not plotted
694 See plot_beachball_mpl for other arguments
695 '''
696 if size_units == 'points':
697 raise BeachballError(
698 'size_units="points" not supported in '
699 'plot_fuzzy_beachball_mpl_pixmap')
701 transform, position, size = choose_transform(
702 axes, size_units, position, size)
704 amps, x, y = mts2amps(
705 mts,
706 grid_resolution=grid_resolution,
707 projection=projection,
708 beachball_type=beachball_type,
709 mask=True,
710 view=view)
712 ncolors = 256
713 cmap = LinearSegmentedColormap.from_list(
714 'dummy', [color_p, color_t], N=ncolors)
716 levels = num.linspace(0, 1., ncolors)
717 if method == 'contourf':
718 axes.contourf(
719 position[0] + y * size, position[1] + x * size, amps.T,
720 levels=levels,
721 cmap=cmap,
722 transform=transform,
723 zorder=zorder,
724 alpha=alpha)
726 elif method == 'imshow':
727 axes.imshow(
728 amps.T,
729 extent=(
730 position[0] + y[0] * size,
731 position[0] + y[-1] * size,
732 position[1] - x[0] * size,
733 position[1] - x[-1] * size),
734 cmap=cmap,
735 transform=transform,
736 zorder=zorder-0.1,
737 alpha=alpha)
738 else:
739 assert False, 'invalid `method` argument'
741 # draw optimum edges
742 if best_mt is not None:
743 best_amps, bx, by = mts2amps(
744 [best_mt],
745 grid_resolution=grid_resolution,
746 projection=projection,
747 beachball_type=beachball_type,
748 mask=False)
750 axes.contour(
751 position[0] + by * size, position[1] + bx * size, best_amps.T,
752 levels=[0.],
753 colors=[best_color],
754 linewidths=linewidth,
755 transform=transform,
756 zorder=zorder,
757 alpha=alpha)
759 phi = num.linspace(0., 2 * PI, 361)
760 x = num.cos(phi)
761 y = num.sin(phi)
762 axes.plot(
763 position[0] + x * size, position[1] + y * size,
764 linewidth=linewidth,
765 color=edgecolor,
766 transform=transform,
767 zorder=zorder,
768 alpha=alpha)
771def plot_beachball_mpl_construction(
772 mt, axes,
773 show='patches',
774 beachball_type='deviatoric',
775 view='top'):
777 mt = deco_part(mt, beachball_type, view)
778 eig = mt.eigensystem()
780 for (group, patches, patches_lower, patches_upper,
781 lines, lines_lower, lines_upper) in eig2gx(eig):
783 if group == 'P':
784 color = 'blue'
785 lw = 1
786 else:
787 color = 'red'
788 lw = 1
790 if show == 'patches':
791 for poly in patches_upper:
792 px, py, pz = poly.T
793 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
795 if show == 'lines':
796 for poly in lines_upper:
797 px, py, pz = poly.T
798 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
801def plot_beachball_mpl_pixmap(
802 mt, axes,
803 beachball_type='deviatoric',
804 position=(0., 0.),
805 size=None,
806 zorder=0,
807 color_t='red',
808 color_p='white',
809 edgecolor='black',
810 linewidth=2,
811 alpha=1.0,
812 projection='lambert',
813 size_units='data',
814 view='top'):
816 if size_units == 'points':
817 raise BeachballError(
818 'size_units="points" not supported in plot_beachball_mpl_pixmap')
820 transform, position, size = choose_transform(
821 axes, size_units, position, size)
823 mt = deco_part(mt, beachball_type, view)
825 ep, en, et, vp, vn, vt = mt.eigensystem()
827 amps, x, y = mts2amps(
828 [mt], projection, beachball_type, grid_resolution=200, mask=False)
830 axes.contourf(
831 position[0] + y * size, position[1] + x * size, amps.T,
832 levels=[-num.inf, 0., num.inf],
833 colors=[color_p, color_t],
834 transform=transform,
835 zorder=zorder,
836 alpha=alpha)
838 axes.contour(
839 position[0] + y * size, position[1] + x * size, amps.T,
840 levels=[0.],
841 colors=[edgecolor],
842 linewidths=linewidth,
843 transform=transform,
844 zorder=zorder,
845 alpha=alpha)
847 phi = num.linspace(0., 2 * PI, 361)
848 x = num.cos(phi)
849 y = num.sin(phi)
850 axes.plot(
851 position[0] + x * size, position[1] + y * size,
852 linewidth=linewidth,
853 color=edgecolor,
854 transform=transform,
855 zorder=zorder,
856 alpha=alpha)
859if __name__ == '__main__':
860 import sys
861 import os
862 import matplotlib.pyplot as plt
863 from pyrocko import model
865 args = sys.argv[1:]
867 data = []
868 for iarg, arg in enumerate(args):
870 if os.path.exists(arg):
871 events = model.load_events(arg)
872 for ev in events:
873 if not ev.moment_tensor:
874 logger.warning('no moment tensor given for event')
875 continue
877 data.append((ev.name, ev.moment_tensor))
878 else:
879 vals = list(map(float, arg.split(',')))
880 mt = mtm.as_mt(vals)
881 data.append(('%i' % (iarg+1), mt))
883 n = len(data)
885 ncols = 1
886 while ncols**2 < n:
887 ncols += 1
889 nrows = ncols
891 fig = plt.figure()
892 axes = fig.add_subplot(1, 1, 1, aspect=1.)
893 axes.axison = False
894 axes.set_xlim(-0.05 - ncols, ncols + 0.05)
895 axes.set_ylim(-0.05 - nrows, nrows + 0.05)
897 for ibeach, (name, mt) in enumerate(data):
898 irow = ibeach // ncols
899 icol = ibeach % ncols
900 plot_beachball_mpl(
901 mt, axes,
902 position=(icol*2-ncols+1, -irow*2+nrows-1),
903 size_units='data')
905 axes.annotate(
906 name,
907 xy=(icol*2-ncols+1, -irow*2+nrows-2),
908 xycoords='data',
909 xytext=(0, 0),
910 textcoords='offset points',
911 verticalalignment='center',
912 horizontalalignment='center',
913 rotation=0.)
915 plt.show()