In [4]:
# Catalogue Homogenization
In [5]:
%matplotlib inline
import os
import sys
import csv
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import eqcat.parsers.isf_catalogue_reader as icr
import eqcat.catalogue_query_tools as cqt
from IPython.utils import io
from eqcat.isc_homogenisor import \
MagnitudeConversionRule, DynamicHomogenisor, HomogenisorPreprocessor
In [6]:
# read the catalogue - why do we bother making an hdf5?
raw_file_name = "Marmaries-catalogue1.txt"
base = os.path.basename(raw_file_name)
db_file_name = os.path.splitext(base)[0] + '.hdf5'
rejection_keywords = ["mining", "geothermal", "explosion", "quarry",
"reservoir", "induced", "rockburst"]
reader = icr.ISFReader(raw_file_name,
rejection_keywords=rejection_keywords)
catalogue = reader.read_file("TUR", "ISC")
In [ ]:
# Mw equivalences considered exact after verification
def equivalent(M):
return M
def zero_sigma(M):
return 0.0
In [ ]:
def from_NIC_MW(M):
"""
Agency-Pairs: (NIC, MW) & (MED_RCMT, MW) returned 85 events
Potential yeild 561 magnitudes
"""
return 0.9 + M
def from_NIC_MW_sigma(M):
return 1.0
def from_CSEM_MW(M):
"""
Agency-Pairs: (CSEM, MW) & (NIC, MW) returned 238 events
Potential yeild 287 magnitudes
"""
return from_NIC_MW(M)
def from_CSEM_MW_sigma(M):
return from_NIC_MW_sigma(M)
def from_IDC_MS(M):
"""
Agency-Pairs: (IDC, MS) & (MED_RCMT, MW) returned 196 events
Potential yeild 1279 magnitudes
"""
return 2.080 + 0.700*M
def from_IDC_MS_sigma(M):
return 0.2
def from_ISC_MS(M):
"""
Agency-Pairs: (ISC, MS) & (MED_RCMT, MW) returned 149 events
Potential yeild 762 magnitudes
"""
return 2.244 + 0.629*M
def from_ISC_MS_sigma(M):
return 0.2
def from_ISCJB_MS(M):
"""
Agency-Pairs: (ISCJB, MS) & (MED_RCMT, MW) returned 79 events
Potential yeild 403 magnitudes
"""
return 2.022 + 0.691*M
def from_ISCJB_MS_sigma(M):
return 0.2
def from_MOS_MS(M):
"""
Agency-Pairs: (MOS, MS) & (MED_RCMT, MW) returned 64 events
Potential yeild 291 magnitudes
"""
return 1.611 + 0.796*M
def from_MOS_MS_sigma(M):
return 0.3
def from_BJI_MS(M):
"""
Agency-Pairs: (BJI, MS) & (MED_RCMT, MW) returned 125 events
Potential yeild 241 magnitudes
"""
return -0.115 + 1.042*M
def from_BJI_MS_sigma(M):
return 0.25
def from_CSEM_MS(M):
"""
Agency-Pairs: (CSEM, MS) & (MED_RCMT, MW) returned 51 events
Potential yeild 128 magnitudes
"""
return 2.844 + 0.532*M
def from_CSEM_MS_sigma(M):
return 0.3
def from_NEIC_MS(M):
"""
Agency-Pairs: (NEIC, MS) & (MED_RCMT, MW) returned 29 events
Potential yeild 65 magnitudes
"""
return 2.345 + 0.612*M
def from_NEIC_MS_sigma(M):
return 0.5
def from_NEIC_MS(M):
"""
Agency-Pairs: (NEIC, MS) & (HRVD, MW) returned 23 events
Potential yeild 65 magnitudes
"""
return 2.389 + 0.605*M
def from_NEIC_MS_sigma(M):
return 0.5
def from_IDC_MS(M):
"""
Agency-Pairs: (IDC, MS) & (ISC, MS) returned 446 events
Potential yeild 1279 magnitudes
"""
return from_ISC_MS(M)
def from_IDC_MS_sigma(M):
return math.sqrt(0.19**2 + from_ISC_MS_sigma(M)**2)
def from_ATH_MD(M):
"""
Agency-Pairs: (ATH, MD) & (MED_RCMT, MW) returned 87 events
Potential yeild 14988 magnitudes
"""
return 0.245 + 1.025*M
def from_ATH_MD_sigma(M):
return 0.246
def from_ISK_MD(M):
"""
Agency-Pairs: (ISK, MD) & (MED_RCMT, MW) returned 55 events
Potential yeild 17210 magnitudes
"""
return 0.009 + 1.070*M
def from_ISK_MD_sigma(M):
return 0.27
def from_HLW_MD(M):
"""
Agency-Pairs: (HLW, MD) & (MED_RCMT, MW) returned 47 events
Potential yeild 510 magnitudes
"""
return 0.2 + M
def from_HLW_MD_sigma(M):
return 0.7
def from_DDA_MD(M):
"""
Agency-Pairs: (DDA, MD) & (MED_RCMT, MW) returned 23 events
Potential yeild 6172 magnitudes
"""
return 0.7 + M
def from_DDA_MD_sigma(M):
return 0.5
def from_GII_MD(M):
"""
Agency-Pairs: (GII, MD) & (MED_RCMT, MW) returned 34 events
Potential yeild 241 magnitudes
"""
return 1.388 + 0.776*M
def from_GII_MD_sigma(M):
return 0.4
def from_CSEM_MD(M):
"""
Agency-Pairs: (CSEM, MD) & (ISK, MD) returned 8512 events
Potential yeild 10863 magnitudes
"""
return from_ISK_MD(M)
def from_CSEM_MD_sigma(M):
return math.sqrt(0.129**2 + from_ISK_MD_sigma(M)**2)
def from_ISC_MB(M):
"""
Agency-Pairs: (ISC, MB) & (MED_RCMT, MW) returned 182 events
Potential yeild 3953 magnitudes
"""
return 0.108 + 1.015*M
def from_ISC_MB_sigma(M):
return 0.224
def from_IDC_MB(M):
"""
Agency-Pairs: (IDC, MB) & (MED_RCMT, MW) returned 205 events
Potential yeild 2455 magnitudes
"""
return 0.197 + 1.056*M
def from_IDC_MB_sigma(M):
return 0.241
def from_NEIC_MB(M):
"""
Agency-Pairs: (NEIC, MB) & (MED_RCMT, MW) returned 182 events
Potential yeild 1849 magnitudes
"""
return 0.113 + 1.005*M
def from_NEIC_MB_sigma(M):
return 0.225
def from_ISCJB_MB(M):
"""
Agency-Pairs: (ISCJB, MB) & (MED_RCMT, MW) returned 91 events
Potential yeild 1161 magnitudes
"""
return 0.110 + 1.017*M
def from_ISCJB_MB_sigma(M):
return 0.2
def from_MOS_MB(M):
"""
Agency-Pairs: (MOS, MB) & (MED_RCMT, MW) returned 189 events
Potential yeild 936 magnitudes
"""
return 0.334 + 0.927*M
def from_MOS_MB_sigma(M):
return 0.225
def from_NIC_MB(M):
"""
Agency-Pairs: (NIC, MB) & (MED_RCMT, MW) returned 128 events
Potential yeild 689 magnitudes
"""
return -1.957 + 1.416*M
def from_NIC_MB_sigma(M):
return 0.276
def from_BJI_MB(M):
"""
Agency-Pairs: (BJI, MB) & (MED_RCMT, MW) returned 161 events
Potential yeild 602 magnitudes
"""
return -1.762 + 1.368*M
def from_BJI_MB_sigma(M):
return 0.283
def from_CSEM_MB(M):
"""
Agency-Pairs: (CSEM, MB) & (MED_RCMT, MW) returned 120 events
Potential yeild 425 magnitudes
"""
return -0.098 + 1.034*M
def from_CSEM_MB_sigma(M):
return 0.218
def from_ISK_ML(M):
"""
Agency-Pairs: (ISK, ML) & (MED_RCMT, MW) returned 162 events
Potential yeild 6545 magnitudes
"""
return 0.578 + 0.893*M
def from_ISK_ML_sigma(M):
return 0.188
def from_ATH_ML(M):
"""
Agency-Pairs: (ATH, ML) & (MED_RCMT, MW) returned 175 events
Potential yeild 6365 magnitudes
"""
return 0.480 + 0.935*M
def from_ATH_ML_sigma(M):
return 0.22
def from_THE_ML(M):
"""
Agency-Pairs: (THE, ML) & (MED_RCMT, MW) returned 172 events
Potential yeild 4637 magnitudes
"""
return 0.676 + 0.899*M
def from_THE_ML_sigma(M):
return 0.275
def from_DDA_ML(M):
"""
Agency-Pairs: (DDA, ML) & (MED_RCMT, MW) returned 73 events
Potential yeild 3350 magnitudes
"""
return 0.483 + 0.910*M
def from_DDA_ML_sigma(M):
return 0.276
def from_IDC_ML(M):
"""
Agency-Pairs: (IDC, ML) & (MED_RCMT, MW) returned 167 events
Potential yeild 2050 magnitudes
"""
return -0.239 + 1.231*M
def from_IDC_ML_sigma(M):
return 0.421
def from_NIC_ML(M):
"""
Agency-Pairs: (NIC, ML) & (MED_RCMT, MW) returned 129 events
Potential yeild 1147 magnitudes
"""
return -0.087 + 1.097*M
def from_NIC_ML_sigma(M):
return 0.261
def from_HLW_ML(M):
"""
Agency-Pairs: (HLW, ML) & (MED_RCMT, MW) returned 55 events
Potential yeild 589 magnitudes
"""
return 0.3 + M
def from_HLW_ML_sigma(M):
return 0.313
def from_CSEM_ML(M):
"""
Agency-Pairs: (CSEM, ML) & (ISK, ML) returned 1969 events
Potential yeild 2999 magnitudes
"""
return from_ISK_ML(M)
def from_CSEM_ML_sigma(M):
return math.sqrt(0.229**2 + from_ISK_ML_sigma(M)**2)
In [ ]:
# origin rules
origin_rules = [
("1900/10/01 - 2015/10/01", [
"IDC", "ZUR_RMT", "THE",
"ISK", "ISCJB", "ISC", "IDC", "ATH", # high resolution + std
"MED_RCMT", "HRVD", "GCMT", # low resolution + std
"NEIC", "HLW", # no error estimate
"MOS", "GII", "DDA", "BJI", "CSEM", "NIC"])
]
In [ ]:
magnitude_rule_set = [
MagnitudeConversionRule("MED_RCMT", "MW", equivalent, zero_sigma),
MagnitudeConversionRule("ZUR_RMT", "MW", equivalent, zero_sigma),
MagnitudeConversionRule("GCMT", "MW", equivalent, zero_sigma),
MagnitudeConversionRule("HRVD", "MW", equivalent, zero_sigma),
MagnitudeConversionRule("NEIC", "MW", equivalent, zero_sigma),
MagnitudeConversionRule("ISK", "ML", from_ISK_ML, from_ISK_ML_sigma), # 0.188
MagnitudeConversionRule("ISCJB", "MB", from_ISCJB_MB, from_ISCJB_MB_sigma), # 0.2
MagnitudeConversionRule("CSEM", "MB", from_CSEM_MB, from_CSEM_MB_sigma), # 0.218
MagnitudeConversionRule("ATH", "ML", from_ATH_ML, from_ATH_ML_sigma), # 0.22
MagnitudeConversionRule("GII", "MB", from_ISC_MB, from_ISC_MB_sigma), # 0.224
MagnitudeConversionRule("NEIC", "MB", from_NEIC_MB, from_NEIC_MB_sigma), # 0.225
MagnitudeConversionRule("MOS", "MB", from_MOS_MB, from_MOS_MB_sigma), # 0.225
MagnitudeConversionRule("GII", "MB", from_IDC_MB, from_IDC_MB_sigma), # 0.241
MagnitudeConversionRule("ATH", "MD", from_ATH_MD, from_ATH_MD_sigma), # 0.246
MagnitudeConversionRule("IDC", "MS", from_IDC_MS, from_IDC_MS_sigma), # 0.25
MagnitudeConversionRule("ISC", "MS", from_ISC_MS, from_ISC_MS_sigma), # 0.25
MagnitudeConversionRule("ISCJB", "MS", from_ISCJB_MS, from_ISCJB_MS_sigma), # 0.25
MagnitudeConversionRule("NIC", "ML", from_NIC_ML, from_NIC_ML_sigma), # 0.261
MagnitudeConversionRule("ISK", "MD", from_ISK_MD, from_ISK_MD_sigma), # 0.27
MagnitudeConversionRule("THE", "ML", from_THE_ML, from_THE_ML_sigma), # 0.275
MagnitudeConversionRule("DDA", "ML", from_DDA_ML, from_DDA_ML_sigma), # 0.276
MagnitudeConversionRule("NIC", "MB", from_NIC_MB, from_NIC_MB_sigma), # 0.276
MagnitudeConversionRule("BJI", "MB", from_BJI_MB, from_BJI_MB_sigma), # 0.283
MagnitudeConversionRule("MOS", "MS", from_MOS_MS, from_MOS_MS_sigma), # 0.3
MagnitudeConversionRule("HLW", "ML", from_HLW_ML, from_HLW_ML_sigma), # 0.313
MagnitudeConversionRule("BJI", "MS", from_BJI_MS, from_BJI_MS_sigma), # 0.35
MagnitudeConversionRule("CSEM", "MS", from_CSEM_MS, from_CSEM_MS_sigma), # 0.35
MagnitudeConversionRule("CSEM", "ML", from_CSEM_ML, from_CSEM_ML_sigma), # 0.4
MagnitudeConversionRule("GII", "MD", from_GII_MD, from_GII_MD_sigma), # 0.4
MagnitudeConversionRule("IDC", "ML", from_IDC_ML, from_IDC_ML_sigma), # 0.421
MagnitudeConversionRule("NEIC", "MS", from_NEIC_MS, from_NEIC_MS_sigma), # 0.5
MagnitudeConversionRule("DDA", "MD", from_DDA_MD, from_DDA_MD_sigma), # 0.5
MagnitudeConversionRule("CSEM", "MD", from_CSEM_MD, from_CSEM_MD_sigma), # 0.5
# MagnitudeConversionRule("HLW", "MD", from_HLW_MD, from_HLW_MD_sigma), # 0.7
# MagnitudeConversionRule("CSEM", "MW", from_CSEM_MW, from_CSEM_MW_sigma), # 1.0
# MagnitudeConversionRule("NIC", "MW", from_CSEM_ML, from_NIC_MW_sigma), # 1.0
]
In [ ]:
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 [ ]:
with io.capture_output() as captured:
output_catalogue = homogenisor.homogenise(
magnitude_rules, origin_rules)
In [ ]:
log_file_name = os.path.splitext(base)[0] + '_exclude_logfile.csv'
homogenisor.dump_log(log_file_name)
homogenised_file_name = os.path.splitext(base)[0] + '_exclude_homogenised.csv'
homogenisor.export_homogenised_to_csv(homogenised_file_name)
In [ ]:
# identify clusters
declusterer = decluster.dec_gardner_knopoff.GardnerKnopoffType1()
decluster_config = {
'time_distance_window': decluster.distance_time_windows.UhrhammerWindow(),
'fs_time_prop': 1.0}
cluster_index, cluster_flag = declusterer.decluster(catalogue, decluster_config)
# purge catalogue
declustered = deepcopy(catalogue)
mainshock_flag = cluster_flag == 0
declustered.purge_catalogue(mainshock_flag)
In [ ]:
# write declustered catalog (doesn't work yet)
base = os.path.basename(catalogue_filename)
declustered_file_name = os.path.splitext(base)[0] + '_declustered.csv'
keys = ['eventID', 'Agency', 'magnitudeType',
'year', 'month', 'day', 'hour', 'minute', 'second', 'timeError',
'longitude', 'latitude', 'SemiMajor90', 'SemiMinor90', 'ErrorStrike',
'depth', 'depthError', 'magnitude', 'sigmaMagnitude']
with open(declustered_file_name, 'wb') as outfile:
writer = csv.writer(outfile, delimiter = ',')
writer.writerow(keys)
writer.writerows(zip(*[declustered.data[key] for key in keys]))