Source code for openquake.hazardlib.gsim.lanzano_2016

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2014-2025 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:`LanzanoEtAl2016`.
"""
import numpy as np
from scipy.constants import g

from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib.gsim import utils
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, PGV, SA
from openquake.baselib.general import CallableDict

_compute_distance = CallableDict()


@_compute_distance.add('rjb')
def _compute_distance_rjb(dist_type, ctx, C):
    """
    Compute the third term of the equation 1:
    FD(Mw,R) = [c1j + c2j(M-Mr)] * log10(R/Rh) con j=1,...4 (eq 4)
    c coeffs are in matrix C
    """
    Mr = 5.0
    Rh = 70

    LATref = -0.33 * ctx.lon + 48.3
    diff = ctx.lat - LATref
    R = np.sqrt(ctx.rjb**2 + C['h']**2)

    dist_term = (diff >= 0) * (C['c11'] + C['c21'] * (ctx.mag - Mr)) *\
                (R <= Rh) * np.log10(R/Rh) +\
                (diff >= 0) * (C['c12'] + C['c22'] * (ctx.mag - Mr)) *\
                (R > Rh) * np.log10(R/Rh) +\
                (diff < 0) * (C['c13'] + C['c23'] * (ctx.mag - Mr)) *\
                (R <= Rh) * np.log10(R/Rh) +\
                (diff < 0) * (C['c14'] + C['c24'] * (ctx.mag - Mr)) *\
                (R > Rh) * np.log10(R/Rh)
    return dist_term


@_compute_distance.add('rhypo')
def _compute_distance_rhypo(dist_type, ctx, C):
    """
    Compute the third term of the equation 1:
    FD(Mw,R) = [c1j + c2j(M-Mr)] * log10(R/Rh) con j=1,...4 (eq 4)
    c coeffs are in matrix C
    """
    Mr = 5.0
    Rh = 70

    LATref = -0.33 * ctx.lon + 48.3
    diff = ctx.lat - LATref
    R = ctx.rhypo

    dist_term = (diff >= 0) * (C['c11'] + C['c21'] * (ctx.mag - Mr)) *\
                (R <= Rh) * np.log10(R/Rh) +\
                (diff >= 0) * (C['c12'] + C['c22'] * (ctx.mag - Mr)) *\
                (R > Rh) * np.log10(R/Rh) +\
                (diff < 0) * (C['c13'] + C['c23'] * (ctx.mag - Mr)) *\
                (R <= Rh) * np.log10(R/Rh) +\
                (diff < 0) * (C['c14'] + C['c24'] * (ctx.mag - Mr)) *\
                (R > Rh) * np.log10(R/Rh)
    return dist_term


def _compute_magnitude(ctx, C):
    """
    Compute the second term of the equation 1:
    Fm(M) = b1(M-Mr) + b2(M-Mr)^2   Eq (5)
    """
    Mr = 5
    return C['a'] + C['b1'] * (ctx.mag - Mr) + C['b2'] * (ctx.mag - Mr)**2


def _get_basin_effect_term(ctx, C):
    """
    Get basin correction for ctx in the Po Plain.
    if ctx.bas == 0 the correction is not necessary,
    otherwise if ctx.bas == 1 the site is in the Po Plain
    and the correction is applied.
    """
    delta = np.zeros(len(ctx.vs30))
    delta[ctx.bas == 1] = 1.0

    return C['dbas'] * delta


def _get_mechanism(ctx, C):
    """
    Compute the part of the second term of the equation 1 (FM(SoF)):
    Get fault type dummy variables
    """
    UN, NF, TF = utils.get_fault_type_dummy_variables(ctx)
    return C['fNF'] * NF + C['fTF'] * TF + C['fUN'] * UN


def _get_site_amplification(ctx, C):
    """
    Compute the fourth term of the equation 1 described on paragraph :
    The functional form Fs in Eq. (1) represents the site amplification and
    it is given by FS = sj Cj , for j = 1,...,3, where sj are the
    coefficients to be determined through the regression analysis,
    while Cj are dummy variables used to denote the five different EC8
    site classes
    """
    ssa, ssb, ssc = _get_site_type_dummy_variables(ctx)
    return (C['sA'] * ssa) + (C['sB'] * ssb) + (C['sC'] * ssc)


def _get_site_type_dummy_variables(ctx):
    """
    Get site type dummy variables, five different EC8 site classes
    The recording ctx are classified into 3 classes,
    based on the shear wave velocity intervals in the uppermost 30 m, Vs30,
    according to the EC8 (CEN 2003):
    class A: Vs30 > 800 m/s
    class B: Vs30 = 360 - 800 m/s
    class C: Vs30 < 360 m/s
    """
    ssa = np.zeros(len(ctx.vs30))
    ssb = np.zeros(len(ctx.vs30))
    ssc = np.zeros(len(ctx.vs30))

    # Class C; 180 m/s <= Vs30 <= 360 m/s.
    idx = ctx.vs30 < 360.0
    ssc[idx] = 1.0
    # Class B; 360 m/s <= Vs30 <= 800 m/s.
    idx = (ctx.vs30 >= 360.0) & (ctx.vs30 < 800)
    ssb[idx] = 1.0
    # Class A; Vs30 > 800 m/s.
    idx = (ctx.vs30 >= 800.0)
    ssa[idx] = 1.0
    return ssa, ssb, ssc


[docs]class LanzanoEtAl2016_RJB(GMPE):
    """
    Implements GMPE developed by G.Lanzano, M. D'Amico, C.Felicetta, R.Puglia,
    L.Luzi, F.Pacor, D.Bindi and published as
    "Ground-Motion Prediction Equations
    for Region-Specific Probabilistic Seismic-Hazard Analysis",
    Bull Seismol. Soc. Am., DOI 10.1785/0120150096
    SA are given up to 4 s.
    The regressions are developed considering the geometrical mean of the
    as-recorded horizontal components
    """
    #: Supported tectonic region type is 'active shallow crust' because the
    #: equations have been derived from data from Italian database ITACA, as
    #: explained in the 'Introduction'.
    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST

    #: Set of :mod:`intensity measure types <openquake.hazardlib.imt>`
    #: this GSIM can calculate. A set should contain classes from module
    #: :mod:`openquake.hazardlib.imt`.
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, PGV, SA}

    #: Supported intensity measure component is the geometric mean of two
    #: horizontal components
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.GEOMETRIC_MEAN

    #: Supported standard deviation types are inter-event, intra-event
    #: and total
    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
        const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT}

    #: Required site parameter
    REQUIRES_SITES_PARAMETERS = {'vs30', 'lon', 'lat', 'bas'}

    #: Required rupture parameters are magnitude and rake (eq. 1).
    REQUIRES_RUPTURE_PARAMETERS = {'rake', 'mag'}

    #: Required distance measure is R Joyner-Boore distance (eq. 1).
    REQUIRES_DISTANCES = {'rjb'}

[docs]    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.compute>`
        for spec of input and result values.
        """
        for m, imt in enumerate(imts):
            C = self.COEFFS[imt]
            dist_type = 'rjb' if "RJB" in self.__class__.__name__ else 'rhypo'

            imean = (_compute_magnitude(ctx, C) +
                     _compute_distance(dist_type, ctx, C) +
                     _get_site_amplification(ctx, C) +
                     _get_basin_effect_term(ctx, C) +
                     _get_mechanism(ctx, C))

            # Convert units to g, but only for PGA and SA (not PGV):
            if imt.string.startswith(("SA", "PGA")):
                mean[m] = np.log((10.0 ** (imean - 2.0)) / g)
            else:
                # PGV
                mean[m] = np.log(10.0 ** imean)

            # Return stddevs in terms of natural log scaling
            sig[m] = np.log(10.0 ** C['SigmaTot'])
            tau[m] = np.log(10.0 ** C['tau'])
            phi[m] = np.log(10.0 ** C['phi'])

    #: Coefficients from SA PGA and PGV from Table S2

    COEFFS = CoeffsTable(sa_damping=5, table="""
    IMT     a       b1      b2      c11     c21     c12     c22     c13     c23     c14     c24     h       fNF     fTF     fUN     sA      sB      sC      dbas    tau     phi     SigmaTot
    0.040   0.122   0.565   -0.015  -1.967  0.305   -0.968  0.026   -1.872  0.567   -2.280  0.443   6.507   0.058   0.199   0.000   0.000   0.028   0.189   -0.094  0.108   0.324   0.342
    0.070   0.258   0.555   -0.017  -2.131  0.299   -0.885  0.087   -1.986  0.593   -2.501  0.486   7.296   0.042   0.182   0.000   0.000   -0.002  0.161   -0.099  0.113   0.340   0.359
    0.100   0.334   0.537   -0.011  -2.125  0.268   -0.742  0.123   -1.996  0.486   -2.573  0.451   7.463   0.032   0.193   0.000   0.000   0.018   0.180   -0.102  0.118   0.354   0.373
    0.150   0.431   0.561   -0.008  -1.979  0.235   -0.804  0.070   -1.890  0.452   -2.443  0.348   7.073   0.024   0.183   0.000   0.000   0.033   0.180   -0.081  0.118   0.353   0.372
    0.200   0.436   0.579   -0.014  -1.847  0.201   -0.810  0.066   -1.820  0.348   -2.321  0.364   6.698   0.022   0.191   0.000   0.000   0.058   0.209   -0.070  0.114   0.342   0.360
    0.250   0.436   0.610   -0.020  -1.790  0.196   -0.905  0.022   -1.683  0.357   -2.205  0.311   6.933   0.035   0.183   0.000   0.000   0.065   0.201   -0.029  0.110   0.331   0.349
    0.300   0.394   0.629   -0.015  -1.759  0.180   -0.963  0.009   -1.639  0.305   -2.065  0.343   7.016   0.038   0.177   0.000   0.000   0.075   0.219   -0.016  0.108   0.323   0.340
    0.350   0.335   0.644   -0.008  -1.724  0.165   -1.009  0.015   -1.631  0.243   -1.979  0.328   7.304   0.044   0.178   0.000   0.000   0.081   0.241   0.012   0.105   0.316   0.333
    0.400   0.284   0.662   -0.010  -1.662  0.155   -1.080  0.004   -1.641  0.182   -1.863  0.282   7.272   0.031   0.176   0.000   0.000   0.088   0.257   0.046   0.103   0.309   0.326
    0.450   0.240   0.683   -0.010  -1.646  0.146   -1.045  -0.017  -1.648  0.169   -1.737  0.266   7.500   0.025   0.169   0.000   0.000   0.078   0.255   0.073   0.102   0.305   0.322
    0.500   0.182   0.699   -0.009  -1.601  0.140   -1.027  -0.027  -1.638  0.166   -1.663  0.270   7.493   0.030   0.168   0.000   0.000   0.086   0.270   0.094   0.101   0.303   0.319
    0.600   0.084   0.727   -0.017  -1.549  0.107   -0.963  0.000   -1.584  0.170   -1.492  0.244   7.135   0.024   0.150   0.000   0.000   0.097   0.280   0.125   0.101   0.303   0.319
    0.700   0.010   0.751   -0.025  -1.509  0.091   -0.956  -0.001  -1.561  0.134   -1.373  0.209   7.044   0.006   0.131   0.000   0.000   0.108   0.292   0.123   0.100   0.301   0.317
    0.800   -0.051  0.771   -0.034  -1.460  0.100   -0.989  0.028   -1.520  0.149   -1.307  0.225   6.943   -0.004  0.122   0.000   0.000   0.104   0.292   0.123   0.101   0.302   0.318
    0.900   -0.114  0.799   -0.036  -1.417  0.114   -0.985  0.016   -1.463  0.180   -1.212  0.259   6.760   -0.004  0.108   0.000   0.000   0.100   0.293   0.131   0.100   0.300   0.316
    1.000   -0.158  0.827   -0.043  -1.373  0.130   -1.009  0.000   -1.411  0.206   -1.189  0.207   6.500   -0.003  0.098   0.000   0.000   0.091   0.289   0.136   0.100   0.300   0.316
    1.200   -0.260  0.879   -0.041  -1.316  0.166   -1.038  -0.020  -1.313  0.267   -1.173  0.121   6.026   0.004   0.091   0.000   0.000   0.087   0.289   0.131   0.100   0.299   0.315
    1.400   -0.323  0.909   -0.052  -1.264  0.157   -1.139  0.053   -1.208  0.260   -1.262  0.116   5.366   -0.004  0.075   0.000   0.000   0.082   0.293   0.119   0.099   0.298   0.314
    1.600   -0.409  0.946   -0.045  -1.238  0.150   -1.214  0.046   -1.128  0.323   -1.278  0.072   5.033   0.000   0.070   0.000   0.000   0.085   0.303   0.115   0.099   0.297   0.313
    1.800   -0.486  0.977   -0.039  -1.218  0.141   -1.239  0.066   -1.103  0.313   -1.315  0.014   4.737   -0.001  0.067   0.000   0.000   0.080   0.303   0.119   0.100   0.299   0.315
    2.000   -0.554  0.997   -0.037  -1.189  0.138   -1.263  0.069   -1.100  0.282   -1.345  0.057   4.241   -0.010  0.056   0.000   0.000   0.081   0.308   0.117   0.099   0.298   0.314
    2.500   -0.742  1.034   -0.027  -1.164  0.151   -1.326  0.045   -1.072  0.299   -1.385  0.060   4.126   0.040   0.054   0.000   0.000   0.085   0.327   0.115   0.100   0.301   0.317
    3.000   -0.881  1.057   -0.019  -1.152  0.165   -1.378  0.018   -1.020  0.339   -1.449  0.084   4.170   0.072   0.045   0.000   0.000   0.089   0.325   0.114   0.101   0.304   0.320
    4.000   -1.084  1.134   0.019   -1.101  0.244   -1.488  -0.153  -0.971  0.414   -1.619  -0.119  4.454   0.073   0.019   0.000   0.000   0.096   0.322   0.131   0.113   0.298   0.318
    pga     0.071   0.603   -0.019  -1.895  0.286   -0.926  0.035   -1.838  0.511   -2.256  0.455   6.701   0.035   0.181   0.000   0.000   0.050   0.203   -0.060  0.106   0.318   0.336
    pgv     -1.142  0.767   -0.005  -1.623  0.230   -1.037  -0.054  -1.596  0.379   -1.741  0.348   5.904   0.022   0.144   0.000   0.000   0.085   0.260   0.037   0.096   0.288   0.304
    """)


