1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5# python 2/3
6from __future__ import absolute_import
8from math import pi as PI
9import logging
10import numpy as num
12from matplotlib.collections import PatchCollection
13from matplotlib.patches import Polygon
14from matplotlib.transforms import Transform
15from matplotlib.colors import LinearSegmentedColormap
17from pyrocko import moment_tensor as mtm
18from pyrocko.util import num_full
20logger = logging.getLogger('pyrocko.plot.beachball')
22NA = num.newaxis
24_view_south = num.array([[0, 0, -1],
25 [0, 1, 0],
26 [1, 0, 0]])
28_view_north = _view_south.T
30_view_east = num.array([[1, 0, 0],
31 [0, 0, -1],
32 [0, 1, 0]])
34_view_west = _view_east.T
37class BeachballError(Exception):
38 pass
41class FixedPointOffsetTransform(Transform):
42 def __init__(self, trans, dpi_scale_trans, fixed_point):
43 Transform.__init__(self)
44 self.input_dims = self.output_dims = 2
45 self.has_inverse = False
46 self.trans = trans
47 self.dpi_scale_trans = dpi_scale_trans
48 self.fixed_point = num.asarray(fixed_point, dtype=num.float64)
50 def transform_non_affine(self, values):
51 fp = self.trans.transform(self.fixed_point)
52 return fp + self.dpi_scale_trans.transform(values)
55def vnorm(points):
56 return num.sqrt(num.sum(points**2, axis=1))
59def clean_poly(points):
60 if not num.all(points[0, :] == points[-1, :]):
61 points = num.vstack((points, points[0:1, :]))
63 dupl = num.concatenate(
64 (num.all(points[1:, :] == points[:-1, :], axis=1), [False]))
65 points = points[num.logical_not(dupl)]
66 return points
69def close_poly(points):
70 if not num.all(points[0, :] == points[-1, :]):
71 points = num.vstack((points, points[0:1, :]))
73 return points
76def circulation(points, axis):
77 # assert num.all(points[:, axis] >= 0.0) or num.all(points[:, axis] <= 0.0)
79 points2 = points[:, ((axis+2) % 3, (axis+1) % 3)].copy()
80 points2 *= 1.0 / num.sqrt(1.0 + num.abs(points[:, axis]))[:, num.newaxis]
82 result = -num.sum(
83 (points2[1:, 0] - points2[:-1, 0]) *
84 (points2[1:, 1] + points2[:-1, 1]))
86 result -= (points2[0, 0] - points2[-1, 0]) \
87 * (points2[0, 1] + points2[-1, 1])
88 return result
91def spoly_cut(l_points, axis=0, nonsimple=True, arcres=181):
92 dphi = 2.*PI / (2*arcres)
94 # cut sub-polygons and gather crossing point information
95 crossings = []
96 snippets = {}
97 for ipath, points in enumerate(l_points):
98 if not num.all(points[0, :] == points[-1, :]):
99 points = num.vstack((points, points[0:1, :]))
101 # get upward crossing points
102 iup = num.where(num.logical_and(points[:-1, axis] <= 0.,
103 points[1:, axis] > 0.))[0]
104 aup = - points[iup, axis] / (points[iup+1, axis] - points[iup, axis])
105 pup = points[iup, :] + aup[:, num.newaxis] * (points[iup+1, :] -
106 points[iup, :])
107 phiup = num.arctan2(pup[:, (axis+2) % 3], pup[:, (axis+1) % 3])
109 for i in range(len(iup)):
110 crossings.append((phiup[i], ipath, iup[i], 1, pup[i], [1, -1]))
112 # get downward crossing points
113 idown = num.where(num.logical_and(points[:-1, axis] > 0.,
114 points[1:, axis] <= 0.))[0]
115 adown = - points[idown+1, axis] / (points[idown, axis] -
116 points[idown+1, axis])
117 pdown = points[idown+1, :] + adown[:, num.newaxis] * (
118 points[idown, :] - points[idown+1, :])
119 phidown = num.arctan2(pdown[:, (axis+2) % 3], pdown[:, (axis+1) % 3])
121 for i in range(idown.size):
122 crossings.append(
123 (phidown[i], ipath, idown[i], -1, pdown[i], [1, -1]))
125 icuts = num.sort(num.concatenate((iup, idown)))
127 for i in range(icuts.size-1):
128 snippets[ipath, icuts[i]] = (
129 ipath, icuts[i+1], points[icuts[i]+1:icuts[i+1]+1])
131 if icuts.size:
132 points_last = num.concatenate((
133 points[icuts[-1]+1:],
134 points[:icuts[0]+1]))
136 snippets[ipath, icuts[-1]] = (ipath, icuts[0], points_last)
137 else:
138 snippets[ipath, 0] = (ipath, 0, points)
140 crossings.sort()
142 # assemble new sub-polygons
143 current = snippets.pop(list(snippets.keys())[0])
144 outs = [[]]
145 while True:
146 outs[-1].append(current[2])
147 for i, c1 in enumerate(crossings):
148 if c1[1:3] == current[:2]:
149 direction = -1 * c1[3]
150 break
151 else:
152 if not snippets:
153 break
154 current = snippets.pop(list(snippets.keys())[0])
155 outs.append([])
156 continue
158 while True:
159 i = (i + direction) % len(crossings)
160 if crossings[i][3] == direction and direction in crossings[i][-1]:
161 break
163 c2 = crossings[i]
164 c2[-1].remove(direction)
166 phi1 = c1[0]
167 phi2 = c2[0]
168 if direction == 1:
169 if phi1 > phi2:
170 phi2 += PI * 2.
172 if direction == -1:
173 if phi1 < phi2:
174 phi2 -= PI * 2.
176 n = int(abs(phi2 - phi1) / dphi) + 2
178 phis = num.linspace(phi1, phi2, n)
179 cpoints = num.zeros((n, 3))
180 cpoints[:, (axis+1) % 3] = num.cos(phis)
181 cpoints[:, (axis+2) % 3] = num.sin(phis)
182 cpoints[:, axis] = 0.0
184 outs[-1].append(cpoints)
186 try:
187 current = snippets[c2[1:3]]
188 del snippets[c2[1:3]]
190 except KeyError:
191 if not snippets:
192 break
194 current = snippets.pop(list(snippets.keys())[0])
195 outs.append([])
197 # separate hemispheres, force polygons closed, remove duplicate points
198 # remove polygons with less than 3 points (4, when counting repeated
199 # endpoint)
201 outs_upper = []
202 outs_lower = []
203 for out in outs:
204 if out:
205 out = clean_poly(num.vstack(out))
206 if out.shape[0] >= 4:
207 if num.sum(out[:, axis]) > 0.0:
208 outs_upper.append(out)
209 else:
210 outs_lower.append(out)
212 if nonsimple and (
213 len(crossings) == 0 or
214 len(outs_upper) == 0 or
215 len(outs_lower) == 0):
217 # check if we are cutting between holes
218 need_divider = False
219 if outs_upper:
220 candis = sorted(
221 outs_upper, key=lambda out: num.min(out[:, axis]))
223 if circulation(candis[0], axis) > 0.0:
224 need_divider = True
226 if outs_lower:
227 candis = sorted(
228 outs_lower, key=lambda out: num.max(out[:, axis]))
230 if circulation(candis[0], axis) < 0.0:
231 need_divider = True
233 if need_divider:
234 phi1 = 0.
235 phi2 = PI*2.
236 n = int(abs(phi2 - phi1) / dphi) + 2
238 phis = num.linspace(phi1, phi2, n)
239 cpoints = num.zeros((n, 3))
240 cpoints[:, (axis+1) % 3] = num.cos(phis)
241 cpoints[:, (axis+2) % 3] = num.sin(phis)
242 cpoints[:, axis] = 0.0
244 outs_upper.append(cpoints)
245 outs_lower.append(cpoints[::-1, :])
247 return outs_lower, outs_upper
250def numpy_rtp2xyz(rtp):
251 r = rtp[:, 0]
252 theta = rtp[:, 1]
253 phi = rtp[:, 2]
254 vecs = num.empty(rtp.shape, dtype=num.float64)
255 vecs[:, 0] = r*num.sin(theta)*num.cos(phi)
256 vecs[:, 1] = r*num.sin(theta)*num.sin(phi)
257 vecs[:, 2] = r*num.cos(theta)
258 return vecs
261def numpy_xyz2rtp(xyz):
262 x, y, z = xyz[:, 0], xyz[:, 1], xyz[:, 2]
263 vecs = num.empty(xyz.shape, dtype=num.float64)
264 vecs[:, 0] = num.sqrt(x**2+y**2+z**2)
265 vecs[:, 1] = num.arctan2(num.sqrt(x**2+y**2), z)
266 vecs[:, 2] = num.arctan2(y, x)
267 return vecs
270def circle_points(aphi, sign=1.0):
271 vecs = num.empty((aphi.size, 3), dtype=num.float64)
272 vecs[:, 0] = num.cos(sign*aphi)
273 vecs[:, 1] = num.sin(sign*aphi)
274 vecs[:, 2] = 0.0
275 return vecs
278def eig2gx(eig, arcres=181):
279 aphi = num.linspace(0., 2.*PI, arcres)
280 ep, en, et, vp, vn, vt = eig
282 mt_sign = num.sign(ep + en + et)
284 groups = []
285 for (pt_name, pt_sign) in [('P', -1.), ('T', 1.)]:
286 patches = []
287 patches_lower = []
288 patches_upper = []
289 lines = []
290 lines_lower = []
291 lines_upper = []
292 for iperm, (va, vb, vc, ea, eb, ec) in enumerate([
293 (vp, vn, vt, ep, en, et),
294 (vt, vp, vn, et, ep, en)]): # (vn, vt, vp, en, et, ep)]):
296 perm_sign = [-1.0, 1.0][iperm]
297 to_e = num.vstack((vb, vc, va))
298 from_e = to_e.T
300 poly_es = []
301 polys = []
302 for sign in (-1., 1.):
303 xphi = perm_sign*pt_sign*sign*aphi
304 denom = eb*num.cos(xphi)**2 + ec*num.sin(xphi)**2
305 if num.any(denom == 0.):
306 continue
308 Y = -ea/denom
309 if num.any(Y < 0.):
310 continue
312 xtheta = num.arctan(num.sqrt(Y))
313 rtp = num.empty(xphi.shape+(3,), dtype=num.float64)
314 rtp[:, 0] = 1.
315 if sign > 0:
316 rtp[:, 1] = xtheta
317 else:
318 rtp[:, 1] = PI - xtheta
320 rtp[:, 2] = xphi
321 poly_e = numpy_rtp2xyz(rtp)
322 poly = num.dot(from_e, poly_e.T).T
323 poly[:, 2] -= 0.001
325 poly_es.append(poly_e)
326 polys.append(poly)
328 if polys:
329 polys_lower, polys_upper = spoly_cut(polys, 2, arcres=arcres)
330 lines.extend(polys)
331 lines_lower.extend(polys_lower)
332 lines_upper.extend(polys_upper)
334 if poly_es:
335 for aa in spoly_cut(poly_es, 0, arcres=arcres):
336 for bb in spoly_cut(aa, 1, arcres=arcres):
337 for cc in spoly_cut(bb, 2, arcres=arcres):
338 for poly_e in cc:
339 poly = num.dot(from_e, poly_e.T).T
340 poly[:, 2] -= 0.001
341 polys_lower, polys_upper = spoly_cut(
342 [poly], 2, nonsimple=False, arcres=arcres)
344 patches.append(poly)
345 patches_lower.extend(polys_lower)
346 patches_upper.extend(polys_upper)
348 if not patches:
349 if mt_sign * pt_sign == 1.:
350 patches_lower.append(circle_points(aphi, -1.0))
351 patches_upper.append(circle_points(aphi, 1.0))
352 lines_lower.append(circle_points(aphi, -1.0))
353 lines_upper.append(circle_points(aphi, 1.0))
355 groups.append((
356 pt_name,
357 patches, patches_lower, patches_upper,
358 lines, lines_lower, lines_upper))
360 return groups
363def extr(points):
364 pmean = num.mean(points, axis=0)
365 return points + pmean*0.05
368def draw_eigenvectors_mpl(eig, axes):
369 vp, vn, vt = eig[3:]
370 for lab, v in [('P', vp), ('N', vn), ('T', vt)]:
371 sign = num.sign(v[2]) + (v[2] == 0.0)
372 axes.plot(sign*v[1], sign*v[0], 'o', color='black')
373 axes.text(sign*v[1], sign*v[0], ' '+lab)
376def project(points, projection='lambert'):
377 points_out = points[:, :2].copy()
378 if projection == 'lambert':
379 factor = 1.0 / num.sqrt(1.0 + points[:, 2])
380 elif projection == 'stereographic':
381 factor = 1.0 / (1.0 + points[:, 2])
382 elif projection == 'orthographic':
383 factor = None
384 else:
385 raise BeachballError(
386 'invalid argument for projection: %s' % projection)
388 if factor is not None:
389 points_out *= factor[:, num.newaxis]
391 return points_out
394def inverse_project(points, projection='lambert'):
395 points_out = num.zeros((points.shape[0], 3))
397 rsqr = points[:, 0]**2 + points[:, 1]**2
398 if projection == 'lambert':
399 points_out[:, 2] = 1.0 - rsqr
400 points_out[:, 1] = num.sqrt(2.0 - rsqr) * points[:, 1]
401 points_out[:, 0] = num.sqrt(2.0 - rsqr) * points[:, 0]
402 elif projection == 'stereographic':
403 points_out[:, 2] = - (rsqr - 1.0) / (rsqr + 1.0)
404 points_out[:, 1] = 2.0 * points[:, 1] / (rsqr + 1.0)
405 points_out[:, 0] = 2.0 * points[:, 0] / (rsqr + 1.0)
406 elif projection == 'orthographic':
407 points_out[:, 2] = num.sqrt(num.maximum(1.0 - rsqr, 0.0))
408 points_out[:, 1] = points[:, 1]
409 points_out[:, 0] = points[:, 0]
410 else:
411 raise BeachballError(
412 'invalid argument for projection: %s' % projection)
414 return points_out
417def deco_part(mt, mt_type='full', view='top'):
418 assert view in ('top', 'north', 'south', 'east', 'west'),\
419 'Allowed views are top, north, south, east and west'
420 mt = mtm.as_mt(mt)
422 if view == 'top':
423 pass
424 elif view == 'north':
425 mt = mt.rotated(_view_north)
426 elif view == 'south':
427 mt = mt.rotated(_view_south)
428 elif view == 'east':
429 mt = mt.rotated(_view_east)
430 elif view == 'west':
431 mt = mt.rotated(_view_west)
433 if mt_type == 'full':
434 return mt
436 res = mt.standard_decomposition()
437 m = dict(
438 dc=res[1][2],
439 deviatoric=res[3][2])[mt_type]
441 return mtm.MomentTensor(m=m)
444def choose_transform(axes, size_units, position, size):
446 if size_units == 'points':
447 transform = FixedPointOffsetTransform(
448 axes.transData,
449 axes.figure.dpi_scale_trans,
450 position)
452 if size is None:
453 size = 12.
455 size = size * 0.5 / 72.
456 position = (0., 0.)
458 elif size_units == 'data':
459 transform = axes.transData
461 if size is None:
462 size = 1.0
464 size = size * 0.5
466 elif size_units == 'axes':
467 transform = axes.transAxes
468 if size is None:
469 size = 1.
471 size = size * .5
473 else:
474 raise BeachballError(
475 'invalid argument for size_units: %s' % size_units)
477 position = num.asarray(position, dtype=num.float64)
479 return transform, position, size
482def mt2beachball(
483 mt,
484 beachball_type='deviatoric',
485 position=(0., 0.),
486 size=None,
487 color_t='red',
488 color_p='white',
489 edgecolor='black',
490 linewidth=2,
491 projection='lambert',
492 view='top'):
494 position = num.asarray(position, dtype=num.float64)
495 size = size or 1
496 mt = deco_part(mt, beachball_type, view)
498 eig = mt.eigensystem()
499 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0:
500 raise BeachballError('eigenvalues are zero')
502 data = []
503 for (group, patches, patches_lower, patches_upper,
504 lines, lines_lower, lines_upper) in eig2gx(eig):
506 if group == 'P':
507 color = color_p
508 else:
509 color = color_t
511 for poly in patches_upper:
512 verts = project(poly, projection)[:, ::-1] * size + \
513 position[NA, :]
514 data.append((verts, color, color, 1.0))
516 for poly in lines_upper:
517 verts = project(poly, projection)[:, ::-1] * size + \
518 position[NA, :]
519 data.append((verts, 'none', edgecolor, linewidth))
520 return data
523def plot_beachball_mpl(
524 mt, axes,
525 beachball_type='deviatoric',
526 position=(0., 0.),
527 size=None,
528 zorder=0,
529 color_t='red',
530 color_p='white',
531 edgecolor='black',
532 linewidth=2,
533 alpha=1.0,
534 arcres=181,
535 decimation=1,
536 projection='lambert',
537 size_units='points',
538 view='top'):
540 '''
541 Plot beachball diagram to a Matplotlib plot
543 :param mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
544 array or sequence which can be converted into an MT object
545 :param beachball_type: ``'deviatoric'`` (default), ``'full'``, or ``'dc'``
546 :param position: position of the beachball in data coordinates
547 :param size: diameter of the beachball either in points or in data
548 coordinates, depending on the ``size_units`` setting
549 :param zorder: (passed through to matplotlib drawing functions)
550 :param color_t: color for compressional quadrants (default: ``'red'``)
551 :param color_p: color for extensive quadrants (default: ``'white'``)
552 :param edgecolor: color for lines (default: ``'black'``)
553 :param linewidth: linewidth in points (default: ``2``)
554 :param alpha: (passed through to matplotlib drawing functions)
555 :param projection: ``'lambert'`` (default), ``'stereographic'``, or
556 ``'orthographic'``
557 :param size_units: ``'points'`` (default) or ``'data'``, where the
558 latter causes the beachball to be projected in the plots data
559 coordinates (axes must have an aspect ratio of 1.0 or the
560 beachball will be shown distorted when using this).
561 :param view: View the beachball from ``top``, ``north``, ``south``,
562 ``east`` or ``west``. Useful for to show beachballs in cross-sections.
563 Default is ``top``.
564 '''
566 transform, position, size = choose_transform(
567 axes, size_units, position, size)
569 mt = deco_part(mt, beachball_type, view)
571 eig = mt.eigensystem()
572 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0:
573 raise BeachballError('eigenvalues are zero')
575 data = []
576 for (group, patches, patches_lower, patches_upper,
577 lines, lines_lower, lines_upper) in eig2gx(eig, arcres):
579 if group == 'P':
580 color = color_p
581 else:
582 color = color_t
584 # plot "upper" features for lower hemisphere, because coordinate system
585 # is NED
587 for poly in patches_upper:
588 verts = project(poly, projection)[:, ::-1] * size + position[NA, :]
589 if alpha == 1.0:
590 data.append(
591 (verts[::decimation], color, color, linewidth))
592 else:
593 data.append(
594 (verts[::decimation], color, 'none', 0.0))
596 for poly in lines_upper:
597 verts = project(poly, projection)[:, ::-1] * size + position[NA, :]
598 data.append(
599 (verts[::decimation], 'none', edgecolor, linewidth))
601 patches = []
602 for (path, facecolor, edgecolor, linewidth) in data:
603 patches.append(Polygon(
604 xy=path, facecolor=facecolor,
605 edgecolor=edgecolor,
606 linewidth=linewidth,
607 alpha=alpha))
609 collection = PatchCollection(
610 patches, zorder=zorder, transform=transform, match_original=True)
612 axes.add_artist(collection)
613 return collection
616def mts2amps(mts, projection, beachball_type, grid_resolution=200, mask=True,
617 view='top'):
619 n_balls = len(mts)
620 nx = grid_resolution
621 ny = grid_resolution
623 x = num.linspace(-1., 1., nx)
624 y = num.linspace(-1., 1., ny)
626 vecs2 = num.zeros((nx * ny, 2), dtype=num.float64)
627 vecs2[:, 0] = num.tile(x, ny)
628 vecs2[:, 1] = num.repeat(y, nx)
630 ii_ok = vecs2[:, 0]**2 + vecs2[:, 1]**2 <= 1.0
631 amps = num_full(nx * ny, num.nan, dtype=num.float64)
633 amps[ii_ok] = 0.
634 for mt in mts:
635 mt = deco_part(mt, beachball_type, view)
637 ep, en, et, vp, vn, vt = mt.eigensystem()
639 vecs3_ok = inverse_project(vecs2[ii_ok, :], projection)
641 to_e = num.vstack((vn, vt, vp))
643 vecs_e = num.dot(to_e, vecs3_ok.T).T
644 rtp = numpy_xyz2rtp(vecs_e)
646 atheta, aphi = rtp[:, 1], rtp[:, 2]
647 amps_ok = ep * num.cos(atheta)**2 + (
648 en * num.cos(aphi)**2 + et * num.sin(aphi)**2) * num.sin(atheta)**2
650 if mask:
651 amps_ok[amps_ok > 0] = 1.
652 amps_ok[amps_ok < 0] = 0.
654 amps[ii_ok] += amps_ok
656 return num.reshape(amps, (ny, nx)) / n_balls, x, y
659def plot_fuzzy_beachball_mpl_pixmap(
660 mts, axes,
661 best_mt=None,
662 beachball_type='deviatoric',
663 position=(0., 0.),
664 size=None,
665 zorder=0,
666 color_t='red',
667 color_p='white',
668 edgecolor='black',
669 best_color='red',
670 linewidth=2,
671 alpha=1.0,
672 projection='lambert',
673 size_units='data',
674 grid_resolution=200,
675 method='imshow',
676 view='top'):
677 '''
678 Plot fuzzy beachball from a list of given MomentTensors
680 :param mts: list of
681 :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
682 array or sequence which can be converted into an MT object
683 :param best_mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or
684 an array or sequence which can be converted into an MT object
685 of most likely or minimum misfit solution to extra highlight
686 :param best_color: mpl color for best MomentTensor edges,
687 polygons are not plotted
689 See plot_beachball_mpl for other arguments
690 '''
691 if size_units == 'points':
692 raise BeachballError(
693 'size_units="points" not supported in '
694 'plot_fuzzy_beachball_mpl_pixmap')
696 transform, position, size = choose_transform(
697 axes, size_units, position, size)
699 amps, x, y = mts2amps(
700 mts,
701 grid_resolution=grid_resolution,
702 projection=projection,
703 beachball_type=beachball_type,
704 mask=True,
705 view=view)
707 ncolors = 256
708 cmap = LinearSegmentedColormap.from_list(
709 'dummy', [color_p, color_t], N=ncolors)
711 levels = num.linspace(0, 1., ncolors)
712 if method == 'contourf':
713 axes.contourf(
714 position[0] + y * size, position[1] + x * size, amps.T,
715 levels=levels,
716 cmap=cmap,
717 transform=transform,
718 zorder=zorder,
719 alpha=alpha)
721 elif method == 'imshow':
722 axes.imshow(
723 amps.T,
724 extent=(
725 position[0] + y[0] * size,
726 position[0] + y[-1] * size,
727 position[1] - x[0] * size,
728 position[1] - x[-1] * size),
729 cmap=cmap,
730 transform=transform,
731 zorder=zorder-0.1,
732 alpha=alpha)
733 else:
734 assert False, 'invalid `method` argument'
736 # draw optimum edges
737 if best_mt is not None:
738 best_amps, bx, by = mts2amps(
739 [best_mt],
740 grid_resolution=grid_resolution,
741 projection=projection,
742 beachball_type=beachball_type,
743 mask=False)
745 axes.contour(
746 position[0] + by * size, position[1] + bx * size, best_amps.T,
747 levels=[0.],
748 colors=[best_color],
749 linewidths=linewidth,
750 transform=transform,
751 zorder=zorder,
752 alpha=alpha)
754 phi = num.linspace(0., 2 * PI, 361)
755 x = num.cos(phi)
756 y = num.sin(phi)
757 axes.plot(
758 position[0] + x * size, position[1] + y * size,
759 linewidth=linewidth,
760 color=edgecolor,
761 transform=transform,
762 zorder=zorder,
763 alpha=alpha)
766def plot_beachball_mpl_construction(
767 mt, axes,
768 show='patches',
769 beachball_type='deviatoric',
770 view='top'):
772 mt = deco_part(mt, beachball_type, view)
773 eig = mt.eigensystem()
775 for (group, patches, patches_lower, patches_upper,
776 lines, lines_lower, lines_upper) in eig2gx(eig):
778 if group == 'P':
779 color = 'blue'
780 lw = 1
781 else:
782 color = 'red'
783 lw = 1
785 if show == 'patches':
786 for poly in patches_upper:
787 px, py, pz = poly.T
788 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
790 if show == 'lines':
791 for poly in lines_upper:
792 px, py, pz = poly.T
793 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
796def plot_beachball_mpl_pixmap(
797 mt, axes,
798 beachball_type='deviatoric',
799 position=(0., 0.),
800 size=None,
801 zorder=0,
802 color_t='red',
803 color_p='white',
804 edgecolor='black',
805 linewidth=2,
806 alpha=1.0,
807 projection='lambert',
808 size_units='data',
809 view='top'):
811 if size_units == 'points':
812 raise BeachballError(
813 'size_units="points" not supported in plot_beachball_mpl_pixmap')
815 transform, position, size = choose_transform(
816 axes, size_units, position, size)
818 mt = deco_part(mt, beachball_type, view)
820 ep, en, et, vp, vn, vt = mt.eigensystem()
822 amps, x, y = mts2amps(
823 [mt], projection, beachball_type, grid_resolution=200, mask=False)
825 axes.contourf(
826 position[0] + y * size, position[1] + x * size, amps.T,
827 levels=[-num.inf, 0., num.inf],
828 colors=[color_p, color_t],
829 transform=transform,
830 zorder=zorder,
831 alpha=alpha)
833 axes.contour(
834 position[0] + y * size, position[1] + x * size, amps.T,
835 levels=[0.],
836 colors=[edgecolor],
837 linewidths=linewidth,
838 transform=transform,
839 zorder=zorder,
840 alpha=alpha)
842 phi = num.linspace(0., 2 * PI, 361)
843 x = num.cos(phi)
844 y = num.sin(phi)
845 axes.plot(
846 position[0] + x * size, position[1] + y * size,
847 linewidth=linewidth,
848 color=edgecolor,
849 transform=transform,
850 zorder=zorder,
851 alpha=alpha)
854if __name__ == '__main__':
855 import sys
856 import os
857 import matplotlib.pyplot as plt
858 from pyrocko import model
860 args = sys.argv[1:]
862 data = []
863 for iarg, arg in enumerate(args):
865 if os.path.exists(arg):
866 events = model.load_events(arg)
867 for ev in events:
868 if not ev.moment_tensor:
869 logger.warning('no moment tensor given for event')
870 continue
872 data.append((ev.name, ev.moment_tensor))
873 else:
874 vals = list(map(float, arg.split(',')))
875 mt = mtm.as_mt(vals)
876 data.append(('%i' % (iarg+1), mt))
878 n = len(data)
880 ncols = 1
881 while ncols**2 < n:
882 ncols += 1
884 nrows = ncols
886 fig = plt.figure()
887 axes = fig.add_subplot(1, 1, 1, aspect=1.)
888 axes.axison = False
889 axes.set_xlim(-0.05 - ncols, ncols + 0.05)
890 axes.set_ylim(-0.05 - nrows, nrows + 0.05)
892 for ibeach, (name, mt) in enumerate(data):
893 irow = ibeach // ncols
894 icol = ibeach % ncols
895 plot_beachball_mpl(
896 mt, axes,
897 position=(icol*2-ncols+1, -irow*2+nrows-1),
898 size_units='data')
900 axes.annotate(
901 name,
902 xy=(icol*2-ncols+1, -irow*2+nrows-2),
903 xycoords='data',
904 xytext=(0, 0),
905 textcoords='offset points',
906 verticalalignment='center',
907 horizontalalignment='center',
908 rotation=0.)
910 plt.show()