Translating lat/lon coordinate on rotated/scaled maps

Question

I need to indicate a lat/lon coordinate in an old map like this: http://www.davidrumsey.com/luna/servlet/detail/RUMSEY~8~1~3017~90020003:Topographical-Map-Of-The-City-and-C

Projected on Google Maps: http://rumsey.geogarage.com/maps/g2784000.html

I know the bounding box of the map as well as the dimensions of the image in pixels, but I couldn't figure out how to locate a lat/lon within this map.

Upper left is (40.698291,-74.079051), Lower Left (40.659855,-73.979638), Lower Right (40.855232,-73.835582) and Upper Right (40.882919,-73.940263).

Is there a standard formula on how to do locate a lot/lon datum in the historic map?

MvG MvG · Accepted Answer · 2013-06-17T13:05:16

As long as you can assume your lat/lon coordinates to form a rectangular grid (which should be a reasonable approximation on a city scale, although it will fail for larger areas or very precise data), you can assume your transformation between coordinates to be a projective transformation. See this post (or this one if you prefer to stay on SO) for details on how to compute the corresponding transformation matrix and use that to convert coordinates.

Based on this post I did a bit of computation with sage:

def pm1(a, b, c, d):
    M = Matrix([a , b, c]).transpose()
    f = M.adjoint() * d
    return M*diagonal_matrix(f)
def pm(a1, a2, b1, b2, c1, c2, d1, d2):
    return pm1(a2, b2, c2, d2)*(pm1(a1, b1, c1, d1).adjoint())
P1 = vector(QQ, [0, 0, 1])
P2 = vector(QQ, [0, 1, 1])
P3 = vector(QQ, [1, 1, 1])
P4 = vector(QQ, [1, 0, 1])
Q1 = vector(QQ, [40698291, -74079051, 1000000]) # Upper left
Q2 = vector(QQ, [40659855, -73979638, 1000000]) # Lower left
Q3 = vector(QQ, [40855232, -73835582, 1000000]) # Lower right
Q4 = vector(QQ, [40882919, -73940263, 1000000]) # Upper right
M = pm(Q1, P1, Q2, P2, Q3, P3, Q4, P4)
M.change_ring(RDF)/1e40

The resulting formula reads like this:

z =   559.910562534*lat - 539.510073656*lon - 64123.7576703
x = (-5629.59680416*lat - 2176.56828347*lon + 67876.8560722)/z
y = ( 8466.61769096*lat - 11263.0392472*lon - 1178932.12938)/z

The resulting coordinates x and y are measuered from zero to one, where zero denotes the left resp. upper edge, and one the right resp. lower edge.

If the assumption of a rectangular grid is not enough, you'll need to know details about the projection used to create the map, so that you can relate map positions to spherical coordinates.

Translating lat/lon coordinate on rotated/scaled maps

1 Answers