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, 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 regridded

• dest_lat (Optional[Any]) – Latitude of destination grid on which the data will be regridded

• extent (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]

• dest_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. 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 numpy as np
>>> import matplotlib as mpl
>>> import matplotlib.pyplot as plt
>>> import cartopy.crs as ccrs
>>> import cartopy.feature as cfeature
>>> import domutils.legs as legs
>>> import domutils.geo_tools as geo_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 = ccrs.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 = ccrs.Geodetic()
>>> geo_cent = proj_ll.as_geocentric()
>>> xyz = geo_cent.transform_points(proj_ll, np.asarray(regular_lons),
...                                          np.asarray(regular_lats))
>>>
>>> #lets rotate points by 45 degrees counter clockwise
>>> theta = np.pi/4
>>> rotation_matrix = geo_tools.rotation_matrix([xyz[1,1,0], #x
...                                              xyz[1,1,1], #y
...                                              xyz[1,1,2]],#z
...                                              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,[xyz[ii,jj,0], #x
...                                                           xyz[ii,jj,1], #y
...                                                           xyz[ii,jj,2]])#z
>>>
>>> #from xyz to lat/lon
>>> rotated_latlon = proj_ll.transform_points(geo_cent, rotated_xyz[:,:,0],
...                                                     rotated_xyz[:,:,1],
...                                                     rotated_xyz[:,:,2])
>>> rotatedLons = rotated_latlon[:,:,0]
>>> rotatedLats = rotated_latlon[:,:,1]
>>> # 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
>>> ProjInds = geo_tools.ProjInds(rotatedLons, 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(map_extent)
>>>
>>> # 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 = ProjInds.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(cfeature.STATES.with_scale('50m'),
...                      linewidth=0.5, edgecolor='0.2',zorder=1)
>>>
>>> #plot border and mask everything outside model domain
>>> ProjInds.plot_border(ax, mask_outside=False, linewidth=2.)
>>>
>>> #uncomment to save figure
>>> plt.savefig('_static/projection_demo.svg')

Example showing how ProjInds can also be used for nearest neighbor interpolation

>>> import numpy as np
>>>
>>> # 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
>>> ProjInds = geo_tools.ProjInds(data_xx=src_lon,   data_yy=src_lat,
...                               dest_lon=dest_lon, dest_lat=dest_lat,
...                               missing=-99.)
>>> #interpolate data with "project_data"
>>> interpolated = ProjInds.project_data(data)
>>> #nearest neighbor output, pts outside the domain are set to missing
>>> #Interpolation with border detection in all directions
>>> print(interpolated)
[[-99. -99. -99. -99.]
 [-99.   3.   1. -99.]
 [-99.   4.   2. -99.]
 [-99. -99. -99. -99.]]
>>>
>>>
>>> #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(data_xx=src_lon,   data_yy=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)
>>> print(interpolated_ext_y)
[[-99.   3.   1. -99.]
 [-99.   3.   1. -99.]
 [-99.   4.   2. -99.]
 [-99.   4.   2. -99.]]
>>> #
>>> # Border detection in X direction (longitudes) only
>>> proj_inds_ext_x = geo_tools.ProjInds(data_xx=src_lon,   data_yy=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)
>>> print(interpolated_ext_x)
[[-99. -99. -99. -99.]
 [  3.   3.   1.   1.]
 [  4.   4.   2.   2.]
 [-99. -99. -99. -99.]]
