Source code for openquake.hmtk.plotting.beachball

# -*- coding: utf-8 -*-
# -------------------------------------------------------------------
# Filename: beachball.py
#  Purpose: Draws a beach ball diagram of an earthquake focal mechanism.
#   Author: Robert Barsch
#    Email: barsch@geophysik.uni-muenchen.de
#
# Copyright (C) 2008-2012 Robert Barsch
# ---------------------------------------------------------------------

"""
Draws a beachball diagram of an earthquake focal mechanism

Most source code provided here are adopted from

1. MatLab script `bb.m`_ written by Andy Michael and Oliver Boyd.
2. ps_meca program from the `Generic Mapping Tools (GMT)`_.

:copyright:
    The ObsPy Development Team (devs@obspy.org)
:license:
    GNU General Public License (GPL)
    (http://www.gnu.org/licenses/gpl.txt)

.. _`Generic Mapping Tools (GMT)`: http://gmt.soest.hawaii.edu
.. _`bb.m`: http://www.ceri.memphis.edu/people/olboyd/Software/Software.html
"""
import io
import matplotlib.pyplot as plt
from matplotlib import patches, collections, transforms, path as mplpath
import numpy as np
from openquake.baselib.general import Param

D2R = np.pi / 180
R2D = 180 / np.pi
EPSILON = 0.00001


