Merge pull request #624 from gem/shahjouei-pezeshk-2016

Implements the Shahjouei & Pezeshk (2016) GMM

debian/changelog

@@ -1,3 +1,6 @@
+  [Abhineet Gupta, Anirudh Rao]
+  * Implements the Shahjouei & Pezeshk (2016) GMPE
+
   [Michele Simionato]
   * Experimental support for the Grid Engine
   * Fixed the serialization to HDF5 of nonparametric ruptures

doc/sphinx/openquake.hazardlib.gsim.rst

@@ -595,6 +595,14 @@ sadigh_1997
     :undoc-members:
     :show-inheritance:
 
+shahjouei_pezeshk_2016
+--------------------------------------------
+
+.. automodule:: openquake.hazardlib.gsim.shahjouei_pezeshk_2016
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
 sharma_2009
 -------------------------------------------
 

openquake/hazardlib/gsim/shahjouei_pezeshk_2016.py

@@ -0,0 +1,209 @@
+# -*- coding: utf-8 -*-
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+#
+# Copyright (C) 2013-2017 GEM Foundation
+#
+# OpenQuake is free software: you can redistribute it and/or modify it
+# under the terms of the GNU Affero General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# OpenQuake is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Module exports :class:'ShahjoueiPezeshk2016'.
+"""
+from __future__ import division
+
+import numpy as np
+
+from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
+from openquake.hazardlib import const
+from openquake.hazardlib.imt import PGA, PGV, SA
+
+
+class ShahjoueiPezeshk2016(GMPE):
+    """
+    Implements GMPE developed by Alireza Shahjouei and Shahram Pezeshk.
+    Published as "Alternative Hybrid Empirical Ground‐Motion Model for
+    Central and Eastern North America Using Hybrid Simulations and
+    NGA‐West2 Models", 2016, Bulletin of the Seismological
+    Society of America, vol. 106, no. 2, 734 - 754.
+    """
+    #: Supported tectonic region type is 'stable continental region'
+    #: equation has been derived from data from Eastern North America (ENA)
+    # 'Introduction', page 735.
+    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.STABLE_CONTINENTAL
+
+    #: Supported intensity measure types are spectral acceleration,
+    #: and peak ground acceleration. See Table 7 on page 743
+    DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
+        PGA,
+        PGV,
+        SA
+    ])
+
+    #: An orientation-independent alternative to :attr:`AVERAGE_HORIZONTAL`.
+    #: Defined at Boore et al. (2006, Bull. Seism. Soc. Am. 96, 1502-1511)
+    #: and is used for all the NGA GMPEs. See page 742.
+    #: :attr:'~openquake.hazardlib.const.IMC.RotD50'.
+    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.RotD50
+
+    #: Supported standard deviation types is total.
+    #: See equation 4 and 5, page 744.
+    #: We use aleatory uncertainty as the total
+    #: standard deviation since page 745 states,
+    #: The epistemic uncertainty for an individual GMM is infrequently employed...
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
+        const.StdDev.TOTAL
+    ])
+
+    #: No site parameters are needed. The GMPE was developed for hard-rock site
+    # with Vs30 >= 3000 m/s (NEHRP site class A) only. Page 734.
+    REQUIRES_SITES_PARAMETERS = set()
+
+    #: Required rupture parameters are magnitude (eq. 4, page 742).
+    REQUIRES_RUPTURE_PARAMETERS = set(('mag', ))
+
+    #: Required distance measure is Rjb (eq. 3 page 742).
+    REQUIRES_DISTANCES = set(('rjb', ))
+
+    #: GMPE not tested against independent implementation so raise
+    #: not verified warning
+    non_verified = True
+
+    def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
+        """
+        See :meth:`superclass method
+        <.base.GroundShakingIntensityModel.get_mean_and_stddevs>`
+        for spec of input and result values.
+        """
+
+        # Extracting dictionary of coefficients specific to required
+        # intensity measure type.
+        C = self.COEFFS[imt]
+
+        imean = (self._compute_magnitude(rup, C) +
+                 self._compute_attenuation(rup, dists, imt, C) +
+                 self._compute_distance(rup, dists, imt, C))
+
+        mean = np.log(10.0 ** (imean))
+
+        istddevs = self._get_stddevs(C, stddev_types, rup, imt,
+                                     num_sites=len(dists.rjb))
+
+        stddevs = np.log(10.0 ** np.array(istddevs))
+
+        return mean, stddevs
+
+    def _get_stddevs(self, C, stddev_types, rup, imt, num_sites):
+        """
+        Return standard deviations as defined in eq. 4 and 5, page 744,
+        based on table 8, page 744.
+        Eq. 5 yields std dev in natural log, so convert to log10
+        """
+        stddevs = []
+        for stddev_type in stddev_types:
+            sigma_mean = self._compute_standard_dev(rup, imt, C)
+            sigma_tot = np.sqrt((sigma_mean ** 2) + (C['SigmaReg'] ** 2))
+            sigma_tot = np.log10(np.exp(sigma_tot))
+            stddevs.append(sigma_tot + np.zeros(num_sites))
+        return stddevs
+
+    def _compute_magnitude(self, rup, C):
+        """
+        Compute the first term of the equation described on p. 742:
+
+        "c1 + (c2 * M) + (c3 * M**2) "
+        """
+        return C['c1'] + (C['c2'] * rup.mag) + (C['c3'] * (rup.mag ** 2))
+
+    def _compute_attenuation(self, rup, dists, imt, C):
+        """
+        Compute the second term of the equation described on p. 742:
+
+        " [(c4 + c5 * M) * min{ log10(R), log10(60.) }] +
+        [(c6 + c7 * M) * max{ min{ log10(R/60.), log10(120./60.) }, 0.}] +
+        [(c8 + c9 * M) * max{ log10(R/120.), 0}] "
+        """
+        vec = np.ones(len(dists.rjb))
+
+        a1 = (np.log10(np.sqrt(dists.rjb ** 2.0 + C['c11'] ** 2.0)),
+              np.log10(60. * vec))
+
+        a = np.column_stack([a1[0], a1[1]])
+
+        b3 = (np.log10(np.sqrt(dists.rjb ** 2.0 + C['c11'] ** 2.0) /
+                      (60. * vec)),
+              np.log10((120. / 60.) * vec))
+
+        b2 = np.column_stack([b3[0], b3[1]])
+        b1 = ([np.min(b2, axis=1), 0. * vec])
+        b = np.column_stack([b1[0], b1[1]])
+
+        c1 = (np.log10(np.sqrt(dists.rjb ** 2.0 + C['c11'] ** 2.0) /
+              (120.) * vec), 0. * vec)
+        c = np.column_stack([c1[0], c1[1]])
+
+        return (((C['c4'] + C['c5'] * rup.mag) * np.min(a, axis=1)) +
+                ((C['c6'] + C['c7'] * rup.mag) * np.max(b, axis=1)) +
+                ((C['c8'] + C['c9'] * rup.mag) * np.max(c, axis=1)))
+
+    def _compute_distance(self, rup, dists, imt, C):
+        """
+        Compute the third term of the equation described on p. 742:
+
+        " c10 * R "
+        """
+        return (C['c10'] * np.sqrt(dists.rjb ** 2.0 + C['c11'] ** 2.0))
+
+    def _compute_standard_dev(self, rup, imt, C):
+        """
+        Compute the the standard deviation in terms of magnitude
+        described on page 744, eq. 4
+        """
+        sigma_mean = 0.
+        if isinstance(imt, (PGA, SA)):
+            psi = -6.898E-3
+        else:
+            psi = -3.054E-5
+        if rup.mag <= 6.5:
+            sigma_mean = (C['c12'] * rup.mag) + C['c13']
+        elif rup.mag > 6.5:
+            sigma_mean = (psi * rup.mag) + C['c14']
+        return sigma_mean
+
+    #: Equation coefficients, described in Table 2 on pp. 1865
+    COEFFS = CoeffsTable(sa_damping=5, table="""
+    IMT    c1      c2      c3       c4      c5     c6       c7        c8      c9      c10        c11    c12     c13   c14   SigmaReg SigmaPar
+    pgv    -2.3891 1.259   -0.07901 -2.9386 0.3034 -0.00929 -0.04605  -2.7548 0.3467  -0.0007623 -4.598 -0.0554 0.978 0.663 0.1      0.288
+    pga    -0.3002 0.5066  -0.04526 -3.224  0.2998 -1.283   0.1045    -3.0856 0.2778  -0.0007711 3.81   -0.041  0.876 0.611 0.194    0.373
+    0.010  -0.3472 0.4838  -0.04093 -3.0832 0.2712 -0.9676  0.04983   -2.9695 0.2693  -0.0006695 -4.434 -0.056  0.982 0.664 0.132    0.281
+    0.020  0.832   0.1934  -0.0206  -3.1134 0.2786 -1.133   0.05994   -3.5023 0.2901  -0.0005857 -4.412 -0.0559 0.983 0.665 0.0928   0.281
+    0.030  1.185   0.1064  -0.01423 -3.1029 0.2792 -1.078   0.05239   -3.5722 0.2865  -0.000622  -4.353 -0.0577 1     0.676 0.0833   0.277
+    0.040  1.246   0.08986 -0.01268 -3.0785 0.2773 -0.9743  0.0416    -3.5083 0.2769  -0.0006818 -4.303 -0.0577 1.01  0.688 0.0798   0.279
+    0.050  1.1793  0.1037  -0.01321 -3.0488 0.2744 -0.8635  0.03077   -3.3986 0.2659  -0.0007439 -4.266 -0.0578 1.03  0.701 0.0776   0.272
+    0.075  0.8045  0.1866  -0.01788 -2.9697 0.266  -0.6122  0.007491  -3.0852 0.2391  -0.0008801 -4.214 -0.0561 1.03  0.721 0.0738   0.252
+    0.100  0.35    0.2871  -0.02381 -2.894  0.2576 -0.4123  -0.01012  -2.7947 0.2163  -0.0009848 4.201  -0.0565 1.05  0.732 0.0717   0.265
+    0.150  -0.5264 0.4782  -0.03519 -2.761  0.2426 -0.1319  -0.03338  -2.3312 0.1818  -0.001125  4.239  -0.0559 1.04  0.724 0.0716   0.276
+    0.200  -1.2884 0.6413  -0.04486 -2.6504 0.2301 0.04637  -0.0469   -1.9927 0.1576  -0.001209  4.325  -0.056  1.03  0.715 0.0743   0.258
+    0.250  -1.9422 0.7789  -0.05295 -2.5573 0.2196 0.1631   -0.05478  -1.7399 0.1398  -0.001258  4.438  -0.0537 1.02  0.712 0.0779   0.268
+    0.300  -2.5071 0.8961  -0.05976 -2.478  0.2107 0.2407   -0.05919  -1.547  0.1265  -0.001286  4.571  -0.0511 1.01  0.718 0.0815   0.284
+    0.400  -3.436  1.085   -0.07059 -2.3495 0.1961 0.3244   -0.06197  -1.2793 0.1085  -0.001304  -4.872 -0.047  0.987 0.725 0.0876   0.34
+    0.500  -4.1699 1.231   -0.07878 -2.251  0.1849 0.3544   -0.06046  -1.1111 0.09757 -0.001294  -5.211 -0.0442 0.981 0.736 0.0923   0.357
+    0.750  -5.4797 1.482   -0.09245 -2.0865 0.1659 0.3284   -0.04979  -0.9131 0.0857  -0.001219  -6.154 -0.0384 0.967 0.76  0.0991   0.374
+    1.000  -6.3464 1.641   -0.1006  -1.9931 0.1546 0.253    -0.03709  -0.8641 0.08405 -0.001123  -7.174 -0.0314 0.933 0.77  0.102    0.392
+    1.500  -7.4087 1.823   -0.1093  -1.9162 0.1438 0.09019  -0.01551  -0.92   0.09103 -0.0009407 -9.253 -0.0227 0.883 0.776 0.105    0.426
+    2.000  -8.0057 1.916   -0.113   -1.9173 0.1418 -0.03828 -0.001252 -1.0327 0.1016  -0.0007926 -11.22 -0.0184 0.857 0.778 0.106    0.44
+    3.000  -8.5793 1.985   -0.1146  -2.0184 0.1499 -0.1744  0.009393  -1.2453 0.1214  -0.0005919 14.38  -0.0189 0.859 0.777 0.107    0.58
+    4.000  -8.8246 1.99    -0.1131  -2.1475 0.1635 -0.1844  0.003919  -1.3849 0.1357  -0.0004855 16.19  -0.016  0.83  0.766 0.107    0.589
+    5.000  -8.9855 1.975   -0.1105  -2.2496 0.1764 -0.1043  -0.01187  -1.4511 0.1446  -0.0004439 16.71  -0.0153 0.826 0.766 0.107    0.631
+    7.500  -9.3927 1.925   -0.1032  -2.3572 0.1973 0.3465   -0.07832  -1.3728 0.149   -0.0005176 14.58  -0.0143 0.815 0.762 0.113    0.721
+    10.000 -9.735  1.879   -0.09666 -2.4139 0.2117 1.01     -0.1678   -1.0631 0.137   -0.000742  11.23  -0.017  0.822 0.752 0.14     0.739
+    """)