This notebook contains:
Read a seismic catalogue and pre-defined seismic sources.
Explore basic methods for qualitative analysis of the catalogue.
Declustering: Aftershocks and foreshocks identification and removal.
Completeness
Estimation of the Gutenberg-Richter parameters.
Exploration of some methods to estimate the maximum magnitude (statistically from the catalogue)
Smoothed seismicity
In this part there are imported the methods from the hmtk library and other dependencies to run the analysis.
In [ ]:
# Dependencies
import os
import numpy as np
from copy import deepcopy
import matplotlib.pyplot as plt
# Classes for importing sources
from hmtk.sources import (simple_fault_source,
complex_fault_source,
area_source)
# Catalogue and sources
from hmtk.parsers.catalogue import CsvCatalogueParser
from hmtk.parsers.catalogue.csv_catalogue_parser import CsvCatalogueWriter
from hmtk.parsers.source_model.nrml04_parser import nrmlSourceModelParser
# Completeness
from hmtk.seismicity.completeness.comp_stepp_1971 import Stepp1971
# Plotting tools
from hmtk.plotting.mapping import HMTKBaseMap
from hmtk.plotting.seismicity.completeness import plot_stepp_1972
from hmtk.plotting.seismicity.catalogue_plots import plot_magnitude_time_scatter
from hmtk.plotting.seismicity.catalogue_plots import plot_depth_histogram
from hmtk.plotting.seismicity.catalogue_plots import plot_magnitude_time_density
from hmtk.plotting.seismicity.max_magnitude.cumulative_moment import plot_cumulative_moment
from hmtk.plotting.seismicity.catalogue_plots import (plot_observed_recurrence,
get_completeness_adjusted_table,
_get_catalogue_bin_limits)
# Seismicity tools: Events and declustering methods
from hmtk.seismicity.selector import CatalogueSelector
from hmtk.seismicity.declusterer.dec_afteran import Afteran
from hmtk.seismicity.declusterer.dec_gardner_knopoff import GardnerKnopoffType1
from hmtk.seismicity.declusterer.distance_time_windows import (GardnerKnopoffWindow,
GruenthalWindow,
UhrhammerWindow)
# Seismicity tools: Recurrence methods
from hmtk.seismicity.occurrence.aki_maximum_likelihood import AkiMaxLikelihood
from hmtk.seismicity.occurrence.b_maximum_likelihood import BMaxLikelihood
from hmtk.seismicity.occurrence.kijko_smit import KijkoSmit
from hmtk.seismicity.occurrence.weichert import Weichert
# Seismicity tools: Recurrence methods
from hmtk.seismicity.max_magnitude.kijko_sellevol_fixed_b import KijkoSellevolFixedb
from hmtk.seismicity.max_magnitude.kijko_sellevol_bayes import KijkoSellevolBayes
from hmtk.seismicity.max_magnitude.kijko_nonparametric_gaussian import KijkoNonParametricGaussian
from hmtk.seismicity.max_magnitude.cumulative_moment_release import CumulativeMoment
# Seismicity tools: Smoothed seismicity
from hmtk.seismicity.smoothing.smoothed_seismicity import SmoothedSeismicity
from hmtk.seismicity.smoothing.kernels.isotropic_gaussian import IsotropicGaussian
print 'Import: ok'
Seismic catalogue format
The hmtk is able to read a catalogue in .csv format (Windows comma separated values). There are a minimum number of fields that are required to import and use the catalogue information into hmtk. The fields require:
eventID* Agency year* month* day* hour* minute* second* longitude* latitude* depth* magnitude*
The fields marked with a star are mandatory.
In [ ]:
#Importing catalogue
catalogue_filename = 'input_data/catalogues/catalogue_example.csv'
parser = CsvCatalogueParser(catalogue_filename) # From .csv to hmtk
# Read and process the catalogue content in a variable called "catalogue"
catalogue = parser.read_file()
print '\n'
print 'The minimum magnitude in the catalogue is Mw =', np.min(catalogue.data['magnitude'])
print 'The maximum magnitude in the catalogue is Mw =', np.max(catalogue.data['magnitude'])
print 'The total number of entries is =', len(catalogue.data['magnitude'])
print '\n'
print 'The catalogue contains:'
print catalogue.data.keys()
print '\n'
In [ ]:
'''
Filtering events based on a minimum magnitude
'''
# Removing from the catalogue events with magnitude smaller than 4.5
catalogue.catalogue_mt_filter([[3000., 4.5]])
print '\n'
print 'The minimum magnitude in the catalogue is Mw =', np.min(catalogue.data['magnitude'])
print 'The maximum magnitude in the catalogue is Mw =', np.max(catalogue.data['magnitude'])
print '\n'
print 'Period covered by the catalogue: %s a %s' % (catalogue.start_year, catalogue.end_year)
In [ ]:
# Extact the limits of the catalogue
min_lon, max_lon, min_lat, max_lat = catalogue.get_bounding_box()
print 'Min Lon = %f, Max Lon = %f, Min Lat = %f, Max Lat = %f' % (min_lon, max_lon, min_lat, max_lat)
In [ ]:
catalogue.sort_catalogue_chronologically()
print 'The catalogue is now arranged cronologically'
Visualization of the catalogue per magnitude vs. time. In case the catalogue contain error magnitude, this information can be used writing "plot_error=True"
plot_magnitude_time_scatter(catalogue, plot_error=True, fmt_string='.')
In [ ]:
plot_magnitude_time_scatter(catalogue, fmt_string='Dr')
Plot showing the density of events magnitude vs time.
It is needed to define the magnitude and time bin width
In [ ]:
mag_bin = 0.1
time_bin = 1.0
plot_magnitude_time_density(catalogue, mag_bin, time_bin)
In [ ]:
# Shows depth histogram every 10 km
plot_depth_histogram(catalogue, 10.)
Visualization of the events in the catalogue by magnitude. Read and visualize seismic sources in nrml format.
The hmtk supports the next type of source modelling:
In [ ]:
# Map configuration
map_config = {'min_lon': -82., 'max_lon': -65., 'min_lat': -5, 'max_lat': 12., 'resolution':'l'}
# Creating a basemap
basemap1 = HMTKBaseMap(map_config, 'Catalogue and seismic sources')
# Adding the catalogue to the basemap
basemap1.add_catalogue(catalogue, overlay=True)
# File containing the seismic sources
source_model_file = 'input_data/source_models/source_model_example.xml'
# Reading the models
parser = nrmlSourceModelParser(source_model_file)
# Parse the seismic sources and save them into a variable called "source_model"
source_model = parser.read_file()
# Adding the seismic sources
basemap1.add_source_model(source_model, area_border='r-', border_width=1.5, alpha=0.5)
Tools used to select events in the catalogue based on:
I. Events inside an area source defined by a polygon
II. Events to a given distance from a fault. The distance metrics allowed by the hmtk are:
- Joyner-Boore distance
- Rupture distance
The hmtk allows to select events based on other definitions or paremeters, such as time period, depth, etc.
In [ ]:
# Creates an instance of the slector method called "selector" and copy the catalogue using "create_copy=True"
selector1 = CatalogueSelector(catalogue, create_copy=True)
In [ ]:
for source in source_model.sources:
if isinstance(source, area_source.mtkAreaSource):
source.select_catalogue(selector1)
print 'Area source %s, name %s, # of events %8.0f' %(source.id, source.name, source.catalogue.get_number_events())
#subcatalogue_area = selector.
else:
continue
In [ ]:
map_config = {'min_lon': -82., 'max_lon': -65., 'min_lat': -5, 'max_lat': 12., 'resolution':'l'}
basemap2 = HMTKBaseMap(map_config, 'Events inside the area source')
basemap2.add_catalogue(source_model.sources[0].catalogue, overlay=True)
basemap2.add_source_model(source_model, area_border='k-')
In [ ]:
# Selecting events from 50 km from the simple fault
selector2 = CatalogueSelector(catalogue, create_copy=True)
for source in source_model.sources:
if isinstance(source, simple_fault_source.mtkSimpleFaultSource):
source.select_catalogue(selector2, 50)
print 'Source number %s, name %s, # of events %8.0f' %(source.id, source.name, source.catalogue.get_number_events())
elif isinstance(source, complex_fault_source.mtkComplexFaultSource):
source.select_catalogue(selector2, 50)
print 'Source number %s, name %s, # of events %8.0f' %(source.id, source.name, source.catalogue.get_number_events())
else:
continue
In [ ]:
map_config = {'min_lon': -82., 'max_lon': -65., 'min_lat': -5, 'max_lat': 12., 'resolution':'l'}
basemap3 = HMTKBaseMap(map_config, 'Events 50 km from the simple fault')
basemap3.add_catalogue(source_model.sources[1].catalogue, overlay=True)
basemap3.add_source_model(source_model, area_border='k-')
In [ ]:
map_config = {'min_lon': -82., 'max_lon': -65., 'min_lat': -5, 'max_lat': 12., 'resolution':'l'}
basemap4 = HMTKBaseMap(map_config, 'Events 50 km from the complex fault')
basemap4.add_catalogue(source_model.sources[2].catalogue, overlay=True)
basemap4.add_source_model(source_model, area_border='k-', alpha=0.2)
The declusterin methods included in the hmtk are:
- Gardner & Knopoff (1974)
- Afteran (Musson, 1999b)
The time-distance windows available are:
- UhrhammerWindow
- GardnerKnopoffWindow
- GruenthalWindow
i.e. from hmtk.seismicity.declusterer.distance_time_windows import UhrhammerWindow
Output:
Cluster index: Flag indicating the cluster or group to which every event is assigned.
Cluster flag: Flag indicating the type of events:
- 0 if is considered a main shock
- 1 if is considered foreshock or aftershock
In [ ]:
# Time-distance window configuration
declust_config_gk = {'time_distance_window': UhrhammerWindow(), 'fs_time_prop': 1.0}
# Instanciate the declustering methode
declustering_gk = GardnerKnopoffType1()
# Declustering
cluster_index_gk, cluster_flag_gk = declustering_gk.decluster(catalogue, declust_config_gk)
# Adding the information to the catalogue: The cluster index and cluter flag
catalogue.data['Cluster_Index_gk'] = cluster_index_gk
catalogue.data['Cluster_Flag_gk'] = cluster_flag_gk
In [ ]:
# Time-distance window configuration
declust_config_af = {'time_distance_window': GardnerKnopoffWindow (), 'time_window': 100.0}
# Instanciate the declustering methode
declustering_af = Afteran()
# Declustering
cluster_index_af, cluster_flag_af = declustering_af.decluster(catalogue, declust_config_af)
# Adding the information to the catalogue: The cluster index and cluter flag
catalogue.data['Cluster_Index_af'] = cluster_index_af
catalogue.data['Cluster_Flag_af'] = cluster_flag_af
In [ ]:
# Print info
data = np.column_stack([catalogue.get_decimal_time(), catalogue.data['magnitude'], catalogue.data['longitude'], catalogue.data['latitude'], cluster_index_gk, cluster_flag_gk])
print ' Time Magnitude Long. Lat. Cluster No. Index'
for row in data:
print '%14.8f %6.2f %8.3f %8.3f %6.0f %6.0f' %(row[0], row[1], row[2], row[3], row[4], row[5])
Removing foreshocks and aftershocks from the catalogue. This example show how to remove the non Poissonian events based on the previous outputs.
Methode Choose
- Gardner and Knopoff choose cluster_flag_gk
- Afteran cluster_flag_af
In [ ]:
# Copying the catalogue and saving it under a new name "catalogue_dec"(declustered catalogue)
catalogue_dec = deepcopy(catalogue)
# Logical indexing: Chossing the outputs for the main events: Cluster_flag = 0
mainshock_flag = cluster_flag_gk == 0
# Filtering the foreshocks and aftershocks in the copy of the catalogue
catalogue_dec.purge_catalogue(mainshock_flag)
# Printing the number of events considered main shocks
print 'Declustering: ok'
print 'Number of main events =', catalogue_dec.get_number_events()
print '\n'
In [ ]:
# Selecting path and name for the output file
output_cat_dec = 'output_data/cat_dec.csv'
# Call the method and save the output file under the name "cat_csv"
cat_csv = CsvCatalogueWriter(output_cat_dec)
# Write the purged catalogue
cat_csv.write_file(catalogue)
print 'Declustered catalogue: ok\n'
Estimates the completeness based on the Stepp methode (1971). The outputs are the completeness table and a plot.
The configuration file contains three variables:
1) 'magnitude_bin': The size of the magnitude bin (float)
2) 'time_bin': The size of the time-step to increment the time analysis (float)
3) 'increment_lock': A boolean ("True"/"False") variable to determine whether the trend in the completeness magnitude
should be fixed ("True") in order to ensure the trend of lower completeness magnitudes for the later part of the
catalogue, or to allow the completeness to vary without constraint ("False").
In [ ]:
# Set up the configuration parameterss
comp_config = {'magnitude_bin': 1.0, 'time_bin': 5.0, 'increment_lock': True}
# Calling the method
completeness_algorithm = Stepp1971()
# Use the catalogue and completeness configuration
completeness_table = completeness_algorithm.completeness(catalogue, comp_config)
print 'Completeness: ok'
# Print the completeness table
print '\n'
print 'Completeness table using Stepp method (1971)'
print completeness_table
print '\n'
# Setting configuration for the completeness plot
completeness_parameters = completeness_algorithm
if os.path.exists('output_data/Stepp_prueba.png'):
os.remove('output_data/Stepp_prueba.png')
plot_stepp_1972.create_stepp_plot(completeness_parameters, 'output_data/Stepp_prueba.png')
In [ ]:
'''
The hmtk allows to use a completeness table proposed by the modeller.
'''
# Table format
completeness_table_a = np.array([[1995, 5.0],
[1983, 6.0],
[1970, 7.0],
[1902, 8.0]])
# Print the table
print '\n'
print 'Completeness table'
print completeness_table_a
print '\n'
Estimation of the Gutenberg-Richter a-value and b-value
I. The inputs needed are a catalogue and a completeness table.
II. In the case where no completeness table is specified, the method assumes
that the catalogue is complete above the minimum magnitude included in
the given catalogue.
III. The methods included in the hmtk are:
- Maximum Likelihood b-value
- Kijko and Smit (2012)
- Weichert (1980)
The average can be estimated using:
- Weighted
- Harmonic
In [ ]:
print 'Completeness properties'
plot_observed_recurrence(catalogue, completeness_table_a, 0.1)
In [ ]:
'''
This method is simply an adaptation of the Aki (1965) approach
in which b-value is determined from the weighted mean of the b-value
for each completeness period.
'''
# Set up the configuration parameters
mle_config = {'magnitude_interval': 1.0, 'Average Type': 'Weighted', 'reference_magnitude': 4.5}
# Calling the method under the name "recurrence"
recurrence = BMaxLikelihood()
# b-value, a-value and unceratinties estimation
bval, sigmab, aval, sigmaa = recurrence.calculate(catalogue_dec, mle_config, completeness = completeness_table)
print '\n'
print 'b-Maximum Likelihood: ok'
print 'bval =', bval, 'sigmab =', sigmab, 'aval =', aval, 'sigmaa =', sigmaa
print '\n'
In [ ]:
'''
Implements the maximum likelihood estimator of Kijko & Smit (2012)
'''
# Set up the configuration parameters
kijko_smit_config = {'magnitude_interval': 1.0}
# Call the method under the name "recurrence"
recurrence = KijkoSmit()
# b-value, a-value and unceratinties estimation
bval, sigmab, aval, sigmaa = recurrence.calculate(catalogue_dec, kijko_smit_config, completeness = completeness_table_a)
print '\n'
print 'Kijko y Smit: ok'
print 'bval =', bval, 'sigmab =', sigmab, 'aval =', aval, 'sigmaa =', sigmaa
print '\n'
In [ ]:
'''
Implements the maximum likelihood estimator of Weichert (1980)
'''
# Set up the configuration parameters
weichert_config = {'magnitude_interval': 1.0, 'bvalue': 1., 'itstab': 1E-5, 'maxiter': 1000}
recurrence = Weichert()
bval, sigmab, rate, sigmarate = recurrence.calculate(catalogue_dec, weichert_config, completeness = completeness_table_a)
print '\n'
print 'Weichert: ok'
print 'bval =', bval, 'sigmab =', sigmab, 'rate =', rate, 'sigmarate =', sigmarate
print '\n'
The hmtk includes methods for the estimation of maximum magnitude, based on statistical analysis of the catalogue.
The methods allocated in the toolkit are:
- Kijko Case I. Fixed b-value
- Kijko Case II. Uncertain b-value
- Kijko Case III. Non-Parametric Gaussian (N-P-G)
- Cumulative Moment Release
In [ ]:
'''
Maximum likelihood estimator of maximum magnitude assuming no uncertainty
in b-value; see Kijko (2004) for more details.
'''
# Set up the configuration parameters
mmax_config = {'input_mmax': 8.0, 'input_mmax_uncertainty': 0.2, 'b-value': 0.8, 'input_mmin': 5.0, 'tolerance': 1.0E-5, 'maximum_iterations': 49}
# Call the medthod under the name "mmax_estimator"
mmax_estimator = KijkoSellevolFixedb()
print 'mmin, mmax, neq, beta'
# Estimation of mmax and uncertainty
mmax, mmax_uncertainty = mmax_estimator.get_mmax(catalogue_dec, mmax_config)
print '\n'
print 'Kijko Case I. Fixed b-value: ok'
print 'mmax =', mmax, 'mmax_uncertainty =', mmax_uncertainty
print '\n'
In [ ]:
'''
Maximum likelihood estimator of maximum magnitude assuming uncertain b-value.
'''
# Set up the configuration parameters
mmax_config = {'input_mmax': 8.9, 'input_mmax_uncertainty': 0.2, 'b-value': 0.8, 'sigma-b': 0.05, 'input_mmin': 5.0, 'tolerance': 1.0E-5, 'maximum_iterations': 1000}
# Call the medthod under the name "mmax_estimator"
mmax_estimator = KijkoSellevolBayes()
# # Estimation of mmax and uncertainty
mmax, mmax_uncertainty = mmax_estimator.get_mmax(catalogue, mmax_config)
print '\n'
print 'Kijko Case II. Uncertain b-value: ok'
print 'mmax =', mmax, 'mmax_uncertainty =', mmax_uncertainty
print '\n'
In [ ]:
'''
Maximum likelihood estimator of maximum magnitude assuming no specified magnitude frequency distribution
'''
# Set up the configuration parameters
mmax_config3 = {'input_mmax': 8.0, 'input_mmax_uncertainty': 0.2, 'number_samples': 51, 'number_earthquakes': 50, 'tolerance': 1.0E-3, 'maximum_iterations': 1000}
# Call the method under the name "mmax_estimator"
mmax_estimator = KijkoNonParametricGaussian()
# # Estimation of mmax and uncertainty
mmax, mmax_uncertainty = mmax_estimator.get_mmax(catalogue_dec, mmax_config3)
print '\n'
print 'mmax =', mmax, 'mmax_uncertainty =', mmax_uncertainty
print '\n'
In [ ]:
'''
Estimator of Maximum Magnitude based on the "Cumulative Moment" method,
an adaptation of the cumulative strain energy estimator of Mmax
(Makropoulos & Burton, 1983), with Mmax uncertainty estimated via Monte Carlo
sampling of the observed magnitude uncertainties.
'''
# Set up the configuration parameters
cm_config = {'number_bootstraps': 500}
# # Call the method under the name "cmax_estimator"
cmmax_estimator = CumulativeMoment()
# Estimation of cmmax and uncertainty
cmmax, cmmax_uncertainty = cmmax_estimator.get_mmax(catalogue_dec, cm_config)
print '\n'
print 'cmax =', cmmax, 'cmmax_uncertainty =', cmmax_uncertainty
print '\n'
In [ ]:
plot_cumulative_moment(catalogue_dec.data['year'], catalogue_dec.data['magnitude'])
Smoothed Seismicty from Frankel (1995)
The method requires as input:
- Grid definition:
Long min, Long max, Lon step, Lat min, Lat max, Lat step, Depth min, Depth max, Depth step
- b-value
- Completeness table
Output:
File in .csv format with a-values in a grid
In [ ]:
print 'Lat min, Lat max, Long min, Lon max'
print min_lat, max_lat, min_lon, max_lon
table = np.array([[1995, 5.0],
[1983, 6.0],
[1970, 7.0],
[1902, 8.0]])
bval = 0.8
# Defining the region for the analysis
grid_limits = [min_lon, max_lon, 0.5, min_lat, max_lat, 0.5, 0., 100., 100.]
# Call the method under the name "smooth_seis"
smooth_seis = SmoothedSeismicity(grid_limits, use_3d=False, bvalue= bval)
# Set up the kernel parameters
config = {'Length_Limit': 1., 'BandWidth': 30., 'increment': False}
# Printing parameters
print 'b-value is:', bval
print '\n'
print 'Completeness table', '\n', table
# Smoothed seismicity estimation
print '\n'
output_data = smooth_seis.run_analysis(catalogue_dec, config, table, smoothing_kernel=IsotropicGaussian())
# Exporting results
smooth_seis.write_to_csv('output_data/smoothed_seismicity.csv')
print 'Smoothed Seismicity: ok'
print '\n'
In [ ]:
from matplotlib.colors import LogNorm
# Map configuration
map_config = {'min_lon': -82., 'max_lon': -65., 'min_lat': -5, 'max_lat': 12., 'resolution':'l'}
# Creating a basemap
basemap2 = HMTKBaseMap(map_config, 'Smoothed Seismicity')
basemap2.add_colour_scaled_points(output_data[:,0], output_data[:,1], output_data[:,4], norm=LogNorm(1E-4, 1))
In [ ]: