# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2014-2018 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:`ToroEtAl1997MblgNSHMP2008`,
:class:`ToroEtAl1997MwNSHMP2008`
"""
import numpy as np
from openquake.hazardlib.gsim.base import CoeffsTable, GMPE
from openquake.hazardlib.gsim.utils import (
mblg_to_mw_johnston_96, mblg_to_mw_atkinson_boore_87, clip_mean)
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, SA
[docs]class ToroEtAl1997MblgNSHMP2008(GMPE):
"""
Implements GMPE developed by G. R. Toro, N. A. Abrahamson, J. F. Sneider
and published in "Model of Strong Ground Motions from Earthquakes in
Central and Eastern North America: Best Estimates and Uncertainties"
(Seismological Research Letters, Volume 68, Number 1, 1997) as utilized
by the National Seismic Hazard Mapping Project (NSHMP) for the 2008 US
hazard model.
This class replicates the algorithm for the Toro et. al. 1997 GMPE as
coded in the subroutine ``getToro`` in the ``hazgridXnga2.f``
Fortran code available at:
http://earthquake.usgs.gov/hazards/products/conterminous/2008/software/
The class assumes rupture magnitude to be in Mblg scale (given that
MFDs for central and eastern US are given in this scale).
The equation implements also the finite-fault correction as given in
"Modification of the Toro et al. 1997 Attenuation Equations for Large
Magnitudes and Short Distances" (available at:
http://www.riskeng.com/downloads/attenuation_equations). The correction
uses Mw. Therefore Mblg is converted to Mw using both the Atkinson & Boore
1987 and Johnston 1996 conversion equations and an average correction term
is computed.
Coefficients are given for the B/C site conditions.
"""
#: Supported tectonic region type is stable continental crust,
#: given that the equations have been derived for central and eastern
#: north America
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.STABLE_CONTINENTAL
#: Supported intensity measure types are spectral acceleration,
#: and peak ground acceleration
DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
PGA,
SA
])
#: Supported intensity measure component is the geometric mean of
#two : horizontal components
#:attr:`~openquake.hazardlib.const.IMC.AVERAGE_HORIZONTAL`,
DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.AVERAGE_HORIZONTAL
#: Supported standard deviation type is only total.
DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
const.StdDev.TOTAL
])
#: No site parameters required
REQUIRES_SITES_PARAMETERS = set()
#: Required rupture parameter is only magnitude (Mblg).
REQUIRES_RUPTURE_PARAMETERS = set(('mag', ))
#: Required distance measure is rjb
REQUIRES_DISTANCES = set(('rjb', ))
[docs] 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.
"""
assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
for stddev_type in stddev_types)
C = self.COEFFS[imt]
mean = self._compute_mean(C, rup.mag, dists.rjb)
stddevs = self._compute_stddevs(C, dists.rjb.size, stddev_types)
mean = clip_mean(imt, mean)
return mean, stddevs
def _compute_mean(self, C, mag, rjb):
"""
Compute ground motion mean value.
"""
# line 1686 in hazgridXnga2.f
ffc = self._compute_finite_fault_correction(mag)
d = np.sqrt(rjb ** 2 + (C['c7'] ** 2) * (ffc ** 2))
# lines 1663, 1694-1696 in hazgridXnga2.f
mean = (
C['c1'] + C['c2'] * (mag - 6.) +
C['c3'] * ((mag - 6.) ** 2) -
C['c4'] * np.log(d) - C['c6'] * d
)
factor = np.log(rjb / 100.)
idx = factor > 0
mean[idx] -= (C['c5'] - C['c4']) * factor[idx]
return mean
def _compute_finite_fault_correction(self, mag):
"""
Compute finite fault correction term as geometric mean of correction
terms obtained from Mw values calculated with Johnston 1996 and
Atkinson and Boore 1987 conversion equations.
Implement equations as in lines 1653 - 1658 in hazgridXnga2.f
"""
mw_j96 = mblg_to_mw_johnston_96(mag)
mw_ab87 = mblg_to_mw_atkinson_boore_87(mag)
t1 = np.exp(-1.25 + 0.227 * mw_j96)
t2 = np.exp(-1.25 + 0.227 * mw_ab87)
return np.sqrt(t1 * t2)
def _compute_stddevs(self, C, num_sites, stddev_types):
"""
Return total standard deviation.
"""
stddevs = []
for _ in stddev_types:
stddevs.append(np.zeros(num_sites) + C['sigma'])
return stddevs
#: Coefficient table obtained from coefficient arrays (tb1, tb2, tb3, tb4,
#: tb5, tb6, tbh) defined from line 1596 - 1614 in hazgridXnga2.f
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT c1 c2 c3 c4 c5 c6 c7 sigma
pga 2.489 1.20 0.0 1.28 1.23 0.0018 9.3 0.7506
0.1 2.91 1.23 0.0 1.12 1.05 0.0043 8.5 0.7506
0.2 2.165 1.24 0.0 0.98 0.74 0.0039 7.5 0.7506
0.3 1.7323 1.51 -0.11 0.96 0.6881 0.0034 7.35 0.7506
0.5 1.109 1.785 -0.2795 0.93 0.6354 0.002732 7.05 0.7506
1.0 0.173 2.05 -0.34 0.90 0.59 0.0019 6.8 0.799
2.0 -0.788 2.52 -0.47 0.93 0.6 0.0012 7.0 0.799
""")
[docs]class ToroEtAl1997MwNSHMP2008(ToroEtAl1997MblgNSHMP2008):
"""
Extend :class:`ToroEtAl1997MblgNSHMP2008` but assumes magnitude to be in
Mw scale.
Coefficients are Mw-specific and no magnitude conversion is considered to
take into account finite-fault correction.
"""
def _compute_finite_fault_correction(self, mag):
"""
Compute finite fault correction term.
"""
return np.exp(-1.25 + 0.227 * mag)
#: Coefficient table obtained from coefficient arrays (tc1, tc2, tc3, tc4,
#: tc5, tc6, th) defined in subroutine getToro in hazgridXnga2.f
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT c1 c2 c3 c4 c5 c6 c7 sigma
pga 2.619 0.81 0.0 1.27 1.16 0.0021 9.3 0.7506
0.1 2.92 0.81 0.0 1.1 1.02 0.004 8.3 0.7506
0.2 2.295 0.84 0.0 0.98 0.66 0.0042 7.5 0.7506
0.3 1.8823 0.964 -0.059 0.951 0.601 0.00367 7.26 0.7506
0.5 1.2887 1.14 -0.1244 0.9227 0.5429 0.00306 7.027 0.7506
1.0 0.383 1.42 -0.2 0.90 0.49 0.0023 6.8 0.799
2.0 -0.558 1.86 -0.31 0.92 0.46 0.0017 6.9 0.799
""")