Catalogue Query Tools

This notebook provides examples on how to use the Catalogue Toolkit to build and explore a catalogue database:



In [ ]:

    
%matplotlib inline
import os, sys
import numpy as np
import matplotlib.pyplot as plt
from eqcat.parsers.isf_catalogue_reader import ISFReader
import eqcat.catalogue_query_tools as cqt

Constructing the Database



In [ ]:

    
# Read in the catalogue
parser = ISFReader("inputs/isc_test_catalogue_isf.txt")
catalogue1 = parser.read_file("ISC_DB1", "ISC Global M >= 5")
print("Catalogue contains: %d events" % catalogue1.get_number_events())



In [ ]:

    
# Build the HDF5 Database
database_file = "outputs/catalogue_db1.hdf5"
if os.path.exists(database_file):
    os.remove(database_file)
_ = catalogue1.build_dataframe(hdf5_file=database_file)

Using the Database



In [ ]:

    
db1 = cqt.CatalogueDB(database_file)

Apply Limiting Selections

By Bounding Box



In [ ]:

    
lower_lon = 15.
upper_lon = 30.
lower_lat = 30.
upper_lat = 45.
bbox = [lower_lon, upper_lon, lower_lat, upper_lat]
selector = cqt.CatalogueSelector(db1)
aegean_cat = selector.select_within_bounding_box(bbox)



In [ ]:

    
number_origins, number_magnitudes = aegean_cat._get_number_origins_magnitudes()
print("Number of Origins = %d, Number of Magnitudes = %d" % (number_origins,
                                                             number_magnitudes))

By Polygon



In [ ]:

    
polygon = np.array([[15.0, 45.0],
                    [30.0, 45.0],
                    [30.0, 30.0],
                    [15.0, 30.0],
                    [15.0, 45.0]])
selector2 = cqt.CatalogueSelector(db1)
aegean_cat_alt = selector2.select_within_polygon(polygon[:, 0], polygon[:, 1])
number_origins, number_magnitudes = aegean_cat_alt._get_number_origins_magnitudes()
print("Number of Origins = %d, Number of Magnitudes = %d" % (number_origins,
                                                             number_magnitudes))

By Magnitude



In [ ]:

    
# Above magnitude 6.0
selector3 = cqt.CatalogueSelector(aegean_cat)
aegean_cat_m6 = selector3.select_within_magnitude_range(lower_mag=6.0, upper_mag=9.0)
number_origins, number_magnitudes = aegean_cat_m6._get_number_origins_magnitudes()
print("Number of Origins = %d, Number of Magnitudes = %d" % (number_origins,
                                                             number_magnitudes))

By Depth



In [ ]:

    
selector4 = cqt.CatalogueSelector(aegean_cat_alt)
aegean_cat_deep = selector4.select_within_depth_range(upper_depth=50.0, lower_depth=200.0)
number_origins, number_magnitudes = aegean_cat_deep._get_number_origins_magnitudes()
print("Number of Origins = %d, Number of Magnitudes = %d" % (number_origins,
                                                             number_magnitudes))

Exploring the Catalogue Database

See Summary of Agencies and Magnitudes in the Catalogue



In [ ]:

    
agency_count = cqt.get_agency_magtype_statistics(db1)

Search for a Specific Agency-Magnitude Combination

Search for body-wave magnitude common to 'BJI' and 'ISC'



In [ ]:

    
query_bji_isc_mb, bji_isc_mb_cat = cqt.get_agency_magnitude_pairs(db1, ("BJI", "mb"), ("ISC", "mb"), no_case=True)



In [ ]:

    
_ = cqt.plot_agency_magnitude_density(query_bji_isc_mb)

Join Query Results

Join together the results of two queries. For example the Global CMT magnitudes are reported under either 'GCMT' or 'HRVD'. So search for both terms



In [ ]:

    
query1, cat1 = cqt.get_agency_magnitude_pairs(db1, ("BJI", "Ms"), ("GCMT", "Mw"), no_case=True)
query2, cat2 = cqt.get_agency_magnitude_pairs(db1, ("BJI", "Ms"), ("HRVD", "Mw"), no_case=True)
query_bji_gcmt_ms = cqt.join_query_results(query1, query2)



In [ ]:

    
_ = cqt.plot_agency_magnitude_density(query_bji_gcmt_ms)

Regression Tools

In this example we compare the $M_S$ scale as recorded by the BJI network with the $M_W$ scale reported by HRVD/GCMT (from the previous query)

Set up the regression



In [ ]:

    
regressor = cqt.CatalogueRegressor(query_bji_gcmt_ms)
regressor.plot_density(overlay=False)

Apply a Linear Model



In [ ]:

    
linear_model = regressor.run_regression("polynomial", # Model Type
                                        [0., 1.]) # Initial guess of parameters
regressor.plot_model_density(False, 0)
# View Results
print("Mw = %.3f + %.3f MS +/- %.3f" % (regressor.results.beta[0],
                                        regressor.results.beta[1],
                                        regressor.standard_deviation))

Overlay another model defined as a Magnitude Conversion Rule



In [ ]:

    
from eqcat.isc_homogenisor import MagnitudeConversionRule
# Define empirical model
def RandomRule1(magnitude):
    return 1.21 + 0.84 * magnitude

def RandomRuleSigma(magnitude):
    return 0.2
# Create Rule
rule1 = MagnitudeConversionRule("BJI", "MS", RandomRule1, RandomRuleSigma,
                                model_name=r"$M_{W_{GCMT}} = 1.2 + 0.767 M_{S_{BJI}}$")
# Plot the model - with overla set to true
regressor.plot_model_density(True, 0)
# Overlay the rule and close the figure (overlay set to False)
regressor.plot_magnitude_conversion_model(rule1, False, line_color="b")

Apply a Piecewise Linear Model



In [ ]:

    
initial_guess = [1.0, 1.0, 0.0]  # [slope 1, slope 2, intercept]

linear_model = regressor.run_regression("2segmentM6.1", # Model Type
                                        initial_guess) # Initial guess of parameters
regressor.plot_model_density(False, 0)
print("Standard Deviation - Segment 1: %.3f, Segment 2: %.3f" % (regressor.standard_deviation[0],
                                                                 regressor.standard_deviation[1]))



In [ ]: