Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/moment_tensor.py: 87%
493 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-23 12:34 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7This module provides various moment tensor related utility functions.
9It can be used to convert between strike-dip-rake and moment tensor
10representations and provides different options to produce random moment
11tensors.
13Moment tensors are represented by :py:class:`MomentTensor` instances. The
14internal representation uses a north-east-down (NED) coordinate system, but it
15can convert from/to the conventions used by the Global CMT catalog
16(up-south-east, USE).
18If not otherwise noted, scalar moment is interpreted as the Frobenius norm
19based scalar moment (see :py:meth:`MomentTensor.scalar_moment`. The scalar
20moment according to the "standard decomposition" can be found in the
21output of :py:meth:`MomentTensor.standard_decomposition`.
22'''
24import math
25import numpy as num
27from .guts import Object, Float
29guts_prefix = 'pf'
31dynecm = 1e-7
34def random_axis(rstate=None):
35 '''
36 Get randomly oriented unit vector.
38 :param rstate: :py:class:`numpy.random.RandomState` object, can be used to
39 create reproducible pseudo-random sequences
40 '''
41 rstate = rstate or num.random
42 while True:
43 axis = rstate.uniform(size=3) * 2.0 - 1.0
44 uabs = math.sqrt(num.sum(axis**2))
45 if 0.001 < uabs < 1.0:
46 return axis / uabs
49def rotation_from_angle_and_axis(angle, axis):
50 '''
51 Build rotation matrix based on axis and angle.
53 :param angle: rotation angle [degrees]
54 :param axis: orientation of rotation axis, either in spherical
55 coordinates ``(theta, phi)`` [degrees], or as a unit vector
56 ``(ux, uy, uz)``.
57 '''
59 if len(axis) == 2:
60 theta, phi = axis
61 ux = math.sin(d2r*theta)*math.cos(d2r*phi)
62 uy = math.sin(d2r*theta)*math.sin(d2r*phi)
63 uz = math.cos(d2r*theta)
65 elif len(axis) == 3:
66 axis = num.asarray(axis)
67 uabs = math.sqrt(num.sum(axis**2))
68 ux, uy, uz = axis / uabs
69 else:
70 assert False
72 ct = math.cos(d2r*angle)
73 st = math.sin(d2r*angle)
74 return num.array([
75 [ct + ux**2*(1.-ct), ux*uy*(1.-ct)-uz*st, ux*uz*(1.-ct)+uy*st],
76 [uy*ux*(1.-ct)+uz*st, ct+uy**2*(1.-ct), uy*uz*(1.-ct)-ux*st],
77 [uz*ux*(1.-ct)-uy*st, uz*uy*(1.-ct)+ux*st, ct+uz**2*(1.-ct)]
78 ])
81def random_rotation(x=None):
82 '''
83 Get random rotation matrix.
85 A random rotation matrix, drawn from a uniform distribution in the space
86 of rotations is returned, after Avro 1992 - "Fast random rotation
87 matrices".
89 :param x: three (uniform random) numbers in the range [0, 1[ used as input
90 to the distribution transformation. If ``None``, random numbers are
91 used. Can be used to create grids of random rotations with uniform
92 density in rotation space.
93 '''
95 if x is not None:
96 x1, x2, x3 = x
97 else:
98 x1, x2, x3 = num.random.random(3)
100 phi = math.pi*2.0*x1
102 zrot = num.array([
103 [math.cos(phi), math.sin(phi), 0.],
104 [-math.sin(phi), math.cos(phi), 0.],
105 [0., 0., 1.]])
107 lam = math.pi*2.0*x2
109 v = num.array([[
110 math.cos(lam)*math.sqrt(x3),
111 math.sin(lam)*math.sqrt(x3),
112 math.sqrt(1.-x3)]]).T
114 house = num.identity(3) - 2.0 * num.dot(v, v.T)
115 return -num.dot(house, zrot)
118def rand(mi, ma):
119 return float(num.random.uniform(mi, ma))
122def randdip(mi, ma):
123 mi_ = 0.5*(math.cos(mi * math.pi/180.)+1.)
124 ma_ = 0.5*(math.cos(ma * math.pi/180.)+1.)
125 return math.acos(rand(mi_, ma_)*2.-1.)*180./math.pi
128def random_strike_dip_rake(
129 strikemin=0., strikemax=360.,
130 dipmin=0., dipmax=90.,
131 rakemin=-180., rakemax=180.):
133 '''
134 Get random strike, dip, rake triplet.
136 .. note::
138 Might not produce a homogeneous distribution of mechanisms. Better use
139 :py:meth:`MomentTensor.random_dc` which is based on
140 :py:func:`random_rotation`.
141 '''
143 strike = rand(strikemin, strikemax)
144 dip = randdip(dipmin, dipmax)
145 rake = rand(rakemin, rakemax)
147 return strike, dip, rake
150def to6(m):
151 '''
152 Get non-redundant components from symmetric 3x3 matrix.
154 :returns: 1D NumPy array with entries ordered like
155 ``(a_xx, a_yy, a_zz, a_xy, a_xz, a_yz)``
156 '''
158 return num.array([m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2]])
161def symmat6(a_xx, a_yy, a_zz, a_xy, a_xz, a_yz):
162 '''
163 Create symmetric 3x3 matrix from its 6 non-redundant values.
164 '''
166 return num.array([[a_xx, a_xy, a_xz],
167 [a_xy, a_yy, a_yz],
168 [a_xz, a_yz, a_zz]], dtype=float)
171def values_to_matrix(values):
172 '''
173 Convert anything to moment tensor represented as a NumPy array.
175 Transforms :py:class:`MomentTensor` objects, tuples, lists and NumPy arrays
176 with 3x3 or 3, 4, 6, or 7 elements into NumPy 3x3 array objects.
178 The ``values`` argument is interpreted depending on shape and type as
179 follows:
181 * ``(strike, dip, rake)``
182 * ``(strike, dip, rake, magnitude)``
183 * ``(mnn, mee, mdd, mne, mnd, med)``
184 * ``(mnn, mee, mdd, mne, mnd, med, magnitude)``
185 * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))``
186 * :py:class:`MomentTensor`
187 '''
189 if isinstance(values, (tuple, list)):
190 values = num.asarray(values, dtype=float)
192 if isinstance(values, MomentTensor):
193 return values.m()
195 elif isinstance(values, num.ndarray):
196 if values.shape == (3,):
197 strike, dip, rake = values
198 rotmat1 = euler_to_matrix(d2r*dip, d2r*strike, -d2r*rake)
199 return num.dot(rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1))
201 elif values.shape == (4,):
202 strike, dip, rake, magnitude = values
203 moment = magnitude_to_moment(magnitude)
204 rotmat1 = euler_to_matrix(d2r*dip, d2r*strike, -d2r*rake)
205 return moment * num.dot(
206 rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1))
208 elif values.shape == (6,):
209 return symmat6(*values)
211 elif values.shape == (7,):
212 magnitude = values[6]
213 moment = magnitude_to_moment(magnitude)
214 mt = symmat6(*values[:6])
215 mt *= moment / (num.linalg.norm(mt) / math.sqrt(2.0))
216 return mt
218 elif values.shape == (3, 3):
219 return num.asarray(values, dtype=float)
221 raise Exception('cannot convert object to 3x3 matrix')
224def moment_to_magnitude(moment):
225 '''
226 Convert scalar moment to moment magnitude Mw.
228 :param moment: scalar moment [Nm]
229 :returns: moment magnitude Mw
231 Moment magnitude is defined as
233 .. math::
235 M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0 \\cdot 10^7) - 10.7
237 where :math:`M_0` is the scalar moment given in [Nm].
239 .. note::
241 Global CMT uses 10.7333333 instead of 10.7, based on [Kanamori 1977],
242 10.7 is from [Hanks and Kanamori 1979].
243 '''
245 return num.log10(moment*1.0e7) / 1.5 - 10.7
248def magnitude_to_moment(magnitude):
249 '''
250 Convert moment magnitude Mw to scalar moment.
252 :param magnitude: moment magnitude
253 :returns: scalar moment [Nm]
255 See :py:func:`moment_to_magnitude`.
256 '''
258 return 10.0**(1.5*(magnitude+10.7))*1.0e-7
261magnitude_1Nm = moment_to_magnitude(1.0)
264def euler_to_matrix(alpha, beta, gamma):
265 '''
266 Given euler angle triplet, create rotation matrix
268 Given coordinate system `(x,y,z)` and rotated system `(xs,ys,zs)`
269 the line of nodes is the intersection between the `x,y` and the `xs,ys`
270 planes.
272 :param alpha: is the angle between the `z`-axis and the `zs`-axis [rad]
273 :param beta: is the angle between the `x`-axis and the line of nodes [rad]
274 :param gamma: is the angle between the line of nodes and the `xs`-axis
275 [rad]
277 Usage for moment tensors::
279 m_unrot = numpy.array([[0,0,-1],[0,0,0],[-1,0,0]])
280 euler_to_matrix(dip,strike,-rake, rotmat)
281 m = num.dot(rotmat.T, num.dot(m_unrot, rotmat))
283 '''
285 ca = math.cos(alpha)
286 cb = math.cos(beta)
287 cg = math.cos(gamma)
288 sa = math.sin(alpha)
289 sb = math.sin(beta)
290 sg = math.sin(gamma)
292 mat = num.array([[cb*cg-ca*sb*sg, sb*cg+ca*cb*sg, sa*sg],
293 [-cb*sg-ca*sb*cg, -sb*sg+ca*cb*cg, sa*cg],
294 [sa*sb, -sa*cb, ca]], dtype=float)
295 return mat
298def matrix_to_euler(rotmat):
299 '''
300 Get Eulerian angle triplet from rotation matrix.
301 '''
303 ex = cvec(1., 0., 0.)
304 ez = cvec(0., 0., 1.)
305 exs = num.dot(rotmat.T, ex)
306 ezs = num.dot(rotmat.T, ez)
307 enodes = num.cross(ez.T, ezs.T).T
308 if num.linalg.norm(enodes) < 1e-10:
309 enodes = exs
310 enodess = num.dot(rotmat, enodes)
311 cos_alpha = float(num.dot(ez.T, ezs)[0, 0])
312 if cos_alpha > 1.:
313 cos_alpha = 1.
315 if cos_alpha < -1.:
316 cos_alpha = -1.
318 alpha = math.acos(cos_alpha)
319 beta = num.mod(math.atan2(enodes[1, 0], enodes[0, 0]), math.pi*2.)
320 gamma = num.mod(-math.atan2(enodess[1, 0], enodess[0, 0]), math.pi*2.)
322 return unique_euler(alpha, beta, gamma)
325def unique_euler(alpha, beta, gamma):
326 '''
327 Uniquify Eulerian angle triplet.
329 Put Eulerian angle triplet into ranges compatible with
330 ``(dip, strike, -rake)`` conventions in seismology::
332 alpha (dip) : [0, pi/2]
333 beta (strike) : [0, 2*pi)
334 gamma (-rake) : [-pi, pi)
336 If ``alpha1`` is near to zero, ``beta`` is replaced by ``beta+gamma`` and
337 ``gamma`` is set to zero, to prevent this additional ambiguity.
339 If ``alpha`` is near to ``pi/2``, ``beta`` is put into the range
340 ``[0,pi)``.
341 '''
343 pi = math.pi
345 alpha = num.mod(alpha, 2.0*pi)
347 if 0.5*pi < alpha and alpha <= pi:
348 alpha = pi - alpha
349 beta = beta + pi
350 gamma = 2.0*pi - gamma
351 elif pi < alpha and alpha <= 1.5*pi:
352 alpha = alpha - pi
353 gamma = pi - gamma
354 elif 1.5*pi < alpha and alpha <= 2.0*pi:
355 alpha = 2.0*pi - alpha
356 beta = beta + pi
357 gamma = pi + gamma
359 alpha = num.mod(alpha, 2.0*pi)
360 beta = num.mod(beta, 2.0*pi)
361 gamma = num.mod(gamma+pi, 2.0*pi)-pi
363 # If dip is exactly 90 degrees, one is still
364 # free to choose between looking at the plane from either side.
365 # Choose to look at such that beta is in the range [0,180)
367 # This should prevent some problems, when dip is close to 90 degrees:
368 if abs(alpha - 0.5*pi) < 1e-10:
369 alpha = 0.5*pi
371 if abs(beta - pi) < 1e-10:
372 beta = pi
374 if abs(beta - 2.*pi) < 1e-10:
375 beta = 0.
377 if abs(beta) < 1e-10:
378 beta = 0.
380 if alpha == 0.5*pi and beta >= pi:
381 gamma = - gamma
382 beta = num.mod(beta-pi, 2.0*pi)
383 gamma = num.mod(gamma+pi, 2.0*pi)-pi
384 assert 0. <= beta < pi
385 assert -pi <= gamma < pi
387 if alpha < 1e-7:
388 beta = num.mod(beta + gamma, 2.0*pi)
389 gamma = 0.
391 return (alpha, beta, gamma)
394def cvec(x, y, z):
395 return num.array([[x, y, z]], dtype=float).T
398def rvec(x, y, z):
399 return num.array([[x, y, z]], dtype=float)
402def eigh_check(a):
403 evals, evecs = num.linalg.eigh(a)
404 assert evals[0] <= evals[1] <= evals[2]
405 return evals, evecs
408r2d = 180. / math.pi
409d2r = 1. / r2d
412def random_mt(x=None, scalar_moment=1.0, magnitude=None):
414 if magnitude is not None:
415 scalar_moment = magnitude_to_moment(magnitude)
417 if x is None:
418 x = num.random.random(6)
420 evals = x[:3] * 2. - 1.0
421 evals /= num.sqrt(num.sum(evals**2)) / math.sqrt(2.0)
422 rotmat = random_rotation(x[3:])
423 return scalar_moment * num.dot(
424 rotmat, num.dot(num.array(num.diag(evals)), rotmat.T))
427def random_m6(*args, **kwargs):
428 return to6(random_mt(*args, **kwargs))
431def random_dc(x=None, scalar_moment=1.0, magnitude=None):
432 if magnitude is not None:
433 scalar_moment = magnitude_to_moment(magnitude)
435 rotmat = random_rotation(x)
436 return scalar_moment * num.dot(
437 rotmat, num.dot(MomentTensor._m_unrot, rotmat.T))
440def sm(m):
441 return '/ %5.2F %5.2F %5.2F \\\n' % (m[0, 0], m[0, 1], m[0, 2]) + \
442 '| %5.2F %5.2F %5.2F |\n' % (m[1, 0], m[1, 1], m[1, 2]) + \
443 '\\ %5.2F %5.2F %5.2F /\n' % (m[2, 0], m[2, 1], m[2, 2])
446def as_mt(mt):
447 '''
448 Convenience function to convert various objects to moment tensor object.
450 Like :py:meth:`MomentTensor.from_values`, but does not create a new
451 :py:class:`MomentTensor` object when ``mt`` already is one.
452 '''
454 if isinstance(mt, MomentTensor):
455 return mt
456 else:
457 return MomentTensor.from_values(mt)
460def pt_axes_to_strike_dip_rake(p_axis, t_axis):
461 n_axis = num.cross(p_axis, t_axis)
462 m_evecs = num.vstack([p_axis, n_axis, t_axis]).T
463 if num.linalg.det(m_evecs) < 0.:
464 m_evecs *= -1.
465 rotmat = num.dot(m_evecs, MomentTensor._u_evecs.T).T
466 if num.linalg.det(rotmat) < 0.:
467 rotmat *= -1.
468 alpha, beta, gamma = [r2d*x for x in matrix_to_euler(rotmat)]
469 return (beta, alpha, -gamma)
472def ned_to_rta(ned):
473 '''
474 Convert north-east-down coordinate vectors to radial-takeoff-azimuth.
476 Output coordinate system has coordinates radial, takeoff angle [deg]
477 (downward is zero), and azimuth [deg] (northward is zero).
479 :param ned:
480 Coordinate vector or array of coordinate vectors (north, east, down).
481 :type ned:
482 :py:class:`numpy.ndarray` of shape ``(N, 3)`` or ``(3,)``
484 :returns:
485 Coordinate vector or array of coordinate vectors (radial, takeoff,
486 azimuth)
487 :rtype:
488 :py:class:`numpy.ndarray` of shape ``(N, 3)`` or ``(3,)``
489 '''
490 if ned.ndim == 1:
491 return ned_to_rta(ned[num.newaxis, :])[0]
493 rta = num.empty_like(ned)
494 rta[:, 0] = num.sqrt(num.sum(ned**2, axis=1))
495 rta[:, 1] = num.arccos(ned[:, 2]) * r2d
496 rta[:, 2] = num.arctan2(ned[:, 1], ned[:, 0]) * r2d
497 return rta
500class MomentTensor(Object):
502 '''
503 Moment tensor object
505 :param m: NumPy array in north-east-down convention
506 :param m_up_south_east: NumPy array in up-south-east convention
507 :param m_east_north_up: NumPy array in east-north-up convention
508 :param strike,dip,rake: fault plane angles in [degrees]
509 :param p_axis,t_axis: initialize double-couple from p and t axes
510 :param scalar_moment: scalar moment in [Nm]
511 :param magnitude: moment magnitude Mw
513 Global CMT catalog moment tensors use the up-south-east (USE) coordinate
514 system convention with :math:`r` (up), :math:`\\theta` (south), and
515 :math:`\\phi` (east).
517 .. math::
518 :nowrap:
520 \\begin{align*}
521 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
522 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
523 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
524 \\end{align*}
526 '''
528 mnn__ = Float.T(default=0.0)
529 mee__ = Float.T(default=0.0)
530 mdd__ = Float.T(default=0.0)
531 mne__ = Float.T(default=0.0)
532 mnd__ = Float.T(default=-1.0)
533 med__ = Float.T(default=0.0)
534 strike1__ = Float.T(default=None, optional=True) # read-only
535 dip1__ = Float.T(default=None, optional=True) # read-only
536 rake1__ = Float.T(default=None, optional=True) # read-only
537 strike2__ = Float.T(default=None, optional=True) # read-only
538 dip2__ = Float.T(default=None, optional=True) # read-only
539 rake2__ = Float.T(default=None, optional=True) # read-only
540 moment__ = Float.T(default=None, optional=True) # read-only
541 magnitude__ = Float.T(default=None, optional=True) # read-only
543 _flip_dc = num.array(
544 [[0., 0., -1.], [0., -1., 0.], [-1., 0., 0.]], dtype=float)
545 _to_up_south_east = num.array(
546 [[0., 0., -1.], [-1., 0., 0.], [0., 1., 0.]], dtype=float).T
547 _to_east_north_up = num.array(
548 [[0., 1., 0.], [1., 0., 0.], [0., 0., -1.]], dtype=float)
549 _m_unrot = num.array(
550 [[0., 0., -1.], [0., 0., 0.], [-1., 0., 0.]], dtype=float)
552 _u_evals, _u_evecs = eigh_check(_m_unrot)
554 @classmethod
555 def random_dc(cls, x=None, scalar_moment=1.0, magnitude=None):
556 '''
557 Create random oriented double-couple moment tensor
559 The rotations used are uniformly distributed in the space of rotations.
560 '''
561 return MomentTensor(
562 m=random_dc(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
564 @classmethod
565 def random_mt(cls, x=None, scalar_moment=1.0, magnitude=None):
566 '''
567 Create random moment tensor
569 Moment tensors produced by this function appear uniformly distributed
570 when shown in a Hudson's diagram. The rotations used are uniformly
571 distributed in the space of rotations.
572 '''
573 return MomentTensor(
574 m=random_mt(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
576 @classmethod
577 def from_values(cls, values):
578 '''
579 Alternative constructor for moment tensor objects
581 This constructor takes a :py:class:`MomentTensor` object, a tuple, list
582 or NumPy array with 3x3 or 3, 4, 6, or 7 elements to build a Moment
583 tensor object.
585 The ``values`` argument is interpreted depending on shape and type as
586 follows:
588 * ``(strike, dip, rake)``
589 * ``(strike, dip, rake, magnitude)``
590 * ``(mnn, mee, mdd, mne, mnd, med)``
591 * ``(mnn, mee, mdd, mne, mnd, med, magnitude)``
592 * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))``
593 * :py:class:`MomentTensor` object
594 '''
596 m = values_to_matrix(values)
597 return MomentTensor(m=m)
599 def __init__(
600 self, m=None, m_up_south_east=None, m_east_north_up=None,
601 strike=0., dip=0., rake=0., scalar_moment=1.,
602 mnn=None, mee=None, mdd=None, mne=None, mnd=None, med=None,
603 strike1=None, dip1=None, rake1=None,
604 strike2=None, dip2=None, rake2=None,
605 p_axis=None, t_axis=None,
606 magnitude=None, moment=None):
608 Object.__init__(self, init_props=False)
610 if any(mxx is not None for mxx in (mnn, mee, mdd, mne, mnd, med)):
611 m = symmat6(mnn, mee, mdd, mne, mnd, med)
613 if p_axis is not None and t_axis is not None:
614 strike, dip, rake = pt_axes_to_strike_dip_rake(p_axis, t_axis)
616 strike = d2r*strike
617 dip = d2r*dip
618 rake = d2r*rake
620 if isinstance(m, num.matrix):
621 m = m.A
623 if isinstance(m_up_south_east, num.matrix):
624 m_up_south_east = m_up_south_east.A
626 if isinstance(m_east_north_up, num.matrix):
627 m_east_north_up = m_east_north_up.A
629 if m_up_south_east is not None:
630 m = num.dot(
631 self._to_up_south_east,
632 num.dot(m_up_south_east, self._to_up_south_east.T))
634 if m_east_north_up is not None:
635 m = num.dot(
636 self._to_east_north_up,
637 num.dot(m_east_north_up, self._to_east_north_up.T))
639 if m is None:
640 if any(x is not None for x in (
641 strike1, dip1, rake1, strike2, dip2, rake2)):
642 raise Exception(
643 'strike1, dip1, rake1, strike2, dip2, rake2 are read-only '
644 'properties')
646 if moment is not None:
647 scalar_moment = moment
649 if magnitude is not None:
650 scalar_moment = magnitude_to_moment(magnitude)
652 rotmat1 = euler_to_matrix(dip, strike, -rake)
653 m = scalar_moment * num.dot(
654 rotmat1.T,
655 num.dot(MomentTensor._m_unrot, rotmat1))
657 self._m = m
658 self._update()
660 def _update(self):
661 m_evals, m_evecs = eigh_check(self._m)
662 if num.linalg.det(m_evecs) < 0.:
663 m_evecs *= -1.
665 rotmat1 = num.dot(m_evecs, MomentTensor._u_evecs.T).T
666 if num.linalg.det(rotmat1) < 0.:
667 rotmat1 *= -1.
669 self._m_eigenvals = m_evals
670 self._m_eigenvecs = m_evecs
672 self._rotmats = sorted(
673 [rotmat1, num.dot(MomentTensor._flip_dc, rotmat1)],
674 key=lambda m: num.abs(m.flat).tolist())
676 @property
677 def mnn(self):
678 return float(self._m[0, 0])
680 @mnn.setter
681 def mnn(self, value):
682 self._m[0, 0] = value
683 self._update()
685 @property
686 def mee(self):
687 return float(self._m[1, 1])
689 @mee.setter
690 def mee(self, value):
691 self._m[1, 1] = value
692 self._update()
694 @property
695 def mdd(self):
696 return float(self._m[2, 2])
698 @mdd.setter
699 def mdd(self, value):
700 self._m[2, 2] = value
701 self._update()
703 @property
704 def mne(self):
705 return float(self._m[0, 1])
707 @mne.setter
708 def mne(self, value):
709 self._m[0, 1] = value
710 self._m[1, 0] = value
711 self._update()
713 @property
714 def mnd(self):
715 return float(self._m[0, 2])
717 @mnd.setter
718 def mnd(self, value):
719 self._m[0, 2] = value
720 self._m[2, 0] = value
721 self._update()
723 @property
724 def med(self):
725 return float(self._m[1, 2])
727 @med.setter
728 def med(self, value):
729 self._m[1, 2] = value
730 self._m[2, 1] = value
731 self._update()
733 @property
734 def strike1(self):
735 return float(self.both_strike_dip_rake()[0][0])
737 @property
738 def dip1(self):
739 return float(self.both_strike_dip_rake()[0][1])
741 @property
742 def rake1(self):
743 return float(self.both_strike_dip_rake()[0][2])
745 @property
746 def strike2(self):
747 return float(self.both_strike_dip_rake()[1][0])
749 @property
750 def dip2(self):
751 return float(self.both_strike_dip_rake()[1][1])
753 @property
754 def rake2(self):
755 return float(self.both_strike_dip_rake()[1][2])
757 def both_strike_dip_rake(self):
758 '''
759 Get both possible (strike,dip,rake) triplets.
760 '''
761 results = []
762 for rotmat in self._rotmats:
763 alpha, beta, gamma = [r2d*x for x in matrix_to_euler(rotmat)]
764 results.append((beta, alpha, -gamma))
766 return results
768 def p_axis(self):
769 '''
770 Get direction of p axis.
772 Use :py:func:`ned_to_rta` to get takeoff angle and azimuth of the
773 returned vectors.
775 :returns:
776 Direction of P (pressure) axis as a (north, east, down) vector.
777 :rtype:
778 :py:class:`numpy.ndarray` of shape ``(3,)``
779 '''
780 return (self._m_eigenvecs.T)[0]
782 def t_axis(self):
783 '''
784 Get direction of t axis.
786 Use :py:func:`ned_to_rta` to get takeoff angle and azimuth of the
787 returned vectors.
789 :returns:
790 Direction of T (tension) axis as a (north, east, down) vector.
791 :rtype:
792 :py:class:`numpy.ndarray` of shape ``(3,)``
793 '''
794 return (self._m_eigenvecs.T)[2]
796 def null_axis(self):
797 '''
798 Get direction of the null axis.
800 Use :py:func:`ned_to_rta` to get takeoff angle and azimuth of the
801 returned vectors.
803 :returns:
804 Direction of null (B) axis as a (north, east, down) vector.
805 :rtype:
806 :py:class:`numpy.ndarray` of shape ``(3,)``
807 '''
808 return self._m_eigenvecs.T[1]
810 def eigenvals(self):
811 '''
812 Get the eigenvalues of the moment tensor in ascending order.
814 :returns: ``(ep, en, et)``
815 '''
817 return self._m_eigenvals
819 def eigensystem(self):
820 '''
821 Get the eigenvalues and eigenvectors of the moment tensor.
823 :returns: ``(ep, en, et, vp, vn, vt)``
824 '''
826 vp = self.p_axis().flatten()
827 vn = self.null_axis().flatten()
828 vt = self.t_axis().flatten()
829 ep, en, et = self._m_eigenvals
830 return ep, en, et, vp, vn, vt
832 def both_slip_vectors(self):
833 '''
834 Get both possible slip directions.
835 '''
836 return [
837 num.dot(rotmat.T, cvec(1., 0., 0.)) for rotmat in self._rotmats]
839 def m(self):
840 '''
841 Get plain moment tensor as 3x3 matrix.
842 '''
843 return self._m.copy()
845 def m6(self):
846 '''
847 Get the moment tensor as a six-element array.
849 :returns: ``(mnn, mee, mdd, mne, mnd, med)``
850 '''
851 return to6(self._m)
853 def m_up_south_east(self):
854 '''
855 Get moment tensor in up-south-east convention as 3x3 matrix.
857 .. math::
858 :nowrap:
860 \\begin{align*}
861 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
862 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
863 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
864 \\end{align*}
865 '''
867 return num.dot(
868 self._to_up_south_east.T,
869 num.dot(self._m, self._to_up_south_east))
871 def m_east_north_up(self):
872 '''
873 Get moment tensor in east-north-up convention as 3x3 matrix.
874 '''
876 return num.dot(
877 self._to_east_north_up.T, num.dot(self._m, self._to_east_north_up))
879 def m6_up_south_east(self):
880 '''
881 Get moment tensor in up-south-east convention as a six-element array.
883 :returns: ``(muu, mss, mee, mus, mue, mse)``
884 '''
885 return to6(self.m_up_south_east())
887 def m6_east_north_up(self):
888 '''
889 Get moment tensor in east-north-up convention as a six-element array.
891 :returns: ``(mee, mnn, muu, men, meu, mnu)``
892 '''
893 return to6(self.m_east_north_up())
895 def m_plain_double_couple(self):
896 '''
897 Get plain double couple with same scalar moment as moment tensor.
898 '''
899 rotmat1 = self._rotmats[0]
900 m = self.scalar_moment() * num.dot(
901 rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1))
902 return m
904 def moment_magnitude(self):
905 '''
906 Get moment magnitude of moment tensor.
907 '''
908 return moment_to_magnitude(self.scalar_moment())
910 def scalar_moment(self):
911 '''
912 Get the scalar moment of the moment tensor (Frobenius norm based)
914 .. math::
916 M_0 = \\frac{1}{\\sqrt{2}}\\sqrt{\\sum_{i,j} |M_{ij}|^2}
918 The scalar moment is calculated based on the Euclidean (Frobenius) norm
919 (Silver and Jordan, 1982). The scalar moment returned by this function
920 differs from the standard decomposition based definition of the scalar
921 moment for non-double-couple moment tensors.
922 '''
923 return num.linalg.norm(self._m_eigenvals) / math.sqrt(2.)
925 @property
926 def moment(self):
927 return float(self.scalar_moment())
929 @moment.setter
930 def moment(self, value):
931 self._m *= value / self.moment
932 self._update()
934 @property
935 def magnitude(self):
936 return float(self.moment_magnitude())
938 @magnitude.setter
939 def magnitude(self, value):
940 self._m *= magnitude_to_moment(value) / self.moment
941 self._update()
943 def __str__(self):
944 mexp = pow(10, math.ceil(num.log10(num.max(num.abs(self._m)))))
945 m = self._m / mexp
946 s = '''Scalar Moment [Nm]: M0 = %g (Mw = %3.1f)
947Moment Tensor [Nm]: Mnn = %6.3f, Mee = %6.3f, Mdd = %6.3f,
948 Mne = %6.3f, Mnd = %6.3f, Med = %6.3f [ x %g ]
949''' % (
950 self.scalar_moment(),
951 self.moment_magnitude(),
952 m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2],
953 mexp)
955 s += self.str_fault_planes()
956 return s
958 def str_fault_planes(self):
959 s = ''
960 for i, sdr in enumerate(self.both_strike_dip_rake()):
961 s += 'Fault plane %i [deg]: ' \
962 'strike = %3.0f, dip = %3.0f, slip-rake = %4.0f\n' \
963 % (i+1, sdr[0], sdr[1], sdr[2])
965 return s
967 def deviatoric(self):
968 '''
969 Get deviatoric part of moment tensor.
971 Returns a new moment tensor object with zero trace.
972 '''
974 m = self.m()
976 trace_m = num.trace(m)
977 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
978 m_devi = m - m_iso
979 mt = MomentTensor(m=m_devi)
980 return mt
982 def standard_decomposition(self):
983 '''
984 Decompose moment tensor into isotropic, DC and CLVD components.
986 Standard decomposition according to e.g. Jost and Herrmann 1989 is
987 returned as::
989 [
990 (moment_iso, ratio_iso, m_iso),
991 (moment_dc, ratio_dc, m_dc),
992 (moment_clvd, ratio_clvd, m_clvd),
993 (moment_devi, ratio_devi, m_devi),
994 (moment, 1.0, m)
995 ]
996 '''
998 epsilon = 1e-6
1000 m = self.m()
1002 trace_m = num.trace(m)
1003 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
1004 moment_iso = abs(trace_m / 3.)
1006 m_devi = m - m_iso
1008 evals, evecs = eigh_check(m_devi)
1010 moment_devi = num.max(num.abs(evals))
1011 moment = moment_iso + moment_devi
1013 iorder = num.argsort(num.abs(evals))
1014 evals_sorted = evals[iorder]
1015 evecs_sorted = (evecs.T[iorder]).T
1017 if moment_devi < epsilon * moment_iso:
1018 signed_moment_dc = 0.
1019 else:
1020 assert -epsilon <= -evals_sorted[0] / evals_sorted[2] <= 0.5
1021 signed_moment_dc = evals_sorted[2] * (1.0 + 2.0 * (
1022 min(0.0, evals_sorted[0] / evals_sorted[2])))
1024 moment_dc = abs(signed_moment_dc)
1025 m_dc_es = num.dot(signed_moment_dc, num.diag([0., -1.0, 1.0]))
1026 m_dc = num.dot(evecs_sorted, num.dot(m_dc_es, evecs_sorted.T))
1028 m_clvd = m_devi - m_dc
1030 moment_clvd = moment_devi - moment_dc
1032 ratio_dc = moment_dc / moment
1033 ratio_clvd = moment_clvd / moment
1034 ratio_iso = moment_iso / moment
1035 ratio_devi = moment_devi / moment
1037 return [
1038 (moment_iso, ratio_iso, m_iso),
1039 (moment_dc, ratio_dc, m_dc),
1040 (moment_clvd, ratio_clvd, m_clvd),
1041 (moment_devi, ratio_devi, m_devi),
1042 (moment, 1.0, m)]
1044 def rotated(self, rot):
1045 '''
1046 Get rotated moment tensor.
1048 :param rot: ratation matrix, coordinate system is NED
1049 :returns: new :py:class:`MomentTensor` object
1050 '''
1052 rotmat = num.array(rot)
1053 return MomentTensor(m=num.dot(rotmat, num.dot(self.m(), rotmat.T)))
1055 def random_rotated(self, angle_std=None, angle=None, rstate=None):
1056 '''
1057 Get distorted MT by rotation around random axis and angle.
1059 :param angle_std: angles are drawn from a normal distribution with
1060 zero mean and given standard deviation [degrees]
1061 :param angle: set angle [degrees], only axis will be random
1062 :param rstate: :py:class:`numpy.random.RandomState` object, can be
1063 used to create reproducible pseudo-random sequences
1064 :returns: new :py:class:`MomentTensor` object
1065 '''
1067 assert (angle_std is None) != (angle is None), \
1068 'either angle or angle_std must be given'
1070 if angle_std is not None:
1071 rstate = rstate or num.random
1072 angle = rstate.normal(scale=angle_std)
1074 axis = random_axis(rstate=rstate)
1075 rot = rotation_from_angle_and_axis(angle, axis)
1076 return self.rotated(rot)
1079def other_plane(strike, dip, rake):
1080 '''
1081 Get the respectively other plane in the double-couple ambiguity.
1082 '''
1084 mt = MomentTensor(strike=strike, dip=dip, rake=rake)
1085 both_sdr = mt.both_strike_dip_rake()
1086 w = [sum([abs(x-y) for x, y in zip(both_sdr[i], (strike, dip, rake))])
1087 for i in (0, 1)]
1089 if w[0] < w[1]:
1090 return both_sdr[1]
1091 else:
1092 return both_sdr[0]
1095def dsdr(sdr1, sdr2):
1096 s1, d1, r1 = sdr1
1097 s2, d2, r2 = sdr2
1099 s1 = s1 % 360.
1100 s2 = s2 % 360.
1101 r1 = r1 % 360.
1102 r2 = r2 % 360.
1104 ds = abs(s1 - s2)
1105 ds = ds if ds <= 180. else 360. - ds
1107 dr = abs(r1 - r2)
1108 dr = dr if dr <= 180. else 360. - dr
1110 dd = abs(d1 - d2)
1112 return math.sqrt(ds**2 + dr**2 + dd**2)
1115def order_like(sdrs, sdrs_ref):
1116 '''
1117 Order strike-dip-rake pair post closely to a given reference pair.
1119 :param sdrs: tuple, ``((strike1, dip1, rake1), (strike2, dip2, rake2))``
1120 :param sdrs_ref: as above but with reference pair
1121 '''
1123 d1 = min(dsdr(sdrs[0], sdrs_ref[0]), dsdr(sdrs[1], sdrs_ref[1]))
1124 d2 = min(dsdr(sdrs[0], sdrs_ref[1]), dsdr(sdrs[1], sdrs_ref[0]))
1125 if d1 < d2:
1126 return sdrs
1127 else:
1128 return sdrs[::-1]
1131def _tpb2q(t, p, b):
1132 eps = 0.001
1133 tqw = 1. + t[0] + p[1] + b[2]
1134 tqx = 1. + t[0] - p[1] - b[2]
1135 tqy = 1. - t[0] + p[1] - b[2]
1136 tqz = 1. - t[0] - p[1] + b[2]
1138 q = num.zeros(4)
1139 if tqw > eps:
1140 q[0] = 0.5 * math.sqrt(tqw)
1141 q[1:] = p[2] - b[1], b[0] - t[2], t[1] - p[0]
1142 elif tqx > eps:
1143 q[0] = 0.5 * math.sqrt(tqx)
1144 q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1145 elif tqy > eps:
1146 q[0] = 0.5 * math.sqrt(tqy)
1147 q[1:] = b[0] - t[2], p[0] + t[1], b[1] + p[2]
1148 elif tqz > eps:
1149 q[0] = 0.5 * math.sqrt(tqz)
1150 q[1:] = t[1] - p[0], b[0] + t[2], b[1] + p[2]
1151 else:
1152 raise Exception('should not reach this line, check theory!')
1153 # q[0] = max(0.5 * math.sqrt(tqx), eps)
1154 # q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1156 q[1:] /= 4.0 * q[0]
1158 q /= math.sqrt(num.sum(q**2))
1160 return q
1163_pbt2tpb = num.array(((0., 0., 1.), (1., 0., 0.), (0., 1., 0.)))
1166def kagan_angle(mt1, mt2):
1167 '''
1168 Given two moment tensors, return the Kagan angle in degrees.
1170 After Kagan (1991) and Tape & Tape (2012).
1171 '''
1173 ai = num.dot(_pbt2tpb, mt1._m_eigenvecs.T)
1174 aj = num.dot(_pbt2tpb, mt2._m_eigenvecs.T)
1176 u = num.dot(ai, aj.T)
1178 tk, pk, bk = u
1180 qk = _tpb2q(tk, pk, bk)
1182 return 2. * r2d * math.acos(num.max(num.abs(qk)))
1185def rand_to_gutenberg_richter(rand, b_value, magnitude_min):
1186 '''
1187 Draw magnitude from Gutenberg Richter distribution.
1188 '''
1189 return magnitude_min + num.log10(1.-rand) / -b_value
1192def magnitude_to_duration_gcmt(magnitudes):
1193 '''
1194 Scaling relation used by Global CMT catalog for most of its events.
1195 '''
1197 mom = magnitude_to_moment(magnitudes)
1198 return (mom / 1.1e16)**(1./3.)
1201if __name__ == '__main__':
1203 import sys
1204 v = list(map(float, sys.argv[1:]))
1205 mt = MomentTensor.from_values(v)
1206 print(mt)