[docs]def Beach(fm, linewidth=2, facecolor="b", bgcolor="w", edgecolor="k", alpha=1.0, xy=(0, 0), width=200, size=100, nofill=False, zorder=100, axes=None): """ Return a beach ball as a collection which can be connected to an current matplotlib axes instance (ax.add_collection). S1, D1, and R1, the strike, dip and rake of one of the focal planes, can be vectors of multiple focal mechanisms. :param fm: Focal mechanism that is either number of mechanisms (NM) by 3 (strike, dip, and rake) or NM x 6 (M11, M22, M33, M12, M13, M23 - the six independent components of the moment tensor, where the coordinate system is 1,2,3 = Up,South,East which equals r,theta,phi). The strike is of the first plane, clockwise relative to north. The dip is of the first plane, defined clockwise and perpendicular to strike, relative to horizontal such that 0 is horizontal and 90 is vertical. The rake is of the first focal plane solution. 90 moves the hanging wall up-dip (thrust), 0 moves it in the strike direction (left-lateral), -90 moves it down-dip (normal), and 180 moves it opposite to strike (right-lateral). :param facecolor: Color to use for quadrants of tension; can be a string, e.g. ``'r'``, ``'b'`` or three component color vector, [R G B]. Defaults to ``'b'`` (blue). :param bgcolor: The background color. Defaults to ``'w'`` (white). :param edgecolor: Color of the edges. Defaults to ``'k'`` (black). :param alpha: The alpha level of the beach ball. Defaults to ``1.0`` (opaque). :param xy: Origin position of the beach ball as tuple. Defaults to ``(0, 0)``. :type width: int or tuple :param width: Symbol size of beach ball, or tuple for elliptically shaped patches. Defaults to size ``200``. :param size: Controls the number of interpolation points for the curves. Minimum is automatically set to ``100``. :param nofill: Do not fill the beach ball, but only plot the planes. :param zorder: Set zorder. Artists with lower zorder values are drawn first. :type axes: :class:`matplotlib.axes.Axes` :param axes: Used to make beach balls circular on non-scaled axes. Also maintains the aspect ratio when resizing the figure. Will not add the returned collection to the axes instance. """ # check if one or two widths are specified (Circle or Ellipse) try: assert len(width) == 2 except TypeError: width = (width, width) mt = None np1 = None if isinstance(fm, MomentTensor): mt = fm np1 = MT2Plane(mt) elif isinstance(fm, NodalPlane): np1 = fm elif len(fm) == 6: mt = MomentTensor(fm[0], fm[1], fm[2], fm[3], fm[4], fm[5], 0) np1 = MT2Plane(mt) elif len(fm) == 3: np1 = NodalPlane(fm[0], fm[1], fm[2]) else: raise TypeError("Wrong input value for 'fm'.") # Only at least size 100, i.e. 100 points in the matrix are allowed if size < 100: size = 100 # Return as collection if mt: (T, N, P) = MT2Axes(mt) if np.fabs(N.val) < EPSILON and np.fabs(T.val + P.val) < EPSILON: colors, p = plotDC(np1, size, xy=xy, width=width) else: colors, p = plotMT( T, N, P, size, plot_zerotrace=True, xy=xy, width=width ) else: colors, p = plotDC(np1, size=size, xy=xy, width=width) if nofill: # XXX: not tested with plotMT col = collections.PatchCollection([p[1]], match_original=False) col.set_facecolor("none") else: col = collections.PatchCollection(p, match_original=False) # Replace color dummies 'b' and 'w' by face and bgcolor fc = [facecolor if c == "b" else bgcolor for c in colors] col.set_facecolors(fc) # Use the given axes to maintain the aspect ratio of beachballs on figure # resize. if axes is not None: # This is what holds the aspect ratio (but breaks the positioning) col.set_transform(transforms.IdentityTransform()) # Next is a dirty hack to fix the positioning: # 1. Need to bring the all patches to the origin (0, 0). for p in col._paths: p.vertices -= xy # 2. Then use the offset property of the collection to position the # patches col.set_offsets(xy) col._transOffset = axes.transData col.set_edgecolor(edgecolor) col.set_alpha(alpha) col.set_linewidth(linewidth) col.set_zorder(zorder) return col
[docs]def Beachball(fm, linewidth=2, facecolor="b", bgcolor="w", edgecolor="k", alpha=1.0, xy=(0, 0), width=200, size=100, nofill=False, zorder=100, outfile=None, format=None, fig=None): """ Draws a beach ball diagram of an earthquake focal mechanism. S1, D1, and R1, the strike, dip and rake of one of the focal planes, can be vectors of multiple focal mechanisms. :param fm: Focal mechanism that is either number of mechanisms (NM) by 3 (strike, dip, and rake) or NM x 6 (M11, M22, M33, M12, M13, M23 - the six independent components of the moment tensor, where the coordinate system is 1,2,3 = Up,South,East which equals r,theta,phi). The strike is of the first plane, clockwise relative to north. The dip is of the first plane, defined clockwise and perpendicular to strike, relative to horizontal such that 0 is horizontal and 90 is vertical. The rake is of the first focal plane solution. 90 moves the hanging wall up-dip (thrust), 0 moves it in the strike direction (left-lateral), -90 moves it down-dip (normal), and 180 moves it opposite to strike (right-lateral). :param facecolor: Color to use for quadrants of tension; can be a string, e.g. ``'r'``, ``'b'`` or three component color vector, [R G B]. Defaults to ``'b'`` (blue). :param bgcolor: The background color. Defaults to ``'w'`` (white). :param edgecolor: Color of the edges. Defaults to ``'k'`` (black). :param alpha: The alpha level of the beach ball. Defaults to ``1.0`` (opaque). :param xy: Origin position of the beach ball as tuple. Defaults to ``(0, 0)``. :type width: int :param width: Symbol size of beach ball. Defaults to ``200``. :param size: Controls the number of interpolation points for the curves. Minimum is automatically set to ``100``. :param nofill: Do not fill the beach ball, but only plot the planes. :param zorder: Set zorder. Artists with lower zorder values are drawn first. :param outfile: Output file string. Also used to automatically determine the output format. Supported file formats depend on your matplotlib backend. Most backends support png, pdf, ps, eps and svg. Defaults to ``None``. :param format: Format of the graph picture. If no format is given the outfile parameter will be used to try to automatically determine the output format. If no format is found it defaults to png output. If no outfile is specified but a format is, than a binary imagestring will be returned. Defaults to ``None``. :param fig: Give an existing figure instance to plot into. New Figure if set to ``None``. """ plot_width = width * 0.95 # plot the figure if not fig: fig = plt.figure(figsize=(3, 3), dpi=100) fig.subplots_adjust(left=0, bottom=0, right=1, top=1) fig.set_figheight(width // 100) fig.set_figwidth(width // 100) ax = fig.add_subplot(111, aspect="equal") # hide axes + ticks ax.axison = False # plot the collection collection = Beach( fm, linewidth=linewidth, facecolor=facecolor, edgecolor=edgecolor, bgcolor=bgcolor, alpha=alpha, nofill=nofill, xy=xy, width=plot_width, size=size, zorder=zorder, ) ax.add_collection(collection) ax.autoscale_view(tight=False, scalex=True, scaley=True) # export if outfile: if format: fig.savefig(outfile, dpi=100, transparent=True, format=format) else: fig.savefig(outfile, dpi=100, transparent=True) elif format and not outfile: imgdata = io.BytesIO() fig.savefig(imgdata, format=format, dpi=100, transparent=True) imgdata.seek(0) return imgdata.read() else: plt.show() return fig
[docs]def plotMT( T, N, P, size=200, plot_zerotrace=True, x0=0, y0=0, xy=(0, 0), width=200 ): """ Uses a principal axis T, N and P to draw a beach ball plot. :param ax: axis object of a matplotlib figure :param T: :class:`~PrincipalAxis` :param N: :class:`~PrincipalAxis` :param P: :class:`~PrincipalAxis` Adapted from ps_tensor / utilmeca.c / `Generic Mapping Tools (GMT)`_. .. _`Generic Mapping Tools (GMT)`: http://gmt.soest.hawaii.edu """ # check if one or two widths are specified (Circle or Ellipse) try: assert len(width) == 2 except TypeError: width = (width, width) collect = [] colors = [] p = Param( # 28 parameters! azi=np.zeros((3, 2)), res=[value / float(size) for value in width], a=np.zeros(3), p=np.zeros(3), v=np.zeros(3), b=1, d=2, m=0, big_iso=0, j=1, j2=0, j3=0, n=0, width=width, xy=xy, radius_size=size * 0.5, x0=x0, y0=y0, x=np.zeros(400), y=np.zeros(400), x2=np.zeros(400), y2=np.zeros(400), x3=np.zeros(400), y3=np.zeros(400), xp1=np.zeros(800), yp1=np.zeros(800), xp2=np.zeros(400), yp2=np.zeros(400), ) p.a[0] = T.strike p.a[1] = N.strike p.a[2] = P.strike p.p[0] = T.dip p.p[1] = N.dip p.p[2] = P.dip p.v[0] = T.val p.v[1] = N.val p.v[2] = P.val vi = (p.v[0] + p.v[1] + p.v[2]) / 3.0 for i in range(3): p.v[i] = p.v[i] - vi if np.fabs(p.v[0] * p.v[0] + p.v[1] * p.v[1] + p.v[2] * p.v[2]) < EPSILON: # pure implosion-explosion if vi > 0.0: cir = patches.Ellipse(xy, width=width[0], height=width[1]) collect.append(cir) colors.append("b") if vi < 0.0: cir = patches.Ellipse(xy, width=width[0], height=width[1]) collect.append(cir) colors.append("w") return colors, collect if np.fabs(p.v[0]) >= np.fabs(p.v[2]): p.d = 0 p.m = 2 if plot_zerotrace: vi = 0.0 f = -p.v[1] / float(p.v[p.d]) iso = vi / float(p.v[p.d]) return _plot(p, f, iso, collect, colors)
def _plot1(p, rgb1, colors, collect): for i in range(p.j): p.xp1[i] = p.x[i] p.yp1[i] = p.y[i] if p.azi[0][0] - p.azi[0][1] > np.pi: p.azi[0][0] -= np.pi * 2.0 elif p.azi[0][1] - p.azi[0][0] > np.pi: p.azi[0][0] += np.pi * 2.0 if p.azi[0][0] < p.azi[0][1]: az = p.azi[0][1] - D2R while az > p.azi[0][0]: si = np.sin(az) co = np.cos(az) p.xp1[i] = p.x0 + p.radius_size * si p.yp1[i] = p.y0 + p.radius_size * co i += 1 az -= D2R else: az = p.azi[0][1] + D2R while az < p.azi[0][0]: si = np.sin(az) co = np.cos(az) p.xp1[i] = p.x0 + p.radius_size * si p.yp1[i] = p.y0 + p.radius_size * co i += 1 az += D2R collect.append(xy2patch(p.xp1[:i], p.yp1[:i], p.res, p.xy)) colors.append(rgb1) for i in range(p.j2): p.xp2[i] = p.x2[i] p.yp2[i] = p.y2[i] if p.azi[1][0] - p.azi[1][1] > np.pi: p.azi[1][0] -= np.pi * 2.0 elif p.azi[1][1] - p.azi[1][0] > np.pi: p.azi[1][0] += np.pi * 2.0 if p.azi[1][0] < p.azi[1][1]: az = p.azi[1][1] - D2R while az > p.azi[1][0]: si = np.sin(az) co = np.cos(az) p.xp2[i] = p.x0 + p.radius_size * si i += 1 p.yp2[i] = p.y0 + p.radius_size * co az -= D2R else: az = p.azi[1][1] + D2R while az < p.azi[1][0]: si = np.sin(az) co = np.cos(az) p.xp2[i] = p.x0 + p.radius_size * si i += 1 p.yp2[i] = p.y0 + p.radius_size * co az += D2R collect.append(xy2patch(p.xp2[:i], p.yp2[:i], p.res, p.xy)) colors.append(rgb1) return colors, collect def _plot2(p, rgb1, colors, collect): for i in range(p.j3): p.xp1[i] = p.x3[i] p.yp1[i] = p.y3[i] for ii in range(p.j): p.xp1[i] = p.x[ii] i += 1 p.yp1[i] = p.y[ii] if p.big_iso: ii = p.j2 - 1 while ii >= 0: p.xp1[i] = p.x2[ii] i += 1 p.yp1[i] = p.y2[ii] ii -= 1 collect.append(xy2patch(p.xp1[:i], p.yp1[:i], p.res, p.xy)) colors.append(rgb1) return colors, collect if p.azi[2][0] - p.azi[0][1] > np.pi: p.azi[2][0] -= np.pi * 2.0 elif p.azi[0][1] - p.azi[2][0] > np.pi: p.azi[2][0] += np.pi * 2.0 if p.azi[2][0] < p.azi[0][1]: az = p.azi[0][1] - D2R while az > p.azi[2][0]: si = np.sin(az) co = np.cos(az) p.xp1[i] = p.x0 + p.radius_size * si i += 1 p.yp1[i] = p.y0 + p.radius_size * co az -= D2R else: az = p.azi[0][1] + D2R while az < p.azi[2][0]: si = np.sin(az) co = np.cos(az) p.xp1[i] = p.x0 + p.radius_size * si i += 1 p.yp1[i] = p.y0 + p.radius_size * co az += D2R collect.append(xy2patch(p.xp1[:i], p.yp1[:i], p.res, p.xy)) colors.append(rgb1) for i in range(p.j2): p.xp2[i] = p.x2[i] p.yp2[i] = p.y2[i] if p.azi[1][0] - p.azi[1][1] > np.pi: p.azi[1][0] -= np.pi * 2.0 elif p.azi[1][1] - p.azi[1][0] > np.pi: p.azi[1][0] += np.pi * 2.0 if p.azi[1][0] < p.azi[1][1]: az = p.azi[1][1] - D2R while az > p.azi[1][0]: si = np.sin(az) co = np.cos(az) p.xp2[i] = p.x0 + p.radius_size * si i += 1 p.yp2[i] = p.y0 + p.radius_size * co az -= D2R else: az = p.azi[1][1] + D2R while az < p.azi[1][0]: si = np.sin(az) co = np.cos(az) p.xp2[i] = p.x0 + p.radius_size * si i += 1 p.yp2[i] = p.y0 + p.radius_size * co az += D2R collect.append(xy2patch(p.xp2[:i], p.yp2[:i], p.res, p.xy)) colors.append(rgb1) return colors, collect def _update_azi(p, f, iso): b, d, m = p.b, p.d, p.m spd = np.sin(p.p[d] * D2R) cpd = np.cos(p.p[d] * D2R) spb = np.sin(p.p[b] * D2R) cpb = np.cos(p.p[b] * D2R) spm = np.sin(p.p[m] * D2R) cpm = np.cos(p.p[m] * D2R) sad = np.sin(p.a[d] * D2R) cad = np.cos(p.a[d] * D2R) sab = np.sin(p.a[b] * D2R) cab = np.cos(p.a[b] * D2R) sam = np.sin(p.a[m] * D2R) cam = np.cos(p.a[m] * D2R) for i in range(360): fir = i * D2R s2alphan = (2.0 + 2.0 * iso) / float( 3.0 + (1.0 - 2.0 * f) * np.cos(2.0 * fir) ) if s2alphan > 1.0: p.big_iso += 1 else: alphan = np.arcsin(np.sqrt(s2alphan)) sfi = np.sin(fir) cfi = np.cos(fir) san = np.sin(alphan) can = np.cos(alphan) xz = can * spd + san * sfi * spb + san * cfi * spm xn = ( can * cpd * cad + san * sfi * cpb * cab + san * cfi * cpm * cam ) xe = ( can * cpd * sad + san * sfi * cpb * sab + san * cfi * cpm * sam ) if np.fabs(xn) < EPSILON and np.fabs(xe) < EPSILON: takeoff = 0.0 az = 0.0 else: az = np.arctan2(xe, xn) if az < 0.0: az += np.pi * 2.0 takeoff = np.arccos( xz / float(np.sqrt(xz * xz + xn * xn + xe * xe)) ) if takeoff > np.pi / 2.0: takeoff = np.pi - takeoff az += np.pi if az > np.pi * 2.0: az -= np.pi * 2.0 r = np.sqrt(2) * np.sin(takeoff / 2.0) si = np.sin(az) co = np.cos(az) if i == 0: p.azi[i][0] = az p.x[i] = p.x0 + p.radius_size * r * si p.y[i] = p.y0 + p.radius_size * r * co azp = az else: if np.fabs(np.fabs(az - azp) - np.pi) < D2R * 10.0: p.azi[p.n][1] = p.azp p.n += 1 p.azi[p.n][0] = p.az if np.fabs(np.fabs(az - azp) - np.pi * 2.0) < D2R * 2.0: if azp < az: p.azi[p.n][0] += np.pi * 2.0 else: p.azi[p.n][0] -= np.pi * 2.0 if p.n == 0: p.x[p.j] = p.x0 + p.radius_size * r * si p.y[p.j] = p.y0 + p.radius_size * r * co p.j += 1 elif p.n == 1: p.x2[p.j2] = p.x0 + p.radius_size * r * si p.y2[p.j2] = p.y0 + p.radius_size * r * co p.j2 += 1 elif p.n == 2: p.x3[p.j3] = p.x0 + p.radius_size * r * si p.y3[p.j3] = p.y0 + p.radius_size * r * co p.j3 += 1 azp = az p.azi[p.n][1] = az def _plot(p, f, iso, collect, colors): # Cliff Frohlich, Seismological Research letters, # Vol 7, Number 1, January-February, 1996 # Unless the isotropic parameter lies in the range # between -1 and 1 - f there will be no nodes whatsoever if iso < -1: cir = patches.Ellipse(p.xy, width=p.width[0], height=p.width[1]) collect.append(cir) colors.append("w") return colors, collect elif iso > 1 - f: cir = patches.Ellipse(p.xy, width=p.width[0], height=p.width[1]) collect.append(cir) colors.append("b") return colors, collect _update_azi(p, f, iso) if p.v[1] < 0.0: rgb1 = "b" rgb2 = "w" else: rgb1 = "w" rgb2 = "b" cir = patches.Ellipse(p.xy, width=p.width[0], height=p.width[1]) collect.append(cir) colors.append(rgb2) if p.n == 0: collect.append(xy2patch(p.x[:360], p.y[:360], p.res, p.xy)) colors.append(rgb1) return colors, collect elif p.n == 1: return _plot1(p, rgb1, colors, collect) elif p.n == 2: return _plot2(p, rgb1, colors, collect)
[docs]def plotDC(np1, size=200, xy=(0, 0), width=200): """ Uses one nodal plane of a double couple to draw a beach ball plot. :param ax: axis object of a matplotlib figure :param np1: :class:`~NodalPlane` Adapted from MATLAB script `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_ written by Andy Michael and Oliver Boyd. """ # check if one or two widths are specified (Circle or Ellipse) try: assert len(width) == 2 except TypeError: width = (width, width) S1 = np1.strike D1 = np1.dip R1 = np1.rake M = 0 if R1 > 180: R1 -= 180 M = 1 if R1 < 0: R1 += 180 M = 1 # Get azimuth and dip of second plane (S2, D2, _R2) = AuxPlane(S1, D1, R1) D = size / 2 if D1 >= 90: D1 = 89.9999 if D2 >= 90: D2 = 89.9999 # arange checked for numerical stablility, np.pi is not multiple of 0.1 phi = np.arange(np.pi, 0.01) l1 = np.sqrt( np.power(90 - D1, 2) / ( np.power(np.sin(phi), 2) + np.power(np.cos(phi), 2) * np.power(90 - D1, 2) / np.power(90, 2) ) ) l2 = np.sqrt( np.power(90 - D2, 2) / ( np.power(np.sin(phi), 2) + np.power(np.cos(phi), 2) * np.power(90 - D2, 2) / np.power(90, 2) ) ) inc = 1 (X1, Y1) = Pol2Cart(phi + S1 * D2R, l1) if M == 1: lo = S1 - 180 hi = S2 if lo > hi: inc = -1 th1 = np.arange(S1 - 180, S2, inc) (Xs1, Ys1) = Pol2Cart(th1 * D2R, 90 * np.ones((1, len(th1)))) (X2, Y2) = Pol2Cart(phi + S2 * D2R, l2) th2 = np.arange(S2 + 180, S1, -inc) else: hi = S1 - 180 lo = S2 - 180 if lo > hi: inc = -1 th1 = np.arange(hi, lo, -inc) (Xs1, Ys1) = Pol2Cart(th1 * D2R, 90 * np.ones((1, len(th1)))) (X2, Y2) = Pol2Cart(phi + S2 * D2R, l2) X2 = X2[::-1] Y2 = Y2[::-1] th2 = np.arange(S2, S1, inc) (Xs2, Ys2) = Pol2Cart(th2 * D2R, 90 * np.ones((1, len(th2)))) X = np.concatenate((X1, Xs1[0], X2, Xs2[0])) Y = np.concatenate((Y1, Ys1[0], Y2, Ys2[0])) X = X * D / 90 Y = Y * D / 90 # calculate resolution res = [value / float(size) for value in width] # construct the patches collect = [patches.Ellipse(xy, width=width[0], height=width[1])] collect.append(xy2patch(Y, X, res, xy)) return ["b", "w"], collect
[docs]def xy2patch(x, y, res, xy): # check if one or two resolutions are specified (Circle or Ellipse) try: assert len(res) == 2 except TypeError: res = (res, res) # transform into the Path coordinate system x = x * res[0] + xy[0] y = y * res[1] + xy[1] verts = list(zip(x.tolist(), y.tolist())) codes = [mplpath.Path.MOVETO] codes.extend([mplpath.Path.LINETO] * (len(x) - 2)) codes.append(mplpath.Path.CLOSEPOLY) path = mplpath.Path(verts, codes) return patches.PathPatch(path)
[docs]def Pol2Cart(th, r): """ """ x = r * np.cos(th) y = r * np.sin(th) return (x, y)
[docs]def StrikeDip(n, e, u): """ Finds strike and dip of plane given normal vector having components n, e, and u. Adapted from MATLAB script `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_ written by Andy Michael and Oliver Boyd. """ r2d = 180 / np.pi if u < 0: n = -n e = -e u = -u strike = np.arctan2(e, n) * r2d strike = strike - 90 while strike >= 360: strike = strike - 360 while strike < 0: strike = strike + 360 x = np.sqrt(np.power(n, 2) + np.power(e, 2)) dip = np.arctan2(x, u) * r2d return (strike, dip)
[docs]def AuxPlane(s1, d1, r1): """ Get Strike and dip of second plane. Adapted from MATLAB script `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_ written by Andy Michael and Oliver Boyd. """ r2d = 180 / np.pi z = (s1 + 90) / r2d z2 = d1 / r2d z3 = r1 / r2d # slick vector in plane 1 sl1 = -np.cos(z3) * np.cos(z) - np.sin(z3) * np.sin(z) * np.cos(z2) sl2 = np.cos(z3) * np.sin(z) - np.sin(z3) * np.cos(z) * np.cos(z2) sl3 = np.sin(z3) * np.sin(z2) (strike, dip) = StrikeDip(sl2, sl1, sl3) n1 = np.sin(z) * np.sin(z2) # normal vector to plane 1 n2 = np.cos(z) * np.sin(z2) h1 = -sl2 # strike vector of plane 2 h2 = sl1 # note h3=0 always so we leave it out # n3 = np.cos(z2) z = h1 * n1 + h2 * n2 z = z / np.sqrt(h1 * h1 + h2 * h2) z = np.arccos(z) rake = 0 if sl3 > 0: rake = z * r2d if sl3 <= 0: rake = -z * r2d return (strike, dip, rake)
[docs]def MT2Plane(mt): """ Calculates a nodal plane of a given moment tensor. :param mt: :class:`~MomentTensor` :return: :class:`~NodalPlane` Adapted from MATLAB script `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_ written by Andy Michael and Oliver Boyd. """ (d, v) = np.linalg.eig(mt.mt) D = np.array([d[1], d[0], d[2]]) V = np.array( [ [v[1, 1], -v[1, 0], -v[1, 2]], [v[2, 1], -v[2, 0], -v[2, 2]], [-v[0, 1], v[0, 0], v[0, 2]], ] ) IMAX = D.argmax() IMIN = D.argmin() AE = (V[:, IMAX] + V[:, IMIN]) / np.sqrt(2.0) AN = (V[:, IMAX] - V[:, IMIN]) / np.sqrt(2.0) AER = np.sqrt(np.power(AE[0], 2) + np.power(AE[1], 2) + np.power(AE[2], 2)) ANR = np.sqrt(np.power(AN[0], 2) + np.power(AN[1], 2) + np.power(AN[2], 2)) AE = AE / AER if not ANR: AN = np.array([np.nan, np.nan, np.nan]) else: AN = AN / ANR if AN[2] <= 0.0: AN1 = AN AE1 = AE else: AN1 = -AN AE1 = -AE (ft, fd, fl) = TDL(AN1, AE1) return NodalPlane(360 - ft, fd, 180 - fl)
[docs]def TDL(AN, BN): """ Helper function for MT2Plane. Adapted from MATLAB script `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_ written by Andy Michael and Oliver Boyd. """ XN = AN[0] YN = AN[1] ZN = AN[2] XE = BN[0] YE = BN[1] ZE = BN[2] AAA = 1.0 / (1000000) CON = 57.2957795 if np.fabs(ZN) < AAA: FD = 90.0 AXN = np.fabs(XN) if AXN > 1.0: AXN = 1.0 FT = np.arcsin(AXN) * CON ST = -XN CT = YN if ST >= 0.0 and CT < 0: FT = 180.0 - FT if ST < 0.0 and CT <= 0: FT = 180.0 + FT if ST < 0.0 and CT > 0: FT = 360.0 - FT FL = np.arcsin(abs(ZE)) * CON SL = -ZE if np.fabs(XN) < AAA: CL = XE / YN else: CL = -YE / XN if SL >= 0.0 and CL < 0: FL = 180.0 - FL if SL < 0.0 and CL <= 0: FL = FL - 180.0 if SL < 0.0 and CL > 0: FL = -FL else: if -ZN > 1.0: ZN = -1.0 FDH = np.arccos(-ZN) FD = FDH * CON SD = np.sin(FDH) if SD == 0: return ST = -XN / SD CT = YN / SD SX = np.fabs(ST) if SX > 1.0: SX = 1.0 FT = np.arcsin(SX) * CON if ST >= 0.0 and CT < 0: FT = 180.0 - FT if ST < 0.0 and CT <= 0: FT = 180.0 + FT if ST < 0.0 and CT > 0: FT = 360.0 - FT SL = -ZE / SD SX = np.fabs(SL) if SX > 1.0: SX = 1.0 FL = np.arcsin(SX) * CON if ST == 0: CL = XE / CT else: XXX = YN * ZN * ZE / SD / SD + YE CL = -SD * XXX / XN if CT == 0: CL = YE / ST if SL >= 0.0 and CL < 0: FL = 180.0 - FL if SL < 0.0 and CL <= 0: FL = FL - 180.0 if SL < 0.0 and CL > 0: FL = -FL return (FT, FD, FL)
[docs]def MT2Axes(mt): """ Calculates the principal axes of a given moment tensor. :param mt: :class:`~MomentTensor` :return: tuple of :class:`~PrincipalAxis` T, N and P Adapted from ps_tensor / utilmeca.c / `Generic Mapping Tools (GMT) <http://gmt.soest.hawaii.edu>`_. """ (D, V) = np.linalg.eigh(mt.mt) pl = np.arcsin(-V[0]) az = np.arctan2(V[2], -V[1]) for i in range(3): if pl[i] <= 0: pl[i] = -pl[i] az[i] += np.pi if az[i] < 0: az[i] += 2 * np.pi if az[i] > 2 * np.pi: az[i] -= 2 * np.pi pl *= R2D az *= R2D T = PrincipalAxis(D[2], az[2], pl[2]) N = PrincipalAxis(D[1], az[1], pl[1]) P = PrincipalAxis(D[0], az[0], pl[0]) return (T, N, P)
[docs]class PrincipalAxis(object): """ A principal axis. Strike and dip values are in degrees. >>> a = PrincipalAxis(1.3, 20, 50) >>> a.dip 50 >>> a.strike 20 >>> a.val 1.3 """ def __init__(self, val=0, strike=0, dip=0): self.val = val self.strike = strike self.dip = dip
[docs]class NodalPlane(object): """ A nodal plane. All values are in degrees. >>> a = NodalPlane(13, 20, 50) >>> a.strike 13 >>> a.dip 20 >>> a.rake 50 """ def __init__(self, strike=0, dip=0, rake=0): self.strike = strike self.dip = dip self.rake = rake
[docs]class MomentTensor(object): """ A moment tensor. >>> a = MomentTensor(1, 1, 0, 0, 0, -1, 26) >>> b = MomentTensor(np.array([1, 1, 0, 0, 0, -1]), 26) >>> c = MomentTensor(np.array([[1, 0, 0], [0, 1, -1], [0, -1, 0]]), 26) >>> a.mt array([[ 1, 0, 0], [ 0, 1, -1], [ 0, -1, 0]]) >>> b.yz -1 >>> a.expo 26 """ def __init__(self, *args): if len(args) == 2: A = args[0] self.expo = args[1] if len(A) == 6: # six independent components self.mt = np.array( [ [A[0], A[3], A[4]], [A[3], A[1], A[5]], [A[4], A[5], A[2]], ] ) elif isinstance(A, np.ndarray) and A.shape == (3, 3): # full matrix self.mt = A else: raise TypeError("Wrong size of input parameter.") elif len(args) == 7: # six independent components self.mt = np.array( [ [args[0], args[3], args[4]], [args[3], args[1], args[5]], [args[4], args[5], args[2]], ] ) self.expo = args[6] else: raise TypeError("Wrong size of input parameter.") @property def xx(self): return self.mt[0][0] @property def xy(self): return self.mt[0][1] @property def xz(self): return self.mt[0][2] @property def yz(self): return self.mt[1][2] @property def yy(self): return self.mt[1][1] @property def zz(self): return self.mt[2][2]
if __name__ == "__main__": import doctest doctest.testmod()