[docs]class LanzanoEtAl2016_Rhypo(LanzanoEtAl2016_RJB):
    """
    Implements GMPE developed by G.Lanzano, M. D'Amico, C.Felicetta, R.Puglia,
    L.Luzi, F.Pacor, D.Bindi and published as "Ground-Motion Prediction Equations
    for Region-Specific Probabilistic Seismic-Hazard Analysis",
    Bull Seismol. Soc. Am., DOI 10.1785/0120150096
    SA are given up to 4 s.
    The regressions are developed considering the geometrical mean of the
    as-recorded horizontal components
    """

    #: Required distance measure is R Hypocentral Distance.
    REQUIRES_DISTANCES = {'rhypo'}

    #: Coefficients from SA PGA and PGV from esupp Table S3
    COEFFS = CoeffsTable(sa_damping=5, table="""
    IMT     a       b1      b2      c11     c21     c12     c22     c13     c23     c14     c24     fNF     fTF     fUN     sA      sB      sC      dbas    tau     phi     SigmaTot
    0.040   0.096   0.586   0.036   -2.119  0.151   -0.970  -0.009  -2.191  0.438   -2.263  0.428   0.029   0.207   0.000   0.000   0.040   0.199   -0.084  0.111   0.333   0.351
    0.070   0.236   0.569   0.030   -2.242  0.137   -0.918  0.055   -2.270  0.437   -2.509  0.497   0.018   0.196   0.000   0.000   0.010   0.173   -0.094  0.116   0.348   0.367
    0.100   0.318   0.553   0.035   -2.221  0.110   -0.785  0.083   -2.281  0.305   -2.586  0.463   0.005   0.205   0.000   0.000   0.027   0.190   -0.097  0.120   0.361   0.380
    0.150   0.411   0.578   0.038   -2.106  0.086   -0.830  0.027   -2.181  0.300   -2.436  0.361   -0.001  0.195   0.000   0.000   0.043   0.189   -0.079  0.120   0.359   0.379
    0.200   0.416   0.594   0.030   -1.992  0.056   -0.833  0.028   -2.110  0.193   -2.312  0.380   0.000   0.203   0.000   0.000   0.067   0.219   -0.071  0.116   0.349   0.368
    0.250   0.421   0.628   0.022   -1.909  0.065   -0.938  -0.019  -1.918  0.238   -2.212  0.313   0.016   0.193   0.000   0.000   0.073   0.209   -0.028  0.113   0.338   0.356
    0.300   0.380   0.648   0.027   -1.868  0.054   -0.997  -0.034  -1.860  0.192   -2.079  0.335   0.020   0.187   0.000   0.000   0.082   0.227   -0.014  0.110   0.330   0.347
    0.350   0.325   0.664   0.032   -1.809  0.043   -1.051  -0.031  -1.831  0.142   -2.003  0.313   0.026   0.188   0.000   0.000   0.088   0.249   0.016   0.108   0.323   0.340
    0.400   0.274   0.682   0.030   -1.745  0.039   -1.119  -0.042  -1.845  0.077   -1.887  0.263   0.012   0.185   0.000   0.000   0.094   0.265   0.049   0.105   0.315   0.332
    0.450   0.232   0.702   0.028   -1.713  0.031   -1.089  -0.064  -1.843  0.063   -1.764  0.247   0.006   0.178   0.000   0.000   0.084   0.263   0.077   0.104   0.311   0.328
    0.500   0.173   0.718   0.028   -1.668  0.027   -1.068  -0.073  -1.835  0.062   -1.689  0.249   0.010   0.177   0.000   0.000   0.092   0.278   0.099   0.103   0.309   0.325
    0.600   0.073   0.746   0.019   -1.639  -0.009  -0.998  -0.052  -1.797  0.073   -1.513  0.215   0.004   0.159   0.000   0.000   0.103   0.287   0.128   0.103   0.308   0.325
    0.700   -0.001  0.771   0.012   -1.602  -0.022  -0.990  -0.060  -1.777  0.045   -1.395  0.169   -0.013  0.140   0.000   0.000   0.113   0.300   0.126   0.102   0.305   0.322
    0.800   -0.062  0.792   0.001   -1.554  -0.006  -1.022  -0.037  -1.736  0.073   -1.329  0.178   -0.024  0.131   0.000   0.000   0.110   0.299   0.126   0.102   0.306   0.323
    0.900   -0.126  0.818   -0.003  -1.521  0.008   -1.010  -0.041  -1.679  0.112   -1.228  0.217   -0.023  0.116   0.000   0.000   0.105   0.300   0.133   0.101   0.304   0.321
    1.000   -0.176  0.844   -0.009  -1.490  0.025   -1.028  -0.051  -1.622  0.140   -1.197  0.175   -0.017  0.109   0.000   0.000   0.097   0.297   0.136   0.101   0.304   0.320
    1.200   -0.278  0.896   -0.008  -1.444  0.073   -1.053  -0.071  -1.526  0.214   -1.175  0.092   -0.008  0.102   0.000   0.000   0.092   0.296   0.131   0.101   0.303   0.320
    1.400   -0.345  0.928   -0.016  -1.417  0.067   -1.141  0.001   -1.437  0.212   -1.255  0.082   -0.014  0.086   0.000   0.000   0.087   0.300   0.118   0.101   0.302   0.319
    1.600   -0.432  0.966   -0.008  -1.402  0.062   -1.211  -0.007  -1.357  0.290   -1.268  0.033   -0.011  0.082   0.000   0.000   0.090   0.309   0.114   0.101   0.302   0.318
    1.800   -0.510  0.999   0.000   -1.393  0.055   -1.235  0.005   -1.339  0.283   -1.300  -0.025  -0.011  0.078   0.000   0.000   0.085   0.309   0.117   0.101   0.304   0.320
    2.000   -0.581  1.020   0.003   -1.380  0.053   -1.255  0.000   -1.350  0.250   -1.324  0.010   -0.017  0.069   0.000   0.000   0.085   0.313   0.116   0.101   0.303   0.319
    2.500   -0.766  1.063   0.015   -1.340  0.079   -1.324  -0.043  -1.315  0.286   -1.365  -0.002  0.032   0.063   0.000   0.000   0.087   0.332   0.118   0.102   0.305   0.322
    3.000   -0.903  1.086   0.023   -1.322  0.096   -1.377  -0.070  -1.242  0.340   -1.430  0.025   0.065   0.053   0.000   0.000   0.091   0.329   0.117   0.102   0.307   0.324
    4.000   -1.105  1.161   0.062   -1.243  0.188   -1.475  -0.209  -1.169  0.418   -1.596  -0.153  0.065   0.027   0.000   0.000   0.098   0.324   0.136   0.101   0.302   0.318
    pga     0.053   0.619   0.023   -2.036  0.150   -0.949  -0.006  -2.129  0.383   -2.257  0.455   0.011   0.192   0.000   0.000   0.059   0.212   -0.059  0.109   0.327   0.344
    pgv     -1.160  0.788   0.033   -1.792  0.124   -1.050  -0.116  -1.868  0.318   -1.740  0.314   0.003   0.155   0.000   0.000   0.092   0.268   0.034   0.099   0.296   0.312
    """)