Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am at the mercy of a govt authority who have provided utility information in a cad format with a 0,0 origin and a nonsense grid. About a dozen files, the grid and origin are entirely arbitrary.

To do an affine on a vector layer in QGIS, I propose to the use the Geoprocessing Plugin - Vector Affine Transformation.

But my math skills betray me....

The matching points for the ground control are

pt1 raw 6.42775 3.61140
pt1 gcp 567742 6450622

pt2 raw 2.7524 4.1539
pt2 gcp 562285 6450858

pt3 raw 6.0293 6.4965
pt3 gcp 566709 6454802

(edited x 2 , there are now 3 points here) A little rough on the correct locations.

The ground control is a UTM format. What would be the form of the equation, so as I can set it up in a spreadsheet and plod away?

The input for the the plugin is Scale X Scale Y Rotation X Y (not used I expect) Translation X Translation Y

Is this the most elegant workflow?

many thanks

share|improve this question
    
Can you give us the local grid coordinates for three points together with their UTM equivalents? Make them as widely-spread as possible. One test point isn't enough :) Nick. –  nhopton Apr 23 '12 at 10:25
    
Duplicate? gis.stackexchange.com/questions/21696/… –  Simbamangu Apr 23 '12 at 16:49
    
Simbamangu, indeed it is a duplicate but I am in another coordinate system. –  Willy Apr 23 '12 at 21:00
add comment

2 Answers 2

up vote 6 down vote accepted

Fitting equations to data is a statistical procedure, so it's best to use statistical software for the purpose. Powerful, easy-to-use software is freely available to all at the R Project.

Because the sample data are too limited to show all the essential capabilities of R, for illustrating its use let's add a fourth link matching the "raw" (source) point (5,5) to the "gcp" (destination) point (565440.0, 6452443.9). Although R can read many file formats, with such a small number of links it's straightforward to type them in, which can be done like this with two lines of code:

source <- t(matrix(c(6.42775,3.61140,  2.7524,4.1539,  6.0293,6.4965, 5.0,5.0), nrow=2))
dest <- t(matrix(c(567742,6450622,  562285,6450858,  566709,6454802, 565440.0,6452443.9), nrow=2))

These create arrays in which each row displays a pair of coordinates. As a check, you can obtain a nicer array-like display of these data by retyping the variable names at the command prompt (>). For example:

> source
        [,1]   [,2]
[1,] 6.42775 3.6114
[2,] 2.75240 4.1539
[3,] 6.02930 6.4965
[4,] 5.00000 5.0000

There they are, one per line. Regardless, once you get the links into R in this form, finding the best (least squares) fit is incredibly simple:

lm(dest ~ source)

That's all there is to it. Here's what the output looks like:

Coefficients:
             [,1]       [,2]     
(Intercept)   558910.8  6444331.5
source1         1461.7      152.8
source2         -156.1     1469.9

This transformation matrix multiplies source coordinates on the right. An example appears below.

It's a good idea to probe the solution more closely. Let's start over, this time saving the fit:

a <- lm(dest ~ source)

Summarize the results to get more information:

> summary(a)
Response Y1 :

Call:
lm(formula = Y1 ~ source)

Residuals:
      1       2       3       4 
-0.3227 -0.4617 -0.5605  1.3449 

Coefficients:
              Estimate Std. Error  t value Pr(>|t|)    
(Intercept)  5.589e+05  4.062e+00 137586.8 4.63e-06 ***
source1      1.462e+03  5.642e-01   2590.7 0.000246 ***
source2     -1.561e+02  7.387e-01   -211.3 0.003013 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.562 on 1 degrees of freedom
Multiple R-squared:     1,      Adjusted R-squared:     1 
F-statistic: 3.447e+06 on 2 and 1 DF,  p-value: 0.0003809 
<SNIP>

