GSCUser85

numpy linalg.lstsq - coordinate translations

Discussion created by GSCUser85 Champion on Sep 26, 2016
Latest reply on Oct 18, 2016 by Dan_Patterson

I have some control points from a local grid to a known grid (a national grid system).

What I am trying to do is come up with the appropriate parameters to define a local prj for the input data.

This involves finding the origin (or rotation point), the amount of rotation, and any shifts or scaling in X & Y.

 

So, I have been messing with this code found here :

how-to-perform-coordinates-affine-transformation-using-python-part-2

 

So, what I have got is (almost the same as the referenced post ):

import sys, os
import numpy as np
localData = [["AJP1", 5671.666, 19466.156], ["GCP CR", 9634.659, 13064.567], \
["MHP2", 10325.093, 7797.802], ["MILE 8", 6508.217, 7511.953]]
refData =[["AJP1", 157171.409, 75402.914], ["GCP CR", 158414.820, 67978.050], \
["MHP2", 157059.513, 62842.546], ["MILE 8", 153418.976, 64023.272]]

inArr = np.array([[l[1], l[2], 0.0] for l in localData])
refArr = np.array([[l[1], l[2], 0.0] for l in refData])
n = inArr.shape[0]
pad = lambda x: np.hstack([x, np.ones((x.shape[0], 1))])
unpad = lambda x: x[:,:-1]
X = pad(inArr)
Y = pad(refArr)
# Solve the least squares problem X * A = Y
# to find our transformation matrix A
A, res, rank, s = np.linalg.lstsq(X, Y)
transform = lambda x: unpad(np.dot(pad(x), A))
print "Target:"
print refArr
print "Result:"
print transform(inArr)
print "Max error:", np.abs(refArr - transform(inArr)).max()

 

which give me this :

>>> 
Target:
[[ 157171.409 75402.914 0. ]
[ 158414.82 67978.05 0. ]
[ 157059.513 62842.546 0. ]
[ 153418.976 64023.272 0. ]]
Result:
[[ 157171.39578228 75402.91459879 0. ]
[ 158414.84986487 67978.04864706 0. ]
[ 157059.48564043 62842.54723945 0. ]
[ 153418.98671243 64023.2715147 0. ]]
Max error: 0.0298648653843

 

So, I can see the fit is quite good.

But how do I use the solution from np.linalg.lstsq to derive the parameters I need for the projection definition of the localData?

In particular, the origin point 0,0 in the target coordinates, and the shifts and rotations that are going on here??

 

Tagging out very own numpy expert and all around math wiz Dan Patterson here.

Or is this not even possible with this approach.

My matrix math is extremely rusty to non-existent. But I am willing to try any methods offered 

 

Thanks in advance.

Neil

Outcomes