# -*- 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:`BergeThierryEtAl2003SIGMA`.
"""
import copy
import numpy as np
from scipy.constants import g
from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, SA
[docs]class BergeThierryEtAl2003SIGMA(GMPE):
"""
Implements GMPE developed by Catherine Berge-Thierry, Fabrice Cotton,
Oona Scoti, Daphne-Anne Griot-Pommera, and Yoshimitsu Fukushima and
published as "New Empirical Response Spectral Attenuation Laws For Moderate
European Earthquakes" (2003, Journal of Earthquake Engineering, 193-222)
The class implements also adjustment of the sigma value as required by
the SIGMA project to make the GMPE usable with Mw (the GMPE was
originally developed for Ms).
"""
#: Supported tectonic region type is active shallow crust, see
#: `Introduction`, page 194.
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST
#: Supported intensity measure types are spectral acceleration, and peak
#: ground acceleration. The original manuscript provide coefficients only
#: SA. For PGA, coefficients are assumed equal to the ones of SA for the
#: smallest period (0.03 s)
DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
PGA,
SA
])
#: Supported intensity measure component is horizontal, see page 196.
DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.HORIZONTAL
#: Supported standard deviation type is total, see table 3, page 203
DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
const.StdDev.TOTAL
])
#: Required site parameters is Vs30, used to distinguish between rock sites
#: (Vs30 >= 800) m/s and alluvium sites (300 < Vs < 800), see section 2.2.3
#: page 201
REQUIRES_SITES_PARAMETERS = set(('vs30', ))
#: Required rupture parameters is magnitude, see equation 1 page 201
REQUIRES_RUPTURE_PARAMETERS = set(('mag', ))
#: Required distance measure is hypocentral distance, see equation 1 page
#: 201
REQUIRES_DISTANCES = set(('rhypo', ))
[docs] def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
return self._get_mean_and_stddevs(sites, rup, dists, imt,
stddev_types)
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.
"""
# extract dictionaries of coefficients specific to required
# intensity measure type
C = self.COEFFS[imt]
# clip distance at 4 km, minimum distance for which the equation is
# valid (see section 2.2.4, page 201). This also avoids singularity
# in the equation
rhypo = dists.rhypo
rhypo[rhypo < 4.] = 4.
mean = C['a'] * rup.mag + C['b'] * rhypo - np.log10(rhypo)
mean[sites.vs30 >= 800] += C['c1']
mean[sites.vs30 < 800] += C['c2']
# convert from log10 to ln, and from cm/s2 to g
mean = mean * np.log(10) - 2 * np.log(10) - np.log(g)
stddevs = self._get_stddevs(C, stddev_types, rhypo.shape[0])
return mean, stddevs
def _get_stddevs(self, C, stddev_types, num_sites):
"""
Return total standard deviation.
"""
assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
for stddev_type in stddev_types)
# adjustement for sigma, as explained in section 3.7.3 in
# Delivrable_SIGMA-V02-2012-D4-18+reviews
sigma = np.zeros(num_sites) + C['sigma'] * np.log(10)
sigma = np.sqrt(sigma ** 2 + (C['a'] ** 2) * (0.2 ** 2))
stddevs = [sigma for _ in stddev_types]
return stddevs
#: Coefficient tables are constructed from the electronic suplements of
#: the original paper. Original coefficients in function of frequency.
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT a b c1 c2 sigma
pga 0.3118000 -0.0009303 1.537000 1.57300 0.2923
0.029412 0.3118000 -0.0009303 1.537000 1.57300 0.2923
0.030000 0.3114000 -0.0009334 1.541000 1.57600 0.2924
0.032258 0.3097000 -0.0009422 1.558000 1.58900 0.2928
0.034000 0.3083000 -0.0009547 1.573000 1.60200 0.2935
0.035714 0.3068000 -0.0009822 1.593000 1.61800 0.2947
0.040000 0.3033000 -0.0011190 1.653000 1.66500 0.2982
0.045455 0.3016000 -0.0012820 1.701000 1.70000 0.3015
0.050000 0.2992000 -0.0013410 1.740000 1.72900 0.3009
0.053000 0.2981000 -0.0014290 1.766000 1.74900 0.3022
0.055556 0.2969000 -0.0014320 1.785000 1.76500 0.3027
0.058824 0.2960000 -0.0014600 1.809000 1.78400 0.3040
0.059999 0.2960000 -0.0014720 1.814000 1.78800 0.3041
0.062500 0.2965000 -0.0015070 1.826000 1.79600 0.3047
0.064998 0.2944000 -0.0015150 1.851000 1.81700 0.3059
0.066667 0.2933000 -0.0015130 1.865000 1.82900 0.3059
0.068966 0.2924000 -0.0015460 1.881000 1.84200 0.3048
0.069999 0.2920000 -0.0015700 1.888000 1.84900 0.3046
0.071429 0.2909000 -0.0015830 1.901000 1.86100 0.3047
0.074074 0.2879000 -0.0015550 1.926000 1.88700 0.3059
0.075002 0.2871000 -0.0015400 1.933000 1.89300 0.3066
0.076923 0.2857000 -0.0015200 1.947000 1.90600 0.3072
0.080000 0.2866000 -0.0015310 1.954000 1.91100 0.3061
0.083333 0.2858000 -0.0015520 1.968000 1.92700 0.3058
0.084998 0.2856000 -0.0015830 1.976000 1.93500 0.3053
0.086957 0.2848000 -0.0015710 1.987000 1.94600 0.3042
0.090001 0.2819000 -0.0015360 2.011000 1.97300 0.3020
0.090909 0.2809000 -0.0015280 2.020000 1.98100 0.3016
0.095238 0.2781000 -0.0015430 2.054000 2.01000 0.3017
0.100000 0.2786000 -0.0014740 2.059000 2.01600 0.3019
0.105260 0.2776000 -0.0014220 2.072000 2.03400 0.3043
0.110000 0.2783000 -0.0014420 2.075000 2.04500 0.3073
0.111110 0.2783000 -0.0014380 2.077000 2.04700 0.3079
0.117650 0.2793000 -0.0014380 2.081000 2.05600 0.3103
0.120000 0.2806000 -0.0014360 2.076000 2.05400 0.3106
0.125000 0.2831000 -0.0014150 2.066000 2.05100 0.3121
0.129030 0.2863000 -0.0014400 2.052000 2.04200 0.3137
0.130010 0.2867000 -0.0014250 2.050000 2.04000 0.3142
0.133330 0.2887000 -0.0013970 2.040000 2.03300 0.3156
0.137930 0.2903000 -0.0013380 2.032000 2.02800 0.3161
0.140000 0.2915000 -0.0013220 2.027000 2.02200 0.3166
0.142860 0.2933000 -0.0013070 2.018000 2.01400 0.3173
0.148150 0.2950000 -0.0012560 2.009000 2.01000 0.3191
0.149990 0.2955000 -0.0012180 2.004000 2.00900 0.3196
0.153850 0.2955000 -0.0010520 1.997000 2.00700 0.3201
0.160000 0.2939000 -0.0008056 1.996000 2.01300 0.3200
0.166670 0.2952000 -0.0007097 1.989000 2.00600 0.3215
0.170010 0.2974000 -0.0006986 1.978000 1.99400 0.3230
0.173910 0.3016000 -0.0007341 1.957000 1.97000 0.3247
0.179990 0.3089000 -0.0007793 1.915000 1.93000 0.3267
0.181820 0.3109000 -0.0007826 1.901000 1.91900 0.3271
0.190010 0.3147000 -0.0007369 1.868000 1.89400 0.3262
0.190480 0.3149000 -0.0007337 1.867000 1.89300 0.3262
0.200000 0.3167000 -0.0006889 1.843000 1.88100 0.3250
0.208330 0.3196000 -0.0006719 1.814000 1.86100 0.3261
0.217390 0.3254000 -0.0006750 1.770000 1.82500 0.3281
0.219780 0.3271000 -0.0006918 1.758000 1.81500 0.3292
0.227270 0.3303000 -0.0006678 1.726000 1.79200 0.3320
0.238100 0.3340000 -0.0006171 1.683000 1.76200 0.3367
0.239980 0.3344000 -0.0005988 1.677000 1.75800 0.3371
0.250000 0.3365000 -0.0005750 1.651000 1.73600 0.3394
0.259740 0.3430000 -0.0007075 1.609000 1.69700 0.3422
0.263160 0.3442000 -0.0007200 1.599000 1.68800 0.3429
0.277780 0.3501000 -0.0007520 1.550000 1.64500 0.3444
0.280030 0.3511000 -0.0007530 1.542000 1.63800 0.3447
0.290020 0.3555000 -0.0007836 1.506000 1.60500 0.3458
0.300030 0.3590000 -0.0008520 1.477000 1.58100 0.3477
0.303030 0.3602000 -0.0008737 1.466000 1.57300 0.3483
0.316960 0.3671000 -0.0009272 1.412000 1.52500 0.3491
0.320000 0.3690000 -0.0009468 1.397000 1.51200 0.3487
0.333330 0.3742000 -0.0010100 1.352000 1.47200 0.3474
0.340020 0.3752000 -0.0010060 1.337000 1.46100 0.3469
0.344830 0.3760000 -0.0009698 1.326000 1.45200 0.3471
0.357140 0.3807000 -0.0009114 1.286000 1.41500 0.3481
0.359970 0.3822000 -0.0009039 1.275000 1.40500 0.3484
0.370370 0.3867000 -0.0008635 1.237000 1.37200 0.3492
0.379940 0.3909000 -0.0008074 1.199000 1.33900 0.3500
0.384620 0.3931000 -0.0007955 1.179000 1.32100 0.3507
0.400000 0.3997000 -0.0007078 1.119000 1.26700 0.3517
0.416670 0.4028000 -0.0006613 1.078000 1.23200 0.3512
0.419990 0.4034000 -0.0006513 1.070000 1.22600 0.3513
0.434780 0.4070000 -0.0006167 1.029000 1.19100 0.3521
0.439950 0.4089000 -0.0006118 1.011000 1.17400 0.3527
0.454550 0.4148000 -0.0005931 0.960100 1.12500 0.3547
0.459980 0.4165000 -0.0005816 0.944000 1.10900 0.3549
0.476190 0.4222000 -0.0005404 0.893100 1.05800 0.3555
0.480080 0.4239000 -0.0005484 0.879500 1.04500 0.3556
0.500000 0.4323000 -0.0005680 0.815000 0.97970 0.3555
0.520020 0.4372000 -0.0005396 0.764200 0.93240 0.3568
0.526320 0.4379000 -0.0005050 0.752200 0.92080 0.3570
0.539960 0.4394000 -0.0004330 0.727100 0.89590 0.3574
0.555560 0.4418000 -0.0003601 0.694100 0.86610 0.3587
0.559910 0.4425000 -0.0003380 0.684400 0.85710 0.3592
0.580050 0.4472000 -0.0002702 0.635700 0.81170 0.3610
0.588240 0.4492000 -0.0002522 0.615700 0.79260 0.3609
0.599880 0.4516000 -0.0002175 0.589500 0.76820 0.3603
0.619960 0.4559000 -0.0001953 0.547800 0.72820 0.3604
0.625000 0.4569000 -0.0001995 0.538300 0.71870 0.3609
0.640200 0.4596000 -0.0001666 0.510600 0.69100 0.3625
0.660070 0.4637000 -0.0001549 0.472500 0.65200 0.3647
0.666670 0.4655000 -0.0001550 0.457300 0.63630 0.3649
0.679810 0.4688000 -0.0001668 0.428400 0.60810 0.3655
0.699790 0.4732000 -0.0001700 0.385700 0.56760 0.3667
0.714290 0.4771000 -0.0002019 0.354800 0.53540 0.3680
0.750190 0.4847000 -0.0003009 0.287100 0.46810 0.3700
0.769230 0.4875000 -0.0003122 0.254500 0.43800 0.3703
0.800000 0.4940000 -0.0002568 0.190600 0.37820 0.3714
0.833330 0.5010000 -0.0001932 0.126400 0.31450 0.3746
0.850340 0.5040000 -0.0001433 0.096150 0.28420 0.3758
0.900090 0.5098000 3.28E-05 0.020060 0.21300 0.3747
0.909090 0.5104000 8.39E-05 0.008195 0.20220 0.3744
1.000000 0.5199000 0.0002516 -0.116200 0.08290 0.3737
1.100100 0.5273000 0.0003908 -0.212300 -0.02900 0.3794
1.111100 0.5278000 0.0004074 -0.220700 -0.03875 0.3804
1.200500 0.5361000 0.0004479 -0.313300 -0.13380 0.3838
1.250000 0.5409000 0.0004860 -0.367900 -0.18910 0.3877
1.300400 0.5444000 0.0005329 -0.411300 -0.23770 0.3920
1.400600 0.5481000 0.0007676 -0.487000 -0.31550 0.3935
1.428600 0.5494000 0.0008272 -0.509500 -0.33850 0.3941
1.499300 0.5527000 0.0009124 -0.560400 -0.39350 0.3932
1.600000 0.5557000 0.0009844 -0.618600 -0.45830 0.3919
1.666700 0.5580000 0.0010850 -0.656400 -0.50260 0.3919
1.798600 0.5620000 0.0012450 -0.725800 -0.58340 0.3942
2.000000 0.5622000 0.0013750 -0.796300 -0.66600 0.4030
2.197800 0.5617000 0.0016520 -0.865600 -0.73950 0.4055
2.398100 0.5641000 0.0018290 -0.940600 -0.81580 0.4093
2.500000 0.5654000 0.0019210 -0.978700 -0.85420 0.4110
2.597400 0.5677000 0.0020060 -1.019000 -0.89800 0.4130
2.801100 0.5666000 0.0022770 -1.071000 -0.94950 0.4202
3.003000 0.5683000 0.0024490 -1.130000 -1.01400 0.4255
3.205100 0.5686000 0.0025360 -1.179000 -1.06900 0.4301
3.333300 0.5705000 0.0025330 -1.220000 -1.11100 0.4329
3.401400 0.5715000 0.0025410 -1.243000 -1.13500 0.4340
3.597100 0.5727000 0.0025730 -1.300000 -1.19400 0.4359
3.802300 0.5712000 0.0026620 -1.350000 -1.24200 0.4365
4.000000 0.5722000 0.0027110 -1.417000 -1.30300 0.4344
4.504500 0.5856000 0.0024490 -1.662000 -1.52000 0.4278
5.000000 0.5990000 0.0021050 -1.886000 -1.72900 0.4233
5.494500 0.6106000 0.0019410 -2.072000 -1.90400 0.4260
5.988000 0.6160000 0.0018800 -2.201000 -2.02800 0.4284
6.993000 0.6175000 0.0017660 -2.373000 -2.19200 0.4299
8.000000 0.6145000 0.0017210 -2.491000 -2.30800 0.4284
9.009000 0.6122000 0.0016370 -2.592000 -2.40800 0.4236
10.00000 0.6086000 0.0015630 -2.668000 -2.48500 0.4183
""")
[docs]class BergeThierryEtAl2003SIGMAMwW(BergeThierryEtAl2003SIGMA):
"""
Mw version of the Berge-Thierry et al. (2003) GMPE. For this conversion
we use the Weatherill et al. (2016) conversion equation between Ms and Mw
"""
[docs] def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
newrup = copy.copy(rup)
if rup.mag <= 6.064:
newrup.mag = (rup.mag - 2.369) / 0.616
else:
newrup.mag = (rup.mag - 0.100) / 0.994
return self._get_mean_and_stddevs(sites, newrup, dists, imt,
stddev_types)