Source code for openquake.hmtk.seismicity.gcmt_catalogue

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4

#
# LICENSE
#
# Copyright (C) 2010-2023 GEM Foundation, G. Weatherill, M. Pagani,
# D. Monelli.
#
# The Hazard Modeller's Toolkit is free software: you can redistribute
# it and/or modify it under the terms of the GNU Affero General Public
# License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>
#
# DISCLAIMER
#
# The software Hazard Modeller's Toolkit (openquake.hmtk) provided herein
# is released as a prototype implementation on behalf of
# scientists and engineers working within the GEM Foundation (Global
# Earthquake Model).
#
# It is distributed for the purpose of open collaboration and in the
# hope that it will be useful to the scientific, engineering, disaster
# risk and software design communities.
#
# The software is NOT distributed as part of GEM’s OpenQuake suite
# (https://www.globalquakemodel.org/tools-products) and must be considered as a
# separate entity. The software provided herein is designed and implemented
# by scientific staff. It is not developed to the design standards, nor
# subject to same level of critical review by professional software
# developers, as GEM’s OpenQuake software suite.
#
# Feedback and contribution to the software is welcome, and can be
# directed to the hazard scientific staff of the GEM Model Facility
# (hazard@globalquakemodel.org).
#
# The Hazard Modeller's Toolkit (openquake.hmtk) is therefore distributed
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
# for more details.
#
# The GEM Foundation, and the authors of the software, assume no
# liability for use of the software.

"""
Implements sets of classes for mapping components of the focal mechanism
"""
import csv
import datetime
import numpy as np
from math import fabs, floor, sqrt, pi
from openquake.hmtk.seismicity import gcmt_utils as utils
from openquake.hmtk.seismicity.catalogue import Catalogue


