Merge remote-tracking branch 'upstream/master'

shakyground2/earthquake.py

@@ -75,48 +75,48 @@ class Earthquake(object):
         date=None,
         time=None,
     ):
-    """
-    Initialize a new Earthquake object
+        """
+        Initialize a new Earthquake object
 
-    Args:
+        Args:
 
-    id:
-        Unique ID for the earthquake
-    lon:
-        Earthquake longitude in decimal degrees (-180 to 180)
-    lat:
-        Earthquake latitude in decimal degrees (-90, to 90)
-    hypo_depth:
-        Hypocentral depth (in km)
-    magnitude:
-        Earthquake magnitude
-    strike:
-        Strike of the rupture (in degrees from north)
-    dip:
-        Dip of the rupture (in degrees from horizontal)
-    rake:
-        Rake of the earthquake rupture in decimal degrees according to the
-        Aki & Richards (1980) convention (-180 to 180)
-    aspect:
-        Aspect ratio (length / width) of the rupture
-    mag_scaling_relation:
-        Magnitude scaling relation as instance of the BaseScalingRelation class or None.
-        Defaults to :class:`shakyground2.magnitude_scaling_relations.PEERScalingRelation`
-        when missing
-    upper_seismogenic_depth:
-        Upper seismogenic depth (in km)
-    lower_seismogenic_depth:
-        Lower seismogenic depth (in km)
-    surface:
-        Rupture surface as instance of :class:`openquake.hazardlib.geo.BaseSurface`, or None
-    hypocenter_position:
-        Hypocentre location within the rupture plane as a tuple of
-        (fraction of along strike length, fraction of down-dip width)
-    date:
-        Date of the earthquake (year, month, day) as datetime.date object
-    time:
-        Time of the earthquake (hour, minute, second) as datetime.time object
-    """
+        id:
+            Unique ID for the earthquake
+        lon:
+            Earthquake longitude in decimal degrees (-180 to 180)
+        lat:
+            Earthquake latitude in decimal degrees (-90, to 90)
+        hypo_depth:
+            Hypocentral depth (in km)
+        magnitude:
+            Earthquake magnitude
+        strike:
+            Strike of the rupture (in degrees from north)
+        dip:
+            Dip of the rupture (in degrees from horizontal)
+        rake:
+            Rake of the earthquake rupture in decimal degrees according to the
+            Aki & Richards (1980) convention (-180 to 180)
+        aspect:
+            Aspect ratio (length / width) of the rupture
+        mag_scaling_relation:
+            Magnitude scaling relation as instance of the BaseScalingRelation class or None.
+            Defaults to :class:`shakyground2.magnitude_scaling_relations.PEERScalingRelation`
+            when missing
+        upper_seismogenic_depth:
+            Upper seismogenic depth (in km)
+        lower_seismogenic_depth:
+            Lower seismogenic depth (in km)
+        surface:
+            Rupture surface as instance of :class:`openquake.hazardlib.geo.BaseSurface`, or None
+        hypocenter_position:
+            Hypocentre location within the rupture plane as a tuple of
+            (fraction of along strike length, fraction of down-dip width)
+        date:
+            Date of the earthquake (year, month, day) as datetime.date object
+        time:
+            Time of the earthquake (hour, minute, second) as datetime.time object
+        """
         self.id = id
         self.lon = valid.longitude(lon)
         self.lat = valid.latitude(lat)

shakyground2/rupture_mechanism.py

