# -*- 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 generating trellis plots, response spectra,
hierarchical clustering plots, Sammon maps and Euclidean distance matrix plots.
"""
import os
import numpy as np
import pandas as pd
from matplotlib import pyplot
from scipy.cluster import hierarchy
from scipy.spatial.distance import pdist, squareform
from scipy import interpolate
from openquake.hazardlib.imt import PGA, SA
from openquake.smt.comparison.sammons import sammon
from openquake.smt.utils import clean_gmm_label, COLORS, GEM_FF_MAPPINGS
from openquake.smt.comparison.utils_gmpes import (get_imtl_unit,
att_curves,
get_rup_pars,
gmpe_check)
# Periods used in spectra plotting (truncated based on max_period in config)
PERIODS = [PGA(), SA(0.025), SA(0.04), SA(0.05), SA(0.07), SA(0.1), SA(0.15),
SA(0.2), SA(0.25), SA(0.3), SA(0.35), SA(0.4), SA(0.45), SA(0.5),
SA(0.6), SA(0.7), SA(0.8), SA(0.9), SA(1.0), SA(1.1), SA(1.2), SA(1.3),
SA(1.4), SA(1.5), SA(1.6), SA(1.7), SA(1.8), SA(1.9), SA(2.0), SA(2.2),
SA(2.4), SA(2.6), SA(2.8), SA(3.0), SA(3.2), SA(3.4), SA(3.6), SA(3.8),
SA(4.0), SA(4.5), SA(5.0), SA(5.5), SA(6.0), SA(6.5), SA(7.0), SA(7.5),
SA(8.0), SA(8.5), SA(9.0), SA(9.5), SA(10.0)]
# Fltering params for plotting observations against GMPEs
MAG_LIM = 0.25 # Mw
DEP_LIM = 15 # km
VS30_LIM = 150 # m/s
DIST_LIM_LOW = 5 # km (near-source distance window)
DIST_LIM_MID = 10 # km (intermediate distance window)
DIST_LIM_MAX = 20 # km (far-source distance window)
### Main Plotting Functions ###
[docs]
def plot_trellis_util(config, output_directory, obs_data_fname):
"""
Generate trellis plots for given run configuration.
"""
# Load observed data if provided
no_obs = True
if obs_data_fname is not None:
data = pd.read_csv(obs_data_fname, low_memory=False)
else:
data = None
# Median, plus sigma, minus sigma per gmc for up to 4 gmc logic trees
gmc_p= {lt: [{}, {}, {}] for lt in config.lt_mapping.keys()}
# Get lt weights
lt_weights = {gmc: getattr(config, config.lt_mapping[gmc]['wei']) for gmc in gmc_p}
# Get config key
cfg_key = f'vs30 = {config.vs30} m/s, GMM sigma epsilon = {config.nstd}'
# Get colours
colors = get_colors(config.custom_color_flag, config.custom_color_list)
# Compute attenuation curves
store_gmm_curves, store_per_imt = {}, {} # For exporting gmm att curves
store_gmm_curves[cfg_key] = {}
store_gmm_curves[cfg_key]['gmm att curves per imt-mag'] = {}
store_gmm_curves[cfg_key]['gmc logic tree curves per imt-mag'] = {}
fig = pyplot.figure(figsize=(len(config.mag_list)*5, len(config.imt_list)*4))
max_pred, min_pred, axs = [], [], []
for i, imt in enumerate(config.imt_list):
store_per_mag = {}
for m, mag in enumerate(config.mag_list):
# Add the axis
ax = fig.add_subplot(len(config.imt_list), len(config.mag_list), m+1+i*len(config.mag_list))
axs.append(ax)
# Get depth params
depth_g = config.depth_list[m]
if config.ztor != -999:
ztor_g = config.ztor[m]
else:
ztor_g = None
# Get rupture params
strike_g, dip_g, aratio_g = get_rup_pars(config.strike,
config.dip,
config.rake,
config.aratio,
config.trt)
# If plotting data get the appropriate subset
if data is not None:
subset = filter_flatfile_trellis(
data, imt, mag, depth_g, config.vs30, config.dist_type)
else:
subset = None
# Per GMPE get attenuation curves
lt_vals_gmc = {lt: {} for lt in lt_weights}
store_per_gmpe = {}
for g, gmpe in enumerate(config.gmpes_list):
# Sub dicts for median, gmm sigma, median +/- nstd * gmm sigma
store_per_gmpe[gmpe] = {}
col = colors[g]
# Perform gmpe check
gmm = gmpe_check(gmpe)
# Get attenuation curves
mean, std, r_vals, tau, phi = att_curves(gmm,
mag,
config.lon,
config.lat,
depth_g,
ztor_g,
aratio_g,
strike_g,
dip_g,
config.rake,
config.trt,
config.rup,
config.vs30,
config.z1pt0,
config.z2pt5,
config.maxR,
1, # Step of 1 km for site spacing
imt,
config.dist_type,
config.up_or_down_dip,
config.volc_back_arc,
config.eshm20_region)
# Get mean, sigma components, mean plus/minus sigma
mean = mean[0][0]
std = std[0][0]
add_sigma = np.exp(mean+config.nstd*std[0])
min_sigma = np.exp(mean-config.nstd*std[0])
# For managing ylim
max_pred.append(np.max([np.exp(mean), add_sigma]))
min_pred.append(np.min([np.exp(mean), min_sigma]))
# Plot predictions and get lt weighted predictions
lt_vals_gmc = trellis_data(gmpe,
r_vals,
mean,
add_sigma,
min_sigma,
col,
config.nstd,
lt_vals_gmc,
lt_weights)
# Get unit of imt for the store
unit = get_imtl_unit(imt)
# Store per gmpe
store_per_gmpe[gmpe]['median (%s)' % unit] = np.exp(mean)
store_per_gmpe[gmpe]['sigma (ln)'] = std
if config.nstd != 0:
store_per_gmpe[gmpe]['median plus sigma (%s)' % unit] = add_sigma
store_per_gmpe[gmpe]['median minus sigma (%s)' % unit] = min_sigma
# Update plots
update_trellis_plots(mag,
imt,
m,
i,
depth_g,
config.vs30,
config.minR,
config.maxR,
r_vals,
config.imt_list,
config.dist_type)
# Plot logic trees if specified and also store
for key_gmc in lt_weights:
store_gmm_curves = trellis_logic_trees(config,
key_gmc,
lt_weights[key_gmc],
lt_vals_gmc[key_gmc],
gmc_p[key_gmc],
store_gmm_curves,
r_vals,
config.nstd,
imt,
mag,
depth_g,
dip_g,
config.rake,
cfg_key,
unit)
# Create key of magnitude and other scenario info
mag_key = f'Mw = {mag}, depth = {depth_g} km, dip = {dip_g} deg, rake = {config.rake} deg'
# Add the distance values to each GMM (avoid's overwrite)
if config.dist_type in ["repi", "rhypo"]:
if r_vals[0] < 1E-09: r_vals[0] = 0 # Precision issue
store_per_gmpe['%s (km)' % config.dist_type] = r_vals
# Store the GMM's info
store_per_mag[mag_key] = store_per_gmpe
# Add grid
pyplot.grid(axis='both', which='both', alpha=0.5)
# Plot data too if required/any retrieved
if subset is not None:
# Set no_obs to False to ensure legend entry added at end of loops
no_obs = False
# NOTE: Units are converted to OQ GSIM units in helper functions
pyplot.scatter(x=subset[GEM_FF_MAPPINGS[config.dist_type]],
y=subset[GEM_FF_MAPPINGS[imt]["col"]],
color="k", marker="x", zorder=0)
# Store per imt
store_per_imt[str(imt)] = store_per_mag
# Store all the curves
store_gmm_curves[cfg_key]['gmm att curves per imt-mag'] = store_per_imt
# Finalise plots
maxy = np.max(max_pred)
miny = np.min(min_pred)
for ax in axs: ax.set_ylim(miny, 2*maxy) # Small buffer in log-space
output = os.path.join(output_directory, 'TrellisPlots.png')
if no_obs is False:
# If any suitable data plotted add to legend
ax.scatter([], [], color='k', marker='x', label='Flatfile Data')
pyplot.legend(loc="center left", bbox_to_anchor=(1.1, 1.05), fontsize='16')
pyplot.savefig(output, bbox_inches='tight', dpi=200, pad_inches=0.2)
pyplot.close()
return store_gmm_curves
[docs]
def plot_spectra_util(config, output_directory, obs_spectra_fname, obs_data_fname):
"""
Plot response spectra for given run configuration. Can also plot an
observed spectrum and the corresponding predictions by the specified
GMPEs.
"""
# Check distances have been provided in the input TOML
if len(config.dist_list) < 1:
raise ValueError("Response spectra have been requested but no distance "
"intervals have been specified in the input toml.")
# If obs spectra csv provided load the data
if obs_spectra_fname is not None:
obs_spectra, max_period, eq_id, st_id = load_obs_spectra(obs_spectra_fname)
else:
max_period = config.max_period
obs_spectra, eq_id, st_id = None, None, None
# Load observed data if provided
no_obs = True
if obs_data_fname is not None:
data = pd.read_csv(obs_data_fname, low_memory=False)
else:
data = None
# Truncate periods to max_period
imt_list = [imt for imt in PERIODS if imt.period <= max_period]
periods = np.array([imt.period for imt in PERIODS if imt.period <= max_period])
# Get gmc lt weights
gmc_weights = {gmc: getattr(config, config.lt_mapping[gmc]['wei']) for gmc in config.lt_mapping.keys()}
# Get colours and make the figure
colors = get_colors(config.custom_color_flag, config.custom_color_list)
fig = pyplot.figure(figsize=(len(config.mag_list)*5, len(config.dist_list)*4))
# Set dicts to store values
lt_vals = {
# Keys for weighted GMM branches to compute LTs with
'med_wei': {ltw: {gmm: {} for gmm in gmc_weights[ltw].keys()}
if gmc_weights[ltw] is not None else {} for ltw in gmc_weights},
'add_wei': {'lt_gmc_1': {gmm: {} for gmm in config.gmpes_list}, # Set for even those without
'lt_gmc_2': {gmm: {} for gmm in config.gmpes_list}, # GMMs as makes assigning vals
'lt_gmc_3': {gmm: {} for gmm in config.gmpes_list}, # later more straightfoward
'lt_gmc_4': {gmm: {} for gmm in config.gmpes_list}},
'min_wei': {'lt_gmc_1': {gmm: {} for gmm in config.gmpes_list},
'lt_gmc_2': {gmm: {} for gmm in config.gmpes_list},
'lt_gmc_3': {gmm: {} for gmm in config.gmpes_list},
'lt_gmc_4': {gmm: {} for gmm in config.gmpes_list}},
# Keys for aggregated gmm LTs
'lt_gmc_1': {},
'lt_gmc_2': {},
'lt_gmc_3': {},
'lt_gmc_4': {},
# Keys for non-weighted individual gmms
"med": {gmm: {} for gmm in config.gmpes_list},
'add': {gmm: {} for gmm in config.gmpes_list},
'min': {gmm: {} for gmm in config.gmpes_list},
# Useful info when exporting
'periods': periods,
'nstd': config.nstd
}
# Plot the data
for d, dist in enumerate(config.dist_list):
for m, mag in enumerate(config.mag_list):
ax = fig.add_subplot(
len(config.dist_list), len(config.mag_list), m+1+d*len(config.mag_list))
# Get depth params
depth_g = config.depth_list[m]
if config.ztor != -999:
ztor_g = config.ztor[m]
else:
ztor_g = None
# Get rupture params
strike_g, dip_g, aratio_g = get_rup_pars(config.strike,
config.dip,
config.rake,
config.aratio,
config.trt)
# If plotting data get the appropriate subset
if data is not None:
subset = filter_flatfile_spectra(
data, imt_list, mag, depth_g, config.vs30, dist, config.dist_type)
else:
subset = None
# Scenario key
sk = f"{config.dist_type}={dist}km, Mw={mag}, depth={depth_g}km, vs30={config.vs30}m/s"
# Iterate over the GMMs
for g, gmpe in enumerate(config.gmpes_list):
# Set stores
rs_50p, sig, rs_ps, rs_ms = [], [], [], []
col = colors[g]
# Perform gmpe check
gmm = gmpe_check(gmpe)
for i, imt in enumerate(imt_list):
# Get mean and sigma
mu, std, r_vals, tau, phi = att_curves(gmm,
mag,
config.lon,
config.lat,
depth_g,
ztor_g,
aratio_g,
strike_g,
dip_g,
config.rake,
config.trt,
config.rup,
config.vs30,
config.z1pt0,
config.z2pt5,
500, # Assume record dist < 500 km
1, # Step of 1 km for site spacing
imt.string,
config.dist_type,
config.up_or_down_dip,
config.volc_back_arc,
config.eshm20_region)
# Interpolate for distances and store
mu = mu[0][0]
f = interpolate.interp1d(r_vals, mu)
try:
rs_50p_dist = np.exp(f(dist))
except:
raise_spectra_dist_error(dist, config.dist_type, r_vals)
rs_50p.append(rs_50p_dist)
f1 = interpolate.interp1d(r_vals, std[0])
sigma_dist = f1(dist)
sig.append(sigma_dist)
if config.nstd != 0:
rs_add_sigma_dist = np.exp(f(dist)+(config.nstd*sigma_dist))[0]
rs_min_sigma_dist = np.exp(f(dist)-(config.nstd*sigma_dist))[0]
rs_ps.append(rs_add_sigma_dist)
rs_ms.append(rs_min_sigma_dist)
# Plot individual GMPEs
if 'plot_lt_only' not in str(gmpe):
ax.plot(periods,
rs_50p,
color=col,
linewidth=2,
linestyle='-',
label=clean_gmm_label(gmpe),
zorder=1)
if config.nstd > 0:
ax.plot(periods, rs_ps, color=col, linewidth=0.75, linestyle='-.')
ax.plot(periods, rs_ms, color=col, linewidth=0.75, linestyle='-.')
# Weight the predictions using logic tree weights
gmc_vals = spectra_data(gmpe,
config.nstd,
gmc_weights,
rs_50p,
rs_ps,
rs_ms,
lt_vals,
sk)
# Plot obs spectra if required
if obs_spectra is not None:
plot_obs_spectra(ax,
obs_spectra,
g,
config.gmpes_list,
config.mag_list,
config.depth_list,
config.dist_list,
config.dist_type,
config.vs30,
eq_id,
st_id)
# Update plots
update_spectra_plots(
ax, mag, depth_g, dist, config.vs30, d, m, config.dist_list, config.dist_type)
# Plot logic trees if required
for key_gmc in gmc_weights:
if gmc_vals[key_gmc][0] != {}: # If none empty LT
lt_vals[key_gmc][sk] = spectra_logic_trees(config,
ax,
config.gmpes_list,
config.nstd,
periods,
key_gmc,
gmc_vals[key_gmc],
sk)
# Plot data too if required/any retrieved
if subset is not None:
# Set no_obs to False to ensure legend entry added at end of loops
no_obs = False
# NOTE: Units are converted to OQ GSIM units in helper functions
for idx_rec, rec in subset.iterrows():
ax.plot(rec["spectra_periods"], rec["spectra_rotD50"],
color='k', linewidth=1.5)
# Add grid and set xlims
ax.set_xlim(min(periods), max(periods))
ax.grid(True)
if config.nstd > 0 or obs_data_fname is not None:
ax.semilogy()
# Finalise the plots and save fig
if len(config.mag_list) * len(config.dist_list) == 1:
bbox_coo = (1.1, 0.5)
fs = '10'
else:
bbox_coo = (1.1, 1.05)
fs = '16'
if no_obs is False:
# If any suitable data plotted add to legend
ax.plot([], [], color='k', linewidth=1.5, label='Flatfile Spectra')
out = os.path.join(output_directory, 'ResponseSpectra.png')
ax.legend(loc="center left", bbox_to_anchor=bbox_coo, fontsize=fs)
fig.savefig(out, bbox_inches='tight', dpi=200, pad_inches=0.2)
pyplot.close()
return lt_vals
[docs]
def plot_ratios_util(config, output_directory):
"""
Generate ratio (GMPE median attenuation/baseline GMPE
median attenuation) plots for given run configuration.
NOTE: The ratios of any specified GMC logic trees against
the baseline GMM are not computed/plotted.
"""
# Get colours
colors = get_colors(config.custom_color_flag, config.custom_color_list)
# Compute ratio curves
fig = pyplot.figure(figsize=(len(config.mag_list)*5, len(config.imt_list)*4))
ratio_store, axs = [], []
for i, imt in enumerate(config.imt_list):
for m, mag in enumerate(config.mag_list):
ax = fig.add_subplot(len(config.imt_list), len(config.mag_list), m+1+i*len(config.mag_list))
axs.append(ax)
# Get depth params
depth_g = config.depth_list[m]
if config.ztor != -999:
ztor_g = config.ztor[m]
else:
ztor_g = None
# Get rupture params
strike_g, dip_g, aratio_g = get_rup_pars(config.strike,
config.dip,
config.rake,
config.aratio,
config.trt)
# Load the baseline GMM and compute baseline
baseline = gmpe_check(config.baseline_gmm)
# Get baseline GMM attenuation curves
results = att_curves(baseline,
mag,
config.lon,
config.lat,
depth_g,
ztor_g,
aratio_g,
strike_g,
dip_g,
config.rake,
config.trt,
config.rup,
config.vs30,
config.z1pt0,
config.z2pt5,
config.maxR,
1, # Step of 1 km for sites
imt,
config.dist_type,
config.up_or_down_dip,
config.volc_back_arc,
config.eshm20_region)
# Get baseline mean
b_mean = results[0][0][0]
if np.all(b_mean) == 0:
# Should only occur in case of using a conditional GMPE
# which also does not support the requested IMT
assert imt not in baseline.params["conditional_gmpe"]
raise ValueError(f"A conditional GMPE which does not "
f"support {imt} has been specified "
f"for as the baseline model in GMPE "
f"ratio plotting.")
# Now compute ratios for each GMM
for g, gmpe in enumerate(config.gmpes_list):
# Perform gmpe check
col = colors[g]
gmm = gmpe_check(gmpe)
# Get attenuation curves for the GMM
results = att_curves(gmm,
mag,
config.lon,
config.lat,
depth_g,
ztor_g,
aratio_g,
strike_g,
dip_g,
config.rake,
config.trt,
config.rup,
config.vs30,
config.z1pt0,
config.z2pt5,
config.maxR,
1, # Step of 1 km for sites
imt,
config.dist_type,
config.up_or_down_dip,
config.volc_back_arc,
config.eshm20_region)
# Get mean and r_vals
mean = results[0][0][0]
r_vals = results[2]
# Compute GMM/baseline
ratio = np.exp(mean)/np.exp(b_mean)
ratio_store.append(ratio)
# Plot ratios
pyplot.semilogy(r_vals,
ratio,
color=col,
linewidth=2,
linestyle='-',
label=clean_gmm_label(gmpe))
# Update plots
update_ratio_plots(mag,
imt,
m,
i,
depth_g,
config.vs30,
config.minR,
config.maxR,
r_vals,
config.imt_list,
config.dist_type)
# Add grid
pyplot.grid(axis='both', which='both', alpha=0.5)
# Finalise plots
for ax in axs: ax.set_ylim(1/2*np.min(ratio_store), 2*np.max(ratio_store)) # Small buffer in log-space
out = os.path.join(output_directory, 'RatioPlots.png')
pyplot.legend(loc="center left", bbox_to_anchor=(1.1, 1.05), fontsize='16')
pyplot.savefig(out, bbox_inches='tight', dpi=200, pad_inches=0.2)
pyplot.close()
[docs]
def compute_matrix_gmpes(config, mtxs_type):
"""
Compute matrix of median ground-motion predictions for each gmpe for the
given run configuration for use within the Sammon maps and hierarchical
clustering dendrograms and Euclidean distance matrix plots.
If any gmpe logic trees are specified in the .toml, then these weights are
used to compute the associated gmpe logic tree (i.e. we can compare not
only gmpes, but which gmpes are most similar to the weighted logic tree
of them too).
:param mtxs_type:
type of predicted ground-motion matrix being computed in
compute_matrix_gmpes (either median, 84th or 16th percentile)
"""
# Get lt weights
lts = {gmc: getattr(config, config.lt_mapping[gmc]['wei'])
for gmc in config.lt_mapping.keys()}
mtxs_median = {}
for i, imt in enumerate(config.imt_list): # Iterate through imt_list
# Dict for storing medians
matrix_medians = np.zeros(
(len(config.gmpes_list),
(len(config.mags_eucl)*int((config.maxR-config.minR)/1))))
# Need to also store GMM LT weighted medians
lt_preds = {
lt: {gm: [] for gm in getattr(config, config.lt_mapping[lt]['wei'])}
for lt in lts if lts[lt] is not None
}
# Iterate over the GMMs
for g, gmpe in enumerate(config.gmpes_list):
# If the GMM is in a logic tree then get weight and LT
if 'lt_weight_gmc' in gmpe:
lt_ini = gmpe.split("lt_weight_gmc")[1]
if 'plot_lt_only' in gmpe:
lt = int(lt_ini.split("_plot_lt_only")[0])
else:
lt = int(lt_ini.split("=")[0])
lt_key = f"lt_gmc_{lt}"
assert lt_key in lt_preds.keys() # Sanity check
wt = getattr(config, f"lt_weight_gmc{lt}")[gmpe]
else:
wt = None
medians = []
for m, mag in enumerate(config.mags_eucl): # Iterate though mags
# Perform gmpe check
gmm = gmpe_check(gmpe)
# Get depth param
depth_g = config.depths_eucl[m]
ztor_g = None # NOTE: No hypo depth constraint used here
# Get rupture params
strike_g, dip_g, aratio_g = get_rup_pars(config.strike,
config.dip,
config.rake,
config.aratio,
config.trt)
mean, std, r_vals, tau, phi = att_curves(gmm,
mag,
config.lon,
config.lat,
depth_g,
ztor_g,
aratio_g,
strike_g,
dip_g,
config.rake,
config.trt,
config.rup,
config.vs30,
config.z1pt0,
config.z2pt5,
config.maxR,
1, # Step of 1 km for site spacing
imt,
config.dist_type,
config.up_or_down_dip,
config.volc_back_arc,
config.eshm20_region)
# Get means further than minR
idx = np.argwhere(r_vals>=config.minR).flatten()
mean = [mean[0][0][idx]]
std = [std[0][0][idx]]
tau = [tau[0][0][idx]]
phi = [phi[0][0][idx]]
# If plotting percentiles check GMM has sigma model
if mtxs_type in ["16th_perc", "18th_perc"] and np.all(std[0]==0):
raise ValueError(f"Cannot perform Euclidean analysis for a "
f"GMPE which lacks a sigma model (GMPE {g+1})")
# Store required percentile of ground-shaking
if mtxs_type == 'median':
preds = (np.exp(mean))
elif mtxs_type == '84th_perc':
nstd = 1 # Median + 1std = ~84th percentile
preds = (np.exp(mean+nstd*std[0]))
else:
assert mtxs_type == '16th_perc'
nstd = 1 # Median - 1std = ~16th percentile
preds = (np.exp(mean-nstd*std[0]))
medians = np.append(medians, preds)
# Store weighted median if gmm in an lt
if wt is not None:
lt_preds[lt_key][gmpe] = np.append(lt_preds[lt_key][gmpe], preds*wt)
# Store medians for gmm for given mag
matrix_medians[:][g] = medians
# Store medians for given imt
mtxs_median[str(imt)] = matrix_medians
# Get any required weighted means now we have medians, for all mags, for each GMM
for gmm_lt in lt_preds.keys():
mtxs_median[f"{imt}_{gmm_lt}"] = pd.DataFrame(lt_preds[gmm_lt].values()).mean(axis=0)
# Store gmpes_list to
mtxs_median['gmpe_list'] = config.gmpes_list.copy()
return mtxs_median
[docs]
def plot_sammons_util(imt_list,
gmpe_list,
mtxs,
namefig,
custom_color_flag,
custom_color_list,
mtxs_type):
"""
Plot Sammon maps for given run configuration. The weighted
mean of the GMPE predictions is plotted if GMM logic tree
weights are specified.
:param imt_list:
A list e.g. ['PGA', 'SA(0.1)', 'SA(1.0)']
:param gmpe_list:
A list e.g. ['BooreEtAl2014', 'CauzziEtAl2014']
:param mtxs:
Matrix of predicted ground-motion for each gmpe per imt
:param namefig:
filename for outputted figure
:param mtxs_type:
type of predicted ground-motion matrix being computed in
compute_matrix_gmpes (either median or 84th or 16th percentile)
"""
# Setup
colors = get_colors(custom_color_flag, custom_color_list)
texts = []
if len(imt_list) < 3:
nrows = 1
else:
nrows = int(np.ceil(len(imt_list)/2))
fig = pyplot.figure()
fig.set_size_inches(12, 6*nrows)
coo_per_imt = {}
for i, imt in enumerate(imt_list):
# Get the data matrix
data = mtxs[imt]
# gmm labels and configs
labels = gmpe_list.copy()
gmm_configs = mtxs['gmpe_list'].copy()
# Add the weighted LTs if any too
for key in mtxs.keys():
check = f"{imt}_lt_gmc"
if check in key:
data = np.vstack((data, mtxs[key]))
labels.append(key.split(f"{imt}_")[1]) # Add label for the gmc
gmm_configs.append(check)
# If only need gmm LT drop the gmms included in it
keep = np.array(['plot_lt_only' not in gmm for gmm in gmm_configs])
data = data[keep]
labels = [gmm for k, gmm in zip(keep, labels) if k]
# Sammon mapping
coo, cost = sammon(data, display=1) # NOTE: each gmm's array in coo has a structure of
coo_per_imt[imt] = coo # of [idx1, idx2, dist, npoints] where idx1 and idx2
fig.add_subplot(nrows, 2, i+1) # are merged at distance of dist into a cluster which
for g, gmpe in enumerate(labels): # containing npoints points
# Get colors and marker
if 'gmcLT' in gmpe:
marker = 'x'
else:
marker = 'o'
col = colors[g]
# Plot data
pyplot.plot(coo[g, 0], coo[g, 1], marker, markersize=9, color=col, label=gmpe)
texts.append(pyplot.text(coo[g, 0]+np.abs(coo[g, 0])*0.02,
coo[g, 1]+np.abs(coo[g, 1])*0.02,
labels[g],
ha='left',
color=col))
# Format plot
pyplot.title(str(imt), fontsize='16')
if mtxs_type == 'median':
pyplot.title(str(imt) + ' (median)', fontsize='14')
elif mtxs_type == '84th_perc':
pyplot.title(str(imt) + ' (84th percentile)', fontsize='14')
else:
assert mtxs_type == '16th_perc'
pyplot.title(str(imt) + ' (16th percentile)', fontsize='14')
pyplot.grid(axis='both', which='both', alpha=0.5)
# Tidy and save
pyplot.savefig(namefig, bbox_inches='tight', dpi=200, pad_inches=0.2)
pyplot.tight_layout()
pyplot.close()
return coo_per_imt
[docs]
def plot_cluster_util(imt_list, gmpe_list, mtxs, namefig, mtxs_type):
"""
Plot hierarchical clusters for given run configuration. The weighted
mean of the GMPE predictions is plotted if GMM logic tree weights are
specified.
:param imt_list:
A list e.g. ['PGA', 'SA(0.1)', 'SA(1.0)']
:param gmpe_list:
A list e.g. ['BooreEtAl2014', 'CauzziEtAl2014']
:param mtxs:
Matrix of predicted ground-motion for each gmpe per imt
:param namefig:
filename for outputted figure
:param mtxs_type:
type of predicted ground-motion matrix being computed in
compute_matrix_gmpes (either median or 84th or 16th percentile)
"""
# Setup
ncols = 2
if len(imt_list) < 3:
nrows = 1
else:
nrows = int(np.ceil(len(imt_list) / 2))
matrix_z = {}
ymax = [0] * len(imt_list)
# Loop over IMTs
for i, imt in enumerate(imt_list):
# Get the data matrix
data = mtxs[imt]
# gmm labels and configs
labels = gmpe_list.copy()
gmm_configs = mtxs['gmpe_list'].copy()
# Add the weighted LTs if any too
for key in mtxs.keys():
check = f"{imt}_lt_gmc"
if check in key:
data = np.vstack((data, mtxs[key]))
labels.append(key.split(f"{imt}_")[1]) # Add label for LT
gmm_configs.append(check)
# If only need gmm LT drop the gmms included in it
keep = np.array(['plot_lt_only' not in gmm for gmm in gmm_configs])
data = data[keep]
labels = [gmm for k, gmm in zip(keep, labels) if k]
# Agglomerative clustering
Z = hierarchy.linkage(
data, method='ward', metric='euclidean', optimal_ordering=True)
matrix_z[imt] = Z
ymax[i] = Z.max(axis=0)[2]
# Create the figure
fig, axs = pyplot.subplots(nrows, ncols)
fig.set_size_inches(12, 6*nrows)
for i, imt in enumerate(imt_list):
if len(imt_list) < 3:
ax = axs[i]
else:
ax = axs[np.unravel_index(i, (nrows, ncols))]
# Plot dendrogram
dn1 = hierarchy.dendrogram(
matrix_z[imt], ax=ax, orientation='right', labels=labels)
ax.set_xlabel('Euclidean Distance', fontsize='12')
if mtxs_type == 'median':
ax.set_title(str(imt) + ' (median)', fontsize='12')
elif mtxs_type == '84th_perc':
ax.set_title(str(imt) + ' (84th percentile)', fontsize='12')
else:
assert mtxs_type == '16th_perc'
ax.set_title(str(imt) + ' (16th percentile)', fontsize='12')
# Remove final plot if not required
if len(imt_list) >= 3 and len(imt_list)/2 != int(len(imt_list)/2):
ax = axs[np.unravel_index(i+1, (nrows, ncols))]
ax.set_visible(False)
if len(imt_list) == 1:
axs[1].set_visible(False)
# Save
pyplot.savefig(namefig, bbox_inches='tight', dpi=200, pad_inches=0.5)
pyplot.tight_layout()
pyplot.close()
return matrix_z
[docs]
def plot_matrix_util(imt_list, gmpe_list, mtxs, namefig, mtxs_type):
"""
Plot Euclidean distance matrices for given run configuration.
:param imt_list:
A list e.g. ['PGA', 'SA(0.1)', 'SA(1.0)']
:param gmpe_list:
A list e.g. ['BooreEtAl2014', 'CauzziEtAl2014']
:param mtxs:
Matrix of predicted ground-motion for each gmpe per imt
:param namefig:
filename for outputted figure.
:param mtxs_type:
type of predicted ground-motion matrix being computed in
compute_matrix_gmpes (either median or 84th or 16th percentile)
"""
# Euclidean
matrix_dist = {}
# Loop over IMTs
for i, imt in enumerate(imt_list):
# Get the data matrix
data = mtxs[imt]
# gmm labels and configs
labels = gmpe_list.copy()
gmm_configs = mtxs['gmpe_list'].copy()
# Add the weighted LTs if any too
for key in mtxs.keys():
check = f"{imt}_lt_gmc"
if check in key:
data = np.vstack((data, mtxs[key]))
labels.append(key.split(f"{imt}_")[1]) # Add label
gmm_configs.append(check)
# If only need gmm LT drop the gmms included in it
keep = np.array(['plot_lt_only' not in gmm for gmm in gmm_configs])
data = data[keep]
labels = [gmm for k, gmm in zip(keep, labels) if k]
# Agglomerative clustering
dist = squareform(pdist(data, 'euclidean'))
matrix_dist[imt] = dist
# Create the figure
ncols = 2
if len(imt_list) < 3:
nrows = 1
else:
nrows = int(np.ceil(len(imt_list) / 2))
fig, axs = pyplot.subplots(nrows, ncols)
fig.set_size_inches(12, 6*nrows)
for i, imt in enumerate(imt_list):
if len(imt_list) < 3:
ax = axs[i]
else:
ax = axs[np.unravel_index(i, (nrows, ncols))]
ax.imshow(matrix_dist[imt], cmap='gray')
# Add title
if mtxs_type == 'median':
ax.set_title(str(imt) + ' (median)', fontsize='14')
elif mtxs_type == '84th_perc':
ax.set_title(str(imt) + ' (84th percentile)', fontsize='14')
else:
assert mtxs_type == '16th_perc'
ax.set_title(str(imt) + ' (16th percentile)', fontsize='14')
# Add axis ticks
ax.xaxis.set_ticks([n for n in range(len(labels))])
ax.xaxis.set_ticklabels(labels, rotation=40)
ax.yaxis.set_ticks([n for n in range(len(labels))])
ax.yaxis.set_ticklabels(labels)
# Remove final plot if not required
if len(imt_list) >= 3 and len(imt_list)/2 != int(len(imt_list)/2):
ax = axs[np.unravel_index(i+1, (nrows, ncols))]
ax.set_visible(False)
# Save
pyplot.savefig(namefig, bbox_inches='tight', dpi=200, pad_inches=0.2)
pyplot.tight_layout()
pyplot.close()
return matrix_dist
### Utils for plots
[docs]
def get_colors(custom_color_flag, custom_color_list):
"""
Get list of colors for plots.
"""
if custom_color_flag is True:
return custom_color_list
else:
return COLORS
### Trellis Utils ###
[docs]
def trellis_data(gmpe,
r_vals,
mean,
add_sigma,
min_sigma,
col,
nstd,
lt_vals_gmc,
lt_weights):
"""
Plot predictions of a single GMPE (if required) and compute weighted
predictions from logic tree(s) (again if required).
"""
# If plotting not only the logic trees, plot each GMPE
if 'plot_lt_only' not in str(gmpe):
pyplot.plot(
r_vals, np.exp(mean), color = col, linewidth=2, linestyle='-', label=clean_gmm_label(gmpe))
# Plot mean with plus/minus sigma too if required
if nstd > 0:
pyplot.plot(r_vals, add_sigma, linewidth=0.75, color=col, linestyle='-.')
pyplot.plot(r_vals, min_sigma, linewidth=0.75, color=col, linestyle='-.')
# Now compute the weighted logic trees
for gmc in lt_vals_gmc.keys():
if lt_weights[gmc] is None:
pass
elif gmpe in lt_weights[gmc]:
if lt_weights[gmc][gmpe] is not None:
if nstd > 0:
lt_vals_gmc[gmc][gmpe] = {
'median': np.exp(mean)*lt_weights[gmc][gmpe],
'add_sigma': add_sigma*lt_weights[gmc][gmpe],
'min_sigma': min_sigma*lt_weights[gmc][gmpe]
}
else:
lt_vals_gmc[gmc][
gmpe] = {'median': np.exp(mean)*lt_weights[gmc][gmpe]}
return lt_vals_gmc
[docs]
def trellis_logic_trees(config,
key_gmc,
gmc,
lt_vals_gmc,
gmc_p,
store_gmm_curves,
r_vals,
nstd,
i,
m,
dep,
dip,
rake,
cfg_key,
unit):
"""
Manages plotting of the logic tree attenuation curves and
adds them to the store of exported attenuation curves.
"""
# If logic tree provided plot and add to attenuation curve store
if gmc is not None:
median, plus_sig, minus_sig = lt_trellis_plot(config,
r_vals,
nstd,
i,
m,
dep,
dip,
rake,
key_gmc,
lt_vals_gmc,
gmc_p[0],
gmc_p[1],
gmc_p[2])
store_gmm_curves[cfg_key][
'gmc logic tree curves per imt-mag'][key_gmc] = {}
store_gmm_curves[cfg_key][
'gmc logic tree curves per imt-mag'][key_gmc]['median (%s)' % unit] = median
if nstd > 0:
store_gmm_curves[
cfg_key]['gmc logic tree curves per imt-mag'][
key_gmc]['median plus sigma (%s)' % unit] = plus_sig
store_gmm_curves[
cfg_key]['gmc logic tree curves per imt-mag'][
key_gmc]['median minus sigma (%s)' % unit] = minus_sig
return store_gmm_curves
[docs]
def lt_trellis_plot(config,
r_vals,
nstd,
i,
m,
dep,
dip,
rake,
key_gmc,
lt_vals_gmc,
median_gmc,
plus_sig_gmc,
minus_sig_gmc):
"""
If required plot trellis from the given GMPE logic tree.
"""
# Get key describing mag-imt combo and some other event info
mk = (f'IMT = {i}, Mw = {m}, depth = {dep} km, dip = {dip} deg, rake = {rake} deg')
# Get logic tree
lt_df_gmc = pd.DataFrame(lt_vals_gmc, index=['median', 'add_sigma', 'min_sigma'])
lt_median = lt_df_gmc.loc['median'].sum()
median_gmc[mk] = lt_median
pyplot.plot(r_vals,
lt_median,
linewidth=2,
color=config.lt_mapping[key_gmc]["col"],
linestyle='--',
label=config.lt_mapping[key_gmc]['label'],
zorder=100)
if nstd > 0:
lt_add = lt_df_gmc.loc['add_sigma'].sum()
lt_min = lt_df_gmc.loc['min_sigma'].sum()
plus_sig_gmc[mk] = lt_add
minus_sig_gmc[mk] = lt_min
# Plot both plus and minus sigma curves
for sigma_val in [lt_add, lt_min]:
pyplot.plot(r_vals,
sigma_val,
linewidth=0.75,
color=config.lt_mapping[key_gmc]["col"],
linestyle='-.',
zorder=100)
return median_gmc, plus_sig_gmc, minus_sig_gmc
[docs]
def update_trellis_plots(mag, imt, m, i, dep, vs30, minR, maxR, r_vals, imt_list, dist_type):
"""
Add titles, axis labels and axis limits to trellis plots.
"""
# Get distance type label
dt_label = get_dist_label(dist_type)
# Bottom row only
if i == len(imt_list)-1:
pyplot.xlabel(dt_label, fontsize='16')
# Top row only
if i == 0:
pyplot.title(f'Mw={mag}, depth={dep}km, vs30={vs30}m/s', fontsize='12')
# Left row only
if m == 0:
if str(imt) in ['PGD', 'SDi']:
pyplot.ylabel(str(imt) + ' (cm)', fontsize='16')
elif str(imt) in ['PGV']:
pyplot.ylabel(str(imt) + ' (cm/s)', fontsize='16')
elif str(imt) in ['IA']:
pyplot.ylabel(str(imt) + ' (m/s)', fontsize='16')
elif str(imt) in ['RSD', 'RSD595', 'RSD575', 'RSD2080', 'DRVT']:
pyplot.ylabel(str(imt) + ' (s)', fontsize='16')
elif str(imt) in ['CAV']:
pyplot.ylabel(str(imt) + ' (g-sec)', fontsize='16')
elif str(imt) in ['MMI']:
pyplot.ylabel(str(imt) + ' (MMI)', fontsize='16')
elif str(imt) in ['FAS', 'EAS']:
pyplot.ylabel(str(imt) + ' (Hz)')
else:
pyplot.ylabel(str(imt) + ' (g)', fontsize='16') # PGA, SA, AvgSA
# xlims (manage because if rrup or rjb will be dependent on finiteness of rupture)
min_r_val = min(r_vals[r_vals>=1])
pyplot.xlim(np.max([min_r_val, minR]), maxR)
# And make loglog
pyplot.loglog()
[docs]
def filter_flatfile_trellis(data, imt, mag, depth, vs30, dist_type):
"""
Filter the dataframe of the provided flatfile for the given imt,
magnitude, focal depth and vs30 for use in trellis plotting.
NOTE: We return RotD50 values which have consistency with OQ units.
"""
# Filter first by magnitude, depth and vs30 first
subset = filter_flatfile(data, mag, depth, vs30, dist_type)
# Check there are values for the given IMT
if imt not in GEM_FF_MAPPINGS.keys():
# Might not be a column with RotD50 values for this IMT
raise ValueError(f'"{imt}" is not an IMT supported in the GEM Global Flatfile.')
imt_col = GEM_FF_MAPPINGS[imt]["col"]
subset = subset.loc[subset[imt_col].notnull()].reset_index(drop=True)
# Convert from flatfile units to those of GMPEs in OQ for given IMT
subset[imt_col] = subset[imt_col] * GEM_FF_MAPPINGS[imt]["conv_factor"]
# End of flatfile filtering
if len(subset) > 0:
return subset
else:
return None
### Spectra Utils ###
[docs]
def spectra_data(gmpe,
nstd,
gmc_weights,
rs_50p,
rs_add_sigma,
rs_min_sigma,
lt_vals,
sk):
"""
Store the spectra for given GMM and if required handle the gmpe
logic trees.
"""
# Store the non-weighted spectra values for given gmm
lt_vals['med'][gmpe][sk] = rs_50p
if nstd > 0:
lt_vals['add'][gmpe][sk] = rs_add_sigma
lt_vals['min'][gmpe][sk] = rs_min_sigma
# Handle the LTs
for gmc in gmc_weights:
if gmc_weights[gmc] is None:
continue
elif gmpe in gmc_weights[gmc]:
if gmc_weights[gmc][gmpe] is not None:
rs_50p_w = np.zeros(len(rs_50p))
rs_add_sigma_w = np.zeros(len(rs_add_sigma))
rs_min_sigma_w = np.zeros(len(rs_min_sigma))
for idx, rs in enumerate(rs_50p):
rs_50p_w[idx] = rs*gmc_weights[gmc][gmpe]
if nstd > 0:
rs_add_sigma_w[idx] = rs_add_sigma[idx]*gmc_weights[gmc][gmpe]
rs_min_sigma_w[idx] = rs_min_sigma[idx]*gmc_weights[gmc][gmpe]
# Store the weighted median for the gmm
lt_vals['med_wei'][gmc][gmpe][sk] = rs_50p_w
# And if nstd > 0 store these weighted branches too
if nstd > 0:
lt_vals['add_wei'][gmc][gmpe][sk] = rs_add_sigma_w
lt_vals['min_wei'][gmc][gmpe][sk] = rs_min_sigma_w
return {
'lt_gmc_1': [lt_vals['med_wei']['lt_gmc_1'],
lt_vals['add_wei']['lt_gmc_1'],
lt_vals['min_wei']['lt_gmc_1']],
'lt_gmc_2': [lt_vals['med_wei']['lt_gmc_2'],
lt_vals['add_wei']['lt_gmc_2'],
lt_vals['min_wei']['lt_gmc_2']],
'lt_gmc_3': [lt_vals['med_wei']['lt_gmc_3'],
lt_vals['add_wei']['lt_gmc_3'],
lt_vals['min_wei']['lt_gmc_3']],
'lt_gmc_4': [lt_vals['med_wei']['lt_gmc_4'],
lt_vals['add_wei']['lt_gmc_4'],
lt_vals['min_wei']['lt_gmc_4']]
}
[docs]
def spectra_logic_trees(config,
ax,
gmpe_list,
nstd,
period,
key_gmc,
ltv,
sk):
"""
Manages plotting and handling of the spectra for each logic tree.
"""
# Get identifier for given GMC in the toml GMMs
check = f'lt_weight_gmc{key_gmc.split("lt_gmc_")[1]}'
# Store medians
wt_per_gmpe_gmc = {
gmpe: ltv[0][gmpe][sk] for gmpe in gmpe_list if check in str(gmpe)
}
lt_median = pd.DataFrame(wt_per_gmpe_gmc, index=period).sum(axis=1).to_dict()
# And plus/minus sigmas too if required
if nstd > 0:
wt_add_sig = {
gmpe: ltv[1][gmpe][sk] for gmpe in gmpe_list if check in str(gmpe)
}
wt_min_sig = {
gmpe: ltv[2][gmpe][sk] for gmpe in gmpe_list if check in str(gmpe)
}
lt_add_sig = pd.DataFrame.from_dict(wt_add_sig, orient='index').sum(axis=0).to_dict()
lt_min_sig = pd.DataFrame.from_dict(wt_min_sig, orient='index').sum(axis=0).to_dict()
else:
lt_add_sig = {}
lt_min_sig = {}
# Plot median logic tree
ax.plot(period,
list(lt_median.values()),
linewidth=2,
color=config.lt_mapping[key_gmc]["col"],
linestyle='--',
label=config.lt_mapping[key_gmc]['label'],
zorder=100)
# Plot plus sigma and minus sigma if required
if nstd > 0:
# Plus sigma
ax.plot(period,
list(lt_add_sig.values()),
linewidth=0.75,
color=config.lt_mapping[key_gmc]["col"],
linestyle='-.',
zorder=100)
# Minus sigma
ax.plot(period,
list(lt_min_sig.values()),
linewidth=0.75,
color=config.lt_mapping[key_gmc]["col"],
linestyle='-.',
zorder=100)
return [lt_median, lt_add_sig, lt_min_sig]
[docs]
def load_obs_spectra(obs_spectra_fname):
"""
If an obs spectra file has been specified get values from the csv
for comparison of observed spectra and spectra computed using GMPE
predictions.
Returns the spectra as a dataframe, the max period of the spectra,
the earthquake ID and the station ID.
"""
# Load the obs spectra
obs_spectra = pd.read_csv(obs_spectra_fname)
# Get values from obs_spectra dataframe...
eq_id = str(obs_spectra['EQ ID'].iloc[0])
st_id = str(obs_spectra['Station Code'].iloc[0])
max_period = obs_spectra['Period (s)'].max()
return obs_spectra, max_period, eq_id, st_id
[docs]
def plot_obs_spectra(ax,
obs_spectra,
g,
gmpe_list,
mag_list,
dep_list,
dist_list,
dist_type,
vs30,
eq_id,
st_id):
"""
Check if an observed spectra must be plotted, and if so plot.
"""
# Plot an observed spectra if inputted...
if obs_spectra is not None and g == len(gmpe_list)-1:
# Get rup params
mw = np.asarray(mag_list, float)[0]
dist = np.asarray(dist_list, float)[0]
depth = np.asarray(dep_list, float)[0]
# Get label for spectra plot
obs_string = (
f"{eq_id}\nrecorded at {st_id} ({dist_type}={dist}km, "
f"\nMw={mw}, depth={depth}km, vs30={vs30}m/s)"
)
# Plot the observed spectra
ax.plot(obs_spectra['Period (s)'],
obs_spectra['SA (g)'],
color='k',
linewidth=3,
linestyle='-',
label=obs_string)
[docs]
def update_spectra_plots(ax, mag, depth_g, dist, vs30, d, m, dist_list, dist_type):
"""
Add titles and axis labels to spectra.
"""
# Title
pyplot.title(
f'Mw={mag}, depth={depth_g}km, {dist_type}={dist}km, vs30={vs30}m/s', fontsize=9.5)
# Bottom row only
if d == len(dist_list)-1:
ax.set_xlabel('Period (s)', fontsize=16)
# Left column only
if m == 0:
ax.set_ylabel('SA (g)', fontsize=16)
[docs]
def raise_spectra_dist_error(dist, dist_type, r_vals):
"""
If there is an interpolation issue when computing the spectra
from the attenuation curves (mean for given distance value, for
given distance metric), then notify the user by raising an error.
"""
rtype = dist_type
assert rtype not in ["repi"] # Should not be interp issues for repi
r_min = int(r_vals.min())
r_max = int(r_vals.max())
raise ValueError(f"Requested spectra distance ({rtype} = {dist} km) is "
f"outside of {rtype} value range for this ground-"
f"shaking scenario (min = {r_min} km, max = {r_max} km)")
[docs]
def filter_flatfile_spectra(data, imts, mag, depth, vs30, dist, dist_type):
"""
Filter the dataframe of the provided flatfile for the given imt,
magnitude, focal depth, vs30 AND distance type for use in spectra
plotting.
NOTE: We return RotD50 values which have consistency with OQ units.
"""
# Filter first by magnitude, depth and vs30 first
subset = filter_flatfile(data, mag, depth, vs30, dist_type)
# Filter by distance (smaller window when closer to source)
if dist <= 50:
dlim = DIST_LIM_LOW
elif dist > 50 and dist <= 100:
dlim = DIST_LIM_MID
else:
dlim = DIST_LIM_MAX
dcol = GEM_FF_MAPPINGS[dist_type]
subset = subset.loc[
subset[dcol].between(dist - dlim, dist + dlim)].reset_index(drop=True)
# Get the column in flatfile corresponding to each period
imt_cols = [GEM_FF_MAPPINGS[imt.string]["col"] for imt in imts
if imt.string in GEM_FF_MAPPINGS]
# Make an array of each record's spectra
subset["spectra_rotD50"] = pd.Series()
subset["spectra_periods"] = pd.Series()
for idx_rec, rec in subset.iterrows():
# Get spectra in correct units
spectra = np.array([rec[col] for col in imt_cols]
) * GEM_FF_MAPPINGS["PGA"]["conv_factor"]
# Get periods for IMTs in flatfile
periods = np.array([imt.period for imt in imts if
imt.string in GEM_FF_MAPPINGS])
# For each record build spectra with conversion to units of g
mask = ~np.isnan(spectra)
subset.at[idx_rec, "spectra_rotD50"] = spectra[mask]
subset.at[idx_rec, "spectra_periods"] = periods[mask]
# End of flatfile filtering
if len(subset) > 0:
return subset
else:
return None
### Utils for Other Plots ###
[docs]
def get_dist_label(dist_type):
"""
Return string representing required distance type.
"""
if dist_type == 'repi':
return 'Repi (km)'
elif dist_type == 'rrup':
return 'Rrup (km)'
elif dist_type == 'rjb':
return 'Rjb (km)'
else:
assert dist_type == 'rhypo'
return 'Rhypo (km)'
[docs]
def update_ratio_plots(mag, imt, m, i, dep, vs30, minR, maxR, r_vals, imt_list, dist_type):
"""
Add titles and axis labels to ratio plots.
"""
# Get distance type label
dt_label = get_dist_label(dist_type)
# Bottom row only
if i == len(imt_list)-1:
pyplot.xlabel(dt_label, fontsize='12')
# Top row only
if i == 0:
pyplot.title(f'Mw={mag}, depth={dep}km, vs30={vs30}m/s', fontsize='12')
# Left row only
if m == 0:
pyplot.ylabel('GMM/baseline for %s' %str(imt), fontsize='14')
# Set xlims
min_r_val = min(r_vals[r_vals>=1])
pyplot.xlim(np.max([min_r_val, minR]), maxR)
[docs]
def filter_flatfile(data, mag, depth, vs30, dist_type):
"""
Filter by mag, depth, vs30 and distance type.
NOTE: Used for filtering of a provided GEM format flatfile for
both trellis and spectra plotting.
"""
# Add rhypo dist
if dist_type == "rhypo":
data["rhypo_dist"] = np.sqrt(
data["epi_dist"]**2 + data["ev_depth_km"]**2)
# Filter by magnitude, depth and vs30
subset = data.loc[
(data.Mw.between(mag - MAG_LIM, mag + MAG_LIM)) &
(data.ev_depth_km.between(depth - DEP_LIM, depth + DEP_LIM)) &
(data.vs30_m_sec.between(vs30 - VS30_LIM, vs30 + VS30_LIM))
].reset_index(drop=True)
return subset