Source code for lasif.rotations

#!/usr/bin/env python# -*- coding: utf-8 -*-"""A collection of functions to rotate vectors, seismograms and moment tensors ona spherical body, e.g. the Earth... note:: **On the used coordinate system** Latitude and longitude are natural geographical coordinates used on Earth. The coordinate system is right handed with the origin at the center of the Earth. The z-axis points directly at the North Pole and the x-axis points at (latitude 0.0/longitude 0.0), e.g. the Greenwich meridian at the equator. The y-axis therefore points at (latitude 0.0/longitue 90.0), e.g. somewhere close to Sumatra. :math:`\\theta` (theta) is the colatitude, e.g. 90.0 - latitude and is the angle from the z-axis. :math:`\phi` (phi) is the longitude and the angle from the x-axis towards the y-axis, a.k.a the azimuth angle. These are also the generally used spherical coordinates. All rotation axes have to be given as [x, y, z] in the just described coordinate system and all rotation angles have to given as degree. A positive rotation will rotate clockwise when looking in the direction of the rotation axis. For convenience reasons, most function in this module work with coordinates given in latitude and longitude.:copyright: Lion Krischer (krischer@geophysik.uni-muenchen.de), 2012-2013:license: GNU General Public License, Version 3 (http://www.gnu.org/copyleft/gpl.html)"""importmathimportnumpyasnpimportsyseps=sys.float_info.epsilondef_get_vector(*args):""" Helper function to make sure vectors always have the same format and dtype. Creates a three component column vector from either a list or three single numbers. If it already is a correct vector, do nothing. >>> vec = _get_vector(1, 2, 3) >>> vec array([ 1., 2., 3.]) >>> print vec.dtype float64 >>> vec = _get_vector([1, 2, 3]) >>> vec array([ 1., 2., 3.]) >>> print vec.dtype float64 >>> vec = _get_vector(np.array([1, 2, 3], dtype="int32")) >>> vec array([ 1., 2., 3.]) >>> print vec.dtype float64 """iflen(args)==1andisinstance(args[0],np.ndarray):returnnp.require(args[0],dtype="float64")eliflen(args)==1andlen(args[0])==3:returnnp.array(args[0],dtype="float64")eliflen(args)==3:returnnp.array(args,dtype="float64")else:raiseNotImplementedError

[docs]defrotate_vector(vector,rotation_axis,angle):""" Takes a vector and rotates it around a rotation axis with a given angle. :param vector: The vector to be rotated given as [x, y, z]. :param rotation_axis: The axis to be rotating around given as [x, y, z]. :param angle: The rotation angle in degree. """rotation_matrix=_get_rotation_matrix(rotation_axis,angle)# Build a column vector.vector=_get_vector(vector)# Rotate the vector.rotated_vector=np.array(rotation_matrix.dot(vector))# Make sure is also works for arrays of vectors.ifrotated_vector.shape[0]>1:returnrotated_vectorelse:returnrotated_vector.ravel()

def_get_rotation_matrix(axis,angle):""" Returns the rotation matrix for the specified axis and angle. """axis=map(float,axis)/np.linalg.norm(axis)angle=np.deg2rad(angle)# Use c1, c2, and c3 as shortcuts for the rotation axis.c1=axis[0]c2=axis[1]c3=axis[2]# Build the rotation matrix.rotation_matrix=np.cos(angle)* \
np.matrix(((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0)))+ \
(1-np.cos(angle))*np.matrix(((c1*c1,c1*c2,c1*c3),(c2*c1,c2*c2,c2*c3),(c3*c1,c3*c2,c3*c3)))+ \
np.sin(angle)*np.matrix(((0,-c3,c2),(c3,0,-c1),(-c2,c1,0)))returnrotation_matrix

[docs]defget_spherical_unit_vectors(lat,lon):""" Returns the spherical unit vectors e_theta, e_phi and e_r for a point on the sphere determined by latitude and longitude which are defined as they are for earth. :param lat: Latitude in degree :param lon: Longitude in degree """colat=lat2colat(lat)# Convert to radian.colat,lon=map(np.deg2rad,[colat,lon])e_theta=_get_vector(np.cos(lon)*np.cos(colat),np.sin(lon)*np.cos(colat),-np.sin(colat))e_phi=_get_vector(-np.sin(lon),np.cos(lon),0.0)e_r=_get_vector(np.cos(lon)*np.sin(colat),np.sin(lon)*np.sin(colat),np.cos(colat))returne_theta,e_phi,e_r

