Merge pull request #363 from GEMScienceTools/tidying_smt

Adding mgmpe unit test + tidying smt

docsrc/contents/smt.rst

@@ -128,7 +128,7 @@ We can specify the inputs to perform a residual analysis with as follows:
         [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]
+        [models.1-KothaEtAl2020]
         sigma_scaling_scalar = 1.05 # scale sigma by factor of 1.05 over all imts
         
         [models.2-KothaEtAl2020]
@@ -150,6 +150,11 @@ We can specify the inputs to perform a residual analysis with as follows:
         d_sigma = 100 # gmpe specific param
         kappa0 = 0.04
         
+        [models.KothaEtAl2020ESHM20] # ESHM20 model
+        sigma_mu_epsilon = 2.85697 
+        c3_epsilon = 1.72    
+        region = 4 # Note that within residuals specify region here, whereas in comparison module toml (below) we specify using the eshm20_region param
+        
         [imts]
         imt_list = ['PGA', 'SA(0.2)', 'SA(0.5)', 'SA(1.0']    
         
@@ -477,9 +482,9 @@ Comparing GMPEs
         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
+        dip =  60
+        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 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/utils_compare_gmpes.py

@@ -78,7 +78,7 @@ def plot_trellis_util(
                     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]
                 std = std[0][0]
@@ -686,6 +686,7 @@ def update_trellis_plots(m, i, n, l, r_vals, imt_list, dist_type):
     """
     Add titles and axis labels to trellis plots
     """
+    # Labels
     if dist_type == 'repi':
         label = 'Repi (km)'
     if dist_type == 'rrup':
@@ -694,15 +695,22 @@ def update_trellis_plots(m, i, n, l, r_vals, imt_list, dist_type):
         label = 'Rjb (km)'
     if dist_type == 'rhypo':
         label = 'Rhypo (km)'
-    if n == 0: #top row only
-        pyplot.title('Mw = ' + str(m), fontsize='16') # Mod if required
-    if n == len(imt_list)-1: #bottom row only
+    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='16')
-    if l == 0: #left row only
+    if l == 0: # Left row only
         if str(i) != 'PGV':
             pyplot.ylabel(str(i) + ' (g)', fontsize='16')
         else:
             pyplot.ylabel('PGV (cm/s)', fontsize='16')
+
+    # xlims
+    if str(i) != 'PGV':
+        pyplot.ylim(1e-03, 2) # g
+    else:
+        pyplot.ylim(0.1, 650) # cm/s
+
     pyplot.loglog()
     pyplot.xlim(1, np.max(r_vals)) # Mod if required
 

openquake/smt/parsers/esm22_flatfile_parser.py

@@ -663,10 +663,6 @@ def _get_ESM18_headers(ESM22,default_string,r_fm_type,r_datetime):
     "U_lp":ESM22.U_lp,
     "V_lp":ESM22.V_lp,
     "W_lp":ESM22.W_lp,
-
-    # GMIM headers --> equivalency assumed for geometric mean and RotD50 within
-    # parser if RotD50 empty (RotD is not provided in ESM22 format flatfile 
-    # from web service)
 
     "U_pga":ESM22.U_pga,
     "V_pga":ESM22.V_pga,

openquake/smt/parsers/esm23_flatfile_parser.py

@@ -141,21 +141,11 @@ def autobuild(cls, dbid, dbname, output_location,
         """
         Quick and dirty full database builder!
         """
-
         # Import ESM 2023 format strong-motion flatfile
         ESM23 = pd.read_csv(ESM23_flatfile_directory)
-
-        # Store for count of records removed during quality checks
-        ESM23_initial_size = len(ESM23)
-
-         # If vs30 is empty use topographic approximation    
-        for idx in range(0,len(ESM23)):
-            if pd.isnull(ESM23.vs30_m_s[idx]):
-                ESM23['vs30_m_s'].iloc[idx]=ESM23['vs30_m_s_wa'].iloc[idx]
-
+
         # Create default values for headers not considered in ESM23 format
-        default_string = pd.Series(np.full(np.size(ESM23.esm_event_id),
-                                           str("")))
+        default_string = pd.Series(np.full(np.size(ESM23.esm_event_id), ""))
 
         # Assign strike-slip to unknown faulting mechanism
         r_fm_type = ESM23.fm_type_code.fillna('SS') 
@@ -164,7 +154,7 @@ def autobuild(cls, dbid, dbname, output_location,
         r_datetime = ESM23.event_time.str.replace('T',' ')
 
         converted_base_data_path=_get_ESM18_headers(
-            ESM23,default_string,r_fm_type,r_datetime)
+            ESM23, default_string, r_fm_type, r_datetime)
 
         if os.path.exists(output_location):
             raise IOError("Target database directory %s already exists!"
@@ -182,9 +172,6 @@ def autobuild(cls, dbid, dbname, output_location,
         print("Storing metadata to file %s" % metadata_file)
         with open(metadata_file, "wb+") as f:
             pickle.dump(database.database, f)
-
-        print(ESM23_initial_size - len(ESM23),
-              'records removed from imported ESM23 flatfile during quality checks.')
 
         return database
 
@@ -598,12 +585,11 @@ def _retreive_ground_motion_from_row(self, row, header_list):
                     scalars["U"][key] * scalars["V"][key])
         return scalars, spectra
 
-def _get_ESM18_headers(ESM23,default_string,r_fm_type,r_datetime):
+def _get_ESM18_headers(ESM23, default_string, r_fm_type, r_datetime):
 
     """
-    Convert first from ESM23 format flatfile to ESM18 format flatfile
+    Convert from ESM23 format flatfile to ESM18 format flatfile
     """
-
     # Construct dataframe with original ESM format 
     ESM_original_headers = pd.DataFrame(
     {
@@ -655,10 +641,10 @@ def _get_ESM18_headers(ESM23,default_string,r_fm_type,r_datetime):
     "st_elevation":ESM23.st_elevation,
 
     "ec8_code":ESM23.ec8_code,
-    "ec8_code_method":ESM23.preferred_estimation_method_vs30_ec8,
-    "ec8_code_ref":ESM23.ec8_ref_programme,
+    "ec8_code_method":default_string,
+    "ec8_code_ref":default_string,
     "vs30_m_sec":ESM23.vs30_m_s,
-    "vs30_ref":ESM23.vs30_ref_auth,
+    "vs30_ref":default_string,
     "vs30_calc_method":default_string, 
     "vs30_meas_type":ESM23.vs30_meas_type,
     "slope_deg":ESM23.slope_deg,
@@ -954,9 +940,9 @@ def _get_ESM18_headers(ESM23,default_string,r_fm_type,r_datetime):
 
     # Output to folder where converted flatfile read into parser   
     DATA = os.path.abspath('')
-    temp_folder=tempfile.mkdtemp()
-    converted_base_data_path = os.path.join(DATA,temp_folder,
+    tmp = tempfile.mkdtemp()
+    converted_base_data_path = os.path.join(DATA, tmp,
                                             'converted_flatfile.csv')
-    ESM_original_headers.to_csv(converted_base_data_path,sep=';')
+    ESM_original_headers.to_csv(converted_base_data_path, sep=';')
 
     return converted_base_data_path