Source code for openquake.hmtk.seismicity.gcmt_utils
# -*- 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."""Set of moment tensor utility functions"""importnumpyasnpfrommathimportfabs,log10,sqrt,acos,atan2,pi
[docs]deftensor_components_to_use(mrr,mtt,mpp,mrt,mrp,mtp):""" Converts components to Up, South, East definition:: USE = [[mrr, mrt, mrp], [mtt, mtt, mtp], [mrp, mtp, mpp]] """returnnp.array([[mrr,mrt,mrp],[mrt,mtt,mtp],[mrp,mtp,mpp]])
[docs]deftensor_components_to_ned(mrr,mtt,mpp,mrt,mrp,mtp):""" Converts components to North, East, Down definition:: NED = [[mtt, -mtp, mrt], [-mtp, mpp, -mrp], [mrt, -mtp, mrr]] """returnnp.array([[mtt,-mtp,mrt],[-mtp,mpp,-mrp],[mrt,-mtp,mrr]])
[docs]defget_azimuth_plunge(vect,degrees=True):""" For a given vector in USE format, retrieve the azimuth and plunge """ifvect[0]>0:vect=-1.0*np.copy(vect)vect_hor=sqrt(vect[1]**2.0+vect[2]**2.0)plunge=atan2(-vect[0],vect_hor)azimuth=atan2(vect[2],-vect[1])ifdegrees:icr=180.0/pireturnicr*azimuth%360.0,icr*plungeelse:returnazimuth%(2.0*pi),plunge
[docs]defuse_to_ned(tensor):""" Converts a tensor in USE coordinate sytem to NED """returnnp.array(ROT_NED_USE.T@tensor@ROT_NED_USE)
[docs]defned_to_use(tensor):""" Converts a tensor in NED coordinate sytem to USE """returnnp.array(ROT_NED_USE@tensor@ROT_NED_USE.T)
[docs]deftensor_to_6component(tensor,frame="USE"):""" Returns a tensor to six component vector [Mrr, Mtt, Mpp, Mrt, Mrp, Mtp] """if"NED"inframe:tensor=ned_to_use(tensor)return[tensor[0,0],tensor[1,1],tensor[2,2],tensor[0,1],tensor[0,2],tensor[1,2],]
[docs]defnormalise_tensor(tensor):""" Normalise the tensor by dividing it by its norm, defined such that np.sqrt(X:X) """tensor_norm=np.linalg.norm(tensor)returntensor/tensor_norm,tensor_norm
[docs]defeigendecompose(tensor,normalise=False):""" Performs and eigendecomposition of the tensor and orders into descending eigenvalues """ifnormalise:tensor,tensor_norm=normalise_tensor(tensor)else:tensor_norm=1.0eigvals,eigvects=np.linalg.eigh(tensor,UPLO="U")isrt=np.argsort(eigvals)eigenvalues=eigvals[isrt]*tensor_normeigenvectors=eigvects[:,isrt]returneigenvalues,eigenvectors
[docs]defmatrix_to_euler(rotmat):"""Inverse of euler_to_matrix()."""defcvec(x,y,z):returnnp.array([[x,y,z]]).Tex=cvec(1.0,0.0,0.0)ez=cvec(0.0,0.0,1.0)exs=rotmat.T@exezs=rotmat.T@ezenodes=np.cross(ez.T,ezs.T).Tifnp.linalg.norm(enodes)<1e-10:enodes=exsenodess=rotmat*enodescos_alpha=float((ez.T@ezs))ifcos_alpha>1.0:cos_alpha=1.0ifcos_alpha<-1.0:cos_alpha=-1.0alpha=acos(cos_alpha)beta=np.mod(atan2(enodes[1,0],enodes[0,0]),pi*2.0)gamma=np.mod(-atan2(enodess[1,0],enodess[0,0]),pi*2.0)returnunique_euler(alpha,beta,gamma)
[docs]defunique_euler(alpha,beta,gamma):""" Uniquify euler angle triplet. Put euler angles into ranges compatible with (dip,strike,-rake) in seismology: alpha (dip) : [0, pi/2] beta (strike) : [0, 2*pi) gamma (-rake) : [-pi, pi) If alpha is near to zero, beta is replaced by beta+gamma and gamma is set to zero, to prevent that additional ambiguity. If alpha is near to pi/2, beta is put into the range [0,pi). """alpha=np.mod(alpha,2.0*pi)if0.5*pi<alphaandalpha<=pi:alpha=pi-alphabeta=beta+pigamma=2.0*pi-gammaelifpi<alphaandalpha<=1.5*pi:alpha=alpha-pigamma=pi-gammaelif1.5*pi<alphaandalpha<=2.0*pi:alpha=2.0*pi-alphabeta=beta+pigamma=pi+gammaalpha=np.mod(alpha,2.0*pi)beta=np.mod(beta,2.0*pi)gamma=np.mod(gamma+pi,2.0*pi)-pi# If dip is exactly 90 degrees, one is still# free to choose between looking at the plane from either side.# Choose to look at such that beta is in the range [0,180)# This should prevent some problems, when dip is close to 90 degrees:iffabs(alpha-0.5*pi)<1e-10:alpha=0.5*piiffabs(beta-pi)<1e-10:beta=piiffabs(beta-2.0*pi)<1e-10:beta=0.0iffabs(beta)<1e-10:beta=0.0ifalpha==0.5*piandbeta>=pi:gamma=-gammabeta=np.mod(beta-pi,2.0*pi)gamma=np.mod(gamma+pi,2.0*pi)-piassert0.0<=beta<piassert-pi<=gamma<piifalpha<1e-7:beta=np.mod(beta+gamma,2.0*pi)gamma=0.0return(alpha,beta,gamma)
[docs]defmoment_magnitude_scalar(moment):""" Uses Hanks & Kanamori formula for calculating moment magnitude from a scalar moment (Nm) """ifisinstance(moment,np.ndarray):return(2.0/3.0)*(np.log10(moment)-9.05)else:return(2.0/3.0)*(log10(moment)-9.05)