In [9]:
%matplotlib inline
import os
import numpy as np
import matplotlib.pyplot as plt
import eqcat.parsers.isf_catalogue_reader as icr
import eqcat.catalogue_query_tools as cqt
In [10]:
raw_file = "Marmaries-catalogue1.txt"
rejection_keywords = ["mining", "geothermal", "explosion", "quarry", "reservoir", "induced", "rockburst"]
reader = icr.ISFReader(raw_file, rejection_keywords=rejection_keywords)
catalogue = reader.read_file("TUR", "ISC")
In [11]:
for event in reader.rejected_catalogue.events[:10]:
print event.id, event.description, event.comment
In [12]:
# Build HDF5 database
CAT_DB_FILE = "catalogue_db.hdf5"
_ = catalogue.build_dataframe(CAT_DB_FILE)
print "Comment this section out after first use!"
In [13]:
CAT_DB_FILE = "catalogue_db.hdf5"
db1 = cqt.CatalogueDB(CAT_DB_FILE)
In [14]:
db1.magnitudes
Out[14]:
In [15]:
temp = db1.origins["eventID"]
In [18]:
# Set up the configuration of the limits
outlier_fraction = 0.001
llon = np.floor(db1.origins.longitude.quantile(outlier_fraction))
ulon = np.ceil(db1.origins.longitude.quantile(1 - outlier_fraction))
llat = np.floor(db1.origins.latitude.quantile(outlier_fraction))
ulat = np.ceil(db1.origins.latitude.quantile(1 - outlier_fraction))
map_config = {"llon": llon, "ulon": ulon, "llat": llat, "ulat": ulat,
"parallel": 5.0, "meridian": 5.0, "resolution": "l"}
cqt.plot_catalogue_map(map_config, db1)
In [19]:
selector = cqt.CatalogueSelector(db1)
turkey_catalogue = selector.select_within_magnitude_range(5,5.5)
cqt.plot_catalogue_map(map_config, turkey_catalogue)
In [20]:
selector = cqt.CatalogueSelector(db1)
turkey_catalogue = selector.select_within_magnitude_range(5,6)
cqt.plot_catalogue_map(map_config, turkey_catalogue)
In [21]:
agency_magnitude_stats = cqt.get_agency_magtype_statistics(db1)
In [22]:
cqt.mine_agency_magnitude_combinations(db1,agency_magnitude_stats,200,no_case=True)
Out[22]:
In [23]:
# place ipynb in same folder as
file_name = "AgencyPairs.txt"
pairs = np.genfromtxt(file_name, dtype="S64,S64,S64,S64,S64,S64,S64,i,S64")
In [ ]:
# have a peek at the data
pairs[:5]
In [ ]:
# parse the data and look at sets
agencies1 = [pair[1][1:-1].replace('_','') for pair in pairs]
mags1 = [pair[2][:-1].upper() for pair in pairs]
agencies2 = [pair[4][1:-1].replace('_','') for pair in pairs]
mags2 = [pair[5][:-1].upper() for pair in pairs]
nums_shared = [pair[7] for pair in pairs]
print set(mags1)
print set(mags2)
print list(set(mags1) - set(mags2))
print list(set(mags2) - set(mags1))
print set(agencies1)
print set(agencies2)
print list(set(agencies1) - set(agencies2))
print list(set(agencies2) - set(agencies1))
In [ ]:
# prioritize magnitudes and print pairs in order of size of overlap
mag_priority = ['MW','Mw']
tex_format = '%s%s/"%s, %s"[input] -> [edge label={%d}] %s%s/"%s, %s"[input];'
data = zip(agencies1, mags1, agencies2, mags2, nums_shared)
data = list(set(data))
data = sorted(data, key=lambda pair: pair[4], reverse=True)
for i, m1 in enumerate(mag_priority):
for m2 in mag_priority[i:]:
for agency1, mag1, agency2, mag2, n_shared in data:
if (mag1 == m1 and mag2 == m2) or (mag1 == m2 and mag2 == m1):
if mag_priority.index(m1) >= mag_priority.index(m2):
print tex_format % (agency1, mag1, agency1, mag1, n_shared,
agency2, mag2, agency2, mag2)
else:
print tex_format % (agency2, mag2, agency2, mag2, n_shared,
agency1, mag1, agency1, mag1)
In [ ]:
ISK_MD_CSEM_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("CSEM", "MW"), ("NIC", "MW"), no_case=True)
#_ = cqt.plot_agency_magnitude_density(ISK_MD_CSEM_MD)
#cqt.plot_catalogue_map(map_config, query3_cat)
In [ ]:
ISK_MD_CSEM_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("ISK", "MD"), ("CSEM", "MD"), no_case=True)
regressor1 = cqt.CatalogueRegressor(ISK_MD_CSEM_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
ISK_MD_ATH_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("NIC", "MW"), ("NIC", "MW"), no_case=True)
_ = cqt.plot_agency_magnitude_density(ISK_MD_ATH_MD)
In [ ]:
ISK_MD_ATH_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("NIC", "MW"), ("NIC", "MW"), no_case=True)
regressor1 = cqt.CatalogueRegressor(ISK_MD_ATH_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
regressor1 = cqt.CatalogueRegressor(ISK_MD_ATH_MD)
results = regressor1.run_regression("2segmentM4.3", [1.0, 1.0, 1.0])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
print regressor1.standard_deviation
In [ ]:
ATH_MD_CSEM_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("ATH", "MD"), ("CSEM", "MD"), no_case=True)
_ = cqt.plot_agency_magnitude_density(ATH_MD_CSEM_MD)
In [ ]:
ATH_MD_CSEM_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("ATH", "MD"), ("CSEM", "MD"), no_case=True)
regressor1 = cqt.CatalogueRegressor(ATH_MD_CSEM_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
CSEM_MD_DDA_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("CSEM", "MD"), ("DDA", "MD"), no_case=True)
_ = cqt.plot_agency_magnitude_density(CSEM_MD_DDA_MD)
In [ ]:
CSEM_MD_DDA_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("CSEM", "MD"), ("DDA", "MD"), no_case=True)
regressor1 = cqt.CatalogueRegressor(CSEM_MD_DDA_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
ISK_MD_DDA_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("ISK", "MD"), ("DDA", "MD"), no_case=True)
_ = cqt.plot_agency_magnitude_density(ISK_MD_DDA_MD)
In [ ]:
ISK_MD_DDA_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("ISK", "MD"), ("DDA", "MD"), no_case=True)
regressor1 = cqt.CatalogueRegressor(ISK_MD_DDA_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
ATH_MD_DDA_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("ATH", "MD"), ("DDA", "MD"), no_case=True)
_ = cqt.plot_agency_magnitude_density(ATH_MD_DDA_MD)
In [ ]:
ATH_MD_DDA_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("DDA", "MD"), ("ATH", "MD"), no_case=True)
regressor1 = cqt.CatalogueRegressor(ATH_MD_DDA_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
regressor1b = cqt.CatalogueRegressor(ATH_MD_DDA_MD)
results = regressor1b.run_regression("exponential", [1.0, 1.0, 1.0])
regressor1b.results.pprint()
regressor1b.plot_model(overlay=False)
print regressor1b.standard_deviation
In [ ]:
DDA_MD_HLW_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("DDA", "MD"), ("HLW", "MD"), no_case=True)
_ = cqt.plot_agency_magnitude_density(DDA_MD_HLW_MD)
In [ ]:
DDA_MD_HLW_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("DDA", "MD"), ("HLW", "MD"), no_case=True)
regressor1 = cqt.CatalogueRegressor(DDA_MD_HLW_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
regressor1 = cqt.CatalogueRegressor(DDA_MD_HLW_MD)
results = regressor1.run_regression("2segmentM4.3", [1.0, 1.0, 1.0])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
print regressor1.standard_deviation
In [ ]:
DDA_MD_HLW_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("ISC", "MB"), ("ATH", "MD"), no_case=True)
In [ ]:
DDA_MD_HLW_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("ISC", "MB"), ("ATH", "MD"), no_case=True)
regressor1 = cqt.CatalogueRegressor(DDA_MD_HLW_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
DDA_MD_HLW_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("ISC", "MB"), ("ISK", "MD"), no_case=True)
In [ ]:
DDA_MD_HLW_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("ISC", "MB"), ("ISK", "MD"), no_case=True)
regressor1 = cqt.CatalogueRegressor(DDA_MD_HLW_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
DDA_MD_HLW_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("NEIC", "MB"), ("ATH", "MD"), no_case=True)
In [ ]:
DDA_MD_HLW_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("NEIC", "MB"), ("ATH", "MD"), no_case=True)
regressor1 = cqt.CatalogueRegressor(DDA_MD_HLW_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
DDA_MD_HLW_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("NEIC", "MB"), ("ISK", "MD"), no_case=True)
In [ ]:
DDA_MD_HLW_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("NEIC", "MB"), ("ISK", "MD"), no_case=True)
regressor1 = cqt.CatalogueRegressor(DDA_MD_HLW_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
DDA_MD_HLW_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("NIC", "MW"), ("NIC", "MW"), no_case=True)
In [ ]:
DDA_MD_HLW_MD, query4_cat = cqt.get_agency_magnitude_pairs(db1, ("NEIC", "MS"), ("CSEM", "MS"), no_case=True)
regressor1 = cqt.CatalogueRegressor(DDA_MD_HLW_MD)
results = regressor1.run_regression("polynomial", [1.02, -0.18])
regressor1.results.pprint()
regressor1.plot_model(overlay=False)
s=regressor1.standard_deviation
print "Model standard deviation is"
print(s)
In [ ]:
DDA_MD_HLW_MD, query3_cat = cqt.get_agency_magnitude_pairs(db1, ("ISK", "MD"), ("NIC", "MW"), no_case=True)
In [ ]:
# origin rules
origin_rules = [
("1900/10/01 - 2015/10/01", ["ZUR_RMT","HRVD","NIC","CSEM","MED_RCMT","NEIC","ATH","ISK","ISC","IDC","DDA","ISCJB","IDC","MOS"])
]
In [ ]:
# Magnitude rules
def CSEM_MW_NIC_MW(magnitude):
return magnitude
def CSEM_MW_NIC_MW_sigma(magnitude):
return 0.0
def ISCJB_MS_CSEM_MS(magnitude):
return 1.135*magnitude-0.765
def ISCJB_MS_CSEM_MS_sigma(magnitude):
return 0.26369784
def NEIC_MS_CSEM_MS(magnitude):
return 1.076*magnitude-0.458
def NEIC_MS_CSEM_MS_sigma(magnitude):
return 0.0648769
def IDC_MS_CSEM_MS(magnitude):
return 1.242*magnitude-1.087
def IDC_MS_CSEM_MS_sigma(magnitude):
return 0.314273
def MOS_MS_CSEM_MS(magnitude):
return 1.276*magnitude-1.328
def MOS_MS_CSEM_MS_sigma(magnitude):
return 0.33142032
def ATH_MD_NIC_MW(magnitude):
return 0.937*magnitude-0.224
def ATH_MD_NIC_MW_sigma(magnitude):
return 0.32
def ISC_MB_IDC_MB(magnitude):
return 0.801*magnitude+0.667
def ISC_MB_IDC_MB_sigma(magnitude):
return 0.1563
def NIC_ML_NIC_MW(magnitude):
return 0.767*magnitude+0.572
def NIC_ML_NIC_MW_sigma(magnitude):
return 0.24669
def DDA_ML_NIC_ML(magnitude):
return 0.789*magnitude+0.754
def DDA_ML_NIC_ML_sigma(magnitude):
return 0.2613
def ISC_MS_CSEM_MS(magnitude):
return 1.092*magnitude-0.586
def ISC_MS_CSEM_MS_sigma(magnitude):
return 0.2254
def CSEM_MS_CSEM_MW(magnitude):
return 0.867*magnitude+1.088
def CSEM_MS_CSEM_MW_sigma(magnitude):
return 0.33265
def IDC_MB_NIC_MW(magnitude):
return 0.751*magnitude+0.705
def IDC_MB_NIC_MW_sigma(magnitude):
return 0.23665
def CSEM_ML_NIC_ML(magnitude):
return 0.699*magnitude+1.234
def CSEM_ML_NIC_ML_sigma(magnitude):
return 0.2567972
def DDA_MD_ATH_MD(magnitude):
return 0.954*magnitude+0.327
def DDA_MD_ATH_MD_sigma(magnitude):
return 0.232736
def ISK_ML_NIC_ML(magnitude):
return 0.755*magnitude+0.931
def ISK_ML_NIC_ML_sigma(magnitude):
return 0.2480
def ISK_MD_ATH_MD(magnitude):
return 1.085*magnitude-0.147
def ISK_MD_ATH_MD_sigma(magnitude):
return 0.2346155
def CSEM_MD_ATH_MD(magnitude):
return 0.929*magnitude+0.286
def CSEM_MD_ATH_MD_sigma(magnitude):
return 0.151265
In [ ]:
from eqcat.isc_homogenisor import MagnitudeConversionRule, DynamicHomogenisor, HomogenisorPreprocessor
In [ ]:
magnitude_rule_set = [
MagnitudeConversionRule("HRVD", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("ZUR_RMT", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("MED_RCMT", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("CSEM", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("GCMT", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("NIC", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("NEIC", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("ISK", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("DDA", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("GII", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("BGS", "MW", CSEM_MW_NIC_MW, CSEM_MW_NIC_MW_sigma),
MagnitudeConversionRule("CSEM", "MS", CSEM_MS_CSEM_MW, CSEM_MS_CSEM_MW_sigma),
MagnitudeConversionRule("ISCJB", "MS", ISCJB_MS_CSEM_MS, ISCJB_MS_CSEM_MS_sigma),
MagnitudeConversionRule("NEIC", "MS", NEIC_MS_CSEM_MS, NEIC_MS_CSEM_MS_sigma),
MagnitudeConversionRule("IDC", "MS", IDC_MS_CSEM_MS, IDC_MS_CSEM_MS_sigma),
MagnitudeConversionRule("ISCJB", "MS", ISCJB_MS_CSEM_MS, ISCJB_MS_CSEM_MS_sigma),
MagnitudeConversionRule("NEIC", "MS", NEIC_MS_CSEM_MS, NEIC_MS_CSEM_MS_sigma),
MagnitudeConversionRule("IDC", "MS", IDC_MS_CSEM_MS, IDC_MS_CSEM_MS_sigma),
MagnitudeConversionRule("MOS", "MS", MOS_MS_CSEM_MS, MOS_MS_CSEM_MS_sigma),
MagnitudeConversionRule("ATH", "MD", ATH_MD_NIC_MW, ATH_MD_NIC_MW_sigma),
MagnitudeConversionRule("NIC", "ML", NIC_ML_NIC_MW, NIC_ML_NIC_MW_sigma),
MagnitudeConversionRule("IDC", "MB", IDC_MB_NIC_MW, IDC_MB_NIC_MW_sigma),
MagnitudeConversionRule("ISC", "MS", ISC_MS_CSEM_MS, ISC_MS_CSEM_MS_sigma),
MagnitudeConversionRule("DDA", "MD", DDA_MD_ATH_MD, DDA_MD_ATH_MD_sigma),
MagnitudeConversionRule("ISK", "MD", ISK_MD_ATH_MD, ISK_MD_ATH_MD_sigma),
MagnitudeConversionRule("CSEM", "MD", CSEM_MD_ATH_MD, CSEM_MD_ATH_MD_sigma),
MagnitudeConversionRule("DDA", "ML", DDA_ML_NIC_ML, DDA_ML_NIC_ML_sigma),
MagnitudeConversionRule("ISC", "MB", ISC_MB_IDC_MB, ISC_MB_IDC_MB_sigma),
MagnitudeConversionRule("CSEM", "ML", CSEM_ML_NIC_ML, CSEM_ML_NIC_ML_sigma),
MagnitudeConversionRule("ISK", "ML", ISK_ML_NIC_ML, ISK_ML_NIC_ML_sigma),
]
magnitude_rules = [
("1900/10/01 - 2015/10/01", magnitude_rule_set)
]
In [ ]:
preprocessor = HomogenisorPreprocessor(rule_type="time")
preprocessed_catalogue = preprocessor.execute(catalogue, origin_rules, magnitude_rules)
homogenisor = DynamicHomogenisor(preprocessed_catalogue, logging=True)
In [ ]:
output_catalogue = homogenisor.homogenise(magnitude_rules, origin_rules)
In [ ]:
homogenisor.dump_log("homogenisation1_logfile.csv")
homogenisor.export_homogenised_to_csv("output_homoge_Marmaris.csv")
In [25]:
(2000 - 500)/12/10
Out[25]:
In [ ]: