numpy linalg.lstsq - coordinate translations

NeilAyres · ‎09-26-2016

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

DanPatterson_Retired · ‎09-26-2016

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

import numpy as np

def cent_pnt(pnts):
    """Return the centre of a point array"""
    pnts = np.asarray(pnts)
    return np.average(pnts, axis=0).squeeze()


def translate(pnts, cent=None):
    """Translate the points about the origin
    :Requires:
    :--------
    :  pnts - 2D sequence of point coordinates (N, 2)
    :  cent - a list, tuple, ndarray of xc, yc
    :
    :Returns:
    :-------
    :  p_trans - input points translated about the origin
    :
    :-------
    """
    pnts = np.asarray(pnts)
    if cent is None:
        p_trans = pnts - cent_pnt(pnts)
    else:
        p_trans = pnts - np.asarray(cent)
    return p_trans    
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

DanPatterson_Retired · ‎09-27-2016

Still busy, but this will delay you some more

Consider the following points:

>>> a = np.array([[1,1], [2,2], [3,3], [4,4]], dtype='float64')
>>> b = affine(a, angle=-22.5)[0]  # rotate to half-way to the x-axis
>>> c = b + [10,10]
>>>
>>> a
array([[ 1.00000000,  1.00000000],
       [ 2.00000000,  2.00000000],
       [ 3.00000000,  3.00000000],
       [ 4.00000000,  4.00000000]])
>>> b
array([[-1.95984445, -0.81179415],
       [-0.65328148, -0.27059805],
       [ 0.65328148,  0.27059805],
       [ 1.95984445,  0.81179415]])
>>> c
array([[ 8.04015555,  9.18820585],
       [ 9.34671852,  9.72940195],
       [ 10.65328148,  10.27059805],
       [ 11.95984445,  10.81179415]])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

This function (although only lightly tested) should translate and rotate the points from one form to another.

def warp(src=None, dest=None):
    """ move the src points based on the dest pnt affine parameters
    """
    a_slope, a_inter, a_cent = reg_params(src)
    b_slope, b_inter, b_cent = reg_params(dest)
    d_slope = b_slope - a_slope
    results = affine(src, d_slope)
    new_pnts = results[0] + b_cent
    return new_pnts, results‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

def reg_params(src=None):
    """ determine the correlation properties of a point set
    """
    a_cent = cent_pnt(src)
    a_t = translate(src, a_cent)
    xt = a_t[:,0]
    yt = a_t[:,1]
    #co_var = np.cov(m=src.T, y=None, rowvar=1, bias=0, ddof=None)
    #sigmX, sigmY = np.var(a, axis=0, dtype=None, out=None, ddof=0, keepdims=False)
    rho = np.corrcoef(x=xt, y=yt, rowvar=1, bias=0, ddof=None)  
    arr = np.vstack([xt, np.ones(len(xt))]).T
    #arr = xt.reshape(len(xt), 1)
    slope, inter = np.linalg.lstsq(arr, yt)[0]
    deg_slope = np.degrees(np.arctan2(slope,1.0))
    return deg_slope, inter, a_cent #, rho‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I left in 3 commented lines for examining outputs of the covariance and x, y variance calculations if one wants to examine them. The other line, just shows an alternate method of producing the arr variable.

So if we put this together (with reg_params duplicated), then we have

def cent_pnt(pnts):
    """Return the centre of a point array"""
    pnts = np.asarray(pnts)
    return np.average(pnts, axis=0).squeeze()


def translate(pnts, cent=None):
    """Translate the points about the origin
    :Requires:
    :--------
    :  pnts - 2D sequence of point coordinates (N, 2)
    :  cent - a list, tuple, ndarray of xc, yc
    :
    :Returns:
    :-------
    :  p_trans - input points translated about the origin
    :
    :-------
    """
    pnts = np.asarray(pnts)
    if cent is None:
        p_trans = pnts - cent_pnt(pnts)
    else:
        p_trans = pnts - np.asarray(cent)
    return p_trans    


