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)

Load in Catalogue - Limit to ISC, GCMT/HRVD, EHB, NEIC, BJI


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())

Define Rule Sets

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'])
]

Magnitude Rules

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

ISC/NEIC


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

BJI

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

Define Magnitude Hierarchy


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)
]

Pre-processing

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)

Harmonise the Catalogue


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)

Export the Homogenised Catalogue to CSV


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)