Python, NumPy and matplotlib - Visualize Matrix Rotations

I am using Python, NumPy and matplotlib to experiment with, and visualize 3-D graphics techniques.  The first step was figuring out how to generate polar plots with matplotlib and to determine the domain and range of spherical to Cartesian coordinate conversions.

The next step is to rotate a point in space using matrices, and then to explore the concept of gimbal lock.  I will probably tackle gimbal lock in a future posting.  Before delving into 3-D rotations, I needed to figure out how to extend my plotting knowledge.

I cleaned up the code from the previous posting, moving the spherical to Cartesian coordinates conversion code to a separate module - sphcoords. sphcoords also includes a function that takes spherical coordinates and returns a 1x3 NumPy matrix that represents a point in space.

import numpy as np
from numpy import matrix

def x_cart(incl_sph_coords, azim_sph_coords):
   return np.sin(incl_sph_coords) * np.cos(azim_sph_coords)

def y_cart(incl_sph_coords, azim_sph_coords):
   return np.sin(incl_sph_coords) * np.sin(azim_sph_coords)

def z_cart(incl_sph_coords):
   return np.cos(incl_sph_coords)

def sph2cart(incl_sph_coord, azim_sph_coord):
   return np.matrix([
   x_cart(incl_sph_coord, azim_sph_coord),
   y_cart(incl_sph_coord, azim_sph_coord),
   z_cart(incl_sph_coord)])

Below is an example of sphcoords function, sph2cart which generates a 1x3 matrix containing Cartesian coordinates. In this case, spherical coordinates of 30 degrees inclination and 180 degrees azimuth are converted.

poi = sphcoords.sph2cart(np.pi/6, -np.pi+(10*azi_incr))

By using matplotlib's subplot capability along with spherical to Cartesian coordinate conversions, I was able to depict a 3-D spherical model of space.  A single red "point of interest" will be rotated using rotational matrices.

The initial position of the red point is arbitrary.  The red point will be initially rotated about the z axis.  The rotated position of the red point will remain along the point's original circle of inclination.  This is what I was hoping  to illustrate with the cyan-colored annotations.

Before going further, let me apologize for my inability to generate circles. The elliptical, squashed figures in my plots will be circular after I figure out how to set fixed aspect ratios.

The Initial Position of the Red Point of Interest

Here are the rotational matrices that were implemented in module rotmat using the NumPy Matrix. 

import numpy as np
from numpy import matrix

def x_rot(rads):
  return np.matrix([
  [1,0,0], 
  [0, np.cos(rads), -np.sin(rads)],
  [0, np.sin(rads), np.cos(rads)]])

def y_rot(rads):
  return np.matrix([
  [np.cos(rads), 0, np.sin(rads)],
  [0, 1, 0], 
  [-np.sin(rads), 0, np.cos(rads)]])

def z_rot(rads):
  return np.matrix([
  [np.cos(rads), -np.sin(rads), 0],
  [np.sin(rads), np.cos(rads), 0],
  [0,0,1]]) 

The 1x3 matrix containing the point of interest's (POI) Cartesian coordinates is rotated about the z axis thusly:

poi = poi*rotmat.z_rot(-5*azi_incr)

The plot below shows the POI rotated 90 degrees about the z axis.  The cyan annotations show the range of the next rotation which is about the x axis.  The POI will be rotated 30 degrees about the x axis.  The most interesting update will be on the upper right plot that will show the the POI jumping from the yellow 30 degrees circle to the black circle at 60 degrees.

After Rotating the POI about the Z axis.

After Rotating the POI about the X axis.

 

The final position of the POI after rotating about each Euler Axis.

 


import numpy as np
import matplotlib.pylab as plt
from numpy import matrix
import rotmat
import sphcoords


