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

Question

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?

Answer 1

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

Answer 2

Oren, your answer is perfect, but somehow the calculations go wrong on low zoom levels, the y part calculates incorrectly.