# -*- 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 with utility functions for gmpes.
"""
import numpy as np
import pandas as pd
import ast
import re
from openquake.hazardlib import valid
from openquake.hazardlib import scalerel
from openquake.hazardlib.geo import Point
from openquake.hazardlib.geo.geodetic import npoints_towards
from openquake.hazardlib.geo import utils as geo_utils
from openquake.hazardlib.site import Site, SiteCollection
from openquake.hazardlib.const import TRT
from openquake.hazardlib.contexts import ContextMaker
from openquake.hazardlib.gsim.mgmpe.modifiable_gmpe import ModifiableGMPE
from openquake.hazardlib.gsim.mgmpe.generic_gmpe_avgsa import (
GmpeIndirectAvgSA, GenericGmpeAvgSA)
from openquake.smt.utils import make_rup, clean_gmm_label
def _get_first_point(sfc, from_point):
"""
Get the first point in the collection of sites from the rupture.
"""
if from_point == 'MP':
return sfc.get_middle_point() # Get midpoint of rup surface
elif from_point == 'TC': # Get the up-dip edge centre point
return sfc.get_top_edge_centroid()
elif from_point == 'BC': # Get the down-dip edge centre point
lon, lat = geo_utils.get_middle_point(sfc.corner_lons[2],
sfc.corner_lats[2],
sfc.corner_lons[3],
sfc.corner_lats[3])
return Point(lon, lat, sfc.corner_depths[2])
elif from_point == 'TL': # Get top left point
idx = 0
elif from_point == 'TR': # Get top right point
idx = 1
elif from_point == 'BR': # Get bottom right
idx = 2
elif from_point == 'BL': # Get bottom left
idx = 3
else:
raise ValueError('Unsupported option from first point')
return Point(
sfc.corner_lons[idx], sfc.corner_lats[idx], sfc.corner_depths[idx])
[docs]
def get_sites_from_rupture(rup,
from_point='TC',
toward_azimuth=90,
direction='positive',
hdist=100,
step=5.,
site_props=''):
"""
Get the sites from the rupture to create the context with.
:param rup:
Rupture object
:param from_point:
A string. Options: 'TC', 'TL', 'TR', 'BR', 'BL'
:return:
A :class:`openquake.hazardlib.site.SiteCollection` instance
"""
from_pnt = _get_first_point(rup.surface, from_point)
r_lon = from_pnt.longitude
r_lat = from_pnt.latitude
r_dep = 0
vdist = 0
npoints = hdist / step
strike = rup.surface.strike
pointsp = []
pointsn = []
if direction == 'positive':
azi = (strike + toward_azimuth) % 360
pointsp = npoints_towards(r_lon,
r_lat,
r_dep,
azi,
hdist,
vdist,
npoints)
if direction == 'negative':
azi = (strike + toward_azimuth + 180) % 360
pointsn = npoints_towards(r_lon,
r_lat,
r_dep,
azi,
hdist,
vdist,
npoints)
sites = []
keys = set(site_props.keys()) - set(['vs30', 'z1pt0', 'z2pt5'])
if len(pointsn):
lons = reversed(pointsn[0][0:])
lats = reversed(pointsn[1][0:])
for lon, lat in zip(lons, lats):
site = Site(Point(lon, lat, 0.0),
vs30=site_props['vs30'],
z1pt0=site_props['z1pt0'],
z2pt5=site_props['z2pt5'])
for key in list(keys):
setattr(site, key, site_props[key])
sites.append(site)
if len(pointsp):
for lon, lat in zip(pointsp[0], pointsp[1]):
site = Site(Point(lon, lat, 0.0),
vs30=site_props['vs30'],
z1pt0=site_props['z1pt0'],
z2pt5=site_props['z2pt5'])
for key in list(keys):
setattr(site, key, site_props[key])
sites.append(site)
return SiteCollection(sites)
[docs]
def get_rup(mag, lon, lat, depth, ztor, aratio, strike, dip, rake, trt):
"""
Create an OQ rupture.
"""
# If TRT specified assign it and an MSR
if trt == 'active_crustal':
rup_trt = TRT.ACTIVE_SHALLOW_CRUST
rup_msr = scalerel.WC1994()
elif trt == 'stable':
rup_trt = TRT.STABLE_CONTINENTAL
rup_msr = scalerel.WC1994()
elif trt == 'slab':
rup_trt = TRT.SUBDUCTION_INTRASLAB
rup_msr = scalerel.strasser2010.StrasserIntraslab()
elif trt == 'interface':
rup_trt = TRT.SUBDUCTION_INTERFACE
rup_msr = scalerel.strasser2010.StrasserInterface()
else:
rup_trt = None
rup_msr = scalerel.WC1994()
# Get rupture
rup = make_rup(lon,
lat,
depth,
msr=rup_msr,
mag=mag,
aratio=aratio,
strike=strike,
dip=dip,
rake=rake,
trt=rup_trt,
ztor=ztor)
return rup
[docs]
def att_curves(gmpe,
mag,
lon,
lat,
depth,
ztor,
aratio,
strike,
dip,
rake,
trt,
oq_rup,
vs30,
z1pt0,
z2pt5,
maxR,
step,
imt,
dist_type,
up_or_down_dip,
volc_back_arc,
eshm20_region):
"""
Compute the ground-motion intensities for the given context created here.
"""
# Make rupture if not provided from XML or CSV
if oq_rup is None:
rup = get_rup(mag, lon, lat, depth, ztor, aratio, strike, dip, rake, trt)
else:
rup = oq_rup
# Set site props
props = {'vs30': vs30,
'z1pt0': z1pt0,
'z2pt5': z2pt5,
'backarc': volc_back_arc,
'vs30measured': False,
'eshm20_region': eshm20_region}
# Check if site up-dip or down-dip of site
if up_or_down_dip == float(1):
direction = 'positive'
elif up_or_down_dip == float(0):
direction = 'negative'
else:
raise ValueError('The site must be specified as up or down dip.')
# Get sites
if dist_type in ['repi', 'rhypo']:
from_pnt = 'MP' # Sites from midpoint of rup surface
else:
from_pnt = 'TC' # Sites from center of top edge
sites = get_sites_from_rupture(rup,
from_point=from_pnt,
toward_azimuth=90,
direction=direction,
hdist=maxR,
step=step,
site_props=props)
# Add main R types to gmpe so can plot against repi, rrup, rjb and rhypo
core_r_types = ['repi', 'rrup', 'rjb', 'rhypo']
orig_r_types = list(gmpe.REQUIRES_DISTANCES)
for core in core_r_types:
if core not in orig_r_types:
orig_r_types.append(core)
gmpe.REQUIRES_DISTANCES = frozenset(orig_r_types)
# Create context
mag_str = [f'{mag:.2f}']
oqp = {'imtls': {k: [] for k in [str(imt)]}, 'mags': mag_str}
ctxm = ContextMaker(rup.tectonic_region_type, [gmpe], oqp)
ctxs = list(ctxm.get_ctxs([rup], sites))
ctxs = ctxs[0]
# Compute ground-motions
mean, std, tau, phi = ctxm.get_mean_stds([ctxs])
if dist_type == 'repi':
distances = ctxs.repi
elif dist_type == 'rrup':
distances = ctxs.rrup
elif dist_type == 'rjb':
distances = ctxs.rjb
elif dist_type == 'rhypo':
distances = ctxs.rhypo
else:
raise ValueError('No valid distance type specified.')
return mean, std, distances, tau, phi
[docs]
def get_rup_pars(strike, dip, rake, aratio, trt):
"""
Get (crude) proxies for strike, dip and aspect ratio if not
provided by the user.
"""
# Strike
if strike == -999:
strike_s = 0
else:
strike_s = strike
# Dip
if dip == -999:
if rake == 0 or rake == 180:
dip_s = 90 # Strike slip
else:
dip_s = 45 # Reverse or normal fault
else:
dip_s = dip
# Prevent assigning neither a trt or an aratio
if trt == -999 and aratio == -999:
msg = ('An aratio must be provided by the user, or alternatively the user '
'must provide a TRT (a trt-dependent aratio is then assigned instead)')
raise ValueError(msg)
# Aspect ratio
if aratio != -999.0 and np.isfinite(aratio):
aratio_s = aratio
else:
if trt in ['slab', 'interface']:
aratio_s = 5
else:
aratio_s = 2 # Crustal
return strike_s, dip_s, aratio_s
[docs]
def build_indirect_avgsa_gmpe(gmpe, avgsa, kw_mgmpe):
"""
Build a GMPE which can be used to predict AvgSA using the
indirect approach. This function can build either the
GmpeIndirectAvgSA class (which requires t_low, t_high, and
the number of periods) OR GenericGmpeAvgSA (specify a list
of averaging periods).
"""
check = list(kw_mgmpe.keys())
if check[0] != "gmpe" or len(check) > 1:
raise ValueError(
"Specification of an indirect approach AvgSA GMPE in combination "
"with additional ModifiableGMPE capabilities is not supported.")
gmm_base = list(kw_mgmpe['gmpe'].keys())[0]
for par in kw_mgmpe['gmpe'][gmm_base].keys():
avgsa[par] = kw_mgmpe['gmpe'][gmm_base][par]
if "GmpeIndirectAvgSA" in gmpe:
return GmpeIndirectAvgSA(gmpe_name=gmm_base, **avgsa)
else:
return GenericGmpeAvgSA(gmpe_name=gmm_base, **avgsa)
[docs]
def construct_gsim_dict(inputs):
"""
Build a dictionary of the arguments for a GMM.
"""
# Build dict
kwargs = {}
parts = re.search(r'\[([^\]]+)\]', inputs) # Square brackets = extra inputs
if parts:
start = parts.group(1) # GMM without square brackets
other = inputs.split(parts.group(0))[1]
other = re.sub(r'\\+n', '\n', other)
other = re.sub(r'[\\\'"]', '', other)
kwargs['gmpe'] = {start: dict(re.findall(r'(\w+)\s*=\s*([^\n]+)', other))}
else:
kwargs['gmpe'] = {inputs: {}} # GMM without any additional arguments
# Force float for appropriate params
for gmpe in kwargs["gmpe"]:
params = kwargs["gmpe"][gmpe]
for param in params:
value = kwargs["gmpe"][gmpe][param]
try:
kwargs["gmpe"][gmpe][param] = float(value)
except:
pass
return kwargs
[docs]
def build_mgmpe(gmpe):
"""
Build a ModifiableGMPE from a string of a GMPE parsed from a Comparison TOML.
NOTE: The way ModifiableGMPEs are built from the comparison TOML makes them
slightly less flexible than when built from an OpenQuake XML. Such limitations
can be seen by inspecting this code. For example, for the Al Atik sigma model,
we always use the "global" tau model. An experienced user can of course modify
below any hard-coded values if they deem necessary.
"""
# All of the inputs for this model
params = pd.Series(gmpe.splitlines(), dtype=object)
# Underlying GMM to modify
base_gsim = re.search(r'gmpe\s*=\s*(.*)', params.iloc[1]).group(1).replace('"','')
# Construct dict of gsim kwargs
kw_mgmpe = construct_gsim_dict(base_gsim)
# Get the mgmpe params
idx_params = []
for idx, par in enumerate(params):
if idx <= 1:
continue
par = str(par).strip()
# Split key and value
if '=' in par:
key, val = par.split('=', 1)
key = key.strip()
val = val.strip().replace('"', '')
else:
key, val = par, None
if any(k in key for k in ['sigma_model', 'site_term', 'basin_term']):
idx_params.append(idx)
elif key == 'fix_total_sigma':
idx_params.append(idx)
fixed_sigma_vector = ast.literal_eval(val)
elif key == 'with_betw_ratio':
idx_params.append(idx)
with_betw_ratio = float(val)
elif key == 'set_between_epsilon':
idx_params.append(idx)
between_epsilon = float(val)
elif key == 'add_delta_to_total_scalar':
idx_params.append(idx)
delta_std_scalar = float(val)
elif key == 'add_delta_to_tau_scalar':
idx_params.append(idx)
delta_tau_scalar = float(val)
elif key == 'add_delta_to_phi_scalar':
idx_params.append(idx)
delta_phi_scalar = float(val)
elif key == 'add_delta_to_total_vector':
idx_params.append(idx)
delta_std_vector = ast.literal_eval(val)
elif key == 'add_delta_to_tau_vector':
idx_params.append(idx)
delta_tau_vector = ast.literal_eval(val)
elif key == 'add_delta_to_phi_vector':
idx_params.append(idx)
delta_phi_vector = ast.literal_eval(val)
elif key == 'set_total_sigma_as_tau_plus_delta':
idx_params.append(idx)
total_set_to_tau_and_delta = float(val)
elif 'scaling' in key:
idx_params.append(idx)
if key == 'median_scaling_scalar':
median_scalar = float(val)
elif key == 'median_scaling_vector':
median_vector = ast.literal_eval(val)
elif key == 'sigma_scaling_scalar':
sigma_scalar = float(val)
elif key == 'sigma_scaling_vector':
sigma_vector = ast.literal_eval(val)
elif key == "conditional_gmpe":
idx_params.append(idx)
re_match = re.search(r'conditional_gmpe\s*=\s*"(.+)"', par, re.DOTALL)
cgmpe_dict = ast.literal_eval(re_match.group(1))
cgmpes = {
imt: construct_gsim_dict(gmpe_str)
for imt, gmpe_str in cgmpe_dict.items()
}
elif key in ["GmpeIndirectAvgSA", "GenericGmpeAvgSA"]:
idx_params.append(idx)
avgsa = ast.literal_eval(val)
# Add the non-gmpe kwargs
for idx_p, param in enumerate(params):
if idx_p > 1 and idx_p not in idx_params:
if 'lt_weight' not in param: # Skip if weight for logic tree
dic_key = param.split('=')[0].strip().replace('"','')
dic_val = param.split('=')[1].strip().replace('"','')
kw_mgmpe['gmpe'][base_gsim][dic_key] = dic_val
# Al Atik 2015 sigma model
if 'al_atik_2015_sigma' in gmpe:
kw_mgmpe['sigma_model_alatik2015'] = {"tau_model": "global", "ergodic": False}
# Fix total sigma per imt
if 'fix_total_sigma' in gmpe:
kw_mgmpe['set_fixed_total_sigma'] = {'total_sigma': fixed_sigma_vector}
# Partition total sigma using a specified ratio of within:between
if 'with_betw_ratio' in gmpe:
kw_mgmpe['add_between_within_stds'] = {'with_betw_ratio': with_betw_ratio}
# Set epsilon for tau and use instead of total sigma
if 'set_between_epsilon' in gmpe:
kw_mgmpe['set_between_epsilon'] = {'epsilon_tau': between_epsilon}
# Add IMT-constant delta to total sigma
if 'add_delta_to_total_scalar' in gmpe:
kw_mgmpe['add_delta_to_total_std_scalar'] = {'delta': delta_std_scalar}
# Add IMT-constant delta to tau
if 'add_delta_to_tau_scalar' in gmpe:
kw_mgmpe['add_delta_to_tau_std_scalar'] = {'delta': delta_tau_scalar}
# Add IMT-constant delta to phi
if 'add_delta_to_phi_scalar' in gmpe:
kw_mgmpe['add_delta_to_phi_std_scalar'] = {'delta': delta_phi_scalar}
# Add IMT-dependent delta to total sigma
if 'add_delta_to_total_vector' in gmpe:
kw_mgmpe['add_delta_to_total_std_vector'] = {'delta': delta_std_vector}
# Add IMT-dependent delta to tau
if 'add_delta_to_tau_vector' in gmpe:
kw_mgmpe['add_delta_to_tau_std_vector'] = {'delta': delta_tau_vector}
# Add IMT-constant delta to phi
if 'add_delta_to_phi_vector' in gmpe:
kw_mgmpe['add_delta_to_phi_std_vector'] = {'delta': delta_phi_vector}
# Set total sigma to sqrt(tau**2 + delta**2)
if 'set_total_sigma_as_tau_plus_delta' in gmpe:
kw_mgmpe['set_total_std_as_tau_plus_delta'] = {'delta': total_set_to_tau_and_delta}
# Scale median by constant factor over all imts
if 'median_scaling_scalar' in gmpe:
kw_mgmpe['set_scale_median_scalar'] = {'scaling_factor': median_scalar}
# Scale median by imt-dependent factor
if 'median_scaling_vector' in gmpe:
kw_mgmpe['set_scale_median_vector'] = {'scaling_factor': median_vector}
# Scale sigma by constant factor over all imts
if 'sigma_scaling_scalar' in gmpe:
kw_mgmpe['set_scale_total_sigma_scalar'] = {'scaling_factor': sigma_scalar}
# Scale sigma by imt-dependent factor
if 'sigma_scaling_vector' in gmpe:
kw_mgmpe['set_scale_total_sigma_vector'] = {'scaling_factor': sigma_vector}
# CY14SiteTerm
if 'CY14SiteTerm' in gmpe:
kw_mgmpe['cy14_site_term'] = {}
# BA08SiteTerm
if 'BA08SiteTerm' in gmpe:
kw_mgmpe['ba08_site_term'] = {}
# BSSA14SiteTerm
if "BSSA14SiteTerm" in gmpe:
kw_mgmpe['bssa14_site_term'] = {}
# NRCan15SiteTerm ("base" kind)
if ('NRCan15SiteTerm' in gmpe and 'NRCan15SiteTermLinear' not in gmpe):
kw_mgmpe['nrcan15_site_term'] = {'kind': 'base'}
# NRCan15SiteTerm ("linear" kind)
if 'NRCan15SiteTermLinear' in gmpe:
kw_mgmpe['nrcan15_site_term'] = {'kind': 'linear'}
# CEUS2020SiteTerm (Stewart et al. 2020)
if 'CEUS2020SiteTerm' in gmpe:
try:
assert "_refVs30=" in gmpe
ref_vs30 = float(gmpe.split("refVs30=")[-1].replace("'",'').replace('"',''))
except:
raise ValueError("If using the CEUS2020SiteTerm the user must also specify a "
"ref vs30 to be used for the non-linear scaling component.")
kw_mgmpe['ceus2020_site_term'] = {"ref_vs30": ref_vs30, 'wimp': None}
# CB14 basin term
if 'CB14BasinTerm' in gmpe:
kw_mgmpe['cb14_basin_term'] = {}
# M9 basin adjustment
if 'M9BasinTerm' in gmpe:
kw_mgmpe['m9_basin_term'] = {}
# Conditional GMPE(s)
if 'conditional_gmpe' in gmpe:
kw_mgmpe['conditional_gmpe'] = cgmpes
# Indirect approach AvgSA GMPE
if "GmpeIndirectAvgSA" in gmpe or "GenericGmpeAvgSA" in gmpe:
return build_indirect_avgsa_gmpe(gmpe, avgsa, kw_mgmpe)
return ModifiableGMPE(**kw_mgmpe)
[docs]
def gmpe_check(gmpe):
"""
This function in effect parses the toml parameters for a GMPE into the
equivalent parameters required for constructing an OQ GSIM object.
:param gmpe: GMM and params parsed from the Comparison toml.
"""
# Modifiable GMPE
if '[ModifiableGMPE]' in gmpe:
return build_mgmpe(gmpe)
# Regular GMPE
else:
# Clean to ensure arguments can be passed (the logic tree weights
# are retained in original GMM strings in utils_compare_gmpes.py)
params = pd.Series(gmpe.splitlines())
idx_to_drop = []
for idx_p, par in enumerate(params):
if 'lt_weight_gmc' in par:
idx_to_drop.append(idx_p)
params = params.drop(idx_to_drop)
gmpe_clean = params.iloc[0].strip()
if len(params) > 1:
for idx_p, par in enumerate(params):
if idx_p > 0:
gmpe_clean = gmpe_clean + '\n' + par
else: # Ensures GSIM aliases work
gmpe_clean = gmpe_clean.replace('[','').replace(']','')
gmm = valid.gsim(gmpe_clean)
return gmm
[docs]
def get_imtl_unit(i):
"""
Return a string of the intensity measure type's physical units of
measurement.
"""
if str(i) in ['PGD', 'SDi']:
unit = 'cm' # PGD, inelastic spectral displacement
elif str(i) in ['PGV']:
unit = 'cm/s' # PGV
elif str(i) in ['IA']:
unit = 'm/s' # Arias intensity
elif str(i) in ['RSD', 'RSD595', 'RSD575', 'RSD2080', 'DRVT']:
unit = 's' # Relative significant duration, DRVT
elif str(i) in ['CAV']:
unit = 'g-sec' # Cumulative absolute velocity
elif str(i) in ['MMI']:
unit = 'MMI' # Modified Mercalli Intensity
elif str(i) in ['FAS', 'EAS']:
unit = str(i) + ' (Hz)' # Fourier/Eff. Amp. Spectrum
else:
if str(i) not in ["PGA", "AvgSA"]:
assert "SA" in str(i)
unit = 'g' # PGA, SA, AvgSA
return unit
[docs]
def matrix_to_df(matrix, ngmms):
"""
Convert matrix of ground-motions to dataframe with
one column per IMT and values being the flattened
array of predictions from each GMPE.
This function also checks that the number of arrays
per IMT is equal to the number of GMPEs specified in
the TOML as a sanity check.
Currently only used in ModifiableGMPE-based unit tests.
"""
store = {}
for imt in matrix.keys():
assert len(matrix[imt]) == ngmms
store[str(imt)] = np.array(matrix[imt]).flatten()
return pd.DataFrame(store)