diff --git a/docsrc/contents/smt.rst b/docsrc/contents/smt.rst index f6b1eb8be..d6bf2fb3c 100644 --- a/docsrc/contents/smt.rst +++ b/docsrc/contents/smt.rst @@ -96,7 +96,7 @@ We can specify the inputs to perform a residual analysis with as follows: 3. We can also specify the GMPEs and intensity measures within a ``.toml`` file. The ``.toml`` file method is required for specifying the inputs of GMPEs with user-specifiable input parameters e.g. regionalisation parameter or logic tree branch parameters. Note that here the GMPEs listed in the first ``.toml`` file are not appropriate for our target region, but have been selected to demonstrate how GMPEs with additional inputs can be specified within a ``.toml`` file. The second ``.toml`` file provides the GMPEs and intensity measures we use for running this demonstration analysis. - The additional input parameters which are specifiable for certain GMPEs are available within their corresponding GSIM files (found in ``oq-engine.openquake.hazardlib.gsim``). Note also that a GMPE sigma model must be provided by the GMPE for the computation of residuals. If a sigma model is not provided by the GMPE, it can be specified as demonstrated below for the YenierAtkinson2015BSSA GMPE. + The additional input parameters which are specifiable for certain GMPEs are available within their corresponding GSIM files (found in ``oq-engine.openquake.hazardlib.gsim``, or for ModifiableGMPE features in ``oq-engine.openquake.hazardlib.gsim.mgmpe.modifiable_gmpe``). Note also that a GMPE sigma model must be provided by the GMPE for the computation of residuals. If a sigma model is not provided by the GMPE, it can be specified as demonstrated below. The ``.toml`` file for specifying GMPEs and intensity measures to consider within a residual analysis should be specified as follows: @@ -104,24 +104,51 @@ We can specify the inputs to perform a residual analysis with as follows: [models] - [models.AbrahamsonGulerce2020SInter] + [models.1-AbrahamsonGulerce2020SInter] region = "GLO" - [models.AbrahamsonGulerce2020SInter] + [models.2-AbrahamsonGulerce2020SInter] region = "CAS" - [models.AbrahamsonGulerce2020SInterCascadia] - [models.YenierAtkinson2015BSSA] - sigma_model = 'al_atik_2015_sigma' + sigma_model = 'al_atik_2015_sigma' # use Al Atik (2015) sigma model + + [models.1-CampbellBozorgnia2014] + fix_total_sigma = "{'PGA': 0.750, 'SA(0.1)': 0.800, 'SA(0.5)': 0.850}" # fix total sigma per imt + + [models.2-CampbellBozorgnia2014] + with_betw_ratio = 1.7 # add between-event and within-event sigma using ratio of 1.7 to partition total sigma + + [models.3-CampbellBozorgnia2014] + set_between_epsilon = 1.5 # set between-event epsilon (i.e. tau epsilon) + + [models.1-AbrahamsonEtAl2014] + median_scaling_scalar = 1.4 # scale median by factor of 1.4 over all imts + + [models.2-AbrahamsonEtAl2014] + median_scaling_vector = "{'PGA': 1.10, 'SA(0.1)': 1.15, 'SA(0.5)': 1.20}" # scale median by imt-dependent factor + + [models.1-KothaEtAl202] + sigma_scaling_scalar = 1.05 # scale sigma by factor of 1.05 over all imts + + [models.2-KothaEtAl2020] + sigma_scaling_vector = "{'PGA': 1.20, 'SA(0.1)': 1.15, 'SA(0.5)': 1.10}" # scale sigma by imt-dependent factor + + [models.1-BooreEtAl2014] + site_term = 'CY14SiteTerm' # use CY14 site term + + [models.2-BooreEtAl2014] + site_term = 'NRCan15SiteTerm' # use NRCan15 non-linear site term + + [models.3-BooreEtAl2014] + site_term = 'NRCan15SiteTermLinear' # use NRCan15 linear site term [models.NGAEastGMPE] - gmpe_table = 'NGAEast_FRANKEL_J15.hdf5' + gmpe_table = 'NGAEast_FRANKEL_J15.hdf5' # use a gmpe table [models.HassaniAtkinson2018] - d_sigma = 100 + d_sigma = 100 # gmpe specific param kappa0 = 0.04 - sigma_model = 'al_atik_2015_sigma' [imts] imt_list = ['PGA', 'SA(0.2)', 'SA(0.5)', 'SA(1.0'] @@ -423,7 +450,7 @@ Comparing GMPEs 2. The tools within the Comparison module include Sammon's Maps, hierarchical clustering and matrix plots of Euclidean distance for both median and 84th percentile of predicted ground-motion per GMPE per intensity measure. Plotting capabilities for response spectra, GMPE sigma with respect to spectral period and trellis plots are also provided in this module. - The inputs for these comparitive tools must be specified within a single ``.toml`` file as specified below. In the ``.toml`` file we have specified the source parameters for earthquakes characteristic of Albania (compressional thrust faulting with magnitudes of interest w.r.t. seismic hazard in the range of Mw 5 to Mw 7), and we have specified a selection of GMPEs which may increase the captured epistemic uncertainty associated with predicting the ground-shaking from earthquakes in/near Albania if implemented in a GMPE logic tree. To plot a GMPE logic tree we must assign model weights using ``lt_weight_gmc1`` or '``lt_weight_gmc2`` in each GMPE depending on if we want to plot the GMPE within GMC logic tree #1 or #2 (up to 2 GMC logic trees can currently be plotted within one trellis or response spectra plot at a time). To plot only the final logic tree and not the individual GMPEs comprising it, we use ``lt_weight_gmc1_plot_lt_only`` or ``lt_weight_gmc2_plot_lt_only`` instead (depending on which GMC we wish to not plot the individual GMPEs for - see the .toml file below for an example of these potential configurations). + The inputs for these comparitive tools must be specified within a single ``.toml`` file as specified below. GMPE parameters can be specified as within the example ``.toml`` file provided above for us in residual analysis. In the ``.toml`` file we have specified the source parameters for earthquakes characteristic of Albania (compressional thrust faulting with magnitudes of interest w.r.t. seismic hazard in the range of Mw 5 to Mw 7), and we have specified a selection of GMPEs which may increase the captured epistemic uncertainty associated with predicting the ground-shaking from earthquakes in/near Albania if implemented in a GMPE logic tree. To plot a GMPE logic tree we must assign model weights using ``lt_weight_gmc1`` or '``lt_weight_gmc2`` in each GMPE depending on if we want to plot the GMPE within GMC logic tree #1 or #2 (up to 2 GMC logic trees can currently be plotted within one trellis or response spectra plot at a time). To plot only the final logic tree and not the individual GMPEs comprising it, we use ``lt_weight_gmc1_plot_lt_only`` or ``lt_weight_gmc2_plot_lt_only`` instead (depending on which GMC we wish to not plot the individual GMPEs for - see the .toml file below for an example of these potential configurations). .. code-block:: ini @@ -447,11 +474,12 @@ Comparing GMPEs # Characterise earthquake for the region of interest as finite rupture [source_properties] - trt = 'ASCR' # Specify a TRT string from ASCR, InSlab, Interface, Stable, Upper_Mantle, Volcanic, Induced, Induced_Geothermal + trt = 'None' # Either string of 'None' to use user-provided aratio OR specify a TRT string from ASCR, InSlab, Interface, Stable, Upper_Mantle, Volcanic, Induced, Induced_Geothermal to assign a trt-dependent proxy aratio ztor = 'None' # Set to string of 'None' to NOT consider otherwise specify as array matching number of mag and depth values strike = -999 dip = 60 # (Albania has predominantly reverse faulting) rake = 90 # (+ 90 for compression, -90 for extension) + aratio = 2 # If set to -999 the user-provided trt string will be used to assign a default trt-dependent aratio trellis_and_rs_mag_list = [5, 6, 7] # mags used only for trellis and response spectra trellis_and_rs_depths = [20, 20, 20] # depth per magnitude for trellis and response spectra diff --git a/openquake/smt/comparison/compare_gmpes.py b/openquake/smt/comparison/compare_gmpes.py index 368cd295f..2fce5b19a 100644 --- a/openquake/smt/comparison/compare_gmpes.py +++ b/openquake/smt/comparison/compare_gmpes.py @@ -65,6 +65,8 @@ def __init__(self, filename): self.up_or_down_dip = float(up_or_down_dip) self.trt = config_file['source_properties']['trt'] + if self.trt == 'None': + self.trt = None self.ztor = config_file['source_properties']['ztor'] if self.ztor == 'None': self.ztor = None @@ -72,7 +74,7 @@ def __init__(self, filename): self.dip = config_file['source_properties']['dip'] self.rake = config_file['source_properties']['rake'] - self.aratio = -999 + self.aratio = config_file['source_properties']['aratio'] # One set of magnitudes for use in trellis plots self.trellis_and_rs_mag_list = config_file['source_properties'][ diff --git a/openquake/smt/comparison/utils_compare_gmpes.py b/openquake/smt/comparison/utils_compare_gmpes.py index 49795d428..fdae64270 100644 --- a/openquake/smt/comparison/utils_compare_gmpes.py +++ b/openquake/smt/comparison/utils_compare_gmpes.py @@ -74,11 +74,10 @@ def plot_trellis_util( strike, dip, depth[l], aratio, rake, trt) # Get attenuation curves - mean, std, r_vals = att_curves(gmm, gmm_orig, depth[l],m, - aratio_g, strike_g, dip_g, - rake,Vs30, Z1, Z25, maxR, - step, i, ztor_m, eshm20_region, - dist_type, trt, up_or_down_dip) + mean, std, r_vals, tau, phi = att_curves( + gmm, gmm_orig, depth[l], m, aratio_g, strike_g, dip_g, + rake,Vs30, Z1, Z25, maxR, step, i, ztor_m, eshm20_region, + dist_type, trt, up_or_down_dip) # Get mean, sigma, mean plus sigma and mean minus sigma mean = mean[0][0] @@ -169,7 +168,8 @@ def plot_spectra_util(trt, ztor, rake, strike, dip, depth, Z1, Z25, Vs30, aratio, rake, trt) - rs_50p, rs_plus_sigma, rs_minus_sigma, sigma = [], [], [], [] + rs_50p, rs_plus_sigma, rs_minus_sigma = [], [], [] + sigma_store, tau_store, phi_store = [], [], [] for k, imt in enumerate(imt_list): if obs_spectra is not None: @@ -188,18 +188,28 @@ def plot_spectra_util(trt, ztor, rake, strike, dip, depth, Z1, Z25, Vs30, ztor_m = None # Get mean and sigma - mu, std, r_vals = att_curves(gmm, gmm_orig, depth[l], m, - aratio_g, strike_g, dip_g, - rake, Vs30, Z1, Z25, dist, - 0.1, imt, ztor_m, eshm20_region, - dist_type, trt, up_or_down_dip) + mu, std, r_vals, tau, phi = att_curves( + gmm, gmm_orig, depth[l], m, aratio_g, strike_g, dip_g, + rake, Vs30, Z1, Z25, dist, 0.1, imt, ztor_m, eshm20_region, + dist_type, trt, up_or_down_dip) + + # Interpolate for distances and store mu = mu[0][0] f = interpolate.interp1d(r_vals, mu) rs_50p_dist = np.exp(f(i)) rs_50p.append(rs_50p_dist) + f1 = interpolate.interp1d(r_vals, std[0]) sigma_dist = f1(i) - sigma.append(sigma_dist) + sigma_store.append(sigma_dist) + + f2 = interpolate.interp1d(r_vals, phi[0]) + phi_dist = f2(i) + phi_store.append(phi_dist) + + f3 = interpolate.interp1d(r_vals, tau[0]) + tau_dist = f3(i) + tau_store.append(tau_dist) if Nstd != 0: rs_plus_sigma_dist = np.exp(f(i)+(Nstd*sigma_dist)) @@ -220,8 +230,12 @@ def plot_spectra_util(trt, ztor, rake, strike, dip, depth, Z1, Z25, Vs30, linewidth=0.75, linestyle='-.') # Plot sigma vs period - ax2.plot(period, sigma, color=col, linewidth=2, linestyle='-', - label=gmpe) + ax2.plot(period, sigma_store, color=col, linewidth=2, + linestyle='-', label='Total sigma \n' + gmpe) + ax2.plot(period, phi_store, color=col, linewidth=2, + linestyle='-.', label='Within-event \n' + gmpe) + ax2.plot(period, tau_store, color=col, linewidth=2, + linestyle='--', label='Between-event \n' + gmpe) # Weight the predictions using logic tree weights lt_vals_gmc1, lt_vals_plus_sig_gmc1, lt_vals_minus_sig_gmc1, \ @@ -241,7 +255,7 @@ def plot_spectra_util(trt, ztor, rake, strike, dip, depth, Z1, Z25, Vs30, # Set axis limits and add grid ax1.set_xlim(min(period), max(period)) - ax2.set_ylim(0.3, 1) + ax2.set_ylim(0, 1) ax1.grid(True) ax2.grid(True) @@ -299,11 +313,10 @@ def compute_matrix_gmpes(trt, ztor, imt_list, mag_list, gmpe_list, rake, strike, else: ztor_m = None - mean, std, r_vals = att_curves(gmm, gmm_orig, depth[l], m, - aratio_g, strike_g, dip_g, - rake, Vs30, Z1, Z25, maxR, - step, i, ztor_m, eshm20_region, - dist_type, trt, up_or_down_dip) + mean, std, r_vals, tau, phi = att_curves( + gmm, gmm_orig, depth[l], m, aratio_g, strike_g, dip_g, + rake, Vs30, Z1, Z25, maxR, step, i, ztor_m, eshm20_region, + dist_type, trt, up_or_down_dip) if mtxs_type == 'median': medians = np.append(medians, (np.exp(mean))) @@ -686,9 +699,11 @@ def update_trellis_plots(m, i, n, l, r_vals, imt_list, dist_type): if n == len(imt_list)-1: #bottom row only pyplot.xlabel(label, fontsize='16') if l == 0: #left row only - pyplot.ylabel(str(i) + ' (g)', fontsize='16') + if str(i) != 'PGV': + pyplot.ylabel(str(i) + ' (g)', fontsize='16') + else: + pyplot.ylabel('PGV (cm/s)', fontsize='16') pyplot.loglog() - pyplot.ylim(0.001, 10) # Mod if required pyplot.xlim(1, np.max(r_vals)) # Mod if required diff --git a/openquake/smt/comparison/utils_gmpes.py b/openquake/smt/comparison/utils_gmpes.py index 379345520..098478ac2 100644 --- a/openquake/smt/comparison/utils_gmpes.py +++ b/openquake/smt/comparison/utils_gmpes.py @@ -20,6 +20,7 @@ """ import numpy as np import pandas as pd +import ast from openquake.hazardlib.geo import Point from openquake.hazardlib.geo.surface import PlanarSurface @@ -152,12 +153,10 @@ def att_curves(gmpe, orig_gmpe, depth, mag, aratio, strike, dip, rake, Vs30, rup_trt = TRT.INDUCED if trt == 'Induced_Geothermal': rup_trt = TRT.GEOTHERMAL - if trt == -999: - rup_trt = gmpe.DEFINED_FOR_TECTONIC_REGION_TYPE - if rup_trt is None: - raise ValueError('Specify a TRT string within the toml file: ASCR, \ - InSlab, Interface, Stable, Upper_Mantle, Volcanic, \ - Induced, Induced_Geothermal') + if rup_trt is None and aratio == -999: + raise ValueError( + 'An ratio must be provided by the user, or alternatively specify a \ + TRT string within the toml file to assign a trt-dependent aratio proxy.') # Get rup rup = get_rupture(0.0, 0.0, depth, WC1994(), mag=mag, aratio=aratio, @@ -212,7 +211,7 @@ def att_curves(gmpe, orig_gmpe, depth, mag, aratio, strike, dip, rake, Vs30, distances[len(distances)-1] = maxR - return mean, std, distances + return mean, std, distances, tau, phi def _get_z1(Vs30, region): @@ -249,7 +248,8 @@ def _get_z25(Vs30, region): def _param_gmpes(strike, dip, depth, aratio, rake, trt): """ - Get proxies for strike, dip, depth and aspect ratio if not provided + Get (crude) proxies for strike, dip, depth and aspect ratio if not provided + by the user. """ # Strike if strike == -999: @@ -295,74 +295,104 @@ def mgmpe_check(gmpe): :param gmpe: gmpe: GMPE to be modified if required (must be a gsim class) """ - # Preserve original GMPE prior and create base version of GMPE orig_gmpe = gmpe - base_gsim = str(gmpe).splitlines()[0].replace('[', '').replace(']', '') + base_gsim = gmpe.__class__.__name__ # Get the additional params if specified inputs = pd.Series(str(gmpe).splitlines()[1:], dtype='object') add_inputs = {} add_as_int = ['eshm20_region'] add_as_str = ['region', 'gmpe_table', 'volc_arc_file'] - + if len(inputs) > 0: # If greater than 0 must add required gsim inputs idx_to_drop = [] for idx, par in enumerate(inputs): + par = str(par) # Drop mgmpe params from the other gsim inputs - if 'al_atik_2015_sigma' in str(par) or 'SiteTerm' in str(par): + if 'sigma_model' in par or 'site_term' in par: + idx_to_drop.append(idx) + if 'fix_total_sigma' in par: + idx_to_drop.append(idx) + base_vector = par.split('=')[1].replace('"','') + fixed_sigma_vector = ast.literal_eval(base_vector) + if 'with_betw_ratio' in par: + idx_to_drop.append(idx) + with_betw_ratio = float(par.split('=')[1]) + if 'set_between_epsilon' in par: idx_to_drop.append(idx) + between_epsilon = float(par.split('=')[1]) + if 'scaling' in par: + idx_to_drop.append(idx) + if 'median_scaling_scalar' in par: + median_scalar = float(par.split('=')[1]) + if 'median_scaling_vector' in par: + base_vector = par.split('=')[1].replace('"','') + median_vector = ast.literal_eval(base_vector) + if 'sigma_scaling_scalar' in par: + sigma_scalar = float(par.split('=')[1]) + if 'sigma_scaling_vector' in par: + base_vector = par.split('=')[1].replace('"','') + sigma_vector = ast.literal_eval(base_vector) + inputs = inputs.drop(np.array(idx_to_drop)) for idx, par in enumerate(inputs): - key = str(par).split('=')[0].strip() + key = par.split('=')[0].strip() if key in add_as_str: val = par.split('=')[1].replace('"', '').strip() elif key in add_as_int: val = int(par.split('=')[1]) else: - val = float(str(par).split('=')[1]) + val = float(par.split('=')[1]) add_inputs[key] = val - # Sigma model implementation - if 'al_atik_2015_sigma' in str(orig_gmpe): - params = {"tau_model": "global", "ergodic": False} - gmpe = mgmpe.ModifiableGMPE(gmpe={base_gsim: add_inputs}, - sigma_model_alatik2015=params) - - # Site term implementations - msg_sigma_and_site_term = 'An alternative sigma model and an alternative \ - site term cannot be specified within a single GMPE implementation.' - msg_multiple_site_terms = 'Two alternative site terms have been specified \ - within the toml for a single GMPE implementation' - - # Check only single site term specified - if ('CY14SiteTerm' in str(orig_gmpe) and - 'NRCan15SiteTerm' in str(orig_gmpe)): - raise ValueError(msg_multiple_site_terms) - - # Check if site term an dsigma model specified - if ('CY14SiteTerm' in str(orig_gmpe) and - 'al_atik_2015_sigma' in str(orig_gmpe) or - 'CY14SiteTerm' in str(orig_gmpe) and - 'al_atik_2015_sigma' in str(orig_gmpe)): - raise ValueError(msg_sigma_and_site_term) + kwargs = {'gmpe':{base_gsim: add_inputs}} # reconstruct the gmpe as kwargs + # Al Atik 2015 sigma model + if 'al_atik_2015_sigma' in str(orig_gmpe): + kwargs['sigma_model_alatik2015'] = {"tau_model": "global", "ergodic": False} + + # Fix total sigma per imt + if 'fix_total_sigma' in str(orig_gmpe): + kwargs['set_fixed_total_sigma'] = {'total_sigma':fixed_sigma_vector} + + # Partition total sigma of gsim using a specified ratio of within:between + if 'with_betw_ratio' in str(orig_gmpe): + kwargs['add_between_within_stds'] = {'with_betw_ratio': with_betw_ratio} + + # Set epsilon for tau + if 'set_between_epsilon' in str(orig_gmpe): + kwargs['set_between_epsilon'] = {'epsilon_tau':between_epsilon} + + # Scale median by constant factor over all imts + if 'median_scaling_scalar' in str(orig_gmpe): + kwargs['set_scale_median_scalar'] = {'scaling_factor':median_scalar} + + # Scale median by imt-dependent factor + if 'median_scaling_vector' in str(orig_gmpe): + kwargs['set_scale_median_vector'] = {'scaling_factor':median_vector} + + # Scale sigma by constant factor over all imts + if 'sigma_scaling_scalar' in str(orig_gmpe): + kwargs['set_scale_total_sigma_scalar'] = {'scaling_factor': sigma_scalar} + + # Scale sigma by imt-dependent factor + if 'sigma_scaling_vector' in str(orig_gmpe): + kwargs['set_scale_total_sigma_vector'] = {'scaling_factor': sigma_vector} + # CY14SiteTerm if 'CY14SiteTerm' in str(orig_gmpe): - params = {} - gmpe = mgmpe.ModifiableGMPE(gmpe={base_gsim: add_inputs}, - cy14_site_term=params) - + kwargs['cy14_site_term'] = {} + # NRCan15SiteTerm (kind = base) if ('NRCan15SiteTerm' in str(orig_gmpe) and - 'NRCan15SiteTermLinear' not in str(orig_gmpe)): - params = {'kind': 'base'} - gmpe = mgmpe.ModifiableGMPE(gmpe={base_gsim: add_inputs}, - nrcan15_site_term=params) + 'NRCan15SiteTermLinear' not in str(orig_gmpe)): + kwargs['nrcan15_site_term'] = {'kind': 'base'} + # NRCan15SiteTerm (kind = linear) if 'NRCan15SiteTermLinear' in str(orig_gmpe): - params = {'kind': 'linear'} - gmpe = mgmpe.ModifiableGMPE(gmpe={base_gsim: add_inputs}, - nrcan15_site_term=params) + kwargs['nrcan15_site_term'] = {'kind': 'linear'} + + gmpe = mgmpe.ModifiableGMPE(**kwargs) # remake gmpe using mgmpe - return gmpe + return gmpe \ No newline at end of file diff --git a/openquake/smt/tests/comparison/data/compare_gmpe_inputs.toml b/openquake/smt/tests/comparison/data/compare_gmpe_inputs.toml index 7ee66c821..edc01c5c6 100644 --- a/openquake/smt/tests/comparison/data/compare_gmpe_inputs.toml +++ b/openquake/smt/tests/comparison/data/compare_gmpe_inputs.toml @@ -31,6 +31,7 @@ ztor = 'None' # Set to string of 'None' to NOT consider strike = -999 dip = -999 rake = 60 +aratio = 2 trellis_and_rs_mag_list = [5,6,7] # mags used only for trellis trellis_and_rs_depths = [20,25,30] # depth per magnitude for trellis diff --git a/openquake/smt/tests/file_samples/example_comparison_module_inputs.toml b/openquake/smt/tests/file_samples/example_comparison_module_inputs.toml index b0b763cfb..05d85fb1a 100644 --- a/openquake/smt/tests/file_samples/example_comparison_module_inputs.toml +++ b/openquake/smt/tests/file_samples/example_comparison_module_inputs.toml @@ -24,8 +24,9 @@ ztor = 'None' # Set to string of 'None' to NOT consider strike = -999 dip = 60 # (Albania has predominantly reverse faulting) rake = 90 # (+ 90 for compression, -90 for extension) -trellis_mag_list = [5, 6, 7] # mags used only for trellis and response spectra -trellis_depths = [20, 20, 20] # depth per magnitude and response spectra +aratio = 2 +trellis_and_rs_mag_list = [5, 6, 7] # mags used only for trellis and response spectra +trellis_and_rs_depths = [20, 20, 20] # depth per magnitude and response spectra # Specify magnitude array for Sammons, Euclidean dist and clustering [mag_values_non_trellis_or_spectra_functions]