[docs]def cmp(a, b): # Python 3 replacement of Python2 cmp return (a > b) - (a < b)
[docs]def cmp_mat(a, b): """ Sorts two matrices returning a positive or zero value """ c = 0 for x, y in zip(a.flat, b.flat): c = cmp(abs(x), abs(y)) if c != 0: return c return c
[docs]class GCMTHypocentre(object): """ Simple representation of the hypocentre """ def __init__(self): """ """ self.source = None self.date = None self.time = None self.longitude = None self.latitude = None self.depth = None self.m_b = None self.m_s = None self.location = None
[docs]class GCMTCentroid(object): """ Representation of a GCMT centroid """ def __init__(self, reference_date, reference_time): """ :param reference_date: Date of hypocentre as instance of :class: datetime.datetime.date :param reference_time: Time of hypocentre as instance of :class: datetime.datetime.time """ self.centroid_type = None self.source = None self.time = reference_time self.time_error = None self.date = reference_date self.longitude = None self.longitude_error = None self.latitude = None self.latitude_error = None self.depth = None self.depth_error = None self.depth_type = None self.centroid_id = None def _get_centroid_time(self, time_diff): """ Calculates the time difference between the date-time classes """ source_time = datetime.datetime.combine(self.date, self.time) second_diff = floor(fabs(time_diff)) microsecond_diff = int(1000.0 * (time_diff - second_diff)) if time_diff < 0.0: source_time = source_time - datetime.timedelta( seconds=int(second_diff), microseconds=microsecond_diff ) else: source_time = source_time + datetime.timedelta( seconds=int(second_diff), microseconds=microsecond_diff ) self.time = source_time.time() self.date = source_time.date()
[docs]class GCMTPrincipalAxes(object): """ Class to represent the eigensystem of the tensor in terms of T-, B- and P- plunge and azimuth """ def __init__(self): """ """ self.t_axis = None self.b_axis = None self.p_axis = None
[docs] def get_moment_tensor_from_principal_axes(self): """ Retreives the moment tensor from the prinicpal axes """ raise NotImplementedError( "Moment tensor from principal axes not yet " "implemented!" )
[docs] def get_azimuthal_projection(self, height=1.0): """ Returns the azimuthal projection of the tensor according to the method of Frohlich (2001) """ raise NotImplementedError( "Get azimuthal projection not yet " "implemented!" )
[docs]class GCMTMomentTensor(object): """ Class to represent a moment tensor and its associated methods """ def __init__(self, reference_frame=None): """ """ self.tensor = None self.tensor_sigma = None self.exponent = None self.eigenvalues = None self.eigenvectors = None if reference_frame: self.ref_frame = reference_frame else: # Default to USE self.ref_frame = "USE"
[docs] def normalise_tensor(self): """ Normalise the tensor by dividing it by its norm, defined such that np.sqrt(X:X) """ self.tensor, tensor_norm = utils.normalise_tensor(self.tensor) return self.tensor / tensor_norm, tensor_norm
def _to_ned(self): """ Switches the reference frame to NED """ if self.ref_frame == "USE": # Rotate return utils.use_to_ned(self.tensor), utils.use_to_ned( self.tensor_sigma ) elif self.ref_frame == "NED": # Alreadt NED return self.tensor, self.tensor_sigma else: raise ValueError( "Reference frame %s not recognised - cannot " "transform to NED!" % self.ref_frame ) def _to_use(self): """ Returns a tensor in the USE reference frame """ if self.ref_frame == "NED": # Rotate return utils.ned_to_use(self.tensor), utils.ned_to_use( self.tensor_sigma ) elif self.ref_frame == "USE": # Already USE return self.tensor, self.tensor_sigma else: raise ValueError( "Reference frame %s not recognised - cannot " "transform to USE!" % self.ref_frame ) def _to_6component(self): """ Returns the unique 6-components of the tensor in USE format [Mrr, Mtt, Mpp, Mrt, Mrp, Mtp] """ return utils.tensor_to_6component(self.tensor, self.ref_frame)
[docs] def eigendecompose(self, normalise=False): """ Performs and eigendecomposition of the tensor and orders into descending eigenvalues """ self.eigenvalues, self.eigenvectors = utils.eigendecompose( self.tensor, normalise ) return self.eigenvalues, self.eigenvectors
[docs] def get_nodal_planes(self): """ Returns the nodal planes by eigendecomposition of the moment tensor """ # Convert reference frame to NED self.tensor, self.tensor_sigma = self._to_ned() self.ref_frame = "NED" # Eigenvalue decomposition # Tensor _, evect = utils.eigendecompose(self.tensor) # Rotation matrix _, rot_vec = utils.eigendecompose( np.array([[0.0, 0.0, -1], [0.0, 0.0, 0.0], [-1.0, 0.0, 0.0]]) ) rotation_matrix = (evect @ rot_vec.T).T if np.linalg.det(rotation_matrix) < 0.0: rotation_matrix @= -1.0 flip_dc = np.array( [[0.0, 0.0, -1.0], [0.0, -1.0, 0.0], [-1.0, 0.0, 0.0]] ) rotation_matrices = sorted( [rotation_matrix, flip_dc @ rotation_matrix], cmp=cmp_mat ) nodal_planes = GCMTNodalPlanes() dip, strike, rake = [ (180.0 / pi) * angle for angle in utils.matrix_to_euler(rotation_matrices[0]) ] # 1st Nodal Plane nodal_planes.nodal_plane_1 = { "strike": strike % 360, "dip": dip, "rake": -rake, } # 2nd Nodal Plane dip, strike, rake = [ (180.0 / pi) * angle for angle in utils.matrix_to_euler(rotation_matrices[1]) ] nodal_planes.nodal_plane_2 = { "strike": strike % 360.0, "dip": dip, "rake": -rake, } return nodal_planes
[docs] def get_principal_axes(self): """ Uses the eigendecomposition to extract the principal axes from the moment tensor - returning an instance of the GCMTPrincipalAxes class """ # Perform eigendecomposition - returns in order P, B, T _ = self.eigendecompose(normalise=True) principal_axes = GCMTPrincipalAxes() # Eigenvalues principal_axes.p_axis = {"eigenvalue": self.eigenvalues[0]} principal_axes.b_axis = {"eigenvalue": self.eigenvalues[1]} principal_axes.t_axis = {"eigenvalue": self.eigenvalues[2]} # Eigen vectors # 1) P axis azim, plun = utils.get_azimuth_plunge(self.eigenvectors[:, 0], True) principal_axes.p_axis["azimuth"] = azim principal_axes.p_axis["plunge"] = plun # 2) B axis azim, plun = utils.get_azimuth_plunge(self.eigenvectors[:, 1], True) principal_axes.b_axis["azimuth"] = azim principal_axes.b_axis["plunge"] = plun # 3) T axis azim, plun = utils.get_azimuth_plunge(self.eigenvectors[:, 2], True) principal_axes.t_axis["azimuth"] = azim principal_axes.t_axis["plunge"] = plun return principal_axes
[docs]class GCMTEvent(object): """ Class to represent full GCMT moment tensor in ndk format """ def __init__(self): """ """ self.identifier = None self.hypocentre = None self.centroid = None self.magnitude = None self.moment = None self.metadata = {} self.moment_tensor = None self.nodal_planes = None self.principal_axes = None self.f_clvd = None self.e_rel = None
[docs] def get_f_clvd(self): """ Returns the statistic f_clvd: the signed ratio of the sizes of the intermediate and largest principal moments:: f_clvd = -b_axis_eigenvalue / max(|t_axis_eigenvalue|,|p_axis_eigenvalue|) """ if not self.principal_axes: # Principal axes not yet defined for moment tensor - raises error raise ValueError("Principal Axes not defined!") denominator = np.max( np.array( [ fabs(self.principal_axes.t_axis["eigenvalue"]), fabs(self.principal_axes.p_axis["eigenvalue"]), ] ) ) self.f_clvd = -self.principal_axes.b_axis["eigenvalue"] / denominator return self.f_clvd
[docs] def get_relative_error(self): """ Returns the relative error statistic (e_rel), defined by Frohlich & Davis (1999): `e_rel = sqrt((U:U) / (M:M))` where M is the moment tensor, U is the uncertainty tensor and : is the tensor dot product """ if not self.moment_tensor: raise ValueError("Moment tensor not defined!") numer = np.tensordot( self.moment_tensor.tensor_sigma, self.moment_tensor.tensor_sigma ) denom = np.tensordot( self.moment_tensor.tensor, self.moment_tensor.tensor ) self.e_rel = sqrt(numer / denom) return self.e_rel
[docs]class GCMTNodalPlanes(object): """ Class to represent the nodal plane distribution of the tensor Each nodal plane is represented as a dictionary of the form: {'strike':, 'dip':, 'rake':} """ def __init__(self): """ """ self.nodal_plane_1 = None self.nodal_plane_2 = None
[docs]class GCMTCatalogue(Catalogue): """ Class to hold a catalogue of moment tensors """ FLOAT_ATTRIBUTE_LIST = [ "second", "timeError", "longitude", "latitude", "SemiMajor90", "SemiMinor90", "ErrorStrike", "depth", "depthError", "magnitude", "sigmaMagnitude", "moment", "strike1", "rake1", "dip1", "strike2", "rake2", "dip2", "eigenvalue_b", "azimuth_b", "plunge_b", "eigenvalue_p", "azimuth_p", "plunge_p", "eigenvalue_t", "azimuth_t", "plunge_t", "f_clvd", "e_rel", ] INT_ATTRIBUTE_LIST = [ "eventID", "year", "month", "day", "hour", "minute", "flag", ] STRING_ATTRIBUTE_LIST = [ "Agency", "magnitudeType", "comment", "centroidID", ] TOTAL_ATTRIBUTE_LIST = list( (set(FLOAT_ATTRIBUTE_LIST).union(set(INT_ATTRIBUTE_LIST))).union( set(STRING_ATTRIBUTE_LIST) ) ) def __init__(self, start_year=None, end_year=None): """ Instantiate catalogue class """ super(GCMTCatalogue, self).__init__() self.gcmts = [] self.number_gcmts = None self.start_year = start_year self.end_year = end_year for attribute in self.TOTAL_ATTRIBUTE_LIST: if attribute in self.FLOAT_ATTRIBUTE_LIST: self.data[attribute] = np.array([], dtype=float) elif attribute in self.INT_ATTRIBUTE_LIST: self.data[attribute] = np.array([], dtype=int)
[docs] def get_number_tensors(self): """ Returns number of CMTs """ return len(self.gcmts)
[docs] def select_catalogue_events(self, id0): """ Orders the events in the catalogue according to an indexing vector :param np.ndarray id0: Pointer array indicating the locations of selected events """ for key in self.data.keys(): if ( isinstance(self.data[key], np.ndarray) and len(self.data[key]) > 0 ): # Dictionary element is numpy array - use logical indexing self.data[key] = self.data[key][id0] elif isinstance(self.data[key], list) and len(self.data[key]) > 0: # Dictionary element is list self.data[key] = [self.data[key][iloc] for iloc in id0] else: continue if len(self.gcmts) > 0: self.gcmts = [self.gcmts[iloc] for iloc in id0] self.number_gcmts = self.get_number_tensors()
[docs] def serialise_to_hmtk_csv(self, filename, centroid_location=True): """ Serialise the catalogue to a simple csv format, designed for comptibility with the GEM Hazard Modeller's Toolkit """ header_list = [ "eventID", "Agency", "year", "month", "day", "hour", "minute", "second", "timeError", "longitude", "latitude", "SemiMajor90", "SemiMinor90", "ErrorStrike", "depth", "depthError", "magnitude", "sigmaMagnitude", ] with open(filename, "wt", newline="") as fid: writer = csv.DictWriter(fid, fieldnames=header_list) headers = dict((header, header) for header in header_list) writer.writerow(headers) print("Writing to simple csv format ...") for iloc, tensor in enumerate(self.gcmts): # Generic Data cmt_dict = { "eventID": iloc + 100000, "Agency": "GCMT", "SemiMajor90": None, "SemiMinor90": None, "ErrorStrike": None, "magnitude": tensor.magnitude, "sigmaMagnitude": None, "depth": None, "depthError": None, } if centroid_location: # Time and location come from centroid cmt_dict["year"] = tensor.centroid.date.year cmt_dict["month"] = tensor.centroid.date.month cmt_dict["day"] = tensor.centroid.date.day cmt_dict["hour"] = tensor.centroid.time.hour cmt_dict["minute"] = tensor.centroid.time.minute cmt_dict["second"] = np.round( float(tensor.centroid.time.second) + float(tensor.centroid.time.microsecond) / 1.0e6, 2, ) cmt_dict["timeError"] = tensor.centroid.time_error cmt_dict["longitude"] = tensor.centroid.longitude cmt_dict["latitude"] = tensor.centroid.latitude cmt_dict["depth"] = tensor.centroid.depth cmt_dict["depthError"] = tensor.centroid.depth_error else: # Time and location come from hypocentre cmt_dict["year"] = tensor.hypocentre.date.year cmt_dict["month"] = tensor.hypocentre.date.month cmt_dict["day"] = tensor.hypocentre.date.day cmt_dict["hour"] = tensor.hypocentre.time.hour cmt_dict["minute"] = tensor.hypocentre.time.minute cmt_dict["second"] = np.round( float(tensor.hypocentre.time.second) + float(tensor.hypocentre.time.microsecond) / 1000000.0, 2, ) cmt_dict["timeError"] = None cmt_dict["longitude"] = tensor.hypocentre.longitude cmt_dict["latitude"] = tensor.hypocentre.latitude cmt_dict["depth"] = tensor.hypocentre.depth cmt_dict["depthError"] = None writer.writerow(cmt_dict) print("done!")
[docs] def gcmt_to_simple_array(self, centroid_location=True): """ Converts the GCMT catalogue to a simple array of [ID, year, month, day, hour, minute, second, long., lat., depth, Mw, strike1, dip1, rake1, strike2, dip2, rake2, b-plunge, b-azimuth, b-eigenvalue, p-plunge, p-azimuth, p-eigenvalue, t-plunge, t-azimuth, t-eigenvalue, moment, f_clvd, erel] """ catalogue = np.zeros([self.get_number_tensors(), 29], dtype=float) for iloc, tensor in enumerate(self.gcmts): catalogue[iloc, 0] = iloc if centroid_location: catalogue[iloc, 1] = float(tensor.centroid.date.year) catalogue[iloc, 2] = float(tensor.centroid.date.month) catalogue[iloc, 3] = float(tensor.centroid.date.day) catalogue[iloc, 4] = float(tensor.centroid.time.hour) catalogue[iloc, 5] = float(tensor.centroid.time.minute) catalogue[iloc, 6] = np.round( float(tensor.centroid.time.second) + float(tensor.centroid.time.microsecond) / 1000000.0, 2, ) catalogue[iloc, 7] = tensor.centroid.longitude catalogue[iloc, 8] = tensor.centroid.latitude catalogue[iloc, 9] = tensor.centroid.depth else: catalogue[iloc, 1] = float(tensor.hypocentre.date.year) catalogue[iloc, 2] = float(tensor.hypocentre.date.month) catalogue[iloc, 3] = float(tensor.hypocentre.date.day) catalogue[iloc, 4] = float(tensor.hypocentre.time.hour) catalogue[iloc, 5] = float(tensor.hypocentre.time.minute) catalogue[iloc, 6] = np.round( float(tensor.centroid.time.second) + float(tensor.centroid.time.microsecond) / 1000000.0, 2, ) catalogue[iloc, 7] = tensor.hypocentre.longitude catalogue[iloc, 8] = tensor.hypocentre.latitude catalogue[iloc, 9] = tensor.hypocentre.depth catalogue[iloc, 10] = tensor.magnitude catalogue[iloc, 11] = tensor.moment catalogue[iloc, 12] = tensor.f_clvd catalogue[iloc, 13] = tensor.e_rel # Nodal planes catalogue[iloc, 14] = tensor.nodal_planes.nodal_plane_1["strike"] catalogue[iloc, 15] = tensor.nodal_planes.nodal_plane_1["dip"] catalogue[iloc, 16] = tensor.nodal_planes.nodal_plane_1["rake"] catalogue[iloc, 17] = tensor.nodal_planes.nodal_plane_2["strike"] catalogue[iloc, 18] = tensor.nodal_planes.nodal_plane_2["dip"] catalogue[iloc, 19] = tensor.nodal_planes.nodal_plane_2["rake"] # Principal axes catalogue[iloc, 20] = tensor.principal_axes.b_axis["eigenvalue"] catalogue[iloc, 21] = tensor.principal_axes.b_axis["azimuth"] catalogue[iloc, 22] = tensor.principal_axes.b_axis["plunge"] catalogue[iloc, 23] = tensor.principal_axes.p_axis["eigenvalue"] catalogue[iloc, 24] = tensor.principal_axes.p_axis["azimuth"] catalogue[iloc, 25] = tensor.principal_axes.p_axis["plunge"] catalogue[iloc, 26] = tensor.principal_axes.t_axis["eigenvalue"] catalogue[iloc, 27] = tensor.principal_axes.t_axis["azimuth"] catalogue[iloc, 28] = tensor.principal_axes.t_axis["plunge"] return catalogue