Skip to content

Commit

Permalink
Merge pull request #8957 from gem/rtgm
Browse files Browse the repository at this point in the history
Initial work for RTGM
  • Loading branch information
micheles authored Aug 28, 2023
2 parents c69e7e9 + 27469f1 commit 7414043
Show file tree
Hide file tree
Showing 13 changed files with 258 additions and 61 deletions.
193 changes: 193 additions & 0 deletions openquake/calculators/postproc/compute_rtgm.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,193 @@
# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2023, 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 to compute the Risk Targeted Ground Motion by using the
rtgmpy module from the USGS.
Useful abbreviations:
- RTGM: Risk Targeted Ground Motion (for geometric mean)
- RTGM_max: Risk Targeted Ground Motion (for maximum component)
- MCE: Maximum Considered Earthquake
- ProbMCE: Probabilistic Maximum Considered Earthquake (UHGM for PGA)
- IMT: the list of IMTs as normalized strings
- UHGM: Uniform Hazard Spectra
- RiskCoeff: RTGM / UHGM_2475y
- DLL: Deterministic Lower Limit
T Site Class A B BC C CD D DE E
0.00 0.50 0.57 0.66 0.73 0.74 0.69 0.61 0.55
0.01 0.50 0.57 0.66 0.73 0.75 0.70 0.62 0.55
0.02 0.52 0.58 0.68 0.74 0.75 0.70 0.62 0.55
0.03 0.60 0.66 0.75 0.79 0.78 0.70 0.62 0.55
0.05 0.81 0.89 0.95 0.96 0.89 0.76 0.62 0.55
0.075 1.04 1.14 1.21 1.19 1.08 0.90 0.71 0.62
0.10 1.12 1.25 1.37 1.37 1.24 1.04 0.82 0.72
0.15 1.12 1.29 1.53 1.61 1.50 1.27 1.00 0.87
0.20 1.01 1.19 1.50 1.71 1.66 1.44 1.15 1.01
0.25 0.90 1.07 1.40 1.71 1.77 1.58 1.30 1.15
0.30 0.81 0.98 1.30 1.66 1.83 1.71 1.44 1.30
0.40 0.69 0.83 1.14 1.53 1.82 1.80 1.61 1.48
0.50 0.60 0.72 1.01 1.38 1.73 1.80 1.68 1.60
0.75 0.46 0.54 0.76 1.07 1.41 1.57 1.60 1.59
1.0 0.37 0.42 0.60 0.86 1.17 1.39 1.51 1.58
1.5 0.26 0.29 0.41 0.60 0.84 1.09 1.35 1.54
2.0 0.21 0.23 0.31 0.45 0.64 0.88 1.19 1.46
3.0 0.15 0.17 0.21 0.31 0.45 0.63 0.89 1.11
4.0 0.12 0.13 0.16 0.24 0.34 0.47 0.66 0.81
5.0 0.10 0.11 0.13 0.19 0.26 0.36 0.49 0.61
7.5 0.063 0.068 0.080 0.11 0.15 0.19 0.26 0.31
10 0.042 0.045 0.052 0.069 0.089 0.11 0.14 0.17
PGAG 0.37 0.43 0.50 0.55 0.56 0.53 0.46 0.42
# PGA_G: PGA for Geometric Mean (no Risk Targeted)
# PGA: PGA for Maximum Component
"""
import logging
import numpy as np
import pandas as pd
from scipy import interpolate
try:
import rtgmpy as rtg
except ImportError:
rtg = None
from openquake.hazardlib.calc.mean_rates import to_rates


def norm_imt(imt):
"""
Normalize the imt string to the USGS format, for instance SA(1.1) -> SA1P1
"""
return imt.replace('(', '').replace(')', '').replace('.', '')

# hard-coded for year 1
IMTs = [norm_imt(im) for im in ['PGA', 'SA(0.2)', 'SA(1.0)']]
Ts = [0, 0.2, 1.0]
LIMITs = np.array([0.5, 1.5, 0.6])


def _find_fact_maxC(T,code):
# find the factor to convert to maximum component based on
# ASCE7-16 and ASCE7-22
if code == 'ASCE7-16':
if T == 0:
fact_maxC = 1.
elif T <= 0.2:
fact_maxC = 1.1
elif T <= 1:
f = interpolate.interp1d([0.2,1],[1.1, 1.3])
fact_maxC = f(T)
elif T <= 5:
f = interpolate.interp1d([1,5],[1.3, 1.5])
fact_maxC = f(T)
else:
fact_maxC = 1.5
elif code == 'ASCE7-22':
if T == 0:
fact_maxC = 1.
elif T <= 0.2:
fact_maxC = 1.2
elif T <= 1:
f = interpolate.interp1d([0.2,1],[1.2, 1.25])
fact_maxC = f(T)
elif T <= 10:
f = interpolate.interp1d([1,5],[1.25, 1.3])
fact_maxC = f(T)
else:
fact_maxC = 1.5
return fact_maxC


def calc_rtgm_df(rtgm_haz, oq):
# Obtaining Risk-Targeted Ground Motions from hazard curves
M = len(IMTs)
export_rtgm = []
riskCoeff = np.zeros(M)
RTGM, UHGM, RTGM_max, MCE = (np.zeros(M), np.zeros(M),
np.zeros(M), np.zeros(M))
results = rtg.BuildingCodeRTGMCalc.calc_rtgm(rtgm_haz, 'ASCE7')
for m, IMT in enumerate(IMTs):
rtgmCalc = results['RTGM'][IMT]['rtgmCalc']
T = Ts[m]
fact = _find_fact_maxC(T, 'ASCE7-16')
RTGM[m] = rtgmCalc['rtgm'] / fact
RTGM_max[m] = rtgmCalc['rtgm']
UHGM[m] = rtgmCalc['uhgm'] / fact # back to GM
riskCoeff[m] = rtgmCalc['riskCoeff']
# note that RTGM_max is the ProbMCEr, while RTGM is used for the
# identification of the sources as the hazard curves are in
# geometric mean
if IMT == 'PGA':
RTGM[m] = UHGM[m]
MCE[m] = RTGM[m] # UHGM in terms of GM: MCEg
else:
MCE[m] = RTGM_max[m]
# this is saved for the next step.
export_rtgm.append(IMT + ': ' + str(RTGM[m]))
print(export_rtgm)
if (MCE < LIMITs).all():
dic = {'IMT': IMTs,
'UHGM_2475yr-GM': UHGM,
'RTGM': RTGM_max,
'ProbMCE': MCE,
'RiskCoeff': riskCoeff,
'DLL': LIMITs,
'MCE>DLL?': RTGM_max > LIMITs,
'GoverningMCE': MCE}
else:
import pdb; pdb.set_trace()
raise NotImplementedError
return pd.DataFrame(dic)


def get_hazdic(hcurves, imtls, invtime, sitecol):
"""
Convert an array of mean hazard curves into a dictionary suitable
for the rtgmpy library
:param hcurves: array of PoEs of shape (N, M, L1)
"""
[site] = sitecol # there must be a single site
hazdic = {
'site': {'name': 'site',
'lon': site.location.x,
'lat': site.location.y,
'Vs30': site.vs30},
'hazCurves': {norm_imt(imt): {'iml': imtls[imt],
'afe': to_rates(hcurves[0, m], invtime)}
for m, imt in enumerate(imtls) if norm_imt(imt) in IMTs}}
return hazdic


def main(dstore):
"""
:param dstore: datastore with the classical calculation
"""
if not rtg:
logging.warning('Missing module rtgmpy: skipping AELO calculation')
return
logging.info('Computing Risk Targeted Ground Motion')
oq = dstore['oqparam']
assert list(oq.hazard_stats()) == ['mean'], oq.hazard_stats()
sitecol = dstore['sitecol']
hcurves = dstore['hcurves-stats'][:, 0] # shape NML1
hazdic = get_hazdic(hcurves, oq.imtls, oq.investigation_time, sitecol)
rtgm_haz = rtg.GroundMotionHazard.from_dict(hazdic)
df = calc_rtgm_df(rtgm_haz, oq)
dstore.create_df('rtgm', df)
4 changes: 2 additions & 2 deletions openquake/calculators/tests/logictree_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -332,7 +332,7 @@ def test_case_19(self):
mean_rates = to_rates(mean_poes)
rates_by_source = self.calc.datastore[
'mean_rates_by_src'][0] # (M, L1, Ns)
aac(mean_rates, rates_by_source.sum(axis=2), atol=2E-7)
aac(mean_rates, rates_by_source.sum(axis=2), atol=5E-7)

def test_case_20(self):
# Source geometry enumeration, apply_to_sources
Expand Down Expand Up @@ -389,7 +389,7 @@ def test_case_20_bis(self):
self.assertEqual(aw.imt, 'PGA')
self.assertEqual(aw.poe, .001)
# the numbers are quite different on macOS, 6.461143e-05 :-(
aac(aw.array['poe'], [6.467104e-05, 0, 0], atol=1E-7)
aac(aw.array['poe'], [6.467104e-05, 0, 0], atol=1E-6)

# testing view_relevant_sources
arr = view('relevant_sources:SA(1.0)', self.calc.datastore)
Expand Down
1 change: 1 addition & 0 deletions openquake/engine/aelo.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ def get_params_from(inputs):
if 'vs30' in inputs:
params['override_vs30'] = '%(vs30)s' % inputs
params['postproc_func'] = 'disagg_by_rel_sources.main'
# params['postproc_func'] = 'compute_rtgm.main'
# params['cachedir'] = datastore.get_datadir()
return params

Expand Down
9 changes: 4 additions & 5 deletions openquake/hazardlib/calc/mean_rates.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,26 +24,25 @@
from openquake.hazardlib.contexts import get_cmakers


def to_rates(probs):
def to_rates(probs, itime=1):
"""
Convert an array of probabilities into an array of rates
>>> numpy.round(to_rates(numpy.array([.8])), 6)
array([1.609438])
"""
pnes = 1. - probs
pnes[pnes == 0] = 1E-45 # minimum float32
return - numpy.log(pnes)
return numpy.clip(- numpy.log(pnes) / itime, 1E-16, 100)


def to_probs(rates):
def to_probs(rates, itime=1):
"""
Convert an array of rates into an array of probabilities
>>> numpy.round(to_probs(numpy.array([1.609438])), 6)
array([0.8])
"""
return 1. - numpy.exp(- rates)
return 1. - numpy.exp(- rates * itime)


def calc_rmap(src_groups, full_lt, sitecol, oq):
Expand Down
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
| poe | imt | distbin | prob |
|---------+-----+---------+---------|
| 0.00211 | PGA | 0 | 0.0 |
| 0.00211 | PGA | 1 | 0.00211 |
| poe | imt | distbin | prob |
|---------+-----+---------+-----------|
| 0.00211 | PGA | 0 | 1.000E-07 |
| 0.00211 | PGA | 1 | 0.00211 |
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
#,,,"lon=36.251, lat=-6.49"
#,,,"generated_by='OpenQuake engine 3.18.0-gitd8caa064d4', start_date='2023-08-28T11:51:11', checksum=2124969217, lon=36.251, lat=-6.49"
src_id,imt,iml,value
1500,PGA,5.00000E-03,6.66609E+00
1500,PGA,7.00146E-03,4.62962E+00
Expand All @@ -17,6 +17,9 @@ src_id,imt,iml,value
1500,PGA,5.57223E-01,4.82794E-05
1500,PGA,7.80274E-01,5.86798E-06
1500,PGA,1.09261E+00,3.09584E-07
1500,PGA,1.52997E+00,1.00000E-07
1500,PGA,2.14241E+00,1.00000E-07
1500,PGA,3.00000E+00,1.00000E-07
1501,PGA,5.00000E-03,6.95542E+00
1501,PGA,7.00146E-03,4.87481E+00
1501,PGA,9.80408E-03,3.28428E+00
Expand All @@ -36,3 +39,4 @@ src_id,imt,iml,value
1501,PGA,1.09261E+00,4.52269E-05
1501,PGA,1.52997E+00,6.58398E-06
1501,PGA,2.14241E+00,2.46933E-07
1501,PGA,3.00000E+00,1.00000E-07
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
#,,,,,,,"generated_by='OpenQuake engine 3.17.0-git5adb6553aa', start_date='2023-05-05T11:16:59', checksum=3344811530, investigation_time=1.0, mag_bin_edges=[5.5, 5.6000000000000005, 5.7, 5.800000000000001, 5.9, 6.0, 6.1000000000000005, 6.2, 6.300000000000001, 6.4, 6.5, 6.6000000000000005], dist_bin_edges=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0], lon_bin_edges=[-3.7204433425622865, -3.540332506921715, -3.3602216712811437, -3.180110835640572, -3.0, -2.8198891643594286, -2.639778328718857, -2.4596674930782854, -2.2795566574377135], lat_bin_edges=[-3.719456, -3.539592, -3.359728, -3.179864, -3.0, -2.8201359999999998, -2.640272, -2.4604079999999997, -2.280544], eps_bin_edges=[-3.0, 3.0], tectonic_region_types=['Active Shallow Crust'], lon=-3.0, lat=-3.0, weights=[0.25, 0.25, 0.25, 0.25], rlz_ids=[1, 0, 2, 3]"
#,,,,,,,"generated_by='OpenQuake engine 3.18.0-git59fd79fead', start_date='2023-08-28T14:18:03', checksum=3344811530, investigation_time=1.0, mag_bin_edges=[5.5, 5.6000000000000005, 5.7, 5.800000000000001, 5.9, 6.0, 6.1000000000000005, 6.2, 6.300000000000001, 6.4, 6.5, 6.6000000000000005], dist_bin_edges=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0], lon_bin_edges=[-3.7204433425622865, -3.540332506921715, -3.3602216712811437, -3.180110835640572, -3.0, -2.8198891643594286, -2.639778328718857, -2.4596674930782854, -2.2795566574377135], lat_bin_edges=[-3.719456, -3.539592, -3.359728, -3.179864, -3.0, -2.8201359999999998, -2.640272, -2.4604079999999997, -2.280544], eps_bin_edges=[-3.0, 3.0], tectonic_region_types=['Active Shallow Crust'], lon=-3.0, lat=-3.0, weights=[0.25, 0.25, 0.25, 0.25], rlz_ids=[1, 0, 2, 3]"
imt,iml,poe,mag,rlz1,rlz0,rlz2,rlz3
SA(0.1),1.00000E-01,3.03098E-03,5.55000E+00,1.51398E-03,5.08386E-03,0.00000E+00,0.00000E+00
SA(0.1),1.00000E-01,3.03098E-03,5.65000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00
Expand Down
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
#,,,,,,,"generated_by='OpenQuake engine 3.17.0-git5adb6553aa', start_date='2023-05-05T11:16:59', checksum=3344811530, investigation_time=1.0, mag_bin_edges=[5.5, 5.6000000000000005, 5.7, 5.800000000000001, 5.9, 6.0, 6.1000000000000005, 6.2, 6.300000000000001, 6.4, 6.5, 6.6000000000000005], dist_bin_edges=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0], lon_bin_edges=[-0.719456, -0.539592, -0.35972799999999994, -0.1798639999999999, 1.1102230246251565e-16, 0.17986400000000013, 0.35972800000000016, 0.5395920000000002, 0.7194560000000002], lat_bin_edges=[-0.719456, -0.539592, -0.35972799999999994, -0.1798639999999999, 1.1102230246251565e-16, 0.17986400000000013, 0.35972800000000016, 0.5395920000000002, 0.7194560000000002], eps_bin_edges=[-3.0, 3.0], tectonic_region_types=['Active Shallow Crust'], lon=0.0, lat=0.0, weights=[0.25, 0.25, 0.25, 0.25], rlz_ids=[3, 2, 0, 1]"
#,,,,,,,"generated_by='OpenQuake engine 3.18.0-git59fd79fead', start_date='2023-08-28T14:18:03', checksum=3344811530, investigation_time=1.0, mag_bin_edges=[5.5, 5.6000000000000005, 5.7, 5.800000000000001, 5.9, 6.0, 6.1000000000000005, 6.2, 6.300000000000001, 6.4, 6.5, 6.6000000000000005], dist_bin_edges=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0], lon_bin_edges=[-0.719456, -0.539592, -0.35972799999999994, -0.1798639999999999, 1.1102230246251565e-16, 0.17986400000000013, 0.35972800000000016, 0.5395920000000002, 0.7194560000000002], lat_bin_edges=[-0.719456, -0.539592, -0.35972799999999994, -0.1798639999999999, 1.1102230246251565e-16, 0.17986400000000013, 0.35972800000000016, 0.5395920000000002, 0.7194560000000002], eps_bin_edges=[-3.0, 3.0], tectonic_region_types=['Active Shallow Crust'], lon=0.0, lat=0.0, weights=[0.25, 0.25, 0.25, 0.25], rlz_ids=[3, 2, 0, 1]"
imt,iml,poe,mag,rlz3,rlz2,rlz0,rlz1
SA(0.1),1.00000E-01,2.45207E-02,5.55000E+00,3.40661E-02,3.84509E-02,5.11046E-03,1.35565E-03
SA(0.1),1.00000E-01,2.45207E-02,5.65000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00
Expand Down
18 changes: 9 additions & 9 deletions openquake/qa_tests_data/disagg/case_6/expected/Dist-0.csv
Original file line number Diff line number Diff line change
@@ -1,13 +1,13 @@
#,,,,,,"generated_by='OpenQuake engine 3.17.0-git5adb6553aa', start_date='2023-05-05T11:01:37', checksum=3577679219, investigation_time=1.0, mag_bin_edges=[6.5, 6.75, 7.0, 7.25], dist_bin_edges=[0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0, 225.0, 250.0, 275.0, 300.0], lon_bin_edges=[175.61128112572965, 176.55755408381978, 177.5038270419099, 178.45010000000002, 179.39637295809015, 180.34264591618026, 181.2889188742704], lat_bin_edges=[-20.822760000000002, -19.92344, -19.02412, -18.1248, -17.22548, -16.32616, -15.426840000000002], eps_bin_edges=[-3.0, 3.0], tectonic_region_types=['Active Shallow Crust'], lon=178.45010000000002, lat=-18.1248, weights=[0.3, 0.3, 0.4], rlz_ids=[2, 1, 0]"
#,,,,,,"generated_by='OpenQuake engine 3.18.0-git59fd79fead', start_date='2023-08-28T14:18:07', checksum=3577679219, investigation_time=1.0, mag_bin_edges=[6.5, 6.75, 7.0, 7.25], dist_bin_edges=[0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0, 225.0, 250.0, 275.0, 300.0], lon_bin_edges=[175.61128112572962, 176.55755408381975, 177.50382704190986, 178.4501, 179.39637295809013, 180.34264591618023, 181.28891887427037], lat_bin_edges=[-20.822760000000002, -19.92344, -19.02412, -18.1248, -17.22548, -16.32616, -15.426840000000002], eps_bin_edges=[-3.0, 3.0], tectonic_region_types=['Active Shallow Crust'], lon=178.4501, lat=-18.1248, weights=[0.3, 0.3, 0.4], rlz_ids=[2, 1, 0]"
imt,iml,poe,dist,rlz2,rlz1,rlz0
PGA,1.39032E-02,2.10500E-03,1.25000E+01,0.00000E+00,0.00000E+00,0.00000E+00
PGA,1.39032E-02,2.10500E-03,3.75000E+01,0.00000E+00,0.00000E+00,0.00000E+00
PGA,1.39032E-02,2.10500E-03,6.25000E+01,0.00000E+00,0.00000E+00,0.00000E+00
PGA,1.39032E-02,2.10500E-03,8.75000E+01,0.00000E+00,0.00000E+00,0.00000E+00
PGA,1.39032E-02,2.10500E-03,1.12500E+02,0.00000E+00,0.00000E+00,0.00000E+00
PGA,1.39032E-02,2.10500E-03,1.37500E+02,0.00000E+00,0.00000E+00,0.00000E+00
PGA,1.39032E-02,2.10500E-03,1.62500E+02,0.00000E+00,0.00000E+00,0.00000E+00
PGA,1.39032E-02,2.10500E-03,1.87500E+02,0.00000E+00,0.00000E+00,0.00000E+00
PGA,1.39032E-02,2.10500E-03,1.25000E+01,1.19904E-14,1.19904E-14,1.19904E-14
PGA,1.39032E-02,2.10500E-03,3.75000E+01,1.19904E-14,1.19904E-14,1.19904E-14
PGA,1.39032E-02,2.10500E-03,6.25000E+01,1.19904E-14,1.19904E-14,1.19904E-14
PGA,1.39032E-02,2.10500E-03,8.75000E+01,1.19904E-14,1.19904E-14,1.19904E-14
PGA,1.39032E-02,2.10500E-03,1.12500E+02,1.19904E-14,1.19904E-14,1.19904E-14
PGA,1.39032E-02,2.10500E-03,1.37500E+02,1.19904E-14,1.19904E-14,1.19904E-14
PGA,1.39032E-02,2.10500E-03,1.62500E+02,1.19904E-14,1.19904E-14,1.19904E-14
PGA,1.39032E-02,2.10500E-03,1.87500E+02,1.19904E-14,1.19904E-14,1.19904E-14
PGA,1.39032E-02,2.10500E-03,2.12500E+02,6.16581E-04,1.08641E-03,2.53621E-03
PGA,1.39032E-02,2.10500E-03,2.37500E+02,9.17975E-05,1.76954E-04,7.13312E-04
PGA,1.39032E-02,2.10500E-03,2.62500E+02,1.92071E-05,4.51343E-05,3.95415E-04
Expand Down
Loading

0 comments on commit 7414043

Please sign in to comment.