Source code for openquake.hazardlib.gsim.sgobba_2020

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2020, 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 :mod:`openquake.hazardlib.gsim.sgobba_2020` implements
:class:`~openquake.hazardlib.gsim.sgobba_2020.SgobbaEtAl2020`
"""

import os
import copy
import numpy as np
import pandas as pd
from scipy.constants import g as gravity_acc
from scipy.spatial import cKDTree
from openquake.hazardlib.geo import Point, Polygon
from openquake.hazardlib import const
from openquake.hazardlib.geo.mesh import Mesh
from openquake.hazardlib.imt import PGA, SA
from openquake.hazardlib.gsim.base import GMPE, CoeffsTable

# From http://www.csgnetwork.com/degreelenllavcalc.html
LEN_1_DEG_LAT_AT_43pt5 = 111.10245
LEN_1_DEG_LON_AT_43pt5 = 80.87665

DATA_FOLDER = os.path.join(os.path.dirname(__file__), 'sgobba_2020')

REGIONS = {
    '1': [[13.37, 42.13], [13.60, 42.24], [13.48, 42.51], [13.19, 42.36]],
    '4': [[13.26, 42.41], [13.43, 42.49], [13.27, 43.02], [12.96, 42.86]],
    '5': [[13.03, 42.90], [13.21, 42.99], [13.10, 43.13], [12.90, 43.06]]}

CONSTS = {'Mh': 5.0,
          'Rref': 1.0,
          'PseudoDepth': 6.0}


def _get_cluster_correction(dat, C, ctx, imt):
    """
    Get cluster correction. The use can specify various options through
    the cluster parameter. The available options are:
    - cluster = None
        In this case the code finds the most appropriate correction using
        the rupture position
    - cluster = 0
        No cluster correction
    - cluser = 1 or 4 or 5
        The code uses the correction for the given cluster
    """
    cluster = dat.cluster
    shape = ctx.sids.shape
    correction = np.zeros_like(shape)
    # st.dev.
    tau_L2L = np.zeros(shape)
    Bp_model = np.zeros(shape)
    phi_P2P = np.zeros(shape)

    # No cluster correction
    if cluster == 0:
        tau_L2L = C['tau_L2L']
        phi_P2P = C['phi_P2P']
        return correction, tau_L2L, Bp_model, phi_P2P
    # the code finds the most appropriate correction
    if cluster is None:
        mesh = Mesh(np.array([ctx.hypo_lon]), np.array([ctx.hypo_lat]))
        # midp = ctx.surface.get_middle_point()
        # mesh = Mesh(np.array([midp.longitude]),np.array([midp.latitude]))

        for key in REGIONS:
            coo = np.array(REGIONS[key])
            pnts = [Point(lo, la) for lo, la in zip(coo[:, 0], coo[:, 1])]
            poly = Polygon(pnts)
            within = poly.intersects(mesh)
            if within.all():
                cluster = int(key)
                break
    # if OUT clusters do not apply corrections
    if cluster is None:
        tau_L2L = C['tau_L2L']
        phi_P2P = C['phi_P2P']
        return correction, tau_L2L, Bp_model, phi_P2P
    else:
        # if IN clusters apply corrections
        # Cluster coefficients
        fname = 'P_model_cluster{:d}.csv'.format(cluster)
        fname = os.path.join(DATA_FOLDER, fname)
        data = np.loadtxt(fname, delimiter=",", skiprows=1)
        # for st.dev.
        fname2 = 'beta_dP2P_cluster{:d}.csv'.format(cluster)
        fname2 = os.path.join(DATA_FOLDER, fname2)
        data2 = np.loadtxt(fname2, delimiter=",", skiprows=1)
        # Compute the coefficients
        correction = np.zeros(shape)
        per = imt.period
        for idx in np.unique(dat.idxs):
            tmp = data[int(idx)]
            correction[dat.idxs == idx] = np.interp(
                per, dat.PERIODS, tmp[0:5])
        # Adding L2L correction
        label = "dL2L_cluster{:d}".format(cluster)
        correction += C[label]
        # compute st.dev.
        for idx in np.unique(dat.idxs):
            tmp2 = data2[int(idx)]
            Bp_model[dat.idxs == idx] = np.interp(per, dat.PERIODS, tmp2[0:5])
        return correction, tau_L2L, Bp_model, phi_P2P


def _get_distance_term(C, mag, ctx):
    """
    Eq.3 - page 3
    """
    term1 = C['c1'] * (mag - C['mref']) + C['c2']
    tmp = np.sqrt(ctx.rjb**2 + CONSTS['PseudoDepth']**2)
    term2 = np.log10(tmp / CONSTS['Rref'])
    term3 = C['c3'] * (tmp - CONSTS['Rref'])
    return term1 * term2 + term3


def _get_magnitude_term(C, mag):
    """
    Eq.2 - page 3
    """
    return np.where(mag <= CONSTS['Mh'],
                    C['b1'] * (mag - CONSTS['Mh']),
                    C['b2'] * (mag - CONSTS['Mh']))


def _get_site_correction(data, shape, imt):
    """
    Get site correction
    """
    correction = np.zeros_like(shape)
    # Compute the coefficients
    correction = np.zeros(shape)
    # stand.dev.
    Bs_model = np.zeros(shape)
    phi_S2Sref = np.zeros(shape)
    per = imt.period
    for idx in np.unique(data.idxs):
        tmp = data.Smodel[int(idx)]
        correction[data.idxs == idx] = np.interp(
            per, data.PERIODS, tmp[0:5])
        tmp2 = data.betaS2S[int(idx)]
        Bs_model[data.idxs == idx] = np.interp(
            per, data.PERIODS, tmp2[0:5])
    return correction, Bs_model, phi_S2Sref


[docs]class Data(object):
    """Helper class"""
    def __init__(self, smodel, cluster, periods, betaS2S, idxs):
        self.Smodel = smodel
        self.cluster = copy.copy(cluster)
        self.PERIODS = periods
        self.betaS2S = betaS2S
        self.idxs = idxs


# NB: the implementation here is HORRIBLE performance-wise,
# because it is using KDTree and pandas!
[docs]class SgobbaEtAl2020(GMPE):
    """
    Implements the GMM proposed by Sgobba et al. (2020).
    Warning:
    This GMM is not meant for national models where it would be too slow,
    it is meant for scenario calculations.

    :param event_id:
        A string identifying an event amongst the ones comprised in the
        list available in the file `event.csv`
    :param directionality:
        A boolean
    :param cluster:
        If set to 'None', the OQ Engine finds the corresponding cluster
        using the rupture epicentral location.
        If cluster=0, no cluster correction applied.
        Otherwise, if an integer ID is provided, that
        corresponds to the cluster id (available cluster indexes are 1, 4
        and 5), the corresponding correction id applied.
    """
    #: Supported tectonic region type is 'active shallow crust'
    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, 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 parameters are lon and lat
    REQUIRES_SITES_PARAMETERS = {'lon', 'lat'}

    #: Required rupture parameters is magnitude, hypo_lat, hypo_lon
    REQUIRES_RUPTURE_PARAMETERS = {'mag', 'hypo_lon', 'hypo_lat'}

    #: Required distance measure is Rjb
    REQUIRES_DISTANCES = {'rjb'}

    PERIODS = np.array([0, 0.2, 0.50251256281407, 1.0, 2.0])

    def __init__(self, event_id=None, directionality=False, cluster=None,
                 site=False, bedrock=False):
        self.event_id = event_id
        self.directionality = directionality
        self.cluster = cluster
        self.site = site
        self.bedrock = bedrock
        # Get site indexes. They are used for the site correction and the
        # cluster (path) correction
        # Load the coordinate of the grid
        fname = os.path.join(DATA_FOLDER, 'grid.csv')
        coo = np.fliplr(np.loadtxt(fname, delimiter=";"))
        self.kdt = cKDTree(coo)

        # Load the table with the between event coefficients
        if event_id is not None:
            self.event_id = event_id
            fname = os.path.join(DATA_FOLDER, 'event.csv')
            # Create dataframe with events
            df = pd.read_csv(fname, dtype={'id': str})
            df.set_index('id', inplace=True)
            self.df = df
            assert event_id in df.index.values
        # Site coefficients
        fname = os.path.join(DATA_FOLDER, "S_model.csv")
        self.Smodel = np.loadtxt(fname, delimiter=",", skiprows=1)
        fname = os.path.join(DATA_FOLDER, "beta_dS2S.csv")
        self.betaS2S = np.loadtxt(fname, delimiter=",", skiprows=1)

[docs]    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        Eq.1 - page 2
        """
        for m, imt in enumerate(imts):
            # Ergodic coeffs
            C = self.COEFFS[imt]
            # between-event
            if self.event_id is not None:
                label = "dBe_{:s}".format(str(imt))
                self.be = self.df.loc[self.event_id][label]
                self.be_std = 0.0
            else:
                self.be_std = C['tau_ev']
                self.be = 0.0
            # Site correction
            points = np.array([ctx.lon, ctx.lat]).T  # shape (N, 2)
            dsts, idxs = self.kdt.query(points)
            dat = Data(self.Smodel, self.cluster, self.PERIODS,
                       self.betaS2S, idxs)
            sc = 0
            phi_S2Sref = C['phi_S2S_ref']
            Bs_model = np.zeros(ctx.sids.shape)
            if self.site and self.bedrock is False:
                sc, Bs_model, phi_S2Sref = _get_site_correction(
                    dat, ctx.sids.shape, imt)

            cc, tau_L2L, Bp_model, phi_P2P = _get_cluster_correction(
                dat, C, ctx, imt)
            # Get mean
            mean[m] = (C['a'] + _get_magnitude_term(C, ctx.mag) +
                       _get_distance_term(C, ctx.mag, ctx) +
                       sc + cc + self.be)
            # To natural logarithm and fraction of g
            mean[m] = np.log(10.0**mean[m] / (gravity_acc*100))
            # Get stds
            std = np.sqrt(C['sigma_0'] ** 2 + self.be_std ** 2 + tau_L2L ** 2 +
                          Bs_model + phi_S2Sref ** 2 + Bp_model + phi_P2P ** 2)
            sig[m] = np.log(10.0 ** std)
            tau[m] = np.log(10.0 ** self.be_std)
            phi[m] = np.log(10.0 ** np.sqrt((std ** 2 - self.be_std ** 2)))

    COEFFS = CoeffsTable(sa_damping=5., table="""\
    IMT                a                   b1                  b2                   c1                   c2                    c3                     mref               tau_ev              tau_L2L               phi_S2S_ref   phi_S2S              phi_P2P             sigma_0            dL2L_cluster1          dL2L_cluster4         dL2L_cluster5
    pga               2.92178299969904    0.549352522898805   0.195787600661646    0.182324348626393    -1.56833817017883     -0.00277072348000775   3.81278167434967   0.148487352282019   0.0244729216674648    0.2011        0.249446934633114    0.129067143622624   0.194457660872416  -0.0144001959308384    -0.0141006684390188   -0.0814912207894038
    0.2               3.23371734550753    0.718110435825521   0.330819511910566    0.101391376086178    -1.47499081134392     -0.00235944669544279   3.52085298608413   0.142257128473462   0.0142903160948886    0.2049        0.25356556670542     0.108002263739343   0.211634711188527  -0.0295313493684611    -0.0242995747709838   -0.0779761977475732
    0.50251256281407  3.16050217205595    0.838494386998919   0.466787811642044    0.105723089676       -1.48056328666322      0                     4.87194107479204   0.118125529818761   0.0131546256192079    0.1419        0.230442366433058    0.0959080118602377  0.191012833298469   0.00929098048437728   -0.00995372305434456  -0.00828100722167989
    1                 2.58227846237728    0.85911311807545    0.519131261495525    0.146088352194266    -1.28019118368202      0                     5.42555199253122   0.124229747688977   2.39299038437967e-08  0.1153        0.212408309867869    0.118568732468557   0.176037658051544   0.000737173026444824  -0.00123578210338215   0.000181351036566464
    2                 1.88792168738756    0.727248116061721   0.47362977053987     0.244695132922949    -1.19816952711971      0                     5.26896508895249   0.127711124129548   8.69064652723658e-09  0.1078        0.189154038588083    0.119572905421336   0.183045950697286   5.60984632803441e-15  -1.18288330352055e-14  9.31778101896791e-15
    """)