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