This is the result for the x coordinates only. The excised material is another panel of output headed by "Response Y2"; it is for the y coordinates. In each panel you want to check that the residual standard error (expressed in world coordinates, usually meters) is acceptably small. Then you want to check that the residuals themselves (listed at the top) are individually within your tolerances for a fit: these are the differences between the data and the predicted values. The residual standard error tells us that the typical discrepancy between a fitted x-coordinate and the intended (gcp) x-coordinate is around 1.56 meters.

As you know, relatively high numerical precision is needed in such work. This summary doesn't have sufficient precision for, say, a world file. Instead, extract the coefficients with

> print(a$coefficients, digits=10)
                      [,1]            [,2]
(Intercept) 558910.7522289 6444331.5258144
source1       1461.6779888     152.7722764
source2       -156.0974128    1469.8705760

As an example of how these coefficients work, we can ask R to perform the matrix multiplication to find where (say) the point (5,5) goes, remembering that the intercept is handled by extending these two-dimensional vectors to three dimensions by inserting a 1 at the beginning:

> print(c(1, 5,5) %*% a$coefficients, digits=10)
            [,1]       [,2]
[1,] 565438.6551 6452444.74

The first coordinate is about 1.35 low and the second one is about 0.84 high: these are the negatives of the reported residuals.

Much, much more can be done here. For example, any array of "raw" points can now be transformed in a single command; more complex fitting can be done (such as polynomial fits of higher order); etc. The power, flexibility, and accuracy of this software recommend it as a general-purpose solution to this problem of transforming points from one coordinate system to another.


Incidentally, repeating these analyses with just the three points in the example produces the affine transformation matrix

                           [,1]            [,2]
(Intercept)      558915.7606071 6444328.3973676
1                  1463.4093542     151.6907917
2                  -158.5614756    1471.4097348

R provides useful facilities for analyzing such matrices. For instance, we can ascertain that this transformation is not a Euclidean motion (that is, a composition of a uniform scaling, rotation, and translation). As an example, its Singular Value Decomposition is obtained via

m <- tail(a$coefficients, 2) # The scale-rotation-skew portion of `a`
a.svd <- svd(m)

The result has three parts: a rotation, two scale factors, and another rotation. The scale factors differ:

> a.svd$d
[1] 1480.859 1470.313

Therefore we can think of this transformation as including a rescaling factor that varies between 1470 and 1481, depending on orientation. You can work out the orientations from the SVD if you like by analyzing a.svd$v and a.svd$u, but in GIS this is rarely done: we just want to get a good fit. What the varying scale factors in the SVD tell us, though, is that your QGIS plugin may not be up to the job. Look for a solution that lets you specify the coefficients of the transformation matrix (a) directly.

share|improve this answer
    
thanks whuber, this explores the solution in detail (an advertisement for R?). I'll have a go at the best fit in a spreadsheet. –  Willy Apr 23 '12 at 21:57
add comment

Oops, you did provided two points. I think the transformation parameters might look like this:

Scale: x: 1470.214334304364 y: 1470.214334304364

Rotation: 5.920175054087165

Translation: x: 567735.57225 y: 6450618.3886

If you can make the affine plug-in do rotation (which I've never been able to do) it might be worth trying these. There definitely appears to be some rotation.

share|improve this answer
    
OK, so previously, as far as translation goes, we were talking Australia vs Antarctica. Now we are within 10 km. As for scale, could be very good. As for rotation, 6 deg too much. On the right track perhaps, what would the form of the equation be? –  Willy Apr 23 '12 at 12:27
    
lights out here.... –  Willy Apr 23 '12 at 12:31
    
Have a look at this, Willy. It includes my preferred solution to affine problems (which didn't seem to be rated too highly). gis.stackexchange.com/questions/21696/… –  nhopton Apr 23 '12 at 13:05
    
Oh, and the third point changes things. Scale: x: 1471.350969541187 y: 1470.2147334304364 Rotation: 6.156702051441602 Trans: x: 562282.2476 y: 6450853.8461 –  nhopton Apr 23 '12 at 13:07
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.