def rot_matrix(rotation=0):
    """Return the rotation matrix given points and rotation angle
    :Requires:
    :--------
    :  pnts - 2D sequence of point coordinates (N, 2)
    :  rotation - angle in degrees
    :
    :Returns:
    :-------
    :  rot_mat - rotation matrix for 2D transform
    :
    :Notes:
    :-----
    :  for homogenous coordinates use
    :    rot_mat = np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])
    :--------
    """
    rad = np.deg2rad(rotation)
    c = np.cos(rad)
    s = np.sin(rad)
    rot_mat = np.array([[c, -s], [s, c]])
    return rot_mat

       
def affine(pnts, angle=0, cent=None):
    """Produce an affine transformation of a point set, given an angle of
    :  rotation and a point to rotate about.
    :Requires:
    :--------
    :  pnts - a series of 2D points a list or ndarray (N,2)
    :  angle - rotation angle, in degrees, relative to the x-axis.
    :  cent - center of rotation, either provided or determined from the
    :         center of the point cloud. [x,y], np.array([x,y])
    :
    :Returns:
    :  p_rot - the rotated points.
    :
    :Notes:
    :-----
    :  The functions that are used are: cent_pnt, translate, rot_matrix.
    :  Once the points are translated, and the rotation matrix is obtained,
    :  the points are transposed to facilite broadcasting, the dot product 
    :  determined and the resultant transposed back to the input shape.
    :
    """
    p = np.asarray(pnts)
    if cent is None:
        cent = cent_pnt(pnts)
    p_trans = translate(p, cent=None)
    rot_mat = rot_matrix(rotation=angle)
    p_rot = np.dot(rot_mat, p_trans.T).T
    return p_rot, cent, rot_mat‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Notice I did a slice on the warp in this example, because I didn't need the rest of the returned values from the warp. Looks pretty good when you compare the from-to results with the inputs.

>>> warp(a, b)[0]
array([[-1.95984445, -0.81179415],
       [-0.65328148, -0.27059805],
       [ 0.65328148,  0.27059805],
       [ 1.95984445,  0.81179415]])
>>> warp(a, c)[0]
array([[ 8.04015555,  9.18820585],
       [ 9.34671852,  9.72940195],
       [ 10.65328148,  10.27059805],
       [ 11.95984445,  10.81179415]])
>>>
>>> # some more ...
>>>
>>> warp(b, c)[0]
array([[ 8.04015555,  9.18820585],
       [ 9.34671852,  9.72940195],
       [ 10.65328148,  10.27059805],
       [ 11.95984445,  10.81179415]])
>>> warp(c, a)[0]
array([[ 1.00000000,  1.00000000],
       [ 2.00000000,  2.00000000],
       [ 3.00000000,  3.00000000],
       [ 4.00000000,  4.00000000]])
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I hope this gets the idea across Neil. This is somewhat simplified and I will still look for the fuller code but the basics entail determining the nature of the parameters of the input points (center, slope etc), then translating, rotating and translating to new center.

NeilAyres · ‎09-27-2016

That's fantastic Dan. Plenty to chew over here, but once i have digested this and hopefully understand what's going on, I will attempt to push this through my problem data.

Will get back to you.

I presume the complicated stuff from WH is in your archive somewhere.

But, this has, I hope, wider application with the broader community here. At least anyone who is faced with the problem of trying to define an ArcGIS style local projection from a set of measured points.

DanPatterson_Retired · ‎09-27-2016

do note, that the approach that was used by the other script uses the transformation matrices of both point clouds since their correlations are not perfect. The examples that I gave have perfect correlations so the full transformation is not needed. My examples were illustrative of how the process works and not a full 'warp'... until I get a chance to get to my files.

NeilAyres · ‎09-27-2016

Thanks Dan.

Looking forward to your digging into your files.

Will have a play with these codes blocks in the next few days.

NeilAyres · ‎09-28-2016

Okay Dan,

Been messing with that code you posted.

It seems that it goes a bit pear shaped if the rotation angle is > 45deg.

Just added a few prints etc... (inside def warp to display the angles and b_cent)

Test 1 :

>>> 
Using Rotation 44.5, Shift [1000, 2000]
Array A 
[[ 1. 1.]
 [ 2. 2.]
 [ 3. 3.]
 [ 4. 4.]]
