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.
 A: 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].
A: 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;
}

A: 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}$. 