[docs]defrotate_lat_lon(lat,lon,rotation_axis,angle):""" Takes a point specified by latitude and longitude and return a new pair of latitude longitude assuming the earth has been rotated around rotation_axis by angle. :param lat: Latitude of original point :param lon: Longitude of original point :param rotation_axis: Rotation axis specified as [x, y, z]. :param angle: Rotation angle in degree. """# Convert to xyz. Do the calculation on the unit sphere as the radius does# not matter.xyz=lat_lon_radius_to_xyz(lat,lon,1.0)# Rotate xyz.new_xyz=rotate_vector(xyz,rotation_axis,angle)new_lat,new_lon,_=xyz_to_lat_lon_radius(new_xyz)returnnew_lat,new_lon

[docs]deflat_lon_radius_to_xyz(lat,lon,r):""" Converts latitude, longitude and radius to x, y, and z. :param lat: The latitude. :param lon: The longitude. :param r: The radius. """colat=lat2colat(lat)# To radiancolat,lon=map(np.deg2rad,[colat,lon])# Do the transformationx=r*np.sin(colat)*np.cos(lon)y=r*np.sin(colat)*np.sin(lon)z=r*np.cos(colat)return_get_vector(x,y,z)

def_get_axis_and_angle_from_rotation_matrix(M):""" Extracts the axis and angle from a rotation matrix. >>> M = _get_rotation_matrix([0.26726124, 0.53452248, 0.80178373], 40) >>> _get_axis_and_angle_from_rotation_matrix(M) \ # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS (array([ 0.26726124, 0.53452248, 0.80178373]), 40.00...) """x=M[2,1]-M[1,2]y=M[0,2]-M[2,0]z=M[1,0]-M[0,1]r=np.sqrt(x**2+y**2+z**2)t=M[0,0]+M[1,1]+M[2,2]theta=np.arctan2(r,t-1)axis=[x,y,z]axis/=np.linalg.norm(axis)returnaxis,np.rad2deg(theta)def_get_rotation_and_base_transfer_matrix(lat,lon,rotation_axis,angle):""" Generates a matrix that rotates a vector/tensor located at lat/lon and given in spherical coordinates for 'angle' degrees around 'rotation_axis'. Furthermore it performs a base change from the spherical unit vectors at lat/lon to the spherical unit vectors at the rotated lat/lon. This should never change the radial component. :param lat: Latitude of the recording point. :param lon: Longitude of the recording point. :param rotation_axis: Rotation axis given as [x, y, z]. :param angle: Rotation angle in degree. >>> mat = _get_rotation_and_base_transfer_matrix(12, 34, [-45, 34, 45], \ -66) >>> mat[2, 2] - 1.0 <= 1E-7 True >>> mat[2, 0] <= 1E-7 True >>> mat[2, 1] <= 1E-7 True >>> mat[0, 2] <= 1E-7 True >>> mat[1, 2] <= 1E-7 True """# Rotate latitude and longitude to obtain the new coordinates after the# rotation.lat_new,lon_new=rotate_lat_lon(lat,lon,rotation_axis,angle)# Get the orthonormal basis vectors at both points. This can be interpreted# as having two sets of basis vectors in the original xyz coordinate# system.e_theta,e_phi,e_r=get_spherical_unit_vectors(lat,lon)e_theta_new,e_phi_new,e_r_new=get_spherical_unit_vectors(lat_new,lon_new)# Rotate the new unit vectors in the opposite direction to simulate a# rotation in the wanted direction.e_theta_new=rotate_vector(e_theta_new,rotation_axis,-angle)e_phi_new=rotate_vector(e_phi_new,rotation_axis,-angle)e_r_new=rotate_vector(e_r_new,rotation_axis,-angle)# Calculate the transfer matrix. This works because both sets of basis# vectors are orthonormal.transfer_matrix=np.matrix(([np.dot(e_theta_new,e_theta),np.dot(e_theta_new,e_phi),np.dot(e_theta_new,e_r)],[np.dot(e_phi_new,e_theta),np.dot(e_phi_new,e_phi),np.dot(e_phi_new,e_r)],[np.dot(e_r_new,e_theta),np.dot(e_r_new,e_phi),np.dot(e_r_new,e_r)]))returntransfer_matrix

