diff --git a/docsrc/contents/smt.rst b/docsrc/contents/smt.rst index 8efa43c0b..88d8ef546 100644 --- a/docsrc/contents/smt.rst +++ b/docsrc/contents/smt.rst @@ -528,11 +528,13 @@ Comparing GMPEs [models.AkkarEtAlRjb2014] lt_weight_gmc2_plot_lt_only = 0.50 + # Also specify a GMM to compute ratios of the attenuation against (GMM/baseline) + [ratios_baseline_gmm.BooreEtAl2020] + [custom_colors] custom_colors_flag = 'False' # Set to "True" for custom colours in plots custom_colors_list = ['lime', 'dodgerblue', 'gold', '0.8'] - 3. Trellis Plots Now that we have defined our inputs for GMPE comparison, we can use each tool within the Comparison module to evaluate how similar the GMPEs predict ground-motion for a given ground-shaking scenario. @@ -573,10 +575,21 @@ Comparing GMPEs > comp.plot_spectra(filename, output_directory, obs_spectra='spectra_chamoli_1991_station_UKHI.csv') Response spectra plots for input parameters specified in toml file: - .. image:: /contents/smt_images/ObsSpectra.png - + .. image:: /contents/smt_images/ObsSpectra.png + +6. Plot of ratios of attenuation curves + + The ratios of the median predictions from each GMM and a baseline GMM (specified in the ``.toml`` - see above) can also be plotted. An example is provided in the demo files: + + .. code-block:: ini + + > # Plot ratios of median attenuation curves for each GMM/median attenuation curves for baseline GMM + > comp.plot_ratios(filename, output_directory) + + Ratio plots for input parameters specified in toml file (note that here the baseline GMM is BooreEtAl2014): + .. image:: /contents/smt_images/RatioPlots.png -6. Sammon's Maps +7. Sammon's Maps We can plot Sammon's Maps to examine how similar the medians (and 16th and 84th percentiles) of predicted ground-motion of each GMPE are (see Sammon, 1969 and Scherbaum et al. 2010 for more details on the Sammon's mapping procedure). @@ -592,7 +605,7 @@ Comparing GMPEs Sammon's Maps (median predicted ground-motion) for input parameters specified in toml file: .. image:: /contents/smt_images/Median_SammonMaps.png -7. Hierarchical Clustering +8. Hierarchical Clustering Dendrograms can be plotted as an alternative tool to evaluate how similarly the predicted ground-motion is by each GMPE. @@ -608,7 +621,7 @@ Comparing GMPEs Dendrograms (median predicted ground-motion) for input parameters specified in toml file: .. image:: /contents/smt_images/Median_Clustering.png -8. Matrix Plots of Euclidean Distance +9. Matrix Plots of Euclidean Distance In addition to Sammon's Maps and hierarchical clustering, we can also plot the Euclidean distance between the predicted ground-motions by each GMPE in a matrix plot. @@ -624,7 +637,7 @@ Comparing GMPEs Matrix plots of Euclidean distance between GMPEs (median predicted ground-motion) for input parameters specified in toml file: .. image:: /contents/smt_images/Median_Euclidean.png -9. Using ModifiableGMPE to modify GMPEs within a ``.toml``. +10. Using ModifiableGMPE to modify GMPEs within a ``.toml``. In addition to specifying predefined arguments for each GMPE, the user can also modify GMPEs using ModifiableGMPE (found in ``oq-engine.openquake.hazardlib.gsim.mgmpe.modifiable_gmpe``). diff --git a/docsrc/contents/smt_images/RatioPlots.png b/docsrc/contents/smt_images/RatioPlots.png new file mode 100644 index 000000000..fa3afbde6 Binary files /dev/null and b/docsrc/contents/smt_images/RatioPlots.png differ diff --git a/openquake/smt/comparison/compare_gmpes.py b/openquake/smt/comparison/compare_gmpes.py index 11d12da7d..1d27c42ff 100644 --- a/openquake/smt/comparison/compare_gmpes.py +++ b/openquake/smt/comparison/compare_gmpes.py @@ -28,8 +28,8 @@ from openquake.hazardlib.imt import from_string from openquake.smt.comparison.utils_compare_gmpes import ( - plot_trellis_util, plot_spectra_util, plot_cluster_util, plot_sammons_util, - plot_euclidean_util, compute_matrix_gmpes) + plot_trellis_util, plot_spectra_util, plot_ratios_util, plot_cluster_util, + plot_sammons_util, plot_euclidean_util, compute_matrix_gmpes) class Configurations(object): @@ -45,7 +45,7 @@ def __init__(self, filename): # Import parameters for comparative plots from .toml file config_file = toml.load(filename) - # Get input params from .toml file + # Get general params self.eshm20_region = config_file['general']['eshm20_region'] self.minR = config_file['general']['minR'] self.maxR = config_file['general']['maxR'] @@ -54,11 +54,7 @@ def __init__(self, filename): self.Nstd = config_file['general']['Nstd'] self.max_period = config_file['general']['max_period'] - self.custom_color_flag = config_file['custom_colors'][ - 'custom_colors_flag'] - self.custom_color_list = config_file['custom_colors'][ - 'custom_colors_list'] - + # Get site params self.Vs30 = config_file['site_properties']['vs30'] self.Z1 = config_file['site_properties']['Z1'] self.Z25 = config_file['site_properties']['Z25'] @@ -66,6 +62,7 @@ def __init__(self, filename): self.up_or_down_dip = float(up_or_down_dip) self.region = config_file['site_properties']['region'] + # Get rupture params self.trt = config_file['source_properties']['trt'] if self.trt == 'None': self.trt = None @@ -75,9 +72,14 @@ def __init__(self, filename): self.strike = config_file['source_properties']['strike'] self.dip = config_file['source_properties']['dip'] self.rake = config_file['source_properties']['rake'] - self.aratio = config_file['source_properties']['aratio'] + # Get custom colors + self.custom_color_flag = config_file['custom_colors'][ + 'custom_colors_flag'] + self.custom_color_list = config_file['custom_colors'][ + 'custom_colors_list'] + # One set of magnitudes for use in trellis plots self.trellis_and_rs_mag_list = config_file['source_properties'][ 'trellis_and_rs_mag_list'] @@ -106,87 +108,148 @@ def __init__(self, filename): for idx, imt in enumerate(self.imt_list): self.imt_list[idx] = from_string(str(self.imt_list[idx])) - # Get model labels - self.gmpe_labels = config_file['gmpe_labels']['gmpes_label'] + # Get models and model labels + (self.gmpe_labels, self.gmpes_list, self.baseline_gmm) = get_gmpes( + config_file) + + # Get lt weights + (self.lt_weights_gmc1, self.lt_weights_gmc2, self.lt_weights_gmc3, + self.lt_weights_gmc4) = get_lt_weights(self.gmpes_list) + + +def get_gmpes(config_file): + """ + Extract strings of the GMPEs from the configuration file. Also get the + labels used in Sammons maps, Clustering (dendrograms) and Euclidean distance + matrix plots. The baseline GMM for computing ratios with is also extracted + if specified within the toml file. + """ + # Get the GMPEs + gmpe_labels = config_file['gmpe_labels']['gmpes_label'] + gmpe_list = [] + config = copy.deepcopy(config_file) + for key in config['models']: + value = get_model(key, config['models']) + gmpe_list.append(value.strip()) + + # Check number of GMPEs matches number of GMPE labels + if len(gmpe_list) != len(gmpe_labels): + raise ValueError("Number of labels must match number of GMPEs.") + + # Get the baseline GMPE used to compute ratios of GMPEs with if required + if 'ratios_baseline_gmm' in config_file.keys(): + if len(config_file['ratios_baseline_gmm']) > 1: + raise ValueError('Only one baseline GMPE should be specified.') + for key in config_file['ratios_baseline_gmm']: + baseline_gmm = get_model(key, config['ratios_baseline_gmm']) + else: + baseline_gmm = None + + return gmpe_labels, gmpe_list, baseline_gmm + + +def get_model(key, models): + """ + Get the model from the toml in the string format required to create an + OpenQuake gsim object from within mgmpe_check (in utils_gmpes.py) + """ + # If the key contains a number we take the second part + if re.search("^\\d+\\-", key): + tmp = re.sub("^\\d+\\-", "", key) + value = f"[{tmp}] " + else: + value = f"[{key}] " + if len(models[key]): + models[key].pop('style', None) + value += '\n' + str(toml.dumps(models[key])) + return value.strip() + + +def get_lt_weights(gmpe_list): + """ + Manage the logic tree weight assigned for each GMPE in the toml (if any) + """ + # If weight is assigned to a GMPE get it + check sum of weights for + # GMPEs with weights allocated is about 1 + weights = [{}, {}, {}, {}] + for gmpe in gmpe_list: + if 'lt_weight' in gmpe: + split_gmpe_str = str(gmpe).splitlines() + for idx, component in enumerate(split_gmpe_str): + if 'lt_weight_gmc1' in component: + weights[0][gmpe] = float(split_gmpe_str[ + idx].split('=')[1]) + if 'lt_weight_gmc2' in component: + weights[1][gmpe] = float(split_gmpe_str[ + idx].split('=')[1]) + if 'lt_weight_gmc3' in component: + weights[2][gmpe] = float(split_gmpe_str[ + idx].split('=')[1]) + if 'lt_weight_gmc4' in component: + weights[3][gmpe] = float(split_gmpe_str[ + idx].split('=')[1]) - # Get models - gmpes_list_initial = [] - config = copy.deepcopy(config_file) - for key in config['models']: - - # If the key contains a number we take the second part - if re.search("^\\d+\\-", key): - tmp = re.sub("^\\d+\\-", "", key) - value = f"[{tmp}] " - else: - value = f"[{key}] " - if len(config['models'][key]): - config['models'][key].pop('style', None) - value += '\n' + str(toml.dumps(config['models'][key])) - value = value.strip() - gmpes_list_initial.append(value.strip()) - - self.gmpes_list = [] - for idx, gmpe in enumerate(gmpes_list_initial): - self.gmpes_list.append(gmpes_list_initial[idx]) - - # Check number of GMPEs matches number of GMPE labels - if len(self.gmpes_list) != len(self.gmpe_labels): - raise ValueError("Number of labels must match number of GMPEs.") - - # If weight is assigned to a GMPE get it + check sum of weights for - # GMPEs with weights allocated = 1 - weights = [{}, {}, {}, {}] - for gmpe in self.gmpes_list: - if 'lt_weight' in gmpe: - split_gmpe_str = str(gmpe).splitlines() - for idx, component in enumerate(split_gmpe_str): - if 'lt_weight_gmc1' in component: - weights[0][gmpe] = float(split_gmpe_str[ - idx].split('=')[1]) - if 'lt_weight_gmc2' in component: - weights[1][gmpe] = float(split_gmpe_str[ - idx].split('=')[1]) - if 'lt_weight_gmc3' in component: - weights[2][gmpe] = float(split_gmpe_str[ - idx].split('=')[1]) - if 'lt_weight_gmc4' in component: - weights[3][gmpe] = float(split_gmpe_str[ - idx].split('=')[1]) - - # Check weights for each logic tree (if present) equal 1.0 - msg = "Sum of GMC logic tree weights must be 1.0" - if weights[0] != {}: - check_weights_gmc1 = np.array(pd.Series(weights[0])) - lt_total_wt_gmc1 = np.sum(check_weights_gmc1, axis=0) - assert abs(lt_total_wt_gmc1-1.0) < 1e-10, msg - self.lt_weights_gmc1 = weights[0] - else: - self.lt_weights_gmc1 = None + # Check weights for each logic tree (if present) equal 1.0 + msg = "Sum of GMC logic tree weights must be 1.0" + if weights[0] != {}: + check_weights_gmc1 = np.array(pd.Series(weights[0])) + lt_total_wt_gmc1 = np.sum(check_weights_gmc1, axis=0) + assert abs(lt_total_wt_gmc1-1.0) < 1e-10, msg + lt_weights_gmc1 = weights[0] + else: + lt_weights_gmc1 = None + + if weights[1] != {}: + check_weights_gmc2 = np.array(pd.Series(weights[1])) + lt_total_wt_gmc2 = np.sum(check_weights_gmc2, axis=0) + assert abs(lt_total_wt_gmc2-1.0) < 1e-10, msg + lt_weights_gmc2 = weights[1] + else: + lt_weights_gmc2 = None + + if weights[2] != {}: + check_weights_gmc3 = np.array(pd.Series(weights[2])) + lt_total_wt_gmc3 = np.sum(check_weights_gmc3, axis=0) + assert abs(lt_total_wt_gmc3-1.0) < 1e-10, msg + lt_weights_gmc3 = weights[2] + else: + lt_weights_gmc3 = None - if weights[1] != {}: - check_weights_gmc2 = np.array(pd.Series(weights[1])) - lt_total_wt_gmc2 = np.sum(check_weights_gmc2, axis=0) - assert abs(lt_total_wt_gmc2-1.0) < 1e-10, msg - self.lt_weights_gmc2 = weights[1] - else: - self.lt_weights_gmc2 = None - - if weights[2] != {}: - check_weights_gmc3 = np.array(pd.Series(weights[2])) - lt_total_wt_gmc3 = np.sum(check_weights_gmc3, axis=0) - assert abs(lt_total_wt_gmc3-1.0) < 1e-10, msg - self.lt_weights_gmc3 = weights[2] - else: - self.lt_weights_gmc3 = None + if weights[3] != {}: + check_weights_gmc4 = np.array(pd.Series(weights[3])) + lt_total_wt_gmc4 = np.sum(check_weights_gmc4, axis=0) + assert abs(lt_total_wt_gmc4-1.0) < 1e-10, msg + lt_weights_gmc4 = weights[3] + else: + lt_weights_gmc4 = None + + return lt_weights_gmc1, lt_weights_gmc2, lt_weights_gmc3, lt_weights_gmc4 + + +def assign_depths_per_mag_bin(config_file, mag_array): + """ + For each magnitude considered within the Sammons Maps, Euclidean distance + and clustering plots assign a depth + """ + # Create dataframe of depth to assign per mag bin + non_trellis_or_spectra_depths = pd.DataFrame(config_file[ + 'mag_values_non_trellis_or_spectra_functions'][ + 'non_trellis_or_spectra_depths'], columns=['mag','depth']) - if weights[3] != {}: - check_weights_gmc4 = np.array(pd.Series(weights[3])) - lt_total_wt_gmc4 = np.sum(check_weights_gmc4, axis=0) - assert abs(lt_total_wt_gmc4-1.0) < 1e-10, msg - self.lt_weights_gmc4 = weights[3] - else: - self.lt_weights_gmc4 = None + # Round each mag in mag_array to closest integer + mag_to_nearest_int = pd.Series(dtype='float') + for mag in mag_array: + mag_to_nearest_int[mag] = np.round(mag+0.001) + + # Assign depth to each mag in mag_array using rounded mags + depth_array_initial = [] + for idx_mag, rounded_mag in enumerate(mag_to_nearest_int): + for idx, val in enumerate(non_trellis_or_spectra_depths['mag']): + if rounded_mag == non_trellis_or_spectra_depths['mag'][idx]: + depth_to_store = non_trellis_or_spectra_depths['depth'][idx] + depth_array_initial.append(depth_to_store) + + return pd.Series(depth_array_initial) def plot_trellis(filename, output_directory): @@ -196,7 +259,6 @@ def plot_trellis(filename, output_directory): toml file providing configuration for use within comparative plotting methods. """ - # Generate config object config = Configurations(filename) store_gmm_curves = plot_trellis_util(config, output_directory) @@ -215,7 +277,6 @@ def plot_spectra(filename, output_directory, obs_spectra=None): csv of an observed spectra to plot and associated event information. An example file can be found in openquake.smt.tests.file_samples. """ - # Generate config object config = Configurations(filename) store_gmc_lts = plot_spectra_util(config, output_directory, obs_spectra) @@ -223,6 +284,23 @@ def plot_spectra(filename, output_directory, obs_spectra=None): return store_gmc_lts +def plot_ratios(filename, output_directory): + """ + Plot ratio (GMPE median attenuation/baseline GMPE median attenuation) for + given run configuration + :param filename: + toml file providing configuration for use within comparative + plotting methods. + """ + config = Configurations(filename) + + if config.baseline_gmm is None: + raise ValueError( + 'User must specify a baseline GMPE to generate ratio plots') + + plot_ratios_util(config, output_directory) + + def plot_cluster(filename, output_directory): """ Plot hierarchical clusters of (1) median, (2) 84th percentile and (3) 16th @@ -231,7 +309,6 @@ def plot_cluster(filename, output_directory): toml file providing configuration for use within comparative plotting methods. """ - # Generate config object with each set of run parameters config = Configurations(filename) if len(config.gmpes_list) < 2: @@ -266,7 +343,6 @@ def plot_sammons(filename, output_directory): toml file providing configuration for use within comparative plotting methods. """ - # Generate config object with each set of run parameters config = Configurations(filename) if len(config.gmpes_list) < 2: @@ -302,7 +378,6 @@ def plot_euclidean(filename, output_directory): toml file providing configuration for use within comparative plotting methods. """ - # Generate config object config = Configurations(filename) if len(config.gmpes_list) < 2: @@ -325,32 +400,4 @@ def plot_euclidean(filename, output_directory): plot_euclidean_util(config.imt_list, config.gmpe_labels, mtxs_16th_perc, os.path.join(output_directory,'16th_perc_Euclidean.png'), - mtxs_type='16th_perc') - - -def assign_depths_per_mag_bin(config_file, mag_array): - """ - For each magnitude considered within the Sammons Maps, Euclidean distance - and clustering plots assign a depth - """ - # Create dataframe of depth to assign per mag bin - non_trellis_or_spectra_depths = pd.DataFrame(config_file[ - 'mag_values_non_trellis_or_spectra_functions'][ - 'non_trellis_or_spectra_depths'], columns=['mag','depth']) - - # Round each mag in mag_array to closest integer - mag_to_nearest_int = pd.Series(dtype='float') - for mag in mag_array: - mag_to_nearest_int[mag] = np.round(mag+0.001) - - # Assign depth to each mag in mag_array using rounded mags - depth_array_initial = [] - for idx_mag, rounded_mag in enumerate(mag_to_nearest_int): - for idx, val in enumerate(non_trellis_or_spectra_depths['mag']): - if rounded_mag == non_trellis_or_spectra_depths['mag'][idx]: - depth_to_store = non_trellis_or_spectra_depths['depth'][idx] - depth_array_initial.append(depth_to_store) - - depths = pd.Series(depth_array_initial) - - return depths \ No newline at end of file + mtxs_type='16th_perc') \ No newline at end of file diff --git a/openquake/smt/comparison/utils_compare_gmpes.py b/openquake/smt/comparison/utils_compare_gmpes.py index 3d87786cc..cad677955 100644 --- a/openquake/smt/comparison/utils_compare_gmpes.py +++ b/openquake/smt/comparison/utils_compare_gmpes.py @@ -45,13 +45,19 @@ def plot_trellis_util(config, output_directory): # Median, plus sigma, minus sigma per gmc for up to 4 gmc logic trees gmc_p= [[{}, {}, {}], [{}, {}, {}], [{}, {}, {}], [{}, {}, {}]] - # Get model weights + basin params + colors + config key + # Get basin params Z1, Z25 = get_z1_z25(config.Z1, config.Z25, config.Vs30, config.region) - colors = get_cols(config.custom_color_flag, config.custom_color_list) + + # Get lt weights lt_weights = [config.lt_weights_gmc1, config.lt_weights_gmc2, config.lt_weights_gmc3, config.lt_weights_gmc4] - cfg_key = 'vs30 = %s m/s, GMM sigma epsilon = %s' % ( - config.Vs30, config.Nstd) + + # Get config key + cfg_key = 'vs30 = %s m/s, GMM sigma epsilon = %s' % (config.Vs30, + 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 @@ -64,6 +70,19 @@ def plot_trellis_util(config, output_directory): for l, m in enumerate(mag_list): fig.add_subplot( len(config.imt_list), len(mag_list), l+1+n*len(mag_list)) + + # get ztor + if config.ztor is not None: + ztor_m = config.ztor[l] + else: + ztor_m = None + + # Get gmpe params + strike_g, dip_g, depth_g, aratio_g = _param_gmpes( + config.strike, config.dip, dep_list[l], config.aratio, + config.rake, config.trt) + + # Per GMPE get attenuation curves lt_vals_gmc = [{}, {}, {}, {}] store_per_gmpe = {} for g, gmpe in enumerate(config.gmpes_list): @@ -75,17 +94,6 @@ def plot_trellis_util(config, output_directory): # Perform mgmpe check gmm = mgmpe_check(gmpe) - # ZTOR value - if config.ztor is not None: - ztor_m = config.ztor[l] - else: - ztor_m = None - - # Get gmpe params - strike_g, dip_g, depth_g, aratio_g = _param_gmpes( - config.strike, config.dip, dep_list[l], config.aratio, - config.rake, config.trt) - # Get attenuation curves mean, std, r_vals, tau, phi = att_curves(gmm, dep_list[l], m, aratio_g, strike_g, @@ -105,27 +113,14 @@ def plot_trellis_util(config, output_directory): minus_sigma = np.exp(mean-config.Nstd*std[0]) # Plot predictions and get lt weighted predictions - lt_vals_gmc = trellis_data(gmpe, r_vals, mean,plus_sigma, - minus_sigma, col, i, m, config.Nstd, + lt_vals_gmc = trellis_data(gmpe, r_vals, mean, plus_sigma, + minus_sigma, col, config.Nstd, lt_vals_gmc, lt_weights) + # Get unit of imt for the store + unit = get_imtl_unit_for_trellis_store(i) + # Store per gmpe - 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']: - pyplot.ylabel(str(i) + ' (Hz)') # Fourier/Eff. Amp. Spectrum - else: - unit = 'g' # PGA, SA, AvgSA store_per_gmpe[gmpe]['median (%s)' % unit] = np.exp(mean) store_per_gmpe[gmpe]['sigma (ln)'] = std if config.Nstd != 0: @@ -185,16 +180,17 @@ def plot_spectra_util(config, output_directory, obs_spectra): obs_spectra = pd.read_csv(obs_spectra) max_period, eq_id, st_id = load_obs_spectra(obs_spectra) else: + max_period = config.max_period obs_spectra, eq_id, st_id = None, None, None # Get gmc lt weights, imts, periods and basin params gmc_weights = [config.lt_weights_gmc1, config.lt_weights_gmc2, config.lt_weights_gmc3, config.lt_weights_gmc4] - imt_list, periods = _get_imts(config.max_period) + imt_list, periods = _get_imts(max_period) Z1, Z25 = get_z1_z25(config.Z1, config.Z25, config.Vs30, config.region) # Get colours and make the figure - colors = get_cols(config.custom_color_flag, config.custom_color_list) + colors = get_colors(config.custom_color_flag, config.custom_color_list) figure = pyplot.figure(figsize=(len(mag_list)*5, len(config.dist_list)*4)) # Set dicts to store values @@ -210,13 +206,18 @@ def plot_spectra_util(config, output_directory, obs_spectra): ax1 = figure.add_subplot( len(config.dist_list), len(mag_list), l+1+n*len(mag_list)) + + strike_g, dip_g, depth_g, aratio_g = _param_gmpes(config.strike, + config.dip, + dep_list[l], + config.aratio, + config.rake, + config.trt) + for g, gmpe in enumerate(config.gmpes_list): rs_50p, sig, rs_ps, rs_ms = [], [], [], [] col = colors[g] gmm = mgmpe_check(gmpe) - strike_g, dip_g, depth_g, aratio_g = _param_gmpes( - config.strike, config.dip, dep_list[l], - config.aratio, config.rake, config.trt) for k, imt in enumerate(imt_list): if config.ztor is not None: @@ -308,6 +309,91 @@ def plot_spectra_util(config, output_directory, obs_spectra): return store_gmc_lts # Returned for unit tests of gmc lt values +def plot_ratios_util(config, output_directory): + """ + Generate ratio (GMPE median attenuation/baseline GMPE median attenuation) + plots for given run configuration + """ + # Get mag and dep lists + mag_list = config.trellis_and_rs_mag_list + dep_list = config.trellis_and_rs_depth_list + + # Get basin params + Z1, Z25 = get_z1_z25(config.Z1, config.Z25, config.Vs30, config.region) + + # Get config key + cfg_key = 'vs30 = %s m/s, GMM sigma epsilon = %s' % (config.Vs30, + config.Nstd) + + # Get colours + colors = get_colors(config.custom_color_flag, config.custom_color_list) + + # Compute ratio curves + fig = pyplot.figure(figsize=(len(mag_list)*5, len(config.imt_list)*4)) + ratio_store = [] + for n, i in enumerate(config.imt_list): + for l, m in enumerate(mag_list): + fig.add_subplot( + len(config.imt_list), len(mag_list), l+1+n*len(mag_list)) + + # ztor value + if config.ztor is not None: + ztor_m = config.ztor[l] + else: + ztor_m = None + + # Get gmpe params + strike_g, dip_g, depth_g, aratio_g = _param_gmpes( + config.strike, config.dip, dep_list[l], config.aratio, + config.rake, config.trt) + + # Load the baseline GMM and compute baseline + baseline = mgmpe_check(config.baseline_gmm) + + # Get baseline GMM attenuation curves + results = att_curves(baseline, dep_list[l], m, aratio_g, + strike_g, dip_g, config.rake, + config.Vs30, Z1, Z25, config.maxR, 1, i, ztor_m, + config.eshm20_region, config.dist_type, + config.trt, config.up_or_down_dip) + b_mean = results[0][0][0] + + # Now compute ratios for each GMM + for g, gmpe in enumerate(config.gmpes_list): + + # Perform mgmpe check + col = colors[g] + gmm = mgmpe_check(gmpe) + + # Get attenuation curves for the GMM + results = att_curves(gmm, dep_list[l], m, aratio_g, strike_g, + dip_g, config.rake, config.Vs30, Z1, Z25, + config.maxR, 1, i, ztor_m, + config.eshm20_region, config.dist_type, + config.trt, config.up_or_down_dip) + + # 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=gmpe) + + # Update plots + update_ratio_plots(config.dist_type, m, i, n, l, config.imt_list, + r_vals, config.minR, config.maxR) + + # Finalise plots + pyplot.legend(loc="center left", bbox_to_anchor=(1.1, 1.05), fontsize='16') + pyplot.savefig(os.path.join(output_directory, 'RatioPlots.png'), + bbox_inches='tight', dpi=200, pad_inches=0.2) + + def compute_matrix_gmpes(config, mtxs_type): """ Compute matrix of median ground-motion predictions for each gmpe for the @@ -335,11 +421,14 @@ def compute_matrix_gmpes(config, mtxs_type): gmm = mgmpe_check(gmpe) - strike_g, dip_g, depth_g, aratio_g = _param_gmpes( - config.strike, config.dip, dep_list[l], config.aratio, - config.rake, config.trt) + strike_g, dip_g, depth_g, aratio_g = _param_gmpes(config.strike, + config.dip, + dep_list[l], + config.aratio, + config.rake, + config.trt) - # ZTOR + # ztor if config.ztor is not None: ztor_m = config.ztor[l] else: @@ -468,7 +557,7 @@ def plot_sammons_util(imt_list, gmpe_list, mtxs, namefig, custom_color_flag, mtxs, gmpe_list = matrix_mean(mtxs, imt_list, gmpe_list) # Setup - colors = get_cols(custom_color_flag, custom_color_list) + colors = get_colors(custom_color_flag, custom_color_list) texts = [] if len(imt_list) < 3: nrows = 1 @@ -593,14 +682,12 @@ def plot_cluster_util(imt_list, gmpe_list, mtxs, namefig, mtxs_type): ### Utils for plots -def get_cols(custom_color_flag, custom_color_list): +def get_colors(custom_color_flag, custom_color_list): """ Get list of colors for plots """ - colors = ['g', 'b', 'y', 'lime', 'dodgerblue', 'gold', '0.8', '0.5', 'r', - 'm', 'mediumseagreen', 'tab:pink', 'tab:orange', 'tab:purple', - 'tab:brown', 'tab:pink', 'tab:red', 'tab:blue', 'tab:cyan', - 'tab:olive', 'aquamarine'] + colors = ['r', 'g', 'b', 'y', 'lime', 'dodgerblue', 'gold', '0.8', 'm', 'k', + 'mediumseagreen', 'tab:orange', 'tab:purple', 'tab:brown', '0.5'] if custom_color_flag == 'True': colors = custom_color_list @@ -622,7 +709,7 @@ def get_z1_z25(Z1, Z25, Vs30, region): ### Trellis utils -def trellis_data(gmpe, r_vals, mean,plus_sigma, minus_sigma, col, i, m, Nstd, +def trellis_data(gmpe, r_vals, mean,plus_sigma, minus_sigma, col, Nstd, lt_vals_gmc, lt_weights): """ Plot predictions of a single GMPE (if required) and compute weighted @@ -644,7 +731,7 @@ def trellis_data(gmpe, r_vals, mean,plus_sigma, minus_sigma, col, i, m, Nstd, pass elif gmpe in lt_weights[idx_gmc]: if lt_weights[idx_gmc][gmpe] is not None: - if Nstd != 0: + if Nstd > 0: lt_vals_gmc[idx_gmc][gmpe] = { 'median': np.exp(mean)*lt_weights[ idx_gmc][gmpe], @@ -680,7 +767,7 @@ def trel_logic_trees(idx_gmc, gmc, lt_vals_gmc, gmc_p, store_gmm_curves, 'gmc logic tree curves per imt-mag'][lt_key][ 'median (%s)' % unit] = median - if Nstd != 0: + if Nstd > 0: store_gmm_curves[ cfg_key]['gmc logic tree curves per imt-mag'][ lt_key]['median plus sigma (%s)' % unit @@ -718,7 +805,7 @@ def lt_trel(r_vals, Nstd, i, m, dep, dip, rake, idx_gmc, lt_vals_gmc, linestyle='--', label=label, zorder=100) - if Nstd != 0: + if Nstd > 0: lt_plus_sigma_gmc = np.sum(lt_df_gmc[:].loc['plus_sigma']) lt_minus_sigma_gmc = np.sum(lt_df_gmc[:].loc['minus_sigma']) @@ -775,7 +862,7 @@ def update_trellis_plots(m, i, n, l, minR, maxR, r_vals, imt_list, dist_type): elif str(i) in ['PGV']: pyplot.ylim(1, 650) # cm/s else: - print('User may want to manually specify y-axis limits') + print('User may want to manually specify y-axis limits for %s' % i) # xlims pyplot.loglog() @@ -783,6 +870,31 @@ def update_trellis_plots(m, i, n, l, minR, maxR, r_vals, imt_list, dist_type): pyplot.xlim(np.max([min_r_val, minR]), maxR) +def get_imtl_unit_for_trellis_store(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']: + pyplot.ylabel(str(i) + ' (Hz)') # Fourier/Eff. Amp. Spectrum + else: + unit = 'g' # PGA, SA, AvgSA + + return unit + + ### Spectra utils def _get_period_values_for_spectra_plots(max_period): """ @@ -865,7 +977,7 @@ def spectra_data(gmpe, Nstd, gmc_weights, rs_50p, rs_plus_sigma, rs_minus_sigma, for idx, rs in enumerate(rs_50p): rs_50p_w[idx] = rs_50p[idx]*gmc_weights[ idx_gmc][gmpe] - if Nstd != 0: + if Nstd > 0: rs_plus_sigma_w[idx] = rs_plus_sigma[ idx]*gmc_weights[idx_gmc][gmpe] rs_minus_sigma_w[idx] = rs_minus_sigma[ @@ -874,8 +986,8 @@ def spectra_data(gmpe, Nstd, gmc_weights, rs_50p, rs_plus_sigma, rs_minus_sigma, # Store the weighted median for the GMPE lt_vals[idx_gmc][gmpe] = {'median': rs_50p_w} - # And if Nstd != 0 store these weighted branches too - if Nstd != 0: + # And if Nstd > 0 store these weighted branches too + if Nstd > 0: lt_vals_plus[idx_gmc][gmpe,'p_sig'] = { 'plus_sigma': rs_plus_sigma_w} lt_vals_minus[idx_gmc][gmpe,'m_sig'] = { @@ -903,7 +1015,7 @@ def lt_spectra(ax1, gmpe, gmpe_list, Nstd, period, idx_gmc, # Plot lt_per_imt_gmc, lt_plus_sig_per_imt, lt_minus_sig_per_imt = {}, {}, {} lt_df_gmc = pd.DataFrame(lt_vals_gmc, index=['median']) - if Nstd != 0: + if Nstd > 0: lt_plus_sig = pd.DataFrame(ltv_plus_sig, index=['plus_sigma']) lt_minus_sig = pd.DataFrame(ltv_minus_sig, index=['minus_sigma']) wt_per_gmpe_gmc, wt_plus_sig, wt_minus_sig = {}, {}, {} @@ -911,7 +1023,7 @@ def lt_spectra(ax1, gmpe, gmpe_list, Nstd, period, idx_gmc, if check in str(gmpe): wt_per_gmpe_gmc[gmpe] = np.array(pd.Series( lt_df_gmc[gmpe].loc['median'])) - if Nstd != 0: + if Nstd > 0: wt_plus_sig[gmpe] = np.array( pd.Series(lt_plus_sig[gmpe,'p_sig'].loc['plus_sigma'])) wt_minus_sig[gmpe] = np.array( @@ -922,7 +1034,7 @@ def lt_spectra(ax1, gmpe, gmpe_list, Nstd, period, idx_gmc, lt_minus_sig = pd.DataFrame(wt_minus_sig, index=period) for idx, imt in enumerate(period): lt_per_imt_gmc[imt] = np.sum(lt_df_gmc.loc[imt]) - if Nstd != 0: + if Nstd > 0: lt_plus_sig_per_imt[imt] = np.sum(lt_plus_sig.loc[imt]) lt_minus_sig_per_imt[imt] = np.sum(lt_minus_sig.loc[imt]) @@ -931,7 +1043,7 @@ def lt_spectra(ax1, gmpe, gmpe_list, Nstd, period, idx_gmc, color=col, linestyle='--', label = label, zorder=100) # Plot mean plus sigma and mean minus sigma if required - if Nstd != 0: + if Nstd > 0: ax1.plot(period, np.array(pd.Series(lt_plus_sig_per_imt)), linewidth=0.75, color=col, linestyle='-.', zorder=100) ax1.plot(period, np.array(pd.Series(lt_minus_sig_per_imt)), @@ -1007,6 +1119,27 @@ def save_spectra_plot(f1, obs_spectra, output_dir, eq_id, st_id): ### Utils for other plots +def update_ratio_plots(dist_type, m, i, n, l, imt_list, r_vals, minR, maxR): + """ + Add titles and axis labels to ratio plots + """ + if dist_type == 'repi': + label = 'Repi (km)' + if dist_type == 'rrup': + label = 'Rrup (km)' + if dist_type == 'rjb': + label = 'Rjb (km)' + if dist_type == 'rhypo': + label = 'Rhypo (km)' + if n == 0: # Top row only + pyplot.title('Mw = ' + str(m), fontsize='16') + if n == len(imt_list)-1: # Bottom row only + pyplot.xlabel(label, fontsize='12') + if l == 0: # Left row only + pyplot.ylabel('GMM/baseline for %s' %str(i), fontsize='14') + min_r_val = min(r_vals[r_vals>=1]) + pyplot.xlim(np.max([min_r_val, minR]), maxR) + def matrix_mean(mtxs, imt_list, gmpe_list): """ For a matrix of predicted ground-motions computed the arithmetic mean per diff --git a/openquake/smt/demos/demo_comparison_analysis_inputs.toml b/openquake/smt/demos/demo_comparison_analysis_inputs.toml index 50991c0c8..bb3cdfe7c 100644 --- a/openquake/smt/demos/demo_comparison_analysis_inputs.toml +++ b/openquake/smt/demos/demo_comparison_analysis_inputs.toml @@ -83,6 +83,9 @@ lt_weight_gmc2_plot_lt_only = 0.50 [models.AkkarEtAlRjb2014] lt_weight_gmc2_plot_lt_only = 0.50 +# Also specify a GMM to compute ratios of the attenuation against (GMM/baseline) +[ratios_baseline_gmm.BooreEtAl2020] + [custom_colors] custom_colors_flag = 'False' # Set to "True" for custom colours in plots) custom_colors_list = ['lime', 'dodgerblue', 'gold', '0.8'] \ No newline at end of file diff --git a/openquake/smt/demos/demo_comparison_script.py b/openquake/smt/demos/demo_comparison_script.py index 85285c2de..7db3e2b83 100644 --- a/openquake/smt/demos/demo_comparison_script.py +++ b/openquake/smt/demos/demo_comparison_script.py @@ -23,12 +23,13 @@ name_analysis = config_file['general']['name_analysis'] output_directory = os.path.join(os.path.abspath(''),name_analysis) - # set the output + # Set the output if not os.path.exists(output_directory): os.makedirs(output_directory) - #Generate plots from config object + # Generate plots from config object attenuation_curve_data[file] = comp.plot_trellis(filename,output_directory) comp.plot_spectra(filename,output_directory) + comp.plot_ratios(filename, output_directory) comp.plot_cluster(filename,output_directory) comp.plot_sammons(filename,output_directory) comp.plot_euclidean(filename,output_directory) diff --git a/openquake/smt/residuals/residual_plotter.py b/openquake/smt/residuals/residual_plotter.py index 93f786786..24480c983 100755 --- a/openquake/smt/residuals/residual_plotter.py +++ b/openquake/smt/residuals/residual_plotter.py @@ -37,9 +37,9 @@ residuals_density_distribution, likelihood, residuals_with_magnitude, residuals_with_vs30, residuals_with_distance, residuals_with_depth) -colors = ['r', 'g', 'b', 'y','lime', 'dodgerblue', 'k', 'gold', '0.8', - 'mediumseagreen', '0.5', 'tab:orange', 'tab:purple','tab:brown', - 'tab:pink'] + +colors = ['r', 'g', 'b', 'y', 'lime', 'dodgerblue', 'gold', '0.8', 'm', 'k', + 'mediumseagreen', 'tab:orange', 'tab:purple', 'tab:brown', '0.5'] class BaseResidualPlot(object): diff --git a/openquake/smt/tests/comparison/comparison_test.py b/openquake/smt/tests/comparison/comparison_test.py index 1a269d697..97c960764 100644 --- a/openquake/smt/tests/comparison/comparison_test.py +++ b/openquake/smt/tests/comparison/comparison_test.py @@ -26,7 +26,7 @@ from openquake.smt.comparison import compare_gmpes as comp from openquake.smt.comparison.utils_compare_gmpes import ( compute_matrix_gmpes, plot_trellis_util, plot_spectra_util, - plot_cluster_util, plot_sammons_util, plot_euclidean_util) + plot_ratios_util, plot_cluster_util, plot_sammons_util, plot_euclidean_util) # Base path @@ -47,6 +47,7 @@ '[CampbellBozorgnia2014] \nlt_weight_gmc1 = 0.5', '[BooreEtAl2014] \nlt_weight_gmc2_plot_lt_only = 0.5', '[KothaEtAl2020] \nlt_weight_gmc2_plot_lt_only = 0.5'] +TARGET_BASELINE_GMPE = '[BooreEtAl2014]' TARGET_TRT = 'ASCR' TARGET_ZTOR = None @@ -76,13 +77,13 @@ def test_configuration_object_check(self): the Configuration object, which stores the inputted parameters for each run. """ - # Check each parameter matches target + # Load config config = comp.Configurations(self.input_file) # Check for target TRT self.assertEqual(config.trt, TARGET_TRT) - # Check for target ZTOR + # Check for target ztor self.assertEqual(config.ztor, TARGET_ZTOR) # Check for target Vs30 @@ -122,12 +123,15 @@ def test_configuration_object_check(self): for imt in range(0, len(config.imt_list)): self.assertEqual(str(config.imt_list[imt]), TARGET_IMTS[imt]) + # Check baseline GMM used to compute ratios + self.assertEqual(config.baseline_gmm, TARGET_BASELINE_GMPE) + def test_mtxs_median_calculation(self): """ Check for matches bewteen the matrix of medians computed using compute_matrix_gmpes and those expected given the input parameters """ - # Check each parameter matches target + # Load config config = comp.Configurations(self.input_file) # Get medians @@ -147,7 +151,7 @@ def test_sammons_and_euclidean_distance_matrix_functions(self): """ TARGET_GMPES.append('mean') # Add mean here to gmpe_list - # Check each parameter matches target + # Load config config = comp.Configurations(self.input_file) # Get medians @@ -188,7 +192,7 @@ def test_clustering_median(self): Check clustering functions for median predicted ground-motion of considered GMPEs in the configuration """ - # Check each parameter matches target + # Load config config = comp.Configurations(self.input_file) # Get medians @@ -213,7 +217,7 @@ def test_clustering_84th_perc(self): Check clustering of 84th percentile of predicted ground-motion of considered GMPEs in the configuration """ - # Check each parameter matches target + # Load config config = comp.Configurations(self.input_file) # Get medians @@ -239,7 +243,7 @@ def test_trellis_and_spectra_functions(self): executed. Also checks correct values are returned for the gmm attenuation curves and spectra. """ - # Check each parameter matches target + # Load config config = comp.Configurations(self.input_file) # Trellis plots @@ -286,9 +290,20 @@ def test_plot_observed_spectra(self): # Check target file created and outputted in expected location self.assertTrue(target_file_spectra) + def test_plot_ratios(self): + """ + Test execution of plotting ratios (median GMM attenuation/median + baseline GMM attenuation). Correctness of values is not examined. + """ + # Load config + config = comp.Configurations(self.input_file) + + # Plot the ratios + plot_ratios_util(config, self.output_directory) + @classmethod def tearDownClass(self): """ Remove the test outputs """ - shutil.rmtree(self.output_directory) + shutil.rmtree(self.output_directory) \ 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 2aaf0a268..0c9ff3e6f 100644 --- a/openquake/smt/tests/comparison/data/compare_gmpe_inputs.toml +++ b/openquake/smt/tests/comparison/data/compare_gmpe_inputs.toml @@ -52,6 +52,9 @@ lt_weight_gmc2_plot_lt_only = 0.5 [models.KothaEtAl2020] lt_weight_gmc2_plot_lt_only = 0.5 +# Baseline GMM to compute ratios of the attenuation against (GMM/baseline) +[ratios_baseline_gmm.BooreEtAl2014] + [custom_colors] custom_colors_flag = 'False' #By default set to "False" (set to "True" for custom colours) diff --git a/openquake/smt/tests/comparison/data/mgmpe_inputs.toml b/openquake/smt/tests/comparison/data/mgmpe_inputs.toml index 59cedb6b4..207790669 100644 --- a/openquake/smt/tests/comparison/data/mgmpe_inputs.toml +++ b/openquake/smt/tests/comparison/data/mgmpe_inputs.toml @@ -97,6 +97,11 @@ gmpe_table = 'NGAEast_FRANKEL_J15.hdf5' # use a gmpe table d_sigma = 100 # gmpe specific param kappa0 = 0.04 +[ratios_baseline_gmm.ModifiableGMPE] # Check mgmpe works for loading baseline model used in ratio plots +gmpe = 'YenierAtkinson2015BSSA' +sigma_model = 'al_atik_2015_sigma' # use Al Atik (2015) sigma model +with_betw_ratio = 1.5 + [custom_colors] custom_colors_flag = 'False' # By default set to "False" (set to "True" for custom colours) custom_colors_list = ['lime','dodgerblue','gold','0.8'] \ No newline at end of file