Comparing GMPEs with Observed Ground Motions

This tutorial illustrates how to use a database of strong ground motions to expected ground motion values from GMPEs with observations.

Loading a Database and Selecting Records

A database (of metadata) is stored as a Python "Pickle" file. To load the database, do the following:



In [ ]:

    
%matplotlib inline
import numpy as np
import cPickle

input_database = "data/demo_records/metadata.pkl"

db1 = cPickle.load(open(input_database, "r"))
print "Number of records: %g" % db1.number_records()

Comparing Records with a Set of GMPEs

In the following example we select a subset of records from a database and compare these observed ground motions against four relevent GMPEs:

1.) Boore & Atkinson (2008) 2.) Akkar & Bommer (2010) 3.) Akkar & Cagnan (2010) 4.) Akkar et al (2014) - Joyner-Boore coefficients

Four intensity measures will also be considered: PGA, Sa (0.2), Sa(1.0), Sa(2.0)



In [ ]:

    
gmpe_list = ["BooreAtkinson2008",
             "AkkarBommer2010",
             "AkkarCagnan2010",
             "AkkarEtAlRjb2014"]

imt_list = ["PGA", "SA(0.2)", "SA(1.0)", "SA(2.0)"]



In [ ]:

    
# Get residuals
from smtk.residuals.gmpe_residuals import Residuals, Likelihood, LLH, EDR
import smtk.residuals.residual_plotter as rspl

# Create an instance of the residual calculator
resid1 = Residuals(gmpe_list, imt_list)
resid1.get_residuals(db1)

Plot the Distribution of Residuals



In [ ]:

    
rspl.ResidualPlot(resid1, "BooreAtkinson2008", "PGA")



In [ ]:

    
rspl.ResidualPlot(resid1, "AkkarEtAlRjb2014", "PGA")

Plots the Ordinary Likelihood Analysis (Scherbaum et al. 2004)



In [ ]:

    
lh_analysis = Likelihood(gmpe_list, imt_list)
lh_analysis.get_residuals(db1)



In [ ]:

    
rspl.LikelihoodPlot(lh_analysis, "BooreAtkinson2008", "PGA")



In [ ]:

    
rspl.LikelihoodPlot(lh_analysis, "AkkarEtAlRjb2014", "PGA")



In [ ]:

    
rspl.LikelihoodPlot(lh_analysis, "AkkarEtAlRjb2014", "SA(1.0)")

Trends in Residuals: 1) Distance



In [ ]:

    
rspl.ResidualWithDistance(resid1, "AkkarEtAlRjb2014", "PGA", plot_type="linear", figure_size=(7, 7))



In [ ]:

    
rspl.ResidualWithDistance(resid1, "BooreAtkinson2008", "PGA", plot_type="linear", figure_size=(7, 7))

Trends in Residuals 2: Magnitude



In [ ]:

    
rspl.ResidualWithMagnitude(resid1, "AkkarEtAlRjb2014", "PGA", figure_size=(7,7))



In [ ]:

    
rspl.ResidualWithMagnitude(resid1, "BooreAtkinson2008", "PGA", figure_size=(7,7))

Trends in Residuals 3: Vs30



In [ ]:

    
rspl.ResidualWithVs30(resid1, "AkkarEtAlRjb2014", "PGA", figure_size=(7,7))



In [ ]:

    
rspl.ResidualWithVs30(resid1, "BooreAtkinson2008", "PGA", figure_size=(7,7))

Trends in Residuals 4: Depths



In [ ]:

    
rspl.ResidualWithDepth(resid1, "AkkarEtAlRjb2014", "PGA", figure_size=(7,7))



In [ ]:

    
rspl.ResidualWithDepth(resid1, "BooreAtkinson2008", "PGA", figure_size=(7,7))

Comparing the GMPEs using the LLH Metric (Scherbaum et al., 2009)



In [ ]:

    
# Create an instance of the residual calculator
llh_1 = LLH(gmpe_list, imt_list)
llh_1.get_residuals(db1)



In [ ]:

    
# Run LLH analysis
llh_values, weights = llh_1.get_loglikelihood_values()

for gmpe in gmpe_list:
    print "LLH(%s) = %12.6f  Weight = %12.6f" % (gmpe, llh_values[gmpe], weights[gmpe])

Comparing the GMPEs using the Euclidian Distance-Based Raking Metric (Kale & Akkar, 2013)



In [ ]:

    
# EDR analysis
edr1 = EDR(gmpe_list, imt_list)
edr1.get_residuals(db1)



In [ ]:

    
edr_values = edr1.get_edr_values()

for gmpe in gmpe_list:
    print "EDR(%s) = %8.4f (sqrt(kappa) = %8.4f)" %(gmpe,
                                                    edr_values[gmpe]["EDR"],
                                                    edr_values[gmpe]["sqrt Kappa"])

Single Station Analysis

This illustrates the application of a set of site specific analyses to determine the single-station phi. Implementation is based on the description and formulae found in Rodgriguez-Marek et al. (2011)

Find the Stations in the Database with the Most Recordings



In [ ]:

    
from smtk.strong_motion_selector import SMRecordSelector, rank_sites_by_record_count
import smtk.residuals.gmpe_residuals as res
import smtk.residuals.residual_plotter as rspl



In [ ]:

    
# Returns a list of the stations with the most records (in descending order)
# cutting off sites with less than 30 records
top_sites = rank_sites_by_record_count(db1, threshold=5)
for site_id in top_sites.keys():
    print "Site ID: %s   Name: %s  Number of Records: %s" %(site_id, top_sites[site_id]["Name"], top_sites[site_id]["Count"])

Calculate the ground motion residuals for each of the sites



In [ ]:

    
# Consider the AkkarEtAlRjb2014 GMPE with PGA, Sa (0.2), Sa(1.0)
gmpe_list = ["AkkarEtAlRjb2014"]
imt_list = ["PGA", "SA(0.2)", "SA(1.0)"]
# Get the residuals for the site
ssa1 = res.SingleStationAnalysis(top_sites.keys(), gmpe_list, imt_list)
ssa1.get_site_residuals(db1)

Plot the ordinary residuals for each site



In [ ]:

    
# Plots the GMPE Residuals for Each Site
rspl.ResidualWithSite(ssa1, "AkkarEtAlRjb2014", "PGA", figure_size=(7,7))



In [ ]:

    
rspl.ResidualWithSite(ssa1, "AkkarEtAlRjb2014", "SA(1.0)", figure_size=(7,7))

Plot the decomposition of the intra-event residuals



In [ ]:

    
# Plots the Intra-event analysis for each site
rspl.IntraEventResidualWithSite(ssa1, "AkkarEtAlRjb2014", "PGA", figure_size=(7,7))



In [ ]:

    
# Plots the Intra-event analysis for each site
rspl.IntraEventResidualWithSite(ssa1, "AkkarEtAlRjb2014", "SA(1.0)", figure_size=(7,7))



In [ ]:

    
for iloc, output in enumerate(ssa1.site_residuals):
    print "Site ID = %s" % (top_sites.keys()[iloc])
    for gmpe in output.site_analysis.keys():
        print "GMPE = %s" % gmpe 
        for imt_type in output.site_analysis[gmpe]:
            site_results = output.site_analysis[gmpe][imt_type]
            print "%10s: dS2ss = %10.8f    phi_ss,s = %10.8f" % (imt_type, site_results["dS2ss"], site_results["phi_ss,s"])
    print "======================================================="