Using numpy linalg to replicate matrix operations: reinventing the wheel

2931
0
11-29-2014 02:38 AM
Labels (1)
DanPatterson_Retired
MVP Emeritus
1 0 2,931

This is simply recorded here so that I don't forget my matrix algebra stuff.
Since SciPy isn't rolled out for 10.2.x and using numpy's matrix operations adds another layer of things that you have to worry about, I thought I would play with numpy's linear algebra module as suggested in this scipy.linalg thread  

The demo here shows to translate and rotate a series of coordinates (represented as an array) about the mean center of the array.  The reverse is performed to transform back to the original.

>>> # working with the linear algebra module to do matrix stuff
>>> import numpy as np
>>> # import the linear algebra module
>>>
>>> from numpy.linalg import linalg as nla     # refer to it as nla
>>> 
>>> # unit square at 1,1 origin with mid-point 1.5,1.5
>>> 
>>> xs = [1,1,1.5,2,2]
>>> ys = [1,2,1.5,2,1]
>>> XYs = [[xs,ys] for i in range(len(xs))]
>>> XYmean = np.mean(XYs, axis=0)              # get the mean
>>> XYmean
array([ 1.5,  1.5])
>>> A = np.array(XYs)
>>> A                          # the array form of XYs list
array([[ 1. ,  1. ],
       [ 1. ,  2. ],
       [ 1.5,  1.5],
       [ 2. ,  2. ],
       [ 2. ,  1. ]])
>>> #   Translate/shift the array by adding the -ve mean
>>> A_t = nla.add(A,-XYmean)   # to subtract, just add -ve values
>>> A_t                        # the array shifted about the mean
array([[-0.5, -0.5],
       [-0.5,  0.5],
       [ 0. ,  0. ],
       [ 0.5,  0.5],
       [ 0.5, -0.5]])
>>> #  form the clockwise and counter-clockwise matrices, eg. 45 deg rotation
>>> a = np.radians(45.0)
>>> r_cw = np.array([[np.cos(a),np.sin(a)],[-np.sin(a),np.cos(a)]]) #rotate clockwise
>>> r_cc = np.array([[np.cos(a),-np.sin(a)],[np.sin(a),np.cos(a)]]) #rotate counter-clockwise 
>>> r_cw
array([[ 0.70710678,  0.70710678],
       [-0.70710678,  0.70710678]])
>>> r_cc
array([[ 0.70710678, -0.70710678],
       [ 0.70710678,  0.70710678]])
>>> A_tr = A_t.dot(r_cw.T)         # rotate the translated matrix by 45 degrees
>>> A_tr
array([[ -7.07106781e-01,  -5.55111512e-17],
       [ -5.55111512e-17,   7.07106781e-01],
       [  0.00000000e+00,   0.00000000e+00],
       [  7.07106781e-01,   5.55111512e-17],
       [  5.55111512e-17,  -7.07106781e-01]])
>>> 
>>> A_rb = A_tr.dot(r_cc.T)               # rotate back
>>> A_rb
array([[-0.5, -0.5],
       [-0.5,  0.5],
       [ 0. ,  0. ],
       [ 0.5,  0.5],
       [ 0.5, -0.5]])
>>> A_final = nla.add(A_rb,XYmean)        # translate back to its original form
>>> A_final
array([[ 1. ,  1. ],
       [ 1. ,  2. ],
       [ 1.5,  1.5],
       [ 2. ,  2. ],
       [ 2. ,  1. ]])
>>> 
>>> #  Done!!!! array translated and rotated, then returned to its original form.

Hope this helps someone else not reinvent the wheel.   Email comments

About the Author
Retired Geomatics Instructor at Carleton University. I am a forum MVP and Moderator. Current interests focus on python-based integration in GIS. See... Py... blog, my GeoNet blog...
Labels