In [ ]:
%matplotlib inline
import os
import numpy as np
import matplotlib.pyplot as plt
# Import the Parsers
from eqcat.parsers.isf_catalogue_reader import ISFReader
from eqcat.isc_homogenisor import (HomogenisorPreprocessor,
DynamicHomogenisor,
MagnitudeConversionRule)
In [ ]:
parser = ISFReader("inputs/isc_test_catalogue_isf.txt",
selected_origin_agencies=["ISC", "GCMT", "HRVD", "NEIC", "EHB", "BJI"],
selected_magnitude_agencies=["ISC", "GCMT", "HRVD", "NEIC", "BJI"])
catalogue = parser.read_file("ISC_DB1", "ISC Global M >= 5")
print("Catalogue contains: %d events" % catalogue.get_number_events())
The catalogue covers the years 2005/06. To illustrate how to apply time variable hierarchies we consider two set of rules:
For the origin the order of preference is:
(For 2005): EHB, ISC, NEIC, GCMT/HRVD, BJI
(For 2006): ISC, EHB, NEIC, BJI, GCMT/HRVD
In [ ]:
origin_rules = [
("2005/01/01 - 2005/12/31", ['EHB', 'ISC', 'NEIC', 'GCMT', 'HRVD', 'BJI']),
("2006/01/01 - 2007/01/01", ['ISC', 'EHB', 'NEIC', 'BJI', 'GCMT', 'HRVD'])
]
In [ ]:
def gcmt_hrvd_mw(magnitude):
"""
For Mw recorded by GCMT take the value with no uncertainty
"""
return magnitude
def gcmt_hrvd_mw_sigma(magnitude):
"""
No additional uncertainty
"""
return 0.0
In [ ]:
def neic_mw(magnitude):
"""
If Mw reported by NEIC,
"""
return magnitude
def neic_mw_sigma(magnitude):
"""
Uncertainty of 0.11 units
"""
return 0.11
def scordillis_ms(magnitude):
"""
Scordilis (2006) indicates ISC and NEIC Ms can treated (almost) equivalently
"""
if magnitude < 6.1:
return 0.67 * magnitude + 2.07
else:
return 0.99 * magnitude + 0.08
def scordillis_ms_sigma(magnitude):
"""
With Magnitude dependent uncertainty
"""
if magnitude < 6.1:
return 0.17
else:
return 0.20
def scordillis_mb(magnitude):
"""
Scordilis (2006) finds NEIC and ISC mb nearly equivalent
"""
return 0.85 * magnitude + 1.03
def scordillis_mb_sigma(magnitude):
"""
"""
return 0.29
For BJI - no analysis has been undertaken. We apply a simple scaling of 0.9 M + 0.15 with uncertainty of 0.2. This is for illustrative purposes only
In [ ]:
def bji_mb(magnitude):
"""
"""
return 0.9 * magnitude + 0.15
def bji_mb_sigma(magnitude):
"""
"""
return 0.2
def bji_ms(magnitude):
"""
"""
return 0.9 * magnitude + 0.15
def bji_ms_sigma(magnitude):
"""
"""
return 0.2
In [ ]:
rule_set_2005 = [
MagnitudeConversionRule("GCMT", "Mw", gcmt_hrvd_mw, gcmt_hrvd_mw_sigma),
MagnitudeConversionRule("HRVD", "Mw", gcmt_hrvd_mw, gcmt_hrvd_mw_sigma),
MagnitudeConversionRule("ISC", "Ms", scordillis_ms, scordillis_ms_sigma),
MagnitudeConversionRule("NEIC", "Ms", scordillis_ms, scordillis_ms_sigma),
MagnitudeConversionRule("ISC", "mb", scordillis_mb, scordillis_mb_sigma),
MagnitudeConversionRule("NEIC", "mb", scordillis_mb, scordillis_mb_sigma),
MagnitudeConversionRule("BJI", "Ms", bji_ms, bji_ms_sigma),
MagnitudeConversionRule("BJI", "mb", bji_mb, bji_mb_sigma)
]
rule_set_2006 = [
MagnitudeConversionRule("GCMT", "Mw", gcmt_hrvd_mw, gcmt_hrvd_mw_sigma),
MagnitudeConversionRule("HRVD", "Mw", gcmt_hrvd_mw, gcmt_hrvd_mw_sigma),
MagnitudeConversionRule("ISC", "Ms", scordillis_ms, scordillis_ms_sigma),
MagnitudeConversionRule("BJI", "Ms", bji_ms, bji_ms_sigma),
MagnitudeConversionRule("NEIC", "Ms", scordillis_ms, scordillis_ms_sigma),
MagnitudeConversionRule("ISC", "mb", scordillis_mb, scordillis_mb_sigma),
MagnitudeConversionRule("BJI", "mb", bji_mb, bji_mb_sigma),
MagnitudeConversionRule("NEIC", "mb", scordillis_mb, scordillis_mb_sigma)
]
magnitude_rules = [
("2005/01/01 - 2005/12/31", rule_set_2005),
("2006/01/01 - 2007/01/01", rule_set_2006)
]
Before executing the homogenisation it is necessary to run a preprocessing step. This searches through the catalogue and identifies which conversion rule to apply:
The preprocessor is instantiated with a string describing the sort of rules to be applied.
"time" - Applies time only
"key" - Applies key rules only
"depth" - Applies depth rules only
"time|key" - Applies joint time and key rules
"time|depth" - Applied joint time and depth rules
"depth|key" - Applies joint depth and key rules
In [ ]:
preprocessor = HomogenisorPreprocessor("time")
catalogue = preprocessor.execute(catalogue, origin_rules, magnitude_rules)
In [ ]:
harmonisor = DynamicHomogenisor(catalogue, logging=True)
homogenised_catalogue = harmonisor.homogenise(magnitude_rules, origin_rules)
As logging was enabled, we can dump the log to a csv file and explore which rules and which hierarchy was applied for each event
In [ ]:
log_file = "outputs/homogenisor_log.csv"
if os.path.exists(log_file):
os.remove(log_file)
harmonisor.dump_log(log_file)
In [ ]:
output_catalogue_file = "outputs/homogeneous_catalogue.csv"
if os.path.exists(output_catalogue_file):
os.remove(output_catalogue_file)
harmonisor.export_homogenised_to_csv(output_catalogue_file)