<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Catalogue Query Tools
This notebook provides examples on how to use the Catalogue Toolkit to build and explore a catalogue database
Step1:  Constructing the Database
Step2:  Using the Database
Step3:  Apply Limiting Selections
By Bounding Box
Step4:  By Polygon
Step5:  By Magnitude
Step6:  By Depth
Step7:  Exploring the Catalogue Database
See Summary of Agencies and Magnitudes in the Catalogue
Step8:  Search for a Specific Agency-Magnitude Combination
Search for body-wave magnitude common to 'BJI' and 'ISC'
Step9:  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
Step10:  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
Step11:  Apply a Linear Model
Step12:  Overlay another model defined as a Magnitude Conversion Rule
Step13:  Apply a Piecewise Linear Model
<ASSISTANT_TASK:>
Python Code:
%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
Explanation: Catalogue Query Tools
This notebook provides examples on how to use the Catalogue Toolkit to build and explore a catalogue database:
End of explanation
# 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())
# 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)
Explanation: Constructing the Database
End of explanation
db1 = cqt.CatalogueDB(database_file)
Explanation: Using the Database
End of explanation
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)
number_origins, number_magnitudes = aegean_cat._get_number_origins_magnitudes()
print("Number of Origins = %d, Number of Magnitudes = %d" % (number_origins,
                                                             number_magnitudes))
Explanation: Apply Limiting Selections
By Bounding Box
End of explanation
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))
Explanation: By Polygon
End of explanation
# 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))
Explanation: By Magnitude
End of explanation
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))
Explanation: By Depth
End of explanation
agency_count = cqt.get_agency_magtype_statistics(db1)
Explanation: Exploring the Catalogue Database
See Summary of Agencies and Magnitudes in the Catalogue
End of explanation
query_bji_isc_mb, bji_isc_mb_cat = cqt.get_agency_magnitude_pairs(db1, ("BJI", "mb"), ("ISC", "mb"), no_case=True)
_ = cqt.plot_agency_magnitude_density(query_bji_isc_mb)
Explanation: Search for a Specific Agency-Magnitude Combination
Search for body-wave magnitude common to 'BJI' and 'ISC'
End of explanation
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)
_ = cqt.plot_agency_magnitude_density(query_bji_gcmt_ms)
Explanation: 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
End of explanation
regressor = cqt.CatalogueRegressor(query_bji_gcmt_ms)
regressor.plot_density(overlay=False)
Explanation: 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
End of explanation
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))
Explanation: Apply a Linear Model
End of explanation
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")
Explanation: Overlay another model defined as a Magnitude Conversion Rule
End of explanation
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]))
Explanation: Apply a Piecewise Linear Model
End of explanation
<END_TASK>