|
|
| GeoCommunity Mailing List |
| |
| Mailing List Archives |
| Subject: | Re: [gislist] Coordinate Geometry assistance needed |
| Date: |
05/23/2007 08:58:47 AM |
| From: |
Quantitative Decisions |
|
|
At 10:58 AM 5/22/2007 -0500, David T. Hughes wrote:
>Of the several hundred points I know the UTM coordinates of THREE. How
>can I rotate the original points into the utm systems and plot them?
There is a straightforward way to compute the unique affine transformation
defined by three such "links." It begins with the elegant trick of
adjoining columns of ones after the y-coordinates, as shown below. This
produces a pair of three-by-three matrices, one formed from the observed
X,Y,1 coordinates and the other formed by the world X,Y,1
coordinates. Multiply the inverse of the first by the second to obtain the
transformation matrix. (This matrix multiplies row vectors on the right.)
In Excel you can do this with a matrix formula, as shown. The example
begins in cell A1. The formula occupies cells F2:H4 (that is, a three by
three array) and is named 'G'. (To enter a matrix formula in Excel, you
select the entire range the result is intended to occupy, type in the
formula, and complete it with a shift-Enter keystroke.)
Grid_X Grid_Y Unit X_Coord Y_Coord Unit Matrix
1008.854 1009.658 1 245476 4114418
1 =(MMULT(MINVERSE(A2:C4), D2:F4))
1012.502 1007.601 1 245488 4114421 1
985.108 1068.04 1 245476 4114442 1
To apply this transformation to any point (x,y), extend it to the
three-vector (x,y,1) and right-multiply it by the matrix, as in this
example (occupying cells A5:E5 of the spreadsheet):
0 0 1 =MMULT(A5:C5, G)
The formula occupies the range D5:E5 (that is, two columns for the world X
and world Y values). Copying the three right columns (C through E) down
all the rows of data will complete the transformation of your points. Plot
them using any software you like.
Any software platform that can invert and multiply three by three matrices
and vectors can perform these calculations. If you have matrix-oriented
software, such as Matlab, Mathematica, or R (which is free), you can do the
whole thing in one swoop. Let the three by three matrix of observed
coordinates (including the column of ones) be 'U', the three by three
matrix of world coordinates be 'V', and the 'n' by three matrix of observed
values you want to transform be 'X' (so that its third column is all
ones). The transformed values are found in the first two columns of the
'n' by three matrix
X * (U^-1 * V).
With these data, the 'G' matrix (equal to U^-1 * V) is
4.26841773841170 1.36788670604437 0.00000000
1.73611468632225 0.96745294319408 0.00000000
239416.907608956 4112061.20542121 1.00000000
The unequal off-diagonal elements (1.367... and 1.736...) indicate this is
not just a rotation and rescaling. The Transform2D extension for ArcView
3.x (available on the ESRI Arcscripts site) describes this matrix as
>A horizontal skew (y --> y + n*x) of size n=0.43472422242978
>Followed by a squeeze (x --> t*x, y --> y/t) of size t=3.3837309329694
>Followed by a rescaling (x --> r*x, y --> r*y) of r=1.324645274864
>Followed by a rotation in the positive (counterclockwise) direction of
>17.76893612854 degrees
>Followed by a translation (x --> x + e, y --> y + f) of (e,f) =
>(239416.90760896, 4112061.2054212).
The tremendous skew and squeeze components, along with a rescaling that
departs appreciably from unity, reveal a lot of distortion. I suspect at
least one of the world coordinates is incorrect. (It is suspicious that
two of the world X values are the same, but even that will not fix the
problem.)
By the way, the transformation computed from any three links will usually
be reasonably accurate within the triangular space they define, but it will
lose accuracy at a quadratic rate as you move away from the triangle's
center. If any of your other points lies outside the triangle, beware!
Cheers,
Bill Huber
_______________________________________________
gislist mailing list
gislist@lists.geocomm.com
http://lists.geocomm.com/mailman/listinfo/gislist
_________________________________
This list is brought to you by
The GeoCommunity
http://www.geocomm.com/
|
|
Sponsored by:
For information
regarding
advertising rates
Click Here!
|
