# How do I calculate where to plot a point on a document, given other reference points?

Imagine a PDF document. This document has been geo-referenced - I have been given two (x,y) coordinates that correspond to two (lng/lat) pairs.

Based on this data, I need to plot a different pair of coordinates on the image.

How can this be done? It is very straightforward if the document is north-up, but if the document it tilted, then it's more complicated. I think I need to calculate and transform a matrix. Any help?

EDIT:

Is there any chance I can give you a sample data, and you could demonstrate how to calculate the point? My intuition tells ne your answer is correct, but I'm not sure how to do the math to calculate the angle between the lines and decompose the matrix equation.

So, data I might ask you to consider is:

The points look like y x (lat lng)
The first point is 307.000000 267.000000 (37.725918 -122.223587)
The second point is 155.000000 400.000000 (37.619473 -122.373886)
The sample point is (37.7 -122.3)


Incidentally, this is for some software I am writing, not for homework.

-
This might be helpful. – Rasmus Jun 23 '11 at 20:55
The answer will depend on the map projection that is used, which could make a significant difference if it covers a large area. – mac Jun 23 '11 at 21:51
Thanks for pointing this out. In this case, the area is very small, and can be assumed to be flat, – Andrew Johnson Jun 23 '11 at 22:03

I'll suppose that the transformation from longitude/latitude $w=\begin{pmatrix}u\\v\end{pmatrix}$ to page coordinates $z=\begin{pmatrix}x\\y\end{pmatrix}$ is of the form $$z=\rho R_\theta w + z_0$$ where $R_\theta=\begin{pmatrix}\cos\theta&-\sin\theta\\ \sin\theta&\cos\theta\end{pmatrix}$ and $z_0=\begin{pmatrix}x_0\\y_0\end{pmatrix}$, for some constants $\rho,\theta,x_0,y_0$. This seems to be what you want, since this transformation is "rotation anticlockwise by $\theta$, then scaling by $\rho$, then translation by $z_0$.

Now suppose for $j=1,2$, we are given $w_j,z_j$ (the two lat/long pairs and their corresponding points) so that $z_j=\rho R_\theta w_j+z_0$.

We wish to compute $\rho, \theta,z_0$ using this data. Since $\rho$ is the scaling factor, we have $\rho=\mathrm{dist}(z_1,z_2)/\mathrm{dist}(w_1,w_2)$, where $\mathrm{dist}$ is the usual Euclidean distance. Similarly, we can compute the angle $\theta$ by calculating the angle that $w_2-w_1$ makes with $z_2-z_1$, and then $z_0=z_1-\rho R_\theta w_1$.

This gives you all of the constants in the general formula, which you can use to find the coordinates corresponding to other lat/long pairs.

Here is some Haskell code which hopefully implements the above scheme and calculates the answer for your example.

neg z = [-x,-y] where [x,y]=z
vplus a b = [a1+b1,a2+b2]
where [a1,a2]=a
[b1,b2]=b
vdiff a b = a vplus neg b
dotprod a b = a1*b1+a2*b2
where [a1,a2]=a
[b1,b2]=b
crossprod a b = a1*b2-a2*b1
where [a1,a2]=a
[b1,b2]=b
scalmult s w = [s*u,s*v] where [u,v]=w
mag z = sqrt $x*x+y*y where [x,y]=z dist a b = mag$ a vdiff b
unitvec z = (1/(mag z)) scalmult z
cos_sin a b = [dotprod a' b', crossprod a' b']
where [a',b']=[unitvec a,unitvec b]
anglefromcs cs = sign * acos c
where [c,s]=cs
sign
| s >= 0 = 1
| otherwise = -1
angle v w = anglefromcs $cos_sin v w rho w1 w2 z1 z2 = (dist z1 z2) / (dist w1 w2) theta w1 w2 z1 z2 = angle (w2 vdiff w1) (z2 vdiff z1) rot theta w = [cos theta * u - sin theta * v, sin theta * u + cos theta * v] where [u,v]=w z0 w1 w2 z1 z2 = z1 vdiff ( rho' scalmult (rot theta' w1 ) ) where theta' = theta w1 w2 z1 z2 rho' = rho w1 w2 z1 z2 z w1 w2 z1 z2 w = rho' scalmult (rot theta' w) vplus z0' where rho' = rho w1 w2 z1 z2 theta' = theta w1 w2 z1 z2 z0' = z0 w1 w2 z1 z2 main = print$  z [37.725918,-122.223587] [37.619473,-122.373886] [307,267] [155,400] [37.7,-122.3]


The output is [226.55468797299545,303.8562915341063].

-
I have edited my question a bit... I think your answer is right, but still beyond me. – Andrew Johnson Jun 23 '11 at 23:20
There should be some way to express thanks other than up votes. If you are ever in Berkeley, I will buy you a beer. – Andrew Johnson Jun 24 '11 at 1:13
I ended up translating to Obj-C, so I posted that as an answer as well. – Andrew Johnson Jun 24 '11 at 4:22
I don't follow the math nor the haskell well, but would just point out that this approach has to fall apart where the difference in x, y, lat, or long is or approaches 0; say where the georeferenced points have x and y values (100, 140) and (100, 160). then you need something like the haversine formula to determine the pixels-per-degree value for the "missing" dimension. – jcomeau_ictx Sep 22 '11 at 7:17

