Get map position when WGS-84 lat/lon when upper left and lower right corners' lat/lon are given

Suppose I have a map, for example from openstreetmaps.org. I know the WGS-84 lat/lon of the upper left and lower right corner of the map. How can I find other positions on the map from given WGS-84 lat/lon coordinates?

Solution

If the map is roughly street/city level, uses a mercator projection (as openstreetmap.org seems to), and isn't too close to the poles, linear interpolation may be accurate enough. Assuming the following:

TL = lat/lon of top left corner
BR = lat/lon of bottom right corner
P = lat/lon of the point you want to locate on the map
(w,h) = width and height of the map you have (pixels?)
the origin of the map image, (0,0), is at its top-left corner

, we could interpolate the (x,y) position corresponding to P as:

x = w * (P.lon - TL.lon) / (BR.lon - TL.lon)
y = h * (P.lat - TL.lat) / (BR.lat - TL.lat)

Common gotcha's:

The lat/lon notation convention lists the latitude first and the longitude second, i.e. "vertical" before "horizontal". This is opposite to the common x,y notation of image coordinates.
Latitude values increase when going in a north-ward direction ("up"), whereas y coordinates in your map image may be increasing when doing down.
If the map covers a larger area, linear interpolation will not be as accurate for latitudes. For a map that spans one degree of latitude and is in the earth's habitable zones (e.g. the bay area), the center latitude will be off by 0.2% or so, which is likely to by less than a pixel (depending on size)

If that's precise enough for your needs, you can stop here!

The more precise math for getting from P's latitude to a pixel y position would start with the mercator math. We know that for a latitude P.lat, the Y position on a projection starting at the equator would be as follows (I'll use a capital Y as unlike the y value we're looking for, Y starts at the equator and increases towards the north):

Y = k * ln((1 + sin(P.lat)) / (1 - sin(P.lat)))

The constant k depends on the vertical scaling of the map, which we may not know. Luckily, it can be deduced observing that y(TL) - y(BR) = h. That gets us:

k = h / (ln((1 + sin(TL.lat)) / (1 - sin(TL.lat))) - ln((1 + sin(BR.lat)) / (1 - sin(BR.lat))))

(yikes! that's four levels of brackets!) With k known, we now have the formula to find out the Y position of any latitude. We just need to correct for: (1) our y value starts at TL.lat, not the equator, and (2) y grows towards the south, rather than to the north. This gets us:

Y(TL.lat) = k * ln((1 + sin(TL.lat)) / (1 - sin(TL.lat)))
Y(P.lat)  = k * ln((1 + sin(P.lat )) / (1 - sin(P.lat )))
y(P.lat)  = -(Y(P.lat) - Y(TL.lat))

So this gets you:

x = w * (P.lon - TL.lon) / (BR.lon - TL.lon) // like before
y = -(Y(P.lat) - Y(TL.lat))                  // where Y(anything) depends just on h, TL.lat and BR.lat