added mgmpe features to smt

docsrc/contents/smt.rst

@@ -96,32 +96,59 @@ 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:
 
     .. code-block:: ini
 
         [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
 

openquake/smt/comparison/compare_gmpes.py

@@ -65,14 +65,16 @@ 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
         self.strike = config_file['source_properties']['strike']
         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'][

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