+"""
+Class to manage the distribution of rupture mechanisms
+"""
+from __future__ import annotations
+import numpy as np
+from enum import Enum
+from typing import Optional, Dict, Tuple
+from openquake.hazardlib.geo import NodalPlane
+from openquake.hazardlib.pmf import PMF
+from shakyground2 import valid
+
+
+class DIP_RANGES(Enum):
+    """
+    For the "global" default rupture mechanism distribution the possible ranges
+    of dip values depend on the style-of-faulting. The possible range of values
+    for each style-of-faulting is provided and they are assumed to be uniformly
+    distributed.
+    """
+    R = (20.0, 45.0)
+    SS = (75.0, 91.0)
+    N = (45.0, 75.0)
+
+
+class RuptureMechanism(object):
+    """
+    General class to handle the mechanism properties of the rupture, which include the
+    strike, dip and rake. The mechanism properties largely constrol the dimensions of the 3D
+    rupture plane that is needed for calculating the finite-source to site distances (e.g.
+    Joyner-Boore distance, Rupture Distance, Rx, Ry0 etc.)
+
+    for the majority of reported earthquakes a single rupture mechanism is not known unless
+    accompanied by a 3D finite fault rupture model. In the absence of any information other
+    than the hypocentre, a "global" distribution of strike, dip and rake values is assumed.
+    In the case that a focal mechanism is available then the RuptureMechanism is described by
+    two equiprobable planes.
+
+    Attributes:
+        mechanism: Distribution of rupture mechanisms (probability mass function) as instance
+                   of :class:`openquake.hazardlib.pmf.PMF`
+
+    """
+    def __init__(self, mechanism: Optional[PMF] = None):
+        if mechanism:
+            assert isinstance(mechanism, PMF)
+            self.mechanism = mechanism
+        else:
+            # Build the default distribution
+            self.mechanism = self.build_mechanism_distribution()
+
+    def __iter__(self):
+        # Iterate over the mechanisms in the data set yielding the probability
+        # of the nodal plane and the nodal plane itself
+        for prob, nodal_plane in self.mechanism.data:
+            yield prob, nodal_plane
+
+    def sample(self, nsamples: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+        """
+        Samples the focal mechanism distribution to return "nsamples" strike,
+        dip and rake values
+
+        Args:
+            nsamples: Number of sample
+        Returns:
+            Tuple of vectors of samples strikes, dips and rakes
+        """
+        assert nsamples > 0
+        sample_planes = self.mechanism.sample(nsamples)
+        strikes = np.empty(nsamples, dtype=float)
+        dips = np.empty(nsamples, dtype=float)
+        rakes = np.empty(nsamples, dtype=float)
+        for i, sample_plane in enumerate(sample_planes):
+            strikes[i] = sample_plane.strike
+            dips[i] = sample_plane.dip
+            rakes[i] = sample_plane.rake
+        return strikes, dips, rakes
+
+    @classmethod
+    def from_strike_dip_rake(cls, strike: Optional[float] = None, dip: Optional[float] = None,
+                             rake: Optional[float] = None) -> RuptureMechanism:
+        """
+        Constructs the rupture mechanism distribution from a simple strike, dip
+        and rake combination, permitting None values if undefined
+
+        Args:
+            strike: Strike of fault (in decimal degrees between 0 and 360)
+            dip: Dip of fault (in decimal degrees between 0 and 90)
+            rake: Rake of fault (in decimal degrees between -180 and 180)
+        """
+        return cls(cls.build_mechanism_distribution(valid.strike(strike),
+                                                    valid.dip(dip),
+                                                    valid.rake(rake)))
+
+    @classmethod
+    def from_nodal_planes(cls, nodal_plane_1: Dict, nodal_plane_2: Dict)\
+            -> RuptureMechanism:
+        """
+        Constructs the rupture mechanism distribution from either one single
+        valid rupture plane or two valud nodal planes
+
+        Args:
+            nodal_plane_1: Dictionary containing strike, dip and rake of the first nodel plane
+            nodal_plane_1: Dictionary containing strike, dip and rake of the second nodel plane
+        """
+        strike_1 = valid.strike(nodal_plane_1["strike"])
+        dip_1 = valid.dip(nodal_plane_1["dip"])
+        rake_1 = valid.rake(nodal_plane_1["rake"])
+        strike_2 = valid.strike(nodal_plane_2["strike"])
+        dip_2 = valid.dip(nodal_plane_2["dip"])
+        rake_2 = valid.rake(nodal_plane_2["rake"])
+        return cls(PMF([(0.5, NodalPlane(strike_1, dip_1, rake_1)),
+                        (0.5, NodalPlane(strike_2, dip_2, rake_2))]))
+
+    @staticmethod
+    def build_mechanism_distribution(strike: Optional[float] = None,
+                                     dip: Optional[float] = None,
+                                     rake: Optional[float] = None) -> PMF:
+        """
+        Builds a mechanism distribution from a partial (or complete) characterisation of the
+        rupture mechanism
+
+        Args:
+            strike: Strike of fault in decimal degrees between 0 and 360, or None
+            dip: Dip of fault in decimal degrees between 0 and 90, or None
+            rake: Rake of fault in decimal degrees between -180 and 180, or None
+
+        Returns:
+            Probability mass function of the mechanism distribution as an instance of the
+            :class:`openquake.hazardlib.pmf.PMF`
+        """
+        if rake is None:
+            # If rake is not defined then describe a complete range of rakes every 15 degrees
+            # from -165.0 to 180.0
+            rakes = np.arange(-165.0, 181.0, 15.0)
+            weight_r = (1.0 / len(rakes)) * np.ones(len(rakes))
+        else:
+            rakes = [rake]
+            weight_r = [1.0]
+
+        if strike is None:
+            # If strike is not defined then describe a complete range of strikes every 15
+            # degrees between 0 and 360 (not included)
+            strikes = np.arange(0.0, 359.0, 15.0)
+            weight_s = (1.0 / len(strikes)) * np.ones(len(strikes))
+        else:
+            strikes = [strike]
+            weight_s = [1.0]
+
+        mechanisms = []
+        for wght1, rake_i in zip(weight_r, rakes):
+            if dip is None:
+                # If dip is undefined then sample uniformly in the range
+                # appropriate to the style of faulting
+                if rake_i >= 45.0 and rake_i < 135.0:
+                    # Reverse
+                    dip_range = DIP_RANGES.R.value
+                elif rake_i < -45.0 and rake_i >= -135.0:
+                    # Normal
+                    dip_range = DIP_RANGES.N.value
+                else:
+                    # Strike-slip
+                    dip_range = DIP_RANGES.SS.value
+                dips = np.arange(dip_range[0], dip_range[1], 5.0)
+                weight_d = (1.0 / len(dips)) * np.ones(len(dips))
+            else:
+                dips = [dip]
+                weight_d = [1.0]
+            for wght2, dip_i in zip(weight_d, dips):
+                for wght3, strike_i in zip(weight_s, strikes):
+                    prob = wght1 * wght2 * wght3
+                    mechanisms.append((prob, NodalPlane(strike_i, dip_i, rake_i)))
+        return PMF(mechanisms)

tests/rupture_mechanism_test.py

+"""
+Unit tests for rupture mechanism
+"""
+import unittest
+import numpy as np
+from openquake.hazardlib.pmf import PMF
+from shakyground2.rupture_mechanism import RuptureMechanism
+
+
+class RuptureMechanismTestCase(unittest.TestCase):
+    """
+    Test the instatiation and operation of the rupture mechanism class
+    """
+    @staticmethod
+    def _mechanism_pmf_to_arrays(mech):
+        """
+        Takes a rupture mechanism defined as a probability mass function and
+        returns the information in respective vectors of probabilities, strikes,
+        dips and rakes
+        """
+        nvals = len(mech.mechanism.data)
+        probs = np.empty(nvals, dtype=float)
+        strikes = np.empty(nvals, dtype=float)
+        dips = np.empty(nvals, dtype=float)
+        rakes = np.empty(nvals, dtype=float)
+        for i, (prob, npl) in enumerate(mech):
+            probs[i] = prob
+            strikes[i] = npl.strike
+            dips[i] = npl.dip
+            rakes[i] = npl.rake
+        return probs, strikes, dips, rakes
+
+    def test_rupture_mechanism_no_prior_information(self):
+        # Checks instantiation of the mechanism when no prior rupture information
+        # is defined - i.e. global distribution of mechanisms
+        mech = RuptureMechanism()
+        probs, strikes, dips, rakes = self._mechanism_pmf_to_arrays(mech)
+        np.testing.assert_array_almost_equal(np.unique(strikes),
+                                             np.arange(0.0, 359.0, 15.0))
+        np.testing.assert_array_almost_equal(np.unique(dips),
+                                             np.arange(20.0, 95.0, 5.0))
+        np.testing.assert_array_almost_equal(np.unique(rakes),
+                                             np.arange(-165.0, 181.0, 15.0))
+
+    def test_rupture_mechanism_single_plane(self):
+        # Checks instantiation when a single combination of strike, dip and rake are input
+        mech = RuptureMechanism.from_strike_dip_rake(strike=20., dip=90.0, rake=10.0)
+        probs, strikes, dips, rakes = self._mechanism_pmf_to_arrays(mech)
+        np.testing.assert_array_almost_equal(probs, np.array([1.0]))
+        np.testing.assert_array_almost_equal(strikes, np.array([20.0]))
+        np.testing.assert_array_almost_equal(dips, np.array([90.0]))
+        np.testing.assert_array_almost_equal(rakes, np.array([10.0]))
+
+    def test_rupture_mechanism_from_nodal_planes(self):
+        # Checks the creation of the class from a pair of nodal planes
+        npl1 = {"strike": 45.0, "dip": 80.0, "rake": 10.}
+        npl2 = {"strike": 135.0, "dip": 90.0, "rake": -170.0}
+        mech = RuptureMechanism.from_nodal_planes(npl1, npl2)
+        probs, strikes, dips, rakes = self._mechanism_pmf_to_arrays(mech)
+        np.testing.assert_array_almost_equal(probs, np.array([0.5, 0.5]))
+        np.testing.assert_array_almost_equal(strikes, np.array([45.0, 135.0]))
+        np.testing.assert_array_almost_equal(dips, np.array([80.0, 90.0]))
+        np.testing.assert_array_almost_equal(rakes, np.array([10.0, -170.0]))
+
+    def test_sample_mechanism(self):
+        # Test the random sampling
+        npl1 = {"strike": 45.0, "dip": 80.0, "rake": 10.}
+        npl2 = {"strike": 135.0, "dip": 90.0, "rake": -170.0}
+        mech = RuptureMechanism.from_nodal_planes(npl1, npl2)
+        np.random.seed(1000)
+        strikes, dips, rakes = mech.sample(5)
+        np.testing.assert_array_almost_equal(strikes,
+                                             np.array([135.0, 45.0, 135.0, 45.0, 135.0]))
+        np.testing.assert_array_almost_equal(dips,
+                                             np.array([90.0, 80.0, 90.0, 80.0, 90.0]))
+        np.testing.assert_array_almost_equal(rakes,
+                                             np.array([-170.0, 10.0, -170.0, 10.0, -170.0]))