def generate_plot(inclin, azi, poi):
# plot x-y plane
plt.subplot(221)
plt.plot(
   sphcoords.x_cart(inclin[0],azi[0]),
   sphcoords.y_cart(inclin[0],azi[0]), "g-o",
   sphcoords.x_cart(inclin[1],azi[1]),
   sphcoords.y_cart(inclin[1],azi[1]), "y-o",
   sphcoords.x_cart(inclin[2],azi[2]),
   sphcoords.y_cart(inclin[2],azi[2]), "y-o",
   sphcoords.x_cart(inclin[3],azi[3]),
   sphcoords.y_cart(inclin[3],azi[3]), "k-o",
   sphcoords.x_cart(inclin[4],azi[4]),
   sphcoords.y_cart(inclin[4],azi[4]), "k-o",
   sphcoords.x_cart(inclin[5],azi[5]),
   sphcoords.y_cart(inclin[5],azi[5]), "b-o",
   sphcoords.x_cart(inclin[6],azi[6]),
   sphcoords.y_cart(inclin[6],azi[6]), "b-o",
   poi[0,0], poi[0,1], "r-o")
plt.grid(True)
plt.ylabel('Y Axis')
plt.xlabel('X Axis')

# plot y-z plane
plt.subplot(222)
plt.plot(
   sphcoords.z_cart(inclin[0]),
   sphcoords.y_cart(inclin[0],azi[0]), "g-o",
   sphcoords.z_cart(inclin[1]),
   sphcoords.y_cart(inclin[1],azi[1]), "y-o",
   sphcoords.z_cart(inclin[2]),
   sphcoords.y_cart(inclin[2],azi[2]), "y-o",
   sphcoords.z_cart(inclin[3]),
   sphcoords.y_cart(inclin[3],azi[3]), "k-o",
   sphcoords.z_cart(inclin[4]),
   sphcoords.y_cart(inclin[4],azi[4]), "k-o",
   sphcoords.z_cart(inclin[5]),
   sphcoords.y_cart(inclin[5],azi[6]), "b-o",
   sphcoords.z_cart(inclin[6]),
   sphcoords.y_cart(inclin[6],azi[6]), "b-o",
   poi[0,2], poi[0,1], "r-o")
plt.xlabel('Z Axis')
plt.grid(True)

#plot x-z plane
plt.subplot(223)
plt.plot(
   sphcoords.x_cart(inclin[0],azi[0]),
   sphcoords.z_cart(inclin[0]),"g-o",
   sphcoords.x_cart(inclin[1],azi[1]),
   sphcoords.z_cart(inclin[1]),"y-o",
   sphcoords.x_cart(inclin[2],azi[2]),
   sphcoords.z_cart(inclin[2]),"y-o",
   sphcoords.x_cart(inclin[3],azi[3]),
   sphcoords.z_cart(inclin[3]),"k-o",
   sphcoords.x_cart(inclin[4],azi[4]),
   sphcoords.z_cart(inclin[4]),"k-o",
   sphcoords.x_cart(inclin[5],azi[5]),
   sphcoords.z_cart(inclin[5]),"b-o",
   sphcoords.x_cart(inclin[6],azi[6]),
   sphcoords.z_cart(inclin[6]),"b-o",
   poi[0,0], poi[0,2], "r-o")
plt.ylabel('Z Axis')
plt.grid(True)
plt.show()

if __name__ == '__main__':
   inclin = []
   azi = []

   inclin_points = range(0,20)
   azi_incr = (2 * np.pi)/float(20)
   azi_points = np.arange(-np.pi, np.pi, azi_incr)

   #90 deg inclination
   inclin.append([np.pi/2 for i in inclin_points])
   #30,150 deg inclinations
   inclin.append([np.pi/6 for i in inclin_points])
   inclin.append([np.pi-np.pi/6 for i in inclin_points])
   #60,120 deg inclinations
   inclin.append([np.pi/3 for i in inclin_points])
   inclin.append([np.pi-np.pi/3 for i in inclin_points])
   # poles
   inclin.append([0])
   inclin.append([np.pi])

   azi.append(azi_points)
   azi.append(azi_points)
   azi.append(azi_points)
   azi.append(azi_points)
   azi.append(azi_points)
   azi.append([np.pi])
   azi.append([np.pi])

   poi = sphcoords.sph2cart(np.pi/6, -np.pi+(10*azi_incr))
   poi = poi*rotmat.z_rot(-5*azi_incr)
   poi = poi*rotmat.x_rot(np.pi/6)
   poi = poi*rotmat.y_rot(np.pi/2)
   generate_plot(inclin, azi, poi)