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

from a cursory look you are trying to solve the matrix, using homogenous coordinates to pick up the rotation terms. I can't deal with it right now, but I will do so at a later time. When WH helped me through some of this some years ago, a translation about the origin was made by subracting the mean X, and Y from the data set. This centers the point cloud about the origin. Solving for the least squares, essentially gives you the correlation coefficient of the points which is related to the axis of rotation of the points from which one gets the angle of rotation. Think of a perfectly straight line of points parallel to the x axis. Translated to the origin would center the points on the x axis with points to the left of the y-axis and some to the right. The deviation in the y direction would be zero, the correlation coefficient would be zero and therefore the axis of rotation would be zero.

I will leave the above paragraph for you to mull over while I attend to business then dig up my complete affine transformation code.

Others may wade in so Neil doesn't have to wait

here some bits