1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6from math import pi as PI
7import logging
8import numpy as num
10from matplotlib.collections import PatchCollection
11from matplotlib.patches import Polygon
12from matplotlib.transforms import Transform
13from matplotlib.colors import LinearSegmentedColormap
15from pyrocko import moment_tensor as mtm
16from pyrocko.util import num_full
18logger = logging.getLogger('pyrocko.plot.beachball')
20NA = num.newaxis
21d2r = num.pi / 180.
24def view_rotation(strike, dip):
25 return mtm.euler_to_matrix(
26 dip*d2r, strike*d2r, -90.*d2r)
29_view_south = view_rotation(90., 90.)
30_view_north = view_rotation(-90., 90.)
31_view_east = view_rotation(0., 90.)
32_view_west = view_rotation(180., 90.)
35class BeachballError(Exception):
36 pass
39class _FixedPointOffsetTransform(Transform):
40 def __init__(self, trans, dpi_scale_trans, fixed_point):
41 Transform.__init__(self)
42 self.input_dims = self.output_dims = 2
43 self.has_inverse = False
44 self.trans = trans
45 self.dpi_scale_trans = dpi_scale_trans
46 self.fixed_point = num.asarray(fixed_point, dtype=num.float64)
48 def transform_non_affine(self, values):
49 fp = self.trans.transform(self.fixed_point)
50 return fp + self.dpi_scale_trans.transform(values)
53def vnorm(points):
54 return num.sqrt(num.sum(points**2, axis=1))
57def clean_poly(points):
58 if not num.all(points[0, :] == points[-1, :]):
59 points = num.vstack((points, points[0:1, :]))
61 dupl = num.concatenate(
62 (num.all(points[1:, :] == points[:-1, :], axis=1), [False]))
63 points = points[num.logical_not(dupl)]
64 return points
67def close_poly(points):
68 if not num.all(points[0, :] == points[-1, :]):
69 points = num.vstack((points, points[0:1, :]))
71 return points
74def circulation(points, axis):
75 # assert num.all(points[:, axis] >= 0.0) or num.all(points[:, axis] <= 0.0)
77 points2 = points[:, ((axis+2) % 3, (axis+1) % 3)].copy()
78 points2 *= 1.0 / num.sqrt(1.0 + num.abs(points[:, axis]))[:, num.newaxis]
80 result = -num.sum(
81 (points2[1:, 0] - points2[:-1, 0]) *
82 (points2[1:, 1] + points2[:-1, 1]))
84 result -= (points2[0, 0] - points2[-1, 0]) \
85 * (points2[0, 1] + points2[-1, 1])
86 return result
89def spoly_cut(l_points, axis=0, nonsimple=True, arcres=181):
90 dphi = 2.*PI / (2*arcres)
92 # cut sub-polygons and gather crossing point information
93 crossings = []
94 snippets = {}
95 for ipath, points in enumerate(l_points):
96 if not num.all(points[0, :] == points[-1, :]):
97 points = num.vstack((points, points[0:1, :]))
99 # get upward crossing points
100 iup = num.where(num.logical_and(points[:-1, axis] <= 0.,
101 points[1:, axis] > 0.))[0]
102 aup = - points[iup, axis] / (points[iup+1, axis] - points[iup, axis])
103 pup = points[iup, :] + aup[:, num.newaxis] * (points[iup+1, :] -
104 points[iup, :])
105 phiup = num.arctan2(pup[:, (axis+2) % 3], pup[:, (axis+1) % 3])
107 for i in range(len(iup)):
108 crossings.append((phiup[i], ipath, iup[i], 1, pup[i], [1, -1]))
110 # get downward crossing points
111 idown = num.where(num.logical_and(points[:-1, axis] > 0.,
112 points[1:, axis] <= 0.))[0]
113 adown = - points[idown+1, axis] / (points[idown, axis] -
114 points[idown+1, axis])
115 pdown = points[idown+1, :] + adown[:, num.newaxis] * (
116 points[idown, :] - points[idown+1, :])
117 phidown = num.arctan2(pdown[:, (axis+2) % 3], pdown[:, (axis+1) % 3])
119 for i in range(idown.size):
120 crossings.append(
121 (phidown[i], ipath, idown[i], -1, pdown[i], [1, -1]))
123 icuts = num.sort(num.concatenate((iup, idown)))
125 for i in range(icuts.size-1):
126 snippets[ipath, icuts[i]] = (
127 ipath, icuts[i+1], points[icuts[i]+1:icuts[i+1]+1])
129 if icuts.size:
130 points_last = num.concatenate((
131 points[icuts[-1]+1:],
132 points[:icuts[0]+1]))
134 snippets[ipath, icuts[-1]] = (ipath, icuts[0], points_last)
135 else:
136 snippets[ipath, 0] = (ipath, 0, points)
138 crossings.sort()
140 # assemble new sub-polygons
141 current = snippets.pop(list(snippets.keys())[0])
142 outs = [[]]
143 while True:
144 outs[-1].append(current[2])
145 for i, c1 in enumerate(crossings):
146 if c1[1:3] == current[:2]:
147 direction = -1 * c1[3]
148 break
149 else:
150 if not snippets:
151 break
152 current = snippets.pop(list(snippets.keys())[0])
153 outs.append([])
154 continue
156 while True:
157 i = (i + direction) % len(crossings)
158 if crossings[i][3] == direction and direction in crossings[i][-1]:
159 break
161 c2 = crossings[i]
162 c2[-1].remove(direction)
164 phi1 = c1[0]
165 phi2 = c2[0]
166 if direction == 1:
167 if phi1 > phi2:
168 phi2 += PI * 2.
170 if direction == -1:
171 if phi1 < phi2:
172 phi2 -= PI * 2.
174 n = int(abs(phi2 - phi1) / dphi) + 2
176 phis = num.linspace(phi1, phi2, n)
177 cpoints = num.zeros((n, 3))
178 cpoints[:, (axis+1) % 3] = num.cos(phis)
179 cpoints[:, (axis+2) % 3] = num.sin(phis)
180 cpoints[:, axis] = 0.0
182 outs[-1].append(cpoints)
184 try:
185 current = snippets[c2[1:3]]
186 del snippets[c2[1:3]]
188 except KeyError:
189 if not snippets:
190 break
192 current = snippets.pop(list(snippets.keys())[0])
193 outs.append([])
195 # separate hemispheres, force polygons closed, remove duplicate points
196 # remove polygons with less than 3 points (4, when counting repeated
197 # endpoint)
199 outs_upper = []
200 outs_lower = []
201 for out in outs:
202 if out:
203 out = clean_poly(num.vstack(out))
204 if out.shape[0] >= 4:
205 if num.sum(out[:, axis]) > 0.0:
206 outs_upper.append(out)
207 else:
208 outs_lower.append(out)
210 if nonsimple and (
211 len(crossings) == 0 or
212 len(outs_upper) == 0 or
213 len(outs_lower) == 0):
215 # check if we are cutting between holes
216 need_divider = False
217 if outs_upper:
218 candis = sorted(
219 outs_upper, key=lambda out: num.min(out[:, axis]))
221 if circulation(candis[0], axis) > 0.0:
222 need_divider = True
224 if outs_lower:
225 candis = sorted(
226 outs_lower, key=lambda out: num.max(out[:, axis]))
228 if circulation(candis[0], axis) < 0.0:
229 need_divider = True
231 if need_divider:
232 phi1 = 0.
233 phi2 = PI*2.
234 n = int(abs(phi2 - phi1) / dphi) + 2
236 phis = num.linspace(phi1, phi2, n)
237 cpoints = num.zeros((n, 3))
238 cpoints[:, (axis+1) % 3] = num.cos(phis)
239 cpoints[:, (axis+2) % 3] = num.sin(phis)
240 cpoints[:, axis] = 0.0
242 outs_upper.append(cpoints)
243 outs_lower.append(cpoints[::-1, :])
245 return outs_lower, outs_upper
248def numpy_rtp2xyz(rtp):
249 r = rtp[:, 0]
250 theta = rtp[:, 1]
251 phi = rtp[:, 2]
252 vecs = num.empty(rtp.shape, dtype=num.float64)
253 vecs[:, 0] = r*num.sin(theta)*num.cos(phi)
254 vecs[:, 1] = r*num.sin(theta)*num.sin(phi)
255 vecs[:, 2] = r*num.cos(theta)
256 return vecs
259def numpy_xyz2rtp(xyz):
260 x, y, z = xyz[:, 0], xyz[:, 1], xyz[:, 2]
261 vecs = num.empty(xyz.shape, dtype=num.float64)
262 vecs[:, 0] = num.sqrt(x**2+y**2+z**2)
263 vecs[:, 1] = num.arctan2(num.sqrt(x**2+y**2), z)
264 vecs[:, 2] = num.arctan2(y, x)
265 return vecs
268def circle_points(aphi, sign=1.0):
269 vecs = num.empty((aphi.size, 3), dtype=num.float64)
270 vecs[:, 0] = num.cos(sign*aphi)
271 vecs[:, 1] = num.sin(sign*aphi)
272 vecs[:, 2] = 0.0
273 return vecs
276def eig2gx(eig, arcres=181):
277 aphi = num.linspace(0., 2.*PI, arcres)
278 ep, en, et, vp, vn, vt = eig
280 mt_sign = num.sign(ep + en + et)
282 groups = []
283 for (pt_name, pt_sign) in [('P', -1.), ('T', 1.)]:
284 patches = []
285 patches_lower = []
286 patches_upper = []
287 lines = []
288 lines_lower = []
289 lines_upper = []
290 for iperm, (va, vb, vc, ea, eb, ec) in enumerate([
291 (vp, vn, vt, ep, en, et),
292 (vt, vp, vn, et, ep, en)]): # (vn, vt, vp, en, et, ep)]):
294 perm_sign = [-1.0, 1.0][iperm]
295 to_e = num.vstack((vb, vc, va))
296 from_e = to_e.T
298 poly_es = []
299 polys = []
300 for sign in (-1., 1.):
301 xphi = perm_sign*pt_sign*sign*aphi
302 denom = eb*num.cos(xphi)**2 + ec*num.sin(xphi)**2
303 if num.any(denom == 0.):
304 continue
306 Y = -ea/denom
307 if num.any(Y < 0.):
308 continue
310 xtheta = num.arctan(num.sqrt(Y))
311 rtp = num.empty(xphi.shape+(3,), dtype=num.float64)
312 rtp[:, 0] = 1.
313 if sign > 0:
314 rtp[:, 1] = xtheta
315 else:
316 rtp[:, 1] = PI - xtheta
318 rtp[:, 2] = xphi
319 poly_e = numpy_rtp2xyz(rtp)
320 poly = num.dot(from_e, poly_e.T).T
321 poly[:, 2] -= 0.001
323 poly_es.append(poly_e)
324 polys.append(poly)
326 if polys:
327 polys_lower, polys_upper = spoly_cut(polys, 2, arcres=arcres)
328 lines.extend(polys)
329 lines_lower.extend(polys_lower)
330 lines_upper.extend(polys_upper)
332 if poly_es:
333 for aa in spoly_cut(poly_es, 0, arcres=arcres):
334 for bb in spoly_cut(aa, 1, arcres=arcres):
335 for cc in spoly_cut(bb, 2, arcres=arcres):
336 for poly_e in cc:
337 poly = num.dot(from_e, poly_e.T).T
338 poly[:, 2] -= 0.001
339 polys_lower, polys_upper = spoly_cut(
340 [poly], 2, nonsimple=False, arcres=arcres)
342 patches.append(poly)
343 patches_lower.extend(polys_lower)
344 patches_upper.extend(polys_upper)
346 if not patches:
347 if mt_sign * pt_sign == 1.:
348 patches_lower.append(circle_points(aphi, -1.0))
349 patches_upper.append(circle_points(aphi, 1.0))
350 lines_lower.append(circle_points(aphi, -1.0))
351 lines_upper.append(circle_points(aphi, 1.0))
353 groups.append((
354 pt_name,
355 patches, patches_lower, patches_upper,
356 lines, lines_lower, lines_upper))
358 return groups
361def extr(points):
362 pmean = num.mean(points, axis=0)
363 return points + pmean*0.05
366def draw_eigenvectors_mpl(eig, axes):
367 vp, vn, vt = eig[3:]
368 for lab, v in [('P', vp), ('N', vn), ('T', vt)]:
369 sign = num.sign(v[2]) + (v[2] == 0.0)
370 axes.plot(sign*v[1], sign*v[0], 'o', color='black')
371 axes.text(sign*v[1], sign*v[0], ' '+lab)
374def project(points, projection='lambert'):
375 points_out = points[:, :2].copy()
376 if projection == 'lambert':
377 factor = 1.0 / num.sqrt(1.0 + points[:, 2])
378 elif projection == 'stereographic':
379 factor = 1.0 / (1.0 + points[:, 2])
380 elif projection == 'orthographic':
381 factor = None
382 else:
383 raise BeachballError(
384 'invalid argument for projection: %s' % projection)
386 if factor is not None:
387 points_out *= factor[:, num.newaxis]
389 return points_out
392def inverse_project(points, projection='lambert'):
393 points_out = num.zeros((points.shape[0], 3))
395 rsqr = points[:, 0]**2 + points[:, 1]**2
396 if projection == 'lambert':
397 points_out[:, 2] = 1.0 - rsqr
398 points_out[:, 1] = num.sqrt(2.0 - rsqr) * points[:, 1]
399 points_out[:, 0] = num.sqrt(2.0 - rsqr) * points[:, 0]
400 elif projection == 'stereographic':
401 points_out[:, 2] = - (rsqr - 1.0) / (rsqr + 1.0)
402 points_out[:, 1] = 2.0 * points[:, 1] / (rsqr + 1.0)
403 points_out[:, 0] = 2.0 * points[:, 0] / (rsqr + 1.0)
404 elif projection == 'orthographic':
405 points_out[:, 2] = num.sqrt(num.maximum(1.0 - rsqr, 0.0))
406 points_out[:, 1] = points[:, 1]
407 points_out[:, 0] = points[:, 0]
408 else:
409 raise BeachballError(
410 'invalid argument for projection: %s' % projection)
412 return points_out
415def deco_part(mt, mt_type='full', view='top'):
416 mt = mtm.as_mt(mt)
418 if isinstance(view, str):
419 if view == 'top':
420 pass
421 elif view == 'north':
422 mt = mt.rotated(_view_north)
423 elif view == 'south':
424 mt = mt.rotated(_view_south)
425 elif view == 'east':
426 mt = mt.rotated(_view_east)
427 elif view == 'west':
428 mt = mt.rotated(_view_west)
429 elif isinstance(view, tuple):
430 mt = mt.rotated(view_rotation(*view))
431 else:
432 raise BeachballError(
433 'Invaild argument for `view`. Allowed values are "top", "north", '
434 '"south", "east", "west" or a tuple of angles `(strike, dip)` '
435 'orienting the view plane.')
437 if mt_type == 'full':
438 return mt
440 res = mt.standard_decomposition()
441 m = dict(
442 dc=res[1][2],
443 deviatoric=res[3][2])[mt_type]
445 return mtm.MomentTensor(m=m)
448def choose_transform(axes, size_units, position, size):
450 if size_units == 'points':
451 transform = _FixedPointOffsetTransform(
452 axes.transData,
453 axes.figure.dpi_scale_trans,
454 position)
456 if size is None:
457 size = 12.
459 size = size * 0.5 / 72.
460 position = (0., 0.)
462 elif size_units == 'data':
463 transform = axes.transData
465 if size is None:
466 size = 1.0
468 size = size * 0.5
470 elif size_units == 'axes':
471 transform = axes.transAxes
472 if size is None:
473 size = 1.
475 size = size * .5
477 else:
478 raise BeachballError(
479 'invalid argument for size_units: %s' % size_units)
481 position = num.asarray(position, dtype=num.float64)
483 return transform, position, size
486def mt2beachball(
487 mt,
488 beachball_type='deviatoric',
489 position=(0., 0.),
490 size=None,
491 color_t='red',
492 color_p='white',
493 edgecolor='black',
494 linewidth=2,
495 projection='lambert',
496 view='top'):
498 position = num.asarray(position, dtype=num.float64)
499 size = size or 1
500 mt = deco_part(mt, beachball_type, view)
502 eig = mt.eigensystem()
503 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0:
504 raise BeachballError('eigenvalues are zero')
506 data = []
507 for (group, patches, patches_lower, patches_upper,
508 lines, lines_lower, lines_upper) in eig2gx(eig):
510 if group == 'P':
511 color = color_p
512 else:
513 color = color_t
515 for poly in patches_upper:
516 verts = project(poly, projection)[:, ::-1] * size + \
517 position[NA, :]
518 data.append((verts, color, color, 1.0))
520 for poly in lines_upper:
521 verts = project(poly, projection)[:, ::-1] * size + \
522 position[NA, :]
523 data.append((verts, 'none', edgecolor, linewidth))
524 return data
527def plot_beachball_mpl(
528 mt, axes,
529 beachball_type='deviatoric',
530 position=(0., 0.),
531 size=None,
532 zorder=0,
533 color_t='red',
534 color_p='white',
535 edgecolor='black',
536 linewidth=2,
537 alpha=1.0,
538 arcres=181,
539 decimation=1,
540 projection='lambert',
541 size_units='points',
542 view='top'):
544 '''
545 Plot beachball diagram to a Matplotlib plot
547 :param mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
548 array or sequence which can be converted into an MT object
549 :param beachball_type: ``'deviatoric'`` (default), ``'full'``, or ``'dc'``
550 :param position: position of the beachball in data coordinates
551 :param size: diameter of the beachball either in points or in data
552 coordinates, depending on the ``size_units`` setting
553 :param zorder: (passed through to matplotlib drawing functions)
554 :param color_t: color for compressional quadrants (default: ``'red'``)
555 :param color_p: color for extensive quadrants (default: ``'white'``)
556 :param edgecolor: color for lines (default: ``'black'``)
557 :param linewidth: linewidth in points (default: ``2``)
558 :param alpha: (passed through to matplotlib drawing functions)
559 :param projection: ``'lambert'`` (default), ``'stereographic'``, or
560 ``'orthographic'``
561 :param size_units: ``'points'`` (default) or ``'data'``, where the
562 latter causes the beachball to be projected in the plots data
563 coordinates (axes must have an aspect ratio of 1.0 or the
564 beachball will be shown distorted when using this).
565 :param view: View the beachball from ``'top'``, ``'north'``, ``'south'``,
566 ``'east'`` or ``'west'``, or project onto plane given by
567 ``(strike, dip)``. Useful to show beachballs in cross-sections.
568 Default is ``'top'``.
569 '''
571 transform, position, size = choose_transform(
572 axes, size_units, position, size)
574 mt = deco_part(mt, beachball_type, view)
576 eig = mt.eigensystem()
577 if eig[0] == 0. and eig[1] == 0. and eig[2] == 0:
578 raise BeachballError('eigenvalues are zero')
580 data = []
581 for (group, patches, patches_lower, patches_upper,
582 lines, lines_lower, lines_upper) in eig2gx(eig, arcres):
584 if group == 'P':
585 color = color_p
586 else:
587 color = color_t
589 # plot "upper" features for lower hemisphere, because coordinate system
590 # is NED
592 for poly in patches_upper:
593 verts = project(poly, projection)[:, ::-1] * size + position[NA, :]
594 if alpha == 1.0:
595 data.append(
596 (verts[::decimation], color, color, linewidth))
597 else:
598 data.append(
599 (verts[::decimation], color, 'none', 0.0))
601 for poly in lines_upper:
602 verts = project(poly, projection)[:, ::-1] * size + position[NA, :]
603 data.append(
604 (verts[::decimation], 'none', edgecolor, linewidth))
606 patches = []
607 for (path, facecolor, edgecolor, linewidth) in data:
608 patches.append(Polygon(
609 xy=path, facecolor=facecolor,
610 edgecolor=edgecolor,
611 linewidth=linewidth,
612 alpha=alpha))
614 collection = PatchCollection(
615 patches, zorder=zorder, transform=transform, match_original=True)
617 axes.add_artist(collection)
618 return collection
621def amplitudes_ned(mt, vecs):
622 ep, en, et, vp, vn, vt = mt.eigensystem()
623 to_e = num.vstack((vn, vt, vp))
624 vecs_e = num.dot(to_e, vecs.T).T
625 rtp = numpy_xyz2rtp(vecs_e)
626 atheta, aphi = rtp[:, 1], rtp[:, 2]
627 return ep * num.cos(atheta)**2 + (
628 en * num.cos(aphi)**2 + et * num.sin(aphi)**2) * num.sin(atheta)**2
631def amplitudes(mt, azimuths, takeoff_angles):
632 azimuths = num.asarray(azimuths, dtype=float)
633 takeoff_angles = num.asarray(takeoff_angles, dtype=float)
634 assert azimuths.size == takeoff_angles.size
635 rtps = num.vstack(
636 (num.ones(azimuths.size), takeoff_angles*d2r, azimuths*d2r)).T
637 vecs = numpy_rtp2xyz(rtps)
638 return amplitudes_ned(mt, vecs)
641def mts2amps(
642 mts,
643 projection,
644 beachball_type,
645 grid_resolution=200,
646 mask=True,
647 view='top'):
649 n_balls = len(mts)
650 nx = grid_resolution
651 ny = grid_resolution
653 x = num.linspace(-1., 1., nx)
654 y = num.linspace(-1., 1., ny)
656 vecs2 = num.zeros((nx * ny, 2), dtype=num.float64)
657 vecs2[:, 0] = num.tile(x, ny)
658 vecs2[:, 1] = num.repeat(y, nx)
660 ii_ok = vecs2[:, 0]**2 + vecs2[:, 1]**2 <= 1.0
661 amps = num_full(nx * ny, num.nan, dtype=num.float64)
663 amps[ii_ok] = 0.
664 vecs3_ok = inverse_project(vecs2[ii_ok, :], projection)
666 for mt in mts:
667 amps_ok = amplitudes_ned(deco_part(mt, beachball_type, view), vecs3_ok)
668 if mask:
669 amps_ok[amps_ok > 0] = 1.
670 amps_ok[amps_ok < 0] = 0.
672 amps[ii_ok] += amps_ok
674 return num.reshape(amps, (ny, nx)) / n_balls, x, y
677def plot_fuzzy_beachball_mpl_pixmap(
678 mts, axes,
679 best_mt=None,
680 beachball_type='deviatoric',
681 position=(0., 0.),
682 size=None,
683 zorder=0,
684 color_t='red',
685 color_p='white',
686 edgecolor='black',
687 best_color='red',
688 linewidth=2,
689 alpha=1.0,
690 projection='lambert',
691 size_units='data',
692 grid_resolution=200,
693 method='imshow',
694 view='top'):
695 '''
696 Plot fuzzy beachball from a list of given MomentTensors
698 :param mts: list of
699 :py:class:`pyrocko.moment_tensor.MomentTensor` object or an
700 array or sequence which can be converted into an MT object
701 :param best_mt: :py:class:`pyrocko.moment_tensor.MomentTensor` object or
702 an array or sequence which can be converted into an MT object
703 of most likely or minimum misfit solution to extra highlight
704 :param best_color: mpl color for best MomentTensor edges,
705 polygons are not plotted
707 See plot_beachball_mpl for other arguments
708 '''
709 if size_units == 'points':
710 raise BeachballError(
711 'size_units="points" not supported in '
712 'plot_fuzzy_beachball_mpl_pixmap')
714 transform, position, size = choose_transform(
715 axes, size_units, position, size)
717 amps, x, y = mts2amps(
718 mts,
719 grid_resolution=grid_resolution,
720 projection=projection,
721 beachball_type=beachball_type,
722 mask=True,
723 view=view)
725 ncolors = 256
726 cmap = LinearSegmentedColormap.from_list(
727 'dummy', [color_p, color_t], N=ncolors)
729 levels = num.linspace(0, 1., ncolors)
730 if method == 'contourf':
731 axes.contourf(
732 position[0] + y * size, position[1] + x * size, amps.T,
733 levels=levels,
734 cmap=cmap,
735 transform=transform,
736 zorder=zorder,
737 alpha=alpha)
739 elif method == 'imshow':
740 axes.imshow(
741 amps.T,
742 extent=(
743 position[0] + y[0] * size,
744 position[0] + y[-1] * size,
745 position[1] - x[0] * size,
746 position[1] - x[-1] * size),
747 cmap=cmap,
748 transform=transform,
749 zorder=zorder-0.1,
750 alpha=alpha)
751 else:
752 assert False, 'invalid `method` argument'
754 # draw optimum edges
755 if best_mt is not None:
756 best_amps, bx, by = mts2amps(
757 [best_mt],
758 grid_resolution=grid_resolution,
759 projection=projection,
760 beachball_type=beachball_type,
761 mask=False)
763 axes.contour(
764 position[0] + by * size, position[1] + bx * size, best_amps.T,
765 levels=[0.],
766 colors=[best_color],
767 linewidths=linewidth,
768 transform=transform,
769 zorder=zorder,
770 alpha=alpha)
772 phi = num.linspace(0., 2 * PI, 361)
773 x = num.cos(phi)
774 y = num.sin(phi)
775 axes.plot(
776 position[0] + x * size, position[1] + y * size,
777 linewidth=linewidth,
778 color=edgecolor,
779 transform=transform,
780 zorder=zorder,
781 alpha=alpha)
784def plot_beachball_mpl_construction(
785 mt, axes,
786 show='patches',
787 beachball_type='deviatoric',
788 view='top'):
790 mt = deco_part(mt, beachball_type, view)
791 eig = mt.eigensystem()
793 for (group, patches, patches_lower, patches_upper,
794 lines, lines_lower, lines_upper) in eig2gx(eig):
796 if group == 'P':
797 color = 'blue'
798 lw = 1
799 else:
800 color = 'red'
801 lw = 1
803 if show == 'patches':
804 for poly in patches_upper:
805 px, py, pz = poly.T
806 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
808 if show == 'lines':
809 for poly in lines_upper:
810 px, py, pz = poly.T
811 axes.plot(*extr(poly).T, color=color, lw=lw, alpha=0.5)
814def plot_beachball_mpl_pixmap(
815 mt, axes,
816 beachball_type='deviatoric',
817 position=(0., 0.),
818 size=None,
819 zorder=0,
820 color_t='red',
821 color_p='white',
822 edgecolor='black',
823 linewidth=2,
824 alpha=1.0,
825 projection='lambert',
826 size_units='data',
827 view='top'):
829 if size_units == 'points':
830 raise BeachballError(
831 'size_units="points" not supported in plot_beachball_mpl_pixmap')
833 transform, position, size = choose_transform(
834 axes, size_units, position, size)
836 mt = deco_part(mt, beachball_type, view)
838 ep, en, et, vp, vn, vt = mt.eigensystem()
840 amps, x, y = mts2amps(
841 [mt], projection, beachball_type, grid_resolution=200, mask=False)
843 axes.contourf(
844 position[0] + y * size, position[1] + x * size, amps.T,
845 levels=[-num.inf, 0., num.inf],
846 colors=[color_p, color_t],
847 transform=transform,
848 zorder=zorder,
849 alpha=alpha)
851 axes.contour(
852 position[0] + y * size, position[1] + x * size, amps.T,
853 levels=[0.],
854 colors=[edgecolor],
855 linewidths=linewidth,
856 transform=transform,
857 zorder=zorder,
858 alpha=alpha)
860 phi = num.linspace(0., 2 * PI, 361)
861 x = num.cos(phi)
862 y = num.sin(phi)
863 axes.plot(
864 position[0] + x * size, position[1] + y * size,
865 linewidth=linewidth,
866 color=edgecolor,
867 transform=transform,
868 zorder=zorder,
869 alpha=alpha)
872if __name__ == '__main__':
873 import sys
874 import os
875 import matplotlib.pyplot as plt
876 from pyrocko import model
878 args = sys.argv[1:]
880 data = []
881 for iarg, arg in enumerate(args):
883 if os.path.exists(arg):
884 events = model.load_events(arg)
885 for ev in events:
886 if not ev.moment_tensor:
887 logger.warning('no moment tensor given for event')
888 continue
890 data.append((ev.name, ev.moment_tensor))
891 else:
892 vals = list(map(float, arg.split(',')))
893 mt = mtm.as_mt(vals)
894 data.append(('%i' % (iarg+1), mt))
896 n = len(data)
898 ncols = 1
899 while ncols**2 < n:
900 ncols += 1
902 nrows = ncols
904 fig = plt.figure()
905 axes = fig.add_subplot(1, 1, 1, aspect=1.)
906 axes.axison = False
907 axes.set_xlim(-0.05 - ncols, ncols + 0.05)
908 axes.set_ylim(-0.05 - nrows, nrows + 0.05)
910 for ibeach, (name, mt) in enumerate(data):
911 irow = ibeach // ncols
912 icol = ibeach % ncols
913 plot_beachball_mpl(
914 mt, axes,
915 position=(icol*2-ncols+1, -irow*2+nrows-1),
916 size_units='data')
918 axes.annotate(
919 name,
920 xy=(icol*2-ncols+1, -irow*2+nrows-2),
921 xycoords='data',
922 xytext=(0, 0),
923 textcoords='offset points',
924 verticalalignment='center',
925 horizontalalignment='center',
926 rotation=0.)
928 plt.show()