Array B 
[[-0.01851178 -2.12123957]
 [-0.00617059 -0.70707986]
 [ 0.00617059 0.70707986]
 [ 0.01851178 2.12123957]]
Array C 
[[ 999.98148822 1997.87876043]
 [ 999.99382941 1999.29292014]
 [ 1000.00617059 2000.70707986]
 [ 1000.01851178 2002.12123957]]
inWarp - Slope A 45.0, B 89.5, D 44.5
Shifts [ 0. 0.]
Warp A -> B
[[-0.01851178 -2.12123957]
 [-0.00617059 -0.70707986]
 [ 0.00617059 0.70707986]
 [ 0.01851178 2.12123957]]
inWarp - Slope A 45.0, B 89.5, D 44.5
Shifts [ 1000. 2000.]
Warp A -> C
[[ 999.98148822 1997.87876043]
 [ 999.99382941 1999.29292014]
 [ 1000.00617059 2000.70707986]
 [ 1000.01851178 2002.12123957]]
inWarp - Slope A 89.5, B 89.5, D 1.67688085639e-12
Shifts [ 1000. 2000.]
Warp B -> C
[[ 999.98148822 1997.87876043]
 [ 999.99382941 1999.29292014]
 [ 1000.00617059 2000.70707986]
 [ 1000.01851178 2002.12123957]]
inWarp - Slope A 89.5, B 45.0, D -44.5
Shifts [ 2.5 2.5]
Warp C -> A
[[ 1. 1.]
 [ 2. 2.]
 [ 3. 3.]
 [ 4. 4.]]
>>>‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Test 2 :

>>> 
Using Rotation 45.5, Shift [1000, 2000]
Array A 
[[ 1.  1.]
 [ 2.  2.]
 [ 3.  3.]
 [ 4.  4.]]
Array B 
[[ 0.01851178 -2.12123957]
 [ 0.00617059 -0.70707986]
 [-0.00617059  0.70707986]
 [-0.01851178  2.12123957]]
Array C 
[[ 1000.01851178  1997.87876043]
 [ 1000.00617059  1999.29292014]
 [  999.99382941  2000.70707986]
 [  999.98148822  2002.12123957]]
inWarp - Slope A 45.0, B -89.5, D -134.5
Shifts [ 0.  0.]
Warp A -> B
[[-0.01851178  2.12123957]
 [-0.00617059  0.70707986]
 [ 0.00617059 -0.70707986]
 [ 0.01851178 -2.12123957]]
inWarp - Slope A 45.0, B -89.5, D -134.5
Shifts [ 1000.  2000.]
Warp A -> C
[[  999.98148822  2002.12123957]
 [  999.99382941  2000.70707986]
 [ 1000.00617059  1999.29292014]
 [ 1000.01851178  1997.87876043]]
inWarp - Slope A -89.5, B -89.5, D -1.67688085639e-12
Shifts [ 1000.  2000.]
Warp B -> C
[[ 1000.01851178  1997.87876043]
 [ 1000.00617059  1999.29292014]
 [  999.99382941  2000.70707986]
 [  999.98148822  2002.12123957]]
inWarp - Slope A -89.5, B 45.0, D 134.5
Shifts [ 2.5  2.5]
Warp C -> A
[[ 4.  4.]
 [ 3.  3.]
 [ 2.  2.]
 [ 1.  1.]]
>>> ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

You see everything gets neatly reversed. And the shifts from C -> A don't look right.

B -> C always works, because that's just a shift and no rotation.

Thanks for your cooperation and work on this so far.

Neil

DanPatterson_Retired · ‎09-28-2016

Neil... as I said... stripped down illustrative examples... not full all-encompassing working code (b-c has a rotation... it is just 0... of course there is 90-angle to look at.)

I have to finish my sabbatical report by Friday, ergo, some of the others that have viewed the thread are free to leap in... (we know who you are... we know you can do it)

NeilAyres · ‎09-28-2016

Sorry Dan, better let you get on with it...

NeilAyres · ‎10-06-2016

Dan_Patterson‌ could I try to revive your interest in this?

Hopefully you are now close to your "archive" and ready to share some of that mathematical wizardry from yourself or Bill Huber.