geo_tools API
A class for interpolating any data to any grids with Cartopy.
In addition to interpolation, this class is usefull to display geographical data. Since projection mappings need to be defined only once for a given dataset/display combination, multi-pannels figures can be made efficiently.
Most of the action happens through the class ProjInds.
The following figure illustrates the convention for storring data in numpy arrays. It is assumed that the first index (rows) represents x-y direction (longitudes):
- class domutils.geo_tools.projinds.ProjInds(data_xx=None, data_yy=None, src_lon=None, src_lat=None, dest_lon=None, dest_lat=None, extent=None, extent_crs=None, dest_crs=None, source_crs=None, fig=None, image_res=(400, 400), extend_x=True, extend_y=True, average=False, smooth_radius=None, min_hits=1, missing=-9999.0)
A class for making, keeping record of, and applying projection indices.
This class handles the projection of gridded data with lat/lon coordinates to cartopy’s own data space in geoAxes.
- Parameters:
src_lon (
Optional[Any]) – (array like) longitude of data. Has to be same dimension as data.src_lat (
Optional[Any]) – (array like) latitude of data. Has to be same dimension as data.dest_lon (
Optional[Any]) – Longitude of destination grid on which the data will be regriddeddest_lat (
Optional[Any]) – Latitude of destination grid on which the data will be regriddedextent (
Optional[Any]) – To be used together with dest_crs, to get destination grid from a cartopy geoAxes (array like) Bounds of geoAxes domain [lon_min, lon_max, lat_min, lat_max]extent_crs (
Optional[Any]) – Crs of extent. By default we assume PlateCarredest_crs (
Optional[Any]) – If provided, cartopy crs instance for the destination grid (necessary for plotting maps)source_crs (
Optional[Any]) – Cartopy crs of the data coordinates. If None, a Geodetic crs is used for using latitude/longitude coordinates.fig (
Optional[Any]) – Instance of figure currently being used. If None (the default) a new figure will be instantiated but never displayed.image_res (
Optional[tuple]) – Pixel density of the figure being generated. eg (300,200) means that the image displayed on the geoAxes will be 300 pixels wide (i.e. longitude in the case of geo_axes) by 200 pixels high (i.e. latitude in the case of geo_axes).extend_x (
Optional[Any]) – Set to True (default) when the destination grid is larger than the source grid. This option allows to mark no data points with Missing rather than extrapolating. Set to False for global grids.extend_y (
Optional[Any]) – Same as extend_x but in the Y direction.average (
Optional[Any]) – (bool) When true, all pts of a source grid falling within a destination grid pt will be averaged. Usefull to display/interpolate hi resolution data to a coarser grid. Weighted averages can be obtained by providing weights to the project_data method.smooth_radius (
Optional[Any]) – Boxcar averaging with a circular area of a given radius given in kilometer. This option allows to perform smoothing at the same time as interpolation. For each point in the destination grid, all source data points within a radius given by smooth_radius will be averaged together.min_hits (
Optional[Any]) – For use in conjunction with average or smooth_radius. This keyword specifies the minimum number of data points for an average to be considered valid.missing (
Optional[float]) – Value which will be assigned to points outside of the data domain.
Example
Simple example showing how to project and display data. Rotation of points on the globe is also demonstrated
import os import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt import cartopy import domutils.legs as legs import domutils.geo_tools as geo_tools import domutils._py_tools as py_tools # make mock data and coordinates # note that there is some regularity to the grid # but that it does not conform to any particular projection regular_lons = [ [-91. , -91 , -91 ], [-90. , -90 , -90 ], [-89. , -89 , -89 ] ] regular_lats = [ [ 44 , 45 , 46 ], [ 44 , 45 , 46 ], [ 44 , 45 , 46 ] ] data_vals = [ [ 6.5, 3.5, .5 ], [ 7.5, 4.5, 1.5 ], [ 8.5, 5.5, 2.5 ] ] missing = -9999. #pixel resolution of image that will be shown in the axes img_res = (800,600) #point density for entire figure mpl.rcParams['figure.dpi'] = 800 #projection and extent of map being displayed proj_aea = cartopy.crs.AlbersEqualArea(central_longitude=-94., central_latitude=35., standard_parallels=(30.,40.)) map_extent=[-94.8,-85.2,43,48.] #------------------------------------------------------------------- #regular lat/lons are boring, lets rotate the coordinate system about # the central data point #use cartopy transforms to get xyz coordinates proj_ll = cartopy.crs.Geodetic() geo_cent = proj_ll.as_geocentric() x, y, z = geo_tools.latlon_to_unit_sphere_xyz(np.asarray(regular_lons), np.asarray(regular_lats), combined=False) # rotate points by 45 degrees counter clockwise theta = np.pi/4 rotation_matrix = geo_tools.rotation_matrix([x[1,1], y[1,1], z[1,1]], theta) rotated_xyz = np.zeros((3,3,3)) for ii, (lat_arr, lon_arr) in enumerate(zip(regular_lats, regular_lons)): for jj, (this_lat, this_lon) in enumerate(zip(lat_arr, lat_arr)): rotated_xyz[ii,jj,:] = np.matmul(rotation_matrix,[x[ii,jj], y[ii,jj], z[ii,jj]]) #from xyz to lat/lon rotatedLons, rotatedLats = geo_tools.unit_sphere_xyz_to_latlon(rotated_xyz[:,:,0], rotated_xyz[:,:,1], rotated_xyz[:,:,2], combined=False) # done rotating #------------------------------------------------------------------- #larger characters mpl.rcParams.update({'font.size': 15}) #instantiate figure fig = plt.figure(figsize=(7.5,6.)) #instantiate object to handle geographical projection of data # onto geoAxes with this specific crs and extent proj_inds = geo_tools.ProjInds(src_lon=rotatedLons, src_lat=rotatedLats, extent=map_extent, dest_crs=proj_aea, image_res=img_res) #axes for this plot ax = fig.add_axes([.01,.1,.8,.8], projection=proj_aea) ax.set_extent(proj_inds.rotated_extent, crs=proj_aea) # Set up colormapping object color_mapping = legs.PalObj(range_arr=[0.,9.], color_arr=['brown','blue','green','orange', 'red','pink','purple','yellow','b_w'], solid='col_dark', excep_val=missing, excep_col='grey_220') #geographical projection of data into axes space proj_data = proj_inds.project_data(data_vals) #plot data & palette color_mapping.plot_data(ax=ax,data=proj_data, palette='right', pal_units='[unitless]', pal_format='{:4.0f}') #palette options #add political boundaries dum = ax.add_feature(cartopy.feature.STATES.with_scale('50m'), linewidth=0.5, edgecolor='0.2',zorder=1) #plot border and mask everything outside model domain proj_inds.plot_border(ax, mask_outside=False, linewidth=2.) # save figure generated_figure = os.path.join(generated_figure_dir, 'test_projinds_simple_example.svg') plt.savefig(generated_figure)
Example showing how ProjInds can also be used for nearest neighbor interpolation
import numpy as np import domutils.geo_tools as geo_tools # Source data on a very simple grid src_lon = [ [-90.1 , -90.1 ], [-89.1 , -89.1 ] ] src_lat = [ [ 44.1 , 45.1 ], [ 44.1 , 45.1 ] ] data = [ [ 3 , 1 ], [ 4 , 2 ] ] # destination grid where we want data # Its larger than the source grid and slightly offset dest_lon = [ [-91., -91, -91, -91 ], [-90., -90, -90, -90 ], [-89., -89, -89, -89 ], [-88., -88, -88, -88 ] ] dest_lat = [ [ 43 , 44, 45, 46 ], [ 43 , 44, 45, 46 ], [ 43 , 44, 45, 46 ], [ 43 , 44, 45, 46 ] ] #instantiate object to handle interpolation proj_inds = geo_tools.ProjInds(src_lon=src_lon, src_lat=src_lat, dest_lon=dest_lon, dest_lat=dest_lat, missing=-99.) #interpolate data with "project_data" interpolated = proj_inds.project_data(data) #nearest neighbor output, pts outside the domain are set to missing #Interpolation with border detection in all directions expected = np.array([[-99., -99., -99., -99.,], [-99., 3., 1., -99.,], [-99., 4., 2., -99.,], [-99., -99., -99., -99.,]]) assert np.allclose(interpolated, expected) #on some domain, border detection is not desirable, it can be turned off # # extend_x here refers to the dimension in data space (longitudes) that are # represented along rows of python array. # for example: # Border detection in Y direction (latitudes) only proj_inds_ext_y = geo_tools.ProjInds(src_lon=src_lon, src_lat=src_lat, dest_lon=dest_lon, dest_lat=dest_lat, missing=-99., extend_x=False) interpolated_ext_y = proj_inds_ext_y.project_data(data) expected = np.array([[-99., 3., 1.,-99.], [-99., 3., 1.,-99.], [-99., 4., 2.,-99.], [-99., 4., 2.,-99.]]) assert np.allclose(interpolated_ext_y, expected) # # Border detection in X direction (longitudes) only proj_inds_ext_x = geo_tools.ProjInds(src_lon=src_lon, src_lat=src_lat, dest_lon=dest_lon, dest_lat=dest_lat, missing=-99., extend_y=False) interpolated_ext_x = proj_inds_ext_x.project_data(data) expected = np.array([[-99., -99., -99., -99.], [ 3., 3., 1., 1.], [ 4., 4., 2., 2.], [-99., -99., -99., -99.]]) assert np.allclose(interpolated_ext_x, expected) # # no border detection proj_inds_no_b = geo_tools.ProjInds(src_lon=src_lon, src_lat=src_lat, dest_lon=dest_lon, dest_lat=dest_lat, missing=-99., extend_x=False, extend_y=False) interpolated_no_b = proj_inds_no_b.project_data(data) expected = np.array([[3., 3., 1., 1.], [3., 3., 1., 1.], [4., 4., 2., 2.], [4., 4., 2., 2.]]) assert np.allclose(interpolated_no_b, expected)
Interpolation to coarser grids can be done with the nearest keyword. With this flag, all high-resoution data falling within a tile of the destination grid will be averaged together. Border detection works as in the example above.
import numpy as np import domutils.geo_tools as geo_tools # source data on a very simple grid src_lon = [ [-88.2 , -88.2 ], [-87.5 , -87.5 ] ] src_lat = [ [ 43.5 , 44.1 ], [ 43.5 , 44.1 ] ] data = [ [ 3 , 1 ], [ 4 , 2 ] ] # coarser destination grid where we want data dest_lon = [ [-92. , -92 , -92 , -92 ], [-90. , -90 , -90 , -90 ], [-88. , -88 , -88 , -88 ], [-86. , -86 , -86 , -86 ] ] dest_lat = [ [ 42 , 44 , 46 , 48 ], [ 42 , 44 , 46 , 48 ], [ 42 , 44 , 46 , 48 ], [ 42 , 44 , 46 , 48 ] ] #instantiate object to handle interpolation #Note the average keyword set to true proj_inds = geo_tools.ProjInds(src_lon=src_lon, src_lat=src_lat, dest_lon=dest_lon, dest_lat=dest_lat, average=True, missing=-99.) #interpolate data with "project_data" interpolated = proj_inds.project_data(data) #Since all high resolution data falls into one of the output #grid tile, they are all aaveraged together: (1+2+3+4)/4 = 2.5 expected = np.array([[-99., -99. ,-99., -99. ], [-99., -99. ,-99., -99. ], [-99., 2.5 ,-99., -99. ], [-99., -99. ,-99., -99. ]]) assert np.allclose(interpolated, expected) #weighted average can be obtained by providing weights for each data pt #being averaged weights = [ [ 0.5 , 1. ], [ 1. , 0.25 ] ] weighted_avg = proj_inds.project_data(data, weights=weights) #result is a weighted average: # (1.*1 + 0.25*2 + 0.5*3 + 1.*4) / (1.+0.25+0.5+1.) = 7.0/2.75 = 2.5454 expected = np.array([[-99., -99. , -99., -99. ], [-99., -99. , -99., -99. ], [-99., 2.54545455, -99., -99. ], [-99., -99. , -99., -99. ]]) assert np.allclose(weighted_avg, expected)Sometimes, it is useful to smooth data during the interpolation process. For example when comparing radar measurement against model output, smoothing can be used to remove the small scales present in the observations but that the model cannot represent.
Use the smoooth_radius to average all source data point within a certain radius (in km) of the destination grid tiles.
import numpy as np import domutils.geo_tools as geo_tools # source data on a very simple grid src_lon = [ [-88.2 , -88.2 ], [-87.5 , -87.5 ] ] src_lat = [ [ 43.5 , 44.1 ], [ 43.5 , 44.1 ] ] data = [ [ 3 , 1 ], [ 4 , 2 ] ] # coarser destination grid where we want data dest_lon = [ [-92. , -92 , -92 , -92 ], [-90. , -90 , -90 , -90 ], [-88. , -88 , -88 , -88 ], [-86. , -86 , -86 , -86 ] ] dest_lat = [ [ 42 , 44 , 46 , 48 ], [ 42 , 44 , 46 , 48 ], [ 42 , 44 , 46 , 48 ], [ 42 , 44 , 46 , 48 ] ] #instantiate object to handle interpolation #All source data points found within 300km of each destination #grid tiles will be averaged together proj_inds = geo_tools.ProjInds(src_lon=src_lon, src_lat=src_lat, dest_lat=dest_lat, dest_lon=dest_lon, smooth_radius=300., missing=-99.) #interpolate and smooth data with "project_data" interpolated = proj_inds.project_data(data) #output is smoother than data source expected = np.array([[-99. , -99. , -99., -99. ], [ 2.66666667, 2.5, 1.5, -99. ], [ 2.5 , 2.5, 2.5, -99. ], [ 2.5 , 2.5, 1.5, -99. ]]) assert np.allclose(interpolated, expected)
- plot_border(ax=None, crs=None, linewidth=0.5, zorder=3, mask_outside=False)
Plot border of a more or less regularly gridded domain
Optionally, mask everything that is outside the domain for simpler figures.
- Parameters:
ax (
Optional[Any]) – Instance of a cartopy geoAxes where border is to be plottedcrs (
Optional[Any]) – Cartopy crs of the data coordinates. If None, a Geodetic crs is used for using latitude/longitude coordinates.linewidth (
Optional[float]) – Width of border being plottedmask_outside (
Optional[bool]) – If True, everything outside the domain will be hiddenzorder (
Optional[int]) – zorder for the mask being applied with maskOutside=True. This number should be larger than any other zorder in the axes.
- Returns:
Nothing
- project_data(data, weights=None, output_avg_weights=False, missing_v=-9999.0)
Project data into geoAxes space
- Parameters:
data (
Any) – (array like), data values to be projected. Must have the same dimension as the lat/lon coordinates used to instantiate the class.weights (
Optional[Any]) – (array like), weights to be applied during averaging (see average and smooth_radius keyword in ProjInds)output_avg_weights (
Optional[bool]) – If True, projected_weights will be outputed along with projected_data “Averaged” and “smooth_radius” interpolation can take a long time on large grids. This option is useful to get radar data + quality index interpolated using one call to project_data instead of two halving computation time. Not available for nearest-neighbor interpolation which is already quite fast anyway.missing_v (
Optional[float]) – Value in input data that will be neglected during the averaging process Has no impact unless average or smooth_radius are used in class instantiation.
- Returns:
projected_data A numpy array of the same shape as destination grid used to instantiate the class
if output_avg_weights=True, returns:
projected_data, projected_weights
- domutils.geo_tools.projinds.geographic_extent_to_rotated(extent_geo, rotated_crs)
Convert a geographic (PlateCarree) extent into rotated-pole coordinates.
- Parameters:
extent_geo (list or tuple) – [lon_min, lon_max, lat_min, lat_max] in true geographic coords.
rotated_crs (cartopy.crs.RotatedPole) – The rotated pole CRS you want to convert into.
- Returns:
extent_rot – [rot_lon_min, rot_lon_max, rot_lat_min, rot_lat_max] in rotated coords.
- Return type:
list
- domutils.geo_tools.projinds.normalize_longitudes(lons)
Normalizes an array of longitudes to the range [-180, 180).
- domutils.geo_tools.lat_lon_extend.lat_lon_extend(lon1_in, lat1_in, lon2_in, lat2_in, half_dist=False, predefined_theta=None)
Extends latitude and longitudes on a sphere
Given two points (pt1 and pt2) on a sphere, this function returns a third point (pt3) at the same distance and in the same direction as between pt1 and pt2. All points are defined with latitudes and longitudes, same altitude is assumed.
The diagram below hopefully make the purpose of this function clearer:
# pt3 output # # ^ # # \ # # \ -> output with half_dist=True # # \ # # pt2 input 2 # # ^ # # \ # # \ # # \ # # pt1 input 1 #
- Parameters:
lon1_in (
Any) – (array like) longitude of point 1lat1_in (
Any) – (array like) latitude of point 1lon2_in (
Any) – (array like) longitude of point 2lat2_in (
Any) – (array like) latitude of point 2half_dist (
bool) – If true, pt 3 will be located at half the distance between pt1 and pt2predefined_theta – (array like) radians. If set, will determine distance between pt2 and pt3
- Returns:
(longitude, latitude) (array like) of extended point 3
Examples
Extend one point.
>>> import numpy as np >>> import domutils.geo_tools as geo_tools >>> # coordinates of pt1 >>> lat1 = 0.; lon1 = 0. >>> # coordinates of pt2 >>> lat2 = 0.; lon2 = 90. >>> # coordinates of extended point >>> lon3, lat3 = geo_tools.lat_lon_extend(lon1,lat1,lon2,lat2) >>> print(lon3, lat3) [180.] [0.]
Extend a number of points all at once
>>> # coordinates for four pt1 >>> lat1 = [ 0., 0., 0., 0. ] >>> lon1 = [ 0., 0., 0., 0. ] >>> # coordinates for four pt2 >>> lat2 = [ 0.,90., 0., 30.] >>> lon2 = [90., 0.,45., 0. ] >>> # coordinates of extended points >>> lon3, lat3 = geo_tools.lat_lon_extend(lon1,lat1,lon2,lat2) >>> with np.printoptions(precision=1,suppress=True): ... print(lon3) ... print(lat3) [180. 180. 90. 0.] [ 0. 0. 0. 60.]
- domutils.geo_tools.lat_lon_range_az.lat_lon_range_az(lon1_in, lat1_in, range_in, azimuth_in, crs=None, earth_radius=None)
Computes lat/lon of a point at a given range and azimuth (meteorological angles)
Given a point (pt1) on a sphere, this function returns the lat/lon of a second point (pt2) at a given range and azimuth from pt1. Same altitude is assumed. Range is a distance following curvature of the earth
The diagram below hopefully make the purpose of this function clearer:
# # # N # # # # ^ - # # | - # # | - azimuth (deg) # # | - # # | - # # pt1 - E # # \ - # # range \ - # # \ - # # \ - # # pt2 output is lat lon of pt2 #
- Internally, two rotations on the sphere are conducted
1st rotation shifts pt1 toward the North and is dictated by range
2nd rotation is around point pt1 and is determined by azimuth
Input arrays will be broadcasted together.
- Parameters:
lon1_in (
Any) – (array like) longitude of point 1lat1_in (
Any) – (array like) latitude of point 1range_in (
Any) – (array like) [km] range (distance) at which pt2 is located from pt1azimuth_in (
Any) – (array like) [degrees] direction (meteorological angle) in thich pt2 is locatedcrs (
Optional[Any]) – Instance of Cartopy geoAxes in which pt coordinates are defined (default is Geodesic for lat/lon specifications)earth_radius (
Optional[Any]) – (km) radius of the earth (default 6371 km)
- Returns:
longitude, latitude (array like) of pt2
Examples
Extend one point.
>>> import numpy as np >>> import domutils.geo_tools as geo_tools >>> # coordinates of pt1 >>> lat1 = 0. >>> lon1 = 0. >>> # range and azimuth >>> C = 2.*np.pi*6371. #circumference of the earth >>> range_km = C/8. >>> azimuth_deg = 0. >>> # coordinates of extended point >>> lon2, lat2 = geo_tools.lat_lon_range_az(lon1,lat1,range_km,azimuth_deg) >>> # print(lon2, lat2) should give approximately >>> # 0.0 45.192099833854435 >>> print(np.allclose((lon2, lat2), (0.0, 45.))) True >>> #Also works for inputs arrays >>> lat1 = [[0., 0], ... [0., 0]] >>> lon1 = [[0., 0], ... [0., 0]] >>> range_km = [[C/8., C/8.], ... [C/8., C/8.]] >>> azimuth_deg = [[0., 90.], ... [-90., 180]] >>> lon2, lat2 = geo_tools.lat_lon_range_az(lon1,lat1,range_km,azimuth_deg) >>> print(np.allclose(lon2, np.array([[ 0.0000000e+00, 4.5000000e+01], ... [-4.5000000e+01, 7.0167093e-15]]))) True
Since arrays are broadcasted together, inputs that consists of repeated values can be passed as floats. Here, we get the lat/lon along a circle 50km from a reference point
>>> #longitudes and latitudes >>> #of points along a circle of 50 km radius around a point at 37 deg N, 88 deg Ouest >>> #note how only azimuths is an array >>> lat0 = 37. >>> lon0 = -88. >>> ranges = 50. >>> azimuths = np.arange(0, 360, 10, dtype='float') >>> circle_lons, circle_lats = geo_tools.lat_lon_range_az(lat0, lon0, ranges, azimuths) >>> print(circle_lons) [37. 38.83132026 40.63355249 42.37690625 44.03018684 45.56010269 46.93060903 48.10234846 49.03230506 49.67387639 49.97769473 49.89369016 49.37504076 48.38468097 46.9046979 44.94791385 42.56907866 39.87096244 37. 34.12903756 31.43092134 29.05208615 27.0953021 25.61531903 24.62495924 24.10630984 24.02230527 24.32612361 24.96769494 25.89765154 27.06939097 28.43989731 29.96981316 31.62309375 33.36644751 35.16867974]
>>> print(circle_lats) [-87.5503392 -87.55592346 -87.57258224 -87.60003229 -87.63779731 -87.68520148 -87.74135996 -87.8051664 -87.87527736 -87.95009446 -88.02774645 -88.10607547 -88.18263607 -88.25472035 -88.31942859 -88.37380702 -88.41506673 -88.44087378 -88.4496608 -88.44087378 -88.41506673 -88.37380702 -88.31942859 -88.25472035 -88.18263607 -88.10607547 -88.02774645 -87.95009446 -87.87527736 -87.8051664 -87.74135996 -87.68520148 -87.63779731 -87.60003229 -87.57258224 -87.55592346]
- domutils.geo_tools.rotation_matrix.rotation_matrix(axis, theta)
Rotation matrix for counterclockwise rotation on a sphere
This formulation was found on SO https://stackoverflow.com/questions/6802577/rotation-of-3d-vector Thanks, “unutbu”
- Parameters:
axis (
Any) – Array-Like containing i,j,k coordinates that defines an axis of rotation starting from the center of the sphere (0,0,0)theta (
float) – The amount of rotation (radians) desired.
- Returns:
A 3x3 matrix (numpy array) for performing the rotation
Example
>>> import numpy as np >>> import domutils.geo_tools as geo_tools >>> #a rotation axis [x,y,z] pointing straight up >>> axis = [0.,0.,1.] >>> # rotation of 45 degrees = pi/4 ratians >>> theta = np.pi/4. >>> #make rotation matrix >>> mm = geo_tools.rotation_matrix(axis, theta) >>> print(mm.shape) (3, 3) >>> # a point on the sphere [x,y,z] which will be rotated >>> origin = [1.,0.,0.] >>> #apply rotation >>> rotated = np.matmul(mm,origin) >>> #rotated coordinates >>> print(rotated) [0.70710678 0.70710678 0. ]
- domutils.geo_tools.rotation_matrix.rotation_matrix_components(axis, theta)
Rotation matrix for counterclockwise rotation on a sphere
This version returns only the components of the rotation matrix(ces) This allows to use numpy array operations to simultaneously perform many application of these rotation matrices
- Parameters:
axis (
Any) – Array-Like containing i,j,k coordinates that defines an axes of rotation starting from the center of the sphere (0,0,0) dimension of array must be (m,3) for m points each defined by i,j,k componentstheta (
float) – The amount of rotation (radians) desired. shape = (m,)
- Returns:
A, B, C, D, E, F, G, H, I
| A B C | | D E F | | G H I |
Each of the components above will be of shape (m,)
Example
>>> import numpy as np >>> import domutils.geo_tools as geo_tools >>> #four rotation axes [x,y,z]; two pointing up, two pointing along y axis >>> axes = [[0.,0.,1.], [0.,0.,1.], [0.,1.,0.], [0.,1.,0.]] >>> # first and third points will be rotated by 45 degrees (pi/4) second and fourth by 90 degrees (pi/2) >>> thetas = [np.pi/4., np.pi/2., np.pi/4., np.pi/2.] >>> #make rotation matrices >>> a, b, c, d, e, f, g, h, i = geo_tools.rotation_matrix_components(axes, thetas) >>> #Four rotation matrices are here generated >>> print(a.shape) (4,) >>> # points on the sphere [x,y,z] which will be rotated >>> # here we use the same point for simplicity but any points can be used >>> oi,oj,ok = np.array([[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]).T #transpose to unpack columns >>> #i, j and k components of each position vector >>> print(oi.shape) (4,) >>> #apply rotation by hand (multiplication of each point (i,j,k vector) by rotation matrix) >>> rotated = np.array([a*oi+b*oj+c*ok, d*oi+e*oj+f*ok, g*oi+h*oj+i*ok]).T >>> print(rotated.shape) (4, 3) >>> print(np.allclose(rotated, np.array([[ 7.07106781e-01, 7.07106781e-01, 0.00000000e+00], ... [ 2.22044605e-16, 1.00000000e+00, 0.00000000e+00], ... [ 7.07106781e-01, 0.00000000e+00, -7.07106781e-01], ... [ 2.22044605e-16, 0.00000000e+00, -1.00000000e+00]]))) True