[docs]defrotate_moment_tensor(Mrr,Mtt,Mpp,Mrt,Mrp,Mtp,lat,lon,rotation_axis,angle):""" Rotates a moment tensor, given in spherical coordinates, located at lat/lon around the rotation axis and simultaneously performs a base change from the spherical unit vectors at lat/lon to the unit vectors at the new coordinates. :param Mrr: A moment tensor component. :param Mtt: A moment tensor component. :param Mpp: A moment tensor component. :param Mrt: A moment tensor component. :param Mrp: A moment tensor component. :param Mtp: A moment tensor component. :param lat: Latitude of the recording point. :param lon: Longitude of the recording point. :param rotation_axis: Rotation axis given as [x, y, z]. :param angle: Rotation angle in degree. """transfer_matrix=_get_rotation_and_base_transfer_matrix(lat,lon,rotation_axis,angle)# Assemble the second order tensor.mt=np.matrix(([Mtt,Mtp,Mrt],[Mtp,Mpp,Mrp],[Mtp,Mrp,Mrr]))# Rotate itrotated_mt=transfer_matrix.dot(mt.dot(transfer_matrix.transpose()))# Return the six independent components in the same order as they were# given.returnrotated_mt[2,2],rotated_mt[0,0],rotated_mt[1,1], \
rotated_mt[0,2],rotated_mt[1,2],rotated_mt[0,1]

[docs]defrotate_data(north_data,east_data,vertical_data,lat,lon,rotation_axis,angle):""" Rotates three component data recorded at lat/lon a certain amount of degrees around a given rotation axis. :param north_data: The north component of the data. :param east_data: The east component of the data. :param vertical_data: The vertical component of the data. Vertical is defined to be up, e.g. radially outwards. :param lat: Latitude of the recording point. :param lon: Longitude of the recording point. :param rotation_axis: Rotation axis given as [x, y, z]. :param angle: Rotation angle in degree. """transfer_matrix=_get_rotation_and_base_transfer_matrix(lat,lon,rotation_axis,angle)# Apply the transfer matrix. Invert north data because they have# to point in the other direction to be consistent with the spherical# coordinates.new_data=np.array(transfer_matrix.dot([-1.0*north_data,east_data,vertical_data]))# Return the transferred data arrays. Again negate north data.north_data=-1.0*new_data[0]east_data=new_data[1]vertical_data=new_data[2]returnnorth_data,east_data,vertical_data

[docs]defget_border_latlng_list(min_lat,max_lat,min_lng,max_lng,number_of_points_per_side=25,rotation_axis=(0,0,1),rotation_angle_in_degree=0):""" Helper function taking a spherical section defined by latitudinal and longitudal extension, rotate it around the given axis and rotation angle and return a list of points outlining the region. Useful for plotting. :param min_lat: The minimum latitude. :param max_lat: The maximum latitude. :param min_lng: The minimum longitude. :param max_lng: The maximum longitude. :param number_of_points_per_side: The number of points per side desired. :param rotation_axis: The rotation axis. Optional. :param rotation_angle_in_degree: The rotation angle in degrees. Optional. """north_border=np.empty((number_of_points_per_side,2))south_border=np.empty((number_of_points_per_side,2))east_border=np.empty((number_of_points_per_side,2))west_border=np.empty((number_of_points_per_side,2))north_border[:,0]=max_latnorth_border[:,1]=np.linspace(min_lng,max_lng,number_of_points_per_side)east_border[:,0]=np.linspace(max_lat,min_lat,number_of_points_per_side)east_border[:,1]=max_lngsouth_border[:,0]=min_latsouth_border[:,1]=np.linspace(max_lng,min_lng,number_of_points_per_side)west_border[:,0]=np.linspace(min_lat,max_lat,number_of_points_per_side)west_border[:,1]=min_lng# Rotate everything.forborderin[north_border,south_border,east_border,west_border]:for_iinxrange(number_of_points_per_side):border[_i,0],border[_i,1]=rotate_lat_lon(border[_i,0],border[_i,1],rotation_axis,rotation_angle_in_degree)# Fix dateline wraparounds.forborderin[north_border,south_border,east_border,west_border]:lngs=border[:,1]lngs[lngs<min_lng]+=360.0# Take care to only use every corner once.borders=np.concatenate([north_border,east_border[1:],south_border[1:],west_border[1:]])borders=list(borders)borders=[tuple(_i)for_iinborders]returnborders