Since I translated Mac's code to Objective-C, I thought I would share it too. This works for me, read from the bottom up:

- (CGPoint) neg:(CGPoint)point {
return CGPointMake(-point.x,-point.y);
}

- (CGPoint) vplus:(CGPoint)a b:(CGPoint)b {
return CGPointMake(a.x+b.x, a.y+b.y);
}

- (CGPoint) vdiff:(CGPoint)a b:(CGPoint)b {
return [self vplus:a b:[self neg:b]];
}

- (float) dotprod:(CGPoint)a b:(CGPoint)b {
return a.x*b.x+a.y*b.y;
}

- (float) crossprod:(CGPoint)a b:(CGPoint)b {
return a.x*b.y - a.y*b.x;
}

- (CGPoint) scalmult:(float)s w:(CGPoint)w {
return CGPointMake(s * w.x, s * w.y);
}

- (float) mag:(CGPoint)z {
return sqrt(z.x*z.x + z.y*z.y);
}

- (float) dist:(CGPoint)a b:(CGPoint)b {
return [self mag:[self vdiff:a b:b]];
}

- (CGPoint) unitvec:(CGPoint)z {
return [self scalmult:1/[self mag:z] w:z];
}

- (CGPoint) cos_sin:(CGPoint)a b:(CGPoint)b {
a = [self unitvec:a];
b = [self unitvec:b];
return CGPointMake([self dotprod:a b:b], [self crossprod:a b:b]);
}

- (float) sign:(float)num {
if (num > 0) return 1;
return -1;
}

- (float) anglefromcs:(CGPoint)cs {
return [self sign:cs.y] * acos(cs.x);
}

- (float) angle:(CGPoint)v w:(CGPoint)w {
return [self anglefromcs:[self cos_sin:v b:w]];
}

- (float) rho:(CGPoint)w1 w2:(CGPoint)w2 z1:(CGPoint)z1 z2:(CGPoint)z2{
return [self dist:z1 b:z2] / [self dist:w1 b:w2];
}

- (float) theta:(CGPoint)w1 w2:(CGPoint)w2
z1:(CGPoint)z1 z2:(CGPoint)z2{
return [self angle:[self vdiff:w2 b:w1] w:[self vdiff:z2 b:z1]];
}

- (CGPoint) rot:(float)theta w:(CGPoint)w {
return CGPointMake(
cos(theta) * w.x - sin(theta)*w.y,
sin(theta) * w.x - cos(theta)*w.y
);
}

- (CGPoint) z0:(CGPoint)w1 w2:(CGPoint)w2 z1:(CGPoint)z1 z2:(CGPoint)z2{

float theta = [self theta:w1 w2:w2 z1:z1 z2:z2];
float rho = [self rho:w1 w2:w2 z1:z1 z2:z2];

CGPoint scalmult = [self scalmult:rho w:[self rot:theta w:w1]];

return [self vdiff:z1 b:scalmult];
}

- (CGPoint) z:(CGPoint)w1 w2:(CGPoint)w2 z1:(CGPoint)z1 z2:(CGPoint)z2 w:(CGPoint)w {

float theta = [self theta:w1 w2:w2 z1:z1 z2:z2];
float rho = [self rho:w1 w2:w2 z1:z1 z2:z2];

CGPoint scalmult = [self scalmult:rho w:[self rot:theta w:w]];
CGPoint z0 = [self z0:w1 w2:w2 z1:z1 z2:z2];
CGPoint vplus = [self vplus:scalmult b:z0];
return vplus;
}

-

Let $${\bf u}=\pmatrix{307\cr267\cr},{\bf v}=\pmatrix{37.7\cr-122.2\cr},{\bf w}=\pmatrix{155\cr400\cr},{\bf z}=\pmatrix{37.6\cr-122.3\cr},{\bf p}=\pmatrix{37.7\cr-122.3\cr}$$ be your first $(y,x)$ pair, your first ${\rm(lat,\ lng)}$ pair, your second $(y,x)$ pair, your second ${\rm(lat,\ lng)}$ pair, and your sample point (where I've truncated some of the numbers to 1 decimal because I'm too lazy to type them out). You want to find a $2\times2$ matrix $A$ such that $A{\bf v}={\bf u}$ and $A{\bf z}={\bf w}$; then your answer will be $A{\bf p}$. Let $A=\pmatrix{a&b\cr c&d\cr}$, where $a,b,c,d$ are unknown. What you want is $$\pmatrix{a&b\cr c&d\cr}\pmatrix{37.7&37.6\cr-122.2&-122.3\cr}=\pmatrix{307&155\cr267&400\cr}$$ where the columns of the first matrix are $\bf v$ and $\bf z$, and the columns of the matrix on the right are $\bf u$ and $\bf w$. Can you solve that equation for $a,b,c,d$? If so, then, as noted, the answer you're looking for is $A{\bf p}$.

-
The transformation for changing coordinates should be affine, not linear, since it need not preserve the origin. – mac Jun 24 '11 at 8:33
@mac, you made some assumptions for your answer, I made some (different) assumptions for mine. I think that shows the original question was underspecified, but it appears that OP knows what he wants, and what he wants is what you've done. – Gerry Myerson Jun 24 '11 at 13:07