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]])
31_view_east = num.array([[0, 0, -1],
32 [1, 0, 0],
33 [0, -1, 0]])
35_view_west = _view_east.T
38class BeachballError(Exception):
39 pass
42class _FixedPointOffsetTransform(Transform):
43 def __init__(self, trans, dpi_scale_trans, fixed_point):
44 Transform.__init__(self)
45 self.input_dims = self.output_dims = 2
46 self.has_inverse = False
47 self.trans = trans
48 self.dpi_scale_trans = dpi_scale_trans
49 self.fixed_point = num.asarray(fixed_point, dtype=num.float64)
51 def transform_non_affine(self, values):
52 fp = self.trans.transform(self.fixed_point)
53 return fp + self.dpi_scale_trans.transform(values)
56def vnorm(points):
57 return num.sqrt(num.sum(points**2, axis=1))
60def clean_poly(points):
61 if not num.all(points[0, :] == points[-1, :]):
62 points = num.vstack((points, points[0:1, :]))
64 dupl = num.concatenate(
65 (num.all(points[1:, :] == points[:-1, :], axis=1), [False]))
66 points = points[num.logical_not(dupl)]
67 return points
70def close_poly(points):
71 if not num.all(points[0, :] == points[-1, :]):
72 points = num.vstack((points, points[0:1, :]))
74 return points
77def circulation(points, axis):
78 # assert num.all(points[:, axis] >= 0.0) or num.all(points[:, axis] <= 0.0)
80 points2 = points[:, ((axis+2) % 3, (axis+1) % 3)].copy()
81 points2 *= 1.0 / num.sqrt(1.0 + num.abs(points[:, axis]))[:, num.newaxis]
83 result = -num.sum(
84 (points2[1:, 0] - points2[:-1, 0]) *
85 (points2[1:, 1] + points2[:-1, 1]))
87 result -= (points2[0, 0] - points2[-1, 0]) \
88 * (points2[0, 1] + points2[-1, 1])
89 return result
92def spoly_cut(l_points, axis=0, nonsimple=True, arcres=181):
93 dphi = 2.*PI / (2*arcres)
95 # cut sub-polygons and gather crossing point information
96 crossings = []
97 snippets = {}
98 for ipath, points in enumerate(l_points):
99 if not num.all(points[0, :] == points[-1, :]):
100 points = num.vstack((points, points[0:1, :]))
102 # get upward crossing points
103 iup = num.where(num.logical_and(points[:-1, axis] <= 0.,
104 points[1:, axis] > 0.))[0]
105 aup = - points[iup, axis] / (points[iup+1, axis] - points[iup, axis])
106 pup = points[iup, :] + aup[:, num.newaxis] * (points[iup+1, :] -
107 points[iup, :])
108 phiup = num.arctan2(pup[:, (axis+2) % 3], pup[:, (axis+1) % 3])
110 for i in range(len(iup)):
111 crossings.append((phiup[i], ipath, iup[i], 1, pup[i], [1, -1]))
113 # get downward crossing points
114 idown = num.where(num.logical_and(points[:-1, axis] > 0.,
115 points[1:, axis] <= 0.))[0]
116 adown = - points[idown+1, axis] / (points[idown, axis] -
117 points[idown+1, axis])
118 pdown = points[idown+1, :] + adown[:, num.newaxis] * (
119 points[idown, :] - points[idown+1, :])
120 phidown = num.arctan2(pdown[:, (axis+2) % 3], pdown[:, (axis+1) % 3])
122 for i in range(idown.size):
123 crossings.append(
124 (phidown[i], ipath, idown[i], -1, pdown[i], [1, -1]))
126 icuts = num.sort(num.concatenate((iup, idown)))
128 for i in range(icuts.size-1):
129 snippets[ipath, icuts[i]] = (
130 ipath, icuts[i+1], points[icuts[i]+1:icuts[i+1]+1])
132 if icuts.size:
133 points_last = num.concatenate((
134 points[icuts[-1]+1:],
135 points[:icuts[0]+1]))
137 snippets[ipath, icuts[-1]] = (ipath, icuts[0], points_last)
138 else:
139 snippets[ipath, 0] = (ipath, 0, points)
141 crossings.sort()
143 # assemble new sub-polygons
144 current = snippets.pop(list(snippets.keys())[0])
145 outs = [[]]
146 while True:
147 outs[-1].append(current[2])
148 for i, c1 in enumerate(crossings):
149 if c1[1:3] == current[:2]:
150 direction = -1 * c1[3]
151 break
152 else:
153 if not snippets:
154 break
155 current = snippets.pop(list(snippets.keys())[0])
156 outs.append([])
157 continue
159 while True:
160 i = (i + direction) % len(crossings)
161 if crossings[i][3] == direction and direction in crossings[i][-1]:
162 break
164 c2 = crossings[i]
165 c2[-1].remove(direction)
167 phi1 = c1[0]
168 phi2 = c2[0]
169 if direction == 1:
170 if phi1 > phi2:
171 phi2 += PI * 2.
173 if direction == -1:
174 if phi1 < phi2:
175 phi2 -= PI * 2.
177 n = int(abs(phi2 - phi1) / dphi) + 2
179 phis = num.linspace(phi1, phi2, n)
180 cpoints = num.zeros((n, 3))
181 cpoints[:, (axis+1) % 3] = num.cos(phis)
182 cpoints[:, (axis+2) % 3] = num.sin(phis)
183 cpoints[:, axis] = 0.0
185 outs[-1].append(cpoints)
187 try:
188 current = snippets[c2[1:3]]
189 del snippets[c2[1:3]]
191 except KeyError:
192 if not snippets:
193 break
195 current = snippets.pop(list(snippets.keys())[0])
196 outs.append([])
198 # separate hemispheres, force polygons closed, remove duplicate points
199 # remove polygons with less than 3 points (4, when counting repeated
200 # endpoint)
202 outs_upper = []
203 outs_lower = []
204 for out in outs:
205 if out:
206 out = clean_poly(num.vstack(out))
207 if out.shape[0] >= 4:
208 if num.sum(out[:, axis]) > 0.0:
209 outs_upper.append(out)
210 else:
211 outs_lower.append(out)
213 if nonsimple and (
214 len(crossings) == 0 or
215 len(outs_upper) == 0 or
216 len(outs_lower) == 0):
218 # check if we are cutting between holes
219 need_divider = False
220 if outs_upper:
221 candis = sorted(
222 outs_upper, key=lambda out: num.min(out[:, axis]))
224 if circulation(candis[0], axis) > 0.0:
225 need_divider = True
227 if outs_lower:
228 candis = sorted(
229 outs_lower, key=lambda out: num.max(out[:, axis]))
231 if circulation(candis[0], axis) < 0.0:
232 need_divider = True
234 if need_divider:
235 phi1 = 0.
236 phi2 = PI*2.
237 n = int(abs(phi2 - phi1) / dphi) + 2
239 phis = num.linspace(phi1, phi2, n)
240 cpoints = num.zeros((n, 3))
241 cpoints[:, (axis+1) % 3] = num.cos(phis)
242 cpoints[:, (axis+2) % 3] = num.sin(phis)
243 cpoints[:, axis] = 0.0
245 outs_upper.append(cpoints)
246 outs_lower.append(cpoints[::-1, :])
248 return outs_lower, outs_upper
251def numpy_rtp2xyz(rtp):
252 r = rtp[:, 0]
253 theta = rtp[:, 1]
254 phi = rtp[:, 2]
255 vecs = num.empty(rtp.shape, dtype=num.float64)
256 vecs[:, 0] = r*num.sin(theta)*num.cos(phi)
257 vecs[:, 1] = r*num.sin(theta)*num.sin(phi)
258 vecs[:, 2] = r*num.cos(theta)
259 return vecs
262def numpy_xyz2rtp(xyz):
263 x, y, z = xyz[:, 0], xyz[:, 1], xyz[:, 2]
264 vecs = num.empty(xyz.shape, dtype=num.float64)
265 vecs[:, 0] = num.sqrt(x**2+y**2+z**2)
266 vecs[:, 1] = num.arctan2(num.sqrt(x**2+y**2), z)
267 vecs[:, 2] = num.arctan2(y, x)
268 return vecs
271def circle_points(aphi, sign=1.0):
272 vecs = num.empty((aphi.size, 3), dtype=num.float64)
273 vecs[:, 0] = num.cos(sign*aphi)
274 vecs[:, 1] = num.sin(sign*aphi)
275 vecs[:, 2] = 0.0
276 return vecs
279def eig2gx(eig, arcres=181):
280 aphi = num.linspace(0., 2.*PI, arcres)
281 ep, en, et, vp, vn, vt = eig
283 mt_sign = num.sign(ep + en + et)
285 groups = []
286 for (pt_name, pt_sign) in [('P', -1.), ('T', 1.)]:
287 patches = []
288 patches_lower = []
289 patches_upper = []
290 lines = []
291 lines_lower = []
292 lines_upper = []
293 for iperm, (va, vb, vc, ea, eb, ec) in enumerate([
294 (vp, vn, vt, ep, en, et),
295 (vt, vp, vn, et, ep, en)]): # (vn, vt, vp, en, et, ep)]):
297 perm_sign = [-1.0, 1.0][iperm]
298 to_e = num.vstack((vb, vc, va))
299 from_e = to_e.T
301 poly_es = []
302 polys = []
303 for sign in (-1., 1.):
304 xphi = perm_sign*pt_sign*sign*aphi
305 denom = eb*num.cos(xphi)**2 + ec*num.sin(xphi)**2
306 if num.any(denom == 0.):
307 continue
309 Y = -ea/denom
310 if num.any(Y < 0.):
311 continue
313 xtheta = num.arctan(num.sqrt(Y))
314 rtp = num.empty(xphi.shape+(3,), dtype=num.float64)
315 rtp[:, 0] = 1.
316 if sign > 0:
317 rtp[:, 1] = xtheta
318 else:
319 rtp[:, 1] = PI - xtheta
321 rtp[:, 2] = xphi
322 poly_e = numpy_rtp2xyz(rtp)
323 poly = num.dot(from_e, poly_e.T).T
324 poly[:, 2] -= 0.001
326 poly_es.append(poly_e)
327 polys.append(poly)
329 if polys:
330 polys_lower, polys_upper = spoly_cut(polys, 2, arcres=arcres)
331 lines.extend(polys)
332 lines_lower.extend(polys_lower)
333 lines_upper.extend(polys_upper)
335 if poly_es:
336 for aa in spoly_cut(poly_es, 0, arcres=arcres):
337 for bb in spoly_cut(aa, 1, arcres=arcres):
338 for cc in spoly_cut(bb, 2, arcres=arcres):
339 for poly_e in cc:
340 poly = num.dot(from_e, poly_e.T).T
341 poly[:, 2] -= 0.001
342 polys_lower, polys_upper = spoly_cut(
343 [poly], 2, nonsimple=False, arcres=arcres)
345 patches.append(poly)
346 patches_lower.extend(polys_lower)
347 patches_upper.extend(polys_upper)
349 if not patches:
350 if mt_sign * pt_sign == 1.:
351 patches_lower.append(circle_points(aphi, -1.0))
352 patches_upper.append(circle_points(aphi, 1.0))
353 lines_lower.append(circle_points(aphi, -1.0))
354 lines_upper.append(circle_points(aphi, 1.0))
356 groups.append((
357 pt_name,
358 patches, patches_lower, patches_upper,
359 lines, lines_lower, lines_upper))
361 return groups
364def extr(points):
365 pmean = num.mean(points, axis=0)
366 return points + pmean*0.05
369def draw_eigenvectors_mpl(eig, axes):
370 vp, vn, vt = eig[3:]
371 for lab, v in [('P', vp), ('N', vn), ('T', vt)]:
372 sign = num.sign(v[2]) + (v[2] == 0.0)
373 axes.plot(sign*v[1], sign*v[0], 'o', color='black')
374 axes.text(sign*v[1], sign*v[0], ' '+lab)
377def project(points, projection='lambert'):
378 points_out = points[:, :2].copy()
379 if projection == 'lambert':
380 factor = 1.0 / num.sqrt(1.0 + points[:, 2])
381 elif projection == 'stereographic':
382 factor = 1.0 / (1.0 + points[:, 2])
383 elif projection == 'orthographic':
384 factor = None
385 else:
386 raise BeachballError(
387 'invalid argument for projection: %s' % projection)
389 if factor is not None:
390 points_out *= factor[:, num.newaxis]
392 return points_out
395def inverse_project(points, projection='lambert'):
396 points_out = num.zeros((points.shape[0], 3))
398 rsqr = points[:, 0]**2 + points[:, 1]**2
399 if projection == 'lambert':
400 points_out[:, 2] = 1.0 - rsqr
401 points_out[:, 1] = num.sqrt(2.0 - rsqr) * points[:, 1]
402 points_out[:, 0] = num.sqrt(2.0 - rsqr) * points[:, 0]
403 elif projection == 'stereographic':
404 points_out[:, 2] = - (rsqr - 1.0) / (rsqr + 1.0)
405 points_out[:, 1] = 2.0 * points[:, 1] / (rsqr + 1.0)
406 points_out[:, 0] = 2.0 * points[:, 0] / (rsqr + 1.0)
407 elif projection == 'orthographic':
408 points_out[:, 2] = num.sqrt(num.maximum(1.0 - rsqr, 0.0))
409 points_out[:, 1] = points[:, 1]
410 points_out[:, 0] = points[:, 0]
411 else:
412 raise BeachballError(
413 'invalid argument for projection: %s' % projection)
415 return points_out
418def deco_part(mt, mt_type='full', view='top'):
419 assert view in ('top', 'north', 'south', 'east', 'west'),\
420 'Allowed views are top, north, south, east and west'
421 mt = mtm.as_mt(mt)
423 if view == 'top':
424 pass
425 elif view == 'north':
426 mt = mt.rotated(_view_north)
427 elif view == 'south':
428 mt = mt.rotated(_view_south)
429 elif view == 'east':
430 mt = mt.rotated(_view_east)
431 elif view == 'west':
432 mt = mt.rotated(_view_west)
434 if mt_type == 'full':
435 return mt
437 res = mt.standard_decomposition()
438 m = dict(
439 dc=res[1][2],
440 deviatoric=res[3][2])[mt_type]
442 return mtm.MomentTensor(m=m)
445def choose_transform(axes, size_units, position, size):
447 if size_units == 'points':
448 transform = _FixedPointOffsetTransform(
449 axes.transData,
450 axes.figure.dpi_scale_trans,
451 position)
453 if size is None:
454 size = 12.
456 size = size * 0.5 / 72.
457 position = (0., 0.)
459 elif size_units == 'data':
460 transform = axes.transData
462 if size is None:
463 size = 1.0
465 size = size * 0.5
467 elif size_units == 'axes':
468 transform = axes.transAxes
469 if size is None:
470 size = 1.
472 size = size * .5
474 else:
475 raise BeachballError(
476 'invalid argument for size_units: %s' % size_units)
478 position = num.asarray(position, dtype=num.float64)
480 return transform, position, size
483def mt2beachball(
484 mt,
485 beachball_type='deviatoric',
486 position=(0., 0.),
487 size=None,
488 color_t='red',
489 color_p='white',
490 edgecolor='black',
491 linewidth=2,
492 projection='lambert',
493 view='top'):
495 position = num.asarray(position, dtype=num.float64)
496 size = size or 1
497 mt = deco_part(mt, beachball_type, view)
499 eig = mt.eigensystem()
500 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0:
501 raise BeachballError('eigenvalues are zero')
503 data = []
504 for (group, patches, patches_lower, patches_upper,
505 lines, lines_lower, lines_upper) in eig2gx(eig):
507 if group == 'P':
508 color = color_p
509 else:
510 color = color_t
512 for poly in patches_upper:
513 verts = project(poly, projection)[:, ::-1] * size + \
514 position[NA, :]
515 data.append((verts, color, color, 1.0))
517 for poly in lines_upper:
518 verts = project(poly, projection)[:, ::-1] * size + \
519 position[NA, :]
520 data.append((verts, 'none', edgecolor, linewidth))
521 return data
524def plot_beachball_mpl(
525 mt, axes,
526 beachball_type='deviatoric',
527 position=(0., 0.),
528 size=None,
529 zorder=0,
530 color_t='red',
531 color_p='white',
532 edgecolor='black',
533 linewidth=2,
534 alpha=1.0,
535 arcres=181,
536 decimation=1,
537 projection='lambert',
538 size_units='points',
539 view='top'):
541 '''
542 Plot beachball diagram to a Matplotlib plot
544 :param mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
545 array or sequence which can be converted into an MT object
546 :param beachball_type: ``'deviatoric'`` (default), ``'full'``, or ``'dc'``
547 :param position: position of the beachball in data coordinates
548 :param size: diameter of the beachball either in points or in data
549 coordinates, depending on the ``size_units`` setting
550 :param zorder: (passed through to matplotlib drawing functions)
551 :param color_t: color for compressional quadrants (default: ``'red'``)
552 :param color_p: color for extensive quadrants (default: ``'white'``)
553 :param edgecolor: color for lines (default: ``'black'``)
554 :param linewidth: linewidth in points (default: ``2``)
555 :param alpha: (passed through to matplotlib drawing functions)
556 :param projection: ``'lambert'`` (default), ``'stereographic'``, or
557 ``'orthographic'``
558 :param size_units: ``'points'`` (default) or ``'data'``, where the
559 latter causes the beachball to be projected in the plots data
560 coordinates (axes must have an aspect ratio of 1.0 or the
561 beachball will be shown distorted when using this).
562 :param view: View the beachball from ``top``, ``north``, ``south``,
563 ``east`` or ``west``. Useful for to show beachballs in cross-sections.
564 Default is ``top``.
565 '''
567 transform, position, size = choose_transform(
568 axes, size_units, position, size)
570 mt = deco_part(mt, beachball_type, view)
572 eig = mt.eigensystem()
573 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0:
574 raise BeachballError('eigenvalues are zero')
576 data = []
577 for (group, patches, patches_lower, patches_upper,
578 lines, lines_lower, lines_upper) in eig2gx(eig, arcres):
580 if group == 'P':
581 color = color_p
582 else:
583 color = color_t
585 # plot "upper" features for lower hemisphere, because coordinate system
586 # is NED
588 for poly in patches_upper:
589 verts = project(poly, projection)[:, ::-1] * size + position[NA, :]
590 if alpha == 1.0:
591 data.append(
592 (verts[::decimation], color, color, linewidth))
593 else:
594 data.append(
595 (verts[::decimation], color, 'none', 0.0))
597 for poly in lines_upper:
598 verts = project(poly, projection)[:, ::-1] * size + position[NA, :]
599 data.append(
600 (verts[::decimation], 'none', edgecolor, linewidth))
602 patches = []
603 for (path, facecolor, edgecolor, linewidth) in data:
604 patches.append(Polygon(
605 xy=path, facecolor=facecolor,
606 edgecolor=edgecolor,
607 linewidth=linewidth,
608 alpha=alpha))
610 collection = PatchCollection(
611 patches, zorder=zorder, transform=transform, match_original=True)
613 axes.add_artist(collection)
614 return collection
617def mts2amps(mts, projection, beachball_type, grid_resolution=200, mask=True,
618 view='top'):
620 n_balls = len(mts)
621 nx = grid_resolution
622 ny = grid_resolution
624 x = num.linspace(-1., 1., nx)
625 y = num.linspace(-1., 1., ny)
627 vecs2 = num.zeros((nx * ny, 2), dtype=num.float64)
628 vecs2[:, 0] = num.tile(x, ny)
629 vecs2[:, 1] = num.repeat(y, nx)
631 ii_ok = vecs2[:, 0]**2 + vecs2[:, 1]**2 <= 1.0
632 amps = num_full(nx * ny, num.nan, dtype=num.float64)
634 amps[ii_ok] = 0.
635 for mt in mts:
636 mt = deco_part(mt, beachball_type, view)
638 ep, en, et, vp, vn, vt = mt.eigensystem()
640 vecs3_ok = inverse_project(vecs2[ii_ok, :], projection)
642 to_e = num.vstack((vn, vt, vp))
644 vecs_e = num.dot(to_e, vecs3_ok.T).T
645 rtp = numpy_xyz2rtp(vecs_e)
647 atheta, aphi = rtp[:, 1], rtp[:, 2]
648 amps_ok = ep * num.cos(atheta)**2 + (
649 en * num.cos(aphi)**2 + et * num.sin(aphi)**2) * num.sin(atheta)**2
651 if mask:
652 amps_ok[amps_ok > 0] = 1.
653 amps_ok[amps_ok < 0] = 0.
655 amps[ii_ok] += amps_ok
657 return num.reshape(amps, (ny, nx)) / n_balls, x, y
660def plot_fuzzy_beachball_mpl_pixmap(
661 mts, axes,
662 best_mt=None,
663 beachball_type='deviatoric',
664 position=(0., 0.),
665 size=None,
666 zorder=0,
667 color_t='red',
668 color_p='white',
669 edgecolor='black',
670 best_color='red',
671 linewidth=2,
672 alpha=1.0,
673 projection='lambert',
674 size_units='data',
675 grid_resolution=200,
676 method='imshow',
677 view='top'):
678 '''
679 Plot fuzzy beachball from a list of given MomentTensors
681 :param mts: list of
682 :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
683 array or sequence which can be converted into an MT object
684 :param best_mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or
685 an array or sequence which can be converted into an MT object
686 of most likely or minimum misfit solution to extra highlight
687 :param best_color: mpl color for best MomentTensor edges,
688 polygons are not plotted
690 See plot_beachball_mpl for other arguments
691 '''
692 if size_units == 'points':
693 raise BeachballError(
694 'size_units="points" not supported in '
695 'plot_fuzzy_beachball_mpl_pixmap')
697 transform, position, size = choose_transform(
698 axes, size_units, position, size)
700 amps, x, y = mts2amps(
701 mts,
702 grid_resolution=grid_resolution,
703 projection=projection,
704 beachball_type=beachball_type,
705 mask=True,
706 view=view)
708 ncolors = 256
709 cmap = LinearSegmentedColormap.from_list(
710 'dummy', [color_p, color_t], N=ncolors)
712 levels = num.linspace(0, 1., ncolors)
713 if method == 'contourf':
714 axes.contourf(
715 position[0] + y * size, position[1] + x * size, amps.T,
716 levels=levels,
717 cmap=cmap,
718 transform=transform,
719 zorder=zorder,
720 alpha=alpha)
722 elif method == 'imshow':
723 axes.imshow(
724 amps.T,
725 extent=(
726 position[0] + y[0] * size,
727 position[0] + y[-1] * size,
728 position[1] - x[0] * size,
729 position[1] - x[-1] * size),
730 cmap=cmap,
731 transform=transform,
732 zorder=zorder-0.1,
733 alpha=alpha)
734 else:
735 assert False, 'invalid `method` argument'
737 # draw optimum edges
738 if best_mt is not None:
739 best_amps, bx, by = mts2amps(
740 [best_mt],
741 grid_resolution=grid_resolution,
742 projection=projection,
743 beachball_type=beachball_type,
744 mask=False)
746 axes.contour(
747 position[0] + by * size, position[1] + bx * size, best_amps.T,
748 levels=[0.],
749 colors=[best_color],
750 linewidths=linewidth,
751 transform=transform,
752 zorder=zorder,
753 alpha=alpha)
755 phi = num.linspace(0., 2 * PI, 361)
756 x = num.cos(phi)
757 y = num.sin(phi)
758 axes.plot(
759 position[0] + x * size, position[1] + y * size,
760 linewidth=linewidth,
761 color=edgecolor,
762 transform=transform,
763 zorder=zorder,
764 alpha=alpha)
767def plot_beachball_mpl_construction(
768 mt, axes,
769 show='patches',
770 beachball_type='deviatoric',
771 view='top'):
773 mt = deco_part(mt, beachball_type, view)
774 eig = mt.eigensystem()
776 for (group, patches, patches_lower, patches_upper,
777 lines, lines_lower, lines_upper) in eig2gx(eig):
779 if group == 'P':
780 color = 'blue'
781 lw = 1
782 else:
783 color = 'red'
784 lw = 1
786 if show == 'patches':
787 for poly in patches_upper:
788 px, py, pz = poly.T
789 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
791 if show == 'lines':
792 for poly in lines_upper:
793 px, py, pz = poly.T
794 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
797def plot_beachball_mpl_pixmap(
798 mt, axes,
799 beachball_type='deviatoric',
800 position=(0., 0.),
801 size=None,
802 zorder=0,
803 color_t='red',
804 color_p='white',
805 edgecolor='black',
806 linewidth=2,
807 alpha=1.0,
808 projection='lambert',
809 size_units='data',
810 view='top'):
812 if size_units == 'points':
813 raise BeachballError(
814 'size_units="points" not supported in plot_beachball_mpl_pixmap')
816 transform, position, size = choose_transform(
817 axes, size_units, position, size)
819 mt = deco_part(mt, beachball_type, view)
821 ep, en, et, vp, vn, vt = mt.eigensystem()
823 amps, x, y = mts2amps(
824 [mt], projection, beachball_type, grid_resolution=200, mask=False)
826 axes.contourf(
827 position[0] + y * size, position[1] + x * size, amps.T,
828 levels=[-num.inf, 0., num.inf],
829 colors=[color_p, color_t],
830 transform=transform,
831 zorder=zorder,
832 alpha=alpha)
834 axes.contour(
835 position[0] + y * size, position[1] + x * size, amps.T,
836 levels=[0.],
837 colors=[edgecolor],
838 linewidths=linewidth,
839 transform=transform,
840 zorder=zorder,
841 alpha=alpha)
843 phi = num.linspace(0., 2 * PI, 361)
844 x = num.cos(phi)
845 y = num.sin(phi)
846 axes.plot(
847 position[0] + x * size, position[1] + y * size,
848 linewidth=linewidth,
849 color=edgecolor,
850 transform=transform,
851 zorder=zorder,
852 alpha=alpha)
855if __name__ == '__main__':
856 import sys
857 import os
858 import matplotlib.pyplot as plt
859 from pyrocko import model
861 args = sys.argv[1:]
863 data = []
864 for iarg, arg in enumerate(args):
866 if os.path.exists(arg):
867 events = model.load_events(arg)
868 for ev in events:
869 if not ev.moment_tensor:
870 logger.warning('no moment tensor given for event')
871 continue
873 data.append((ev.name, ev.moment_tensor))
874 else:
875 vals = list(map(float, arg.split(',')))
876 mt = mtm.as_mt(vals)
877 data.append(('%i' % (iarg+1), mt))
879 n = len(data)
881 ncols = 1
882 while ncols**2 < n:
883 ncols += 1
885 nrows = ncols
887 fig = plt.figure()
888 axes = fig.add_subplot(1, 1, 1, aspect=1.)
889 axes.axison = False
890 axes.set_xlim(-0.05 - ncols, ncols + 0.05)
891 axes.set_ylim(-0.05 - nrows, nrows + 0.05)
893 for ibeach, (name, mt) in enumerate(data):
894 irow = ibeach // ncols
895 icol = ibeach % ncols
896 plot_beachball_mpl(
897 mt, axes,
898 position=(icol*2-ncols+1, -irow*2+nrows-1),
899 size_units='data')
901 axes.annotate(
902 name,
903 xy=(icol*2-ncols+1, -irow*2+nrows-2),
904 xycoords='data',
905 xytext=(0, 0),
906 textcoords='offset points',
907 verticalalignment='center',
908 horizontalalignment='center',
909 rotation=0.)
911 plt.show()