>>> #
>>> # no border detection
>>> proj_inds_no_b = geo_tools.ProjInds(data_xx=src_lon,   data_yy=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)
>>> print(interpolated_no_b)
[[3. 3. 1. 1.]
 [3. 3. 1. 1.]
 [4. 4. 2. 2.]
 [4. 4. 2. 2.]]

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
>>> # 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
>>> ProjInds = geo_tools.ProjInds(data_xx=src_lon,   data_yy=src_lat,
...                               dest_lon=dest_lon, dest_lat=dest_lat,
...                               average=True, missing=-99.)
>>>
>>> #interpolate data with "project_data"
>>> interpolated = ProjInds.project_data(data)
>>>
>>> #Since all high resolution data falls into one of the output
>>> #grid tile, they are all aaveraged togetherL:  (1+2+3+4)/4 = 2.5
>>> print(interpolated)
[[-99.  -99.  -99.  -99. ]
 [-99.  -99.  -99.  -99. ]
 [-99.    2.5 -99.  -99. ]
 [-99.  -99.  -99.  -99. ]]
>>>
>>> #weighted average can be obtained by providing weights for each data pt
>>> #being averaged
>>> weights   =     [ [  0.5   ,  1.    ],
...                   [  1.    ,  0.25  ] ]
>>>
>>> weighted_avg = ProjInds.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
>>> print(weighted_avg)
[[-99.         -99.         -99.         -99.        ]
 [-99.         -99.         -99.         -99.        ]
 [-99.           2.54545455 -99.         -99.        ]
 [-99.         -99.         -99.         -99.        ]]

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.

>>> # 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
>>> ProjInds = geo_tools.ProjInds(data_xx=src_lon,    data_yy=src_lat,
...                               dest_lat=dest_lat,  dest_lon=dest_lon,
...                               smooth_radius=300., missing=-99.)
>>>
>>> #interpolate and smooth data with "project_data"
>>> interpolated = ProjInds.project_data(data)
>>>
>>> #output is smoother than data source
>>> print(interpolated)
[[-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.        ]]

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 plotted

• crs (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 plotted

• mask_outside (Optional[bool]) – If True, everything outside the domain will be hidden

• zorder (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.lat_lon_extend.lat_lon_extend(lon1_in, lat1_in, lon2_in, lat2_in, crs=None, 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 1

• lat1_in (Any) – (array like) latitude of point 1

• lon2_in (Any) – (array like) longitude of point 2

• lat2_in (Any) – (array like) latitude of point 2

• crs (Optional[Any]) – Instance of Cartopy geoAxes in which pt coordinates are defined (default is PlateCarree)

• half_dist (bool) – If true, pt 3 will be located at half the distance between pt1 and pt2

• predefined_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.  59.8]
>>> #           ^
>>> #           Presumably not =60. because cartopy uses a spheroid

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 1

• lat1_in (Any) – (array like) latitude of point 1

• range_in (Any) – (array like) [km] range (distance) at which pt2 is located from pt1

• azimuth_in (Any) – (array like) [degrees] direction (meteorological angle) in thich pt2 is located

• crs (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.192099833854435)))
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.82129197 40.61352908 42.34696182 43.99045278 45.51079424
 46.8720648  48.03508603 48.9570966  49.59184727 49.89044518 49.80343238
 49.28472885 48.29808711 46.82635859 44.88285245 42.52223827 39.8463912
 37.         34.1536088  31.47776173 29.11714755 27.17364141 25.70191289
 24.71527115 24.19656762 24.10955482 24.40815273 25.0429034  25.96491397
 27.1279352  28.48920576 30.00954722 31.65303818 33.38647092 35.17870803]
>>> print(circle_lats)
[-87.55334031 -87.55889412 -87.57546173 -87.60276046 -87.64031502
 -87.68745105 -87.74328587 -87.80671606 -87.87640222 -87.95075157
 -88.02790057 -88.10570197 -88.18172472 -88.25328003 -88.31749246
 -88.37143711 -88.4123562  -88.43794478 -88.44665642 -88.43794478
 -88.4123562  -88.37143711 -88.31749246 -88.25328003 -88.18172472
 -88.10570197 -88.02790057 -87.95075157 -87.87640222 -87.80671606
 -87.74328587 -87.68745105 -87.64031502 -87.60276046 -87.57546173
 -87.55889412]

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 components

• theta (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