Source code for obspy.core.event.source

# -*- coding: utf-8 -*-
"""
obspy.core.event.source - Classes for handling seismic source characteristics
=============================================================================
This module provides a class hierarchy to consistently handle event metadata.
This class hierarchy is closely modelled after the de-facto standard format
`QuakeML <https://quake.ethz.ch/quakeml/>`_.

.. figure:: /_images/Event.png

.. note::

    For handling additional information not covered by the QuakeML standard and
    how to output it to QuakeML see the :ref:`ObsPy Tutorial <quakeml-extra>`.

:copyright:
    The ObsPy Development Team (devs@obspy.org)
:license:
    GNU Lesser General Public License, Version 3
    (https://www.gnu.org/copyleft/lesser.html)
"""
from __future__ import (absolute_import, division, print_function,
                        unicode_literals)
from future.builtins import *  # NOQA

import numpy as np

from obspy.core.event.base import (
    _event_type_class_factory, ResourceIdentifier, CreationInfo)
from obspy.core.event.header import (
    EvaluationMode, EvaluationStatus, MomentTensorCategory, MTInversionType,
    SourceTimeFunctionType, ATTRIBUTE_HAS_ERRORS)


__Axis = _event_type_class_factory(
    "__Axis",
    class_attributes=[("azimuth", float, ATTRIBUTE_HAS_ERRORS),
                      ("plunge", float, ATTRIBUTE_HAS_ERRORS),
                      ("length", float, ATTRIBUTE_HAS_ERRORS)])


[docs]class Axis(__Axis): """ This class describes an eigenvector of a moment tensor expressed in its principal-axes system. It uses the angles azimuth, plunge, and the eigenvalue length. :type azimuth: float :param azimuth: Azimuth of eigenvector of moment tensor expressed in principal-axes system. Measured clockwise from South-North direction at epicenter. Unit: deg :type azimuth_errors: :class:`~obspy.core.event.base.QuantityError` :param azimuth_errors: AttribDict containing error quantities. :type plunge: float :param plunge: Plunge of eigenvector of moment tensor expressed in principal-axes system. Measured against downward vertical direction at epicenter. Unit: deg :type plunge_errors: :class:`~obspy.core.event.base.QuantityError` :param plunge_errors: AttribDict containing error quantities. :type length: float :param length: Eigenvalue of moment tensor expressed in principal-axes system. Unit: Nm :type length_errors: :class:`~obspy.core.event.base.QuantityError` :param length_errors: AttribDict containing error quantities. .. note:: For handling additional information not covered by the QuakeML standard and how to output it to QuakeML see the :ref:`ObsPy Tutorial <quakeml-extra>`. """
__NodalPlane = _event_type_class_factory( "__NodalPlane", class_attributes=[("strike", float, ATTRIBUTE_HAS_ERRORS), ("dip", float, ATTRIBUTE_HAS_ERRORS), ("rake", float, ATTRIBUTE_HAS_ERRORS)])
[docs]class NodalPlane(__NodalPlane): """ This class describes a nodal plane using the attributes strike, dip, and rake. For a definition of the angles see Aki & Richards (1980). :type strike: float :param strike: Strike angle of nodal plane. Unit: deg :type strike_errors: :class:`~obspy.core.event.base.QuantityError` :param strike_errors: AttribDict containing error quantities. :type dip: float :param dip: Dip angle of nodal plane. Unit: deg :type dip_errors: :class:`~obspy.core.event.base.QuantityError` :param dip_errors: AttribDict containing error quantities. :type rake: float :param rake: Rake angle of nodal plane. Unit: deg :type rake_errors: :class:`~obspy.core.event.base.QuantityError` :param rake_errors: AttribDict containing error quantities. .. note:: For handling additional information not covered by the QuakeML standard and how to output it to QuakeML see the :ref:`ObsPy Tutorial <quakeml-extra>`. """
__NodalPlanes = _event_type_class_factory( "__NodalPlanes", class_attributes=[("nodal_plane_1", NodalPlane), ("nodal_plane_2", NodalPlane), ("preferred_plane", int)])
[docs]class NodalPlanes(__NodalPlanes): """ This class describes the nodal planes of a double-couple moment-tensor solution. The attribute ``preferred_plane`` can be used to define which plane is the preferred one. :type nodal_plane_1: :class:`~obspy.core.event.source.NodalPlane`, optional :param nodal_plane_1: First nodal plane of double-couple moment tensor solution. :type nodal_plane_2: :class:`~obspy.core.event.source.NodalPlane`, optional :param nodal_plane_2: Second nodal plane of double-couple moment tensor solution. :type preferred_plane: int, optional :param preferred_plane: Indicator for preferred nodal plane of moment tensor solution. It can take integer values ``1`` or ``2``. .. note:: For handling additional information not covered by the QuakeML standard and how to output it to QuakeML see the :ref:`ObsPy Tutorial <quakeml-extra>`. """
__PrincipalAxes = _event_type_class_factory( "__PrincipalAxes", class_attributes=[("t_axis", Axis), ("p_axis", Axis), ("n_axis", Axis)])
[docs]class PrincipalAxes(__PrincipalAxes): """ This class describes the principal axes of a double-couple moment tensor solution. t_axis and p_axis are required, while n_axis is optional. :type t_axis: :class:`~obspy.core.event.source.Axis` :param t_axis: T (tension) axis of a double-couple moment tensor solution. :type p_axis: :class:`~obspy.core.event.source.Axis` :param p_axis: P (pressure) axis of a double-couple moment tensor solution. :type n_axis: :class:`~obspy.core.event.source.Axis`, optional :param n_axis: N (neutral) axis of a double-couple moment tensor solution. .. note:: For handling additional information not covered by the QuakeML standard and how to output it to QuakeML see the :ref:`ObsPy Tutorial <quakeml-extra>`. """
__Tensor = _event_type_class_factory( "__Tensor", class_attributes=[("m_rr", float, ATTRIBUTE_HAS_ERRORS), ("m_tt", float, ATTRIBUTE_HAS_ERRORS), ("m_pp", float, ATTRIBUTE_HAS_ERRORS), ("m_rt", float, ATTRIBUTE_HAS_ERRORS), ("m_rp", float, ATTRIBUTE_HAS_ERRORS), ("m_tp", float, ATTRIBUTE_HAS_ERRORS)])
[docs]class Tensor(__Tensor): """ The Tensor class represents the six moment-tensor elements Mrr, Mtt, Mpp, Mrt, Mrp, Mtp in the spherical coordinate system defined by local upward vertical (r), North-South (t), and West-East (p) directions. :type m_rr: float :param m_rr: Moment-tensor element Mrr. Unit: Nm :type m_rr_errors: :class:`~obspy.core.event.base.QuantityError` :param m_rr_errors: AttribDict containing error quantities. :type m_tt: float :param m_tt: Moment-tensor element Mtt. Unit: Nm :type m_tt_errors: :class:`~obspy.core.event.base.QuantityError` :param m_tt_errors: AttribDict containing error quantities. :type m_pp: float :param m_pp: Moment-tensor element Mpp. Unit: Nm :type m_pp_errors: :class:`~obspy.core.event.base.QuantityError` :param m_pp_errors: AttribDict containing error quantities. :type m_rt: float :param m_rt: Moment-tensor element Mrt. Unit: Nm :type m_rt_errors: :class:`~obspy.core.event.base.QuantityError` :param m_rt_errors: AttribDict containing error quantities. :type m_rp: float :param m_rp: Moment-tensor element Mrp. Unit: Nm :type m_rp_errors: :class:`~obspy.core.event.base.QuantityError` :param m_rp_errors: AttribDict containing error quantities. :type m_tp: float :param m_tp: Moment-tensor element Mtp. Unit: Nm :type m_tp_errors: :class:`~obspy.core.event.base.QuantityError` :param m_tp_errors: AttribDict containing error quantities. .. note:: For handling additional information not covered by the QuakeML standard and how to output it to QuakeML see the :ref:`ObsPy Tutorial <quakeml-extra>`. """
__SourceTimeFunction = _event_type_class_factory( "__SourceTimeFunction", class_attributes=[("type", SourceTimeFunctionType), ("duration", float), ("rise_time", float), ("decay_time", float)])
[docs]class SourceTimeFunction(__SourceTimeFunction): """ Source time function used in moment-tensor inversion. :type type: str :param type: Type of source time function. See :class:`~obspy.core.event.header.SourceTimeFunctionType` for allowed values. :type duration: float :param duration: Source time function duration. Unit: s :type rise_time: float, optional :param rise_time: Source time function rise time. Unit: s :type decay_time: float, optional :param decay_time: Source time function decay time. Unit: s .. note:: For handling additional information not covered by the QuakeML standard and how to output it to QuakeML see the :ref:`ObsPy Tutorial <quakeml-extra>`. """
__MomentTensor = _event_type_class_factory( "__MomentTensor", class_attributes=[("resource_id", ResourceIdentifier), ("derived_origin_id", ResourceIdentifier), ("moment_magnitude_id", ResourceIdentifier), ("scalar_moment", float, ATTRIBUTE_HAS_ERRORS), ("tensor", Tensor), ("variance", float), ("variance_reduction", float), ("double_couple", float), ("clvd", float), ("iso", float), ("greens_function_id", ResourceIdentifier), ("filter_id", ResourceIdentifier), ("source_time_function", SourceTimeFunction), ("method_id", ResourceIdentifier), ("category", MomentTensorCategory), ("inversion_type", MTInversionType), ("creation_info", CreationInfo)], class_contains=["comments", "data_used"])
[docs]class MomentTensor(__MomentTensor): """ This class represents a moment tensor solution for an event. It is an optional part of a FocalMechanism description. :type resource_id: :class:`~obspy.core.event.base.ResourceIdentifier` :param resource_id: Resource identifier of MomentTensor. :type force_resource_id: bool, optional :param force_resource_id: If set to False, the automatic initialization of `resource_id` attribute in case it is not specified will be skipped. :type derived_origin_id: :class:`~obspy.core.event.base.ResourceIdentifier` :param derived_origin_id: Refers to the resource_id of the Origin derived in the moment tensor inversion. :type moment_magnitude_id: :class:`~obspy.core.event.base.ResourceIdentifier`, optional :param moment_magnitude_id: Refers to the publicID of the Magnitude object which represents the derived moment magnitude. :type scalar_moment: float, optional :param scalar_moment: Scalar moment as derived in moment tensor inversion. Unit: Nm :type scalar_moment_errors: :class:`~obspy.core.event.base.QuantityError` :param scalar_moment_errors: AttribDict containing error quantities. :type tensor: :class:`~obspy.core.event.source.Tensor`, optional :param tensor: Tensor object holding the moment tensor elements. :type variance: float, optional :param variance: Variance of moment tensor inversion. :type variance_reduction: float, optional :param variance_reduction: Variance reduction of moment tensor inversion, given in percent (Dreger 2003). This is a goodness-of-fit measure. :type double_couple: float, optional :param double_couple: Double couple parameter obtained from moment tensor inversion (decimal fraction between 0 and 1). :type clvd: float, optional :param clvd: CLVD (compensated linear vector dipole) parameter obtained from moment tensor inversion (decimal fraction between 0 and 1). :type iso: float, optional :param iso: Isotropic part obtained from moment tensor inversion (decimal fraction between 0 and 1). :type greens_function_id: :class:`~obspy.core.event.base.ResourceIdentifier`, optional :param greens_function_id: Resource identifier of the Green’s function used in moment tensor inversion. :type filter_id: :class:`~obspy.core.event.base.ResourceIdentifier`, optional :param filter_id: Resource identifier of the filter setup used in moment tensor inversion. :type source_time_function: :class:`~obspy.core.event.source.SourceTimeFunction`, optional :param source_time_function: Source time function used in moment-tensor inversion. :type data_used: list of :class:`~obspy.core.event.base.DataUsed`, optional :param data_used: Describes waveform data used for moment-tensor inversion. :type method_id: :class:`~obspy.core.event.base.ResourceIdentifier`, optional :param method_id: Resource identifier of the method used for moment-tensor inversion. :type category: str, optional :param category: Moment tensor category. See :class:`~obspy.core.event.header.MomentTensorCategory` for allowed values. :type inversion_type: str, optional :param inversion_type: Moment tensor inversion type. Users should avoid to give contradictory information in inversion_type and method_id. See :class:`~obspy.core.event.header.MTInversionType` for allowed values. :type comments: list of :class:`~obspy.core.event.base.Comment`, optional :param comments: Additional comments. :type creation_info: :class:`~obspy.core.event.base.CreationInfo`, optional :param creation_info: Creation information used to describe author, version, and creation time. .. note:: For handling additional information not covered by the QuakeML standard and how to output it to QuakeML see the :ref:`ObsPy Tutorial <quakeml-extra>`. """
__FocalMechanism = _event_type_class_factory( "__FocalMechanism", class_attributes=[("resource_id", ResourceIdentifier), ("triggering_origin_id", ResourceIdentifier), ("nodal_planes", NodalPlanes), ("principal_axes", PrincipalAxes), ("azimuthal_gap", float), ("station_polarity_count", int), ("misfit", float), ("station_distribution_ratio", float), ("method_id", ResourceIdentifier), ("evaluation_mode", EvaluationMode), ("evaluation_status", EvaluationStatus), ("moment_tensor", MomentTensor), ("creation_info", CreationInfo)], class_contains=['waveform_id', 'comments'])
[docs]class FocalMechanism(__FocalMechanism): """ This class describes the focal mechanism of an event. It includes different descriptions like nodal planes, principal axes, and a moment tensor. The moment tensor description is provided by objects of the class MomentTensor which can be specified as child elements of FocalMechanism. :type resource_id: :class:`~obspy.core.event.base.ResourceIdentifier` :param resource_id: Resource identifier of FocalMechanism. :type force_resource_id: bool, optional :param force_resource_id: If set to False, the automatic initialization of `resource_id` attribute in case it is not specified will be skipped. :type triggering_origin_id: :class:`~obspy.core.event.base.ResourceIdentifier`, optional :param triggering_origin_id: Refers to the resource_id of the triggering origin. :type nodal_planes: :class:`~obspy.core.event.source.NodalPlanes`, optional :param nodal_planes: Nodal planes of the focal mechanism. :type principal_axes: :class:`~obspy.core.event.source.PrincipalAxes`, optional :param principal_axes: Principal axes of the focal mechanism. :type azimuthal_gap: float, optional :param azimuthal_gap: Largest azimuthal gap in distribution of stations used for determination of focal mechanism. Unit: deg :type station_polarity_count: int, optional :param station_polarity_count: :type misfit: float, optional :param misfit: Fraction of misfit polarities in a first-motion focal mechanism determination. Decimal fraction between 0 and 1. :type station_distribution_ratio: float, optional :param station_distribution_ratio: Station distribution ratio (STDR) parameter. Indicates how the stations are distributed about the focal sphere (Reasenberg and Oppenheimer 1985). Decimal fraction between 0 and 1. :type method_id: :class:`~obspy.core.event.base.ResourceIdentifier`, optional :param method_id: Resource identifier of the method used for determination of the focal mechanism. :type waveform_id: list of :class:`~obspy.core.event.base.WaveformStreamID`, optional :param waveform_id: Refers to a set of waveform streams from which the focal mechanism was derived. :type evaluation_mode: str, optional :param evaluation_mode: Evaluation mode of FocalMechanism. See :class:`~obspy.core.event.header.EvaluationMode` for allowed values. :type evaluation_status: str, optional :param evaluation_status: Evaluation status of FocalMechanism. See :class:`~obspy.core.event.header.EvaluationStatus` for allowed values. :type moment_tensor: :class:`~obspy.core.event.source.MomentTensor`, optional :param moment_tensor: Moment tensor description for this focal mechanism. :type comments: list of :class:`~obspy.core.event.base.Comment`, optional :param comments: Additional comments. :type creation_info: :class:`~obspy.core.event.base.CreationInfo`, optional :param creation_info: Creation information used to describe author, version, and creation time. .. note:: For handling additional information not covered by the QuakeML standard and how to output it to QuakeML see the :ref:`ObsPy Tutorial <quakeml-extra>`. """
[docs]def farfield(mt, points, type): """ Returns the P/S farfield radiation pattern based on [Aki1980]_ eq. 4.29. :param mt: Focal mechanism NM x 6 (Mxx, Myy, Mzz, Mxy, Mxz, Myz - the six independent components of the moment tensor) :param points: 3D vector array with shape [3,npts] (x,y,z) or [2,npts] (theta,phi) The normalized displacement of the moment tensor source is computed at these points. :type type: str :param type: 'P' or 'S' (P or S wave). :return: 3D vector array with shape [3,npts] that contains the displacement vector for each grid point """ type = type.upper() if type not in ("P", "S"): msg = ("type must be 'P' or 'S'") raise ValueError(msg) is_p_wave = type == "P" ndim, npoints = points.shape if ndim == 2: # points are given as theta,phi _points = np.empty((3, npoints)) _points[0] = np.sin(points[0]) * np.cos(points[1]) _points[1] = np.sin(points[0]) * np.sin(points[1]) _points[2] = np.cos(points[0]) points = _points ndim = 3 elif ndim == 3: # points are given as x,y,z, (same system as the moment tensor) pass else: raise ValueError('points should have shape 2 x npoints or 3 x npoints') m_pq = _fullmt(mt) # precompute directional cosine array dists = np.sqrt(points[0] * points[0] + points[1] * points[1] + points[2] * points[2]) gammas = points / dists # initialize displacement array disp = np.empty((ndim, npoints)) # loop through points if is_p_wave: for ipoint in range(npoints): # loop through displacement component [n index] gamma = gammas[:, ipoint] gammapq = np.outer(gamma, gamma) gammatimesmt = gammapq * m_pq for n in range(ndim): disp[n, ipoint] = gamma[n] * np.sum(gammatimesmt.flatten()) else: for ipoint in range(npoints): # loop through displacement component [n index] gamma = gammas[:, ipoint] m_p = np.dot(m_pq, gamma) for n in range(ndim): psum = 0.0 for p in range(ndim): deltanp = int(n == p) psum += (gamma[n] * gamma[p] - deltanp) * m_p[p] disp[n, ipoint] = psum return disp
[docs]def _fullmt(mt): """takes 6 comp moment tensor and returns full 3x3 moment tensor""" mt_full = np.array(([[mt[0], mt[3], mt[4]], [mt[3], mt[1], mt[5]], [mt[4], mt[5], mt[2]]])) return mt_full
if __name__ == '__main__': import doctest doctest.testmod(exclude_empty=True)