Sunday, January 23, 2011

Mercator

There's a ship anchored in Oostende harbour named Mercator, after the Flemish cartographer Gerardus Mercator. I've seen it a few times and registered that there must be something vaguely interesting about it, but the light bulb in my brain never flicked on. It turns out that it did a lot of boring stuff (with the exception of bringing back two Moai from Easter Island) but the cartogopher after whom it was named is of much more interest to me lately: Google Maps uses a variant of the Mercator projection to produce a flat rectangular map of the Earth - an ellipsoid. The mathematics behind this projection can be seen at Wolfram.

I set myself a task. I have a Garmin GPSMap 62st: a great little device (if a little limited), but the utility I can derive from it is tightly coupled to the availability of maps. What's the point in knowing you're at 13°09′47″S 72°32′44″W if the map around is just empty black space? Armed with a 1024 x 1024 pixel map (the maximum tile size for the 62st) in Mercator projection, I isolated two points, in diagonally opposite corners, to get the most distance between them. Google Earth allowed me to drop markers and get the corresponding latitude and longitude co-ordinates for each pair of (x,y) pixels.

Degrees are easily converted into radians (and back again), and radians of latitude can easily be converted into the unitless measures of height used in the Mercator projection, while degrees of longitude don't need any conversion at all. If you look at a complete map of the earth, you'll see the areas surrounding the poles are most distorted. Anyway, long story short:

Given a pair of (N,E) coordinates and their respective (X,Y) positions on the raster map, it's pretty simple to work out what the left, top, right and bottom boundaries of the entire map are, and hence easy to create the KML and KMZ file required to correctly position a raster map (of the Mercator projection) on the 62st.

Longitude first because it's easiest: X2-X1 gives the pixel width between the two selected points on the raster map; X1-X0 gives the pixel padding to the left edge, and X3 - X2 gives the pixel padding to the right. We know the angle of longitude at E1 and E2, so we can work out a ratio of degrees to pixels using (X2-X1)/(E2-E1). Substituting the ratio into X3-X2 and X1-X0, and adding the correct starting point, we get the longitudinal boundaries.

Latitude has some tricky hyperbolic and trigonmetric functions (as seen on Wolfram): an angle of latitude N (in radians) is represented at height H (no units) above the equator. N1's height would be H1 = atan(sinh(N1)) and N2's height H2 = atan(sinh(N2)). The formula (Y2-Y1)/(N2-N1) gives us the ratio of pixels to "heights". It's simple to get to Y0 and Y3 from Y2 (and/or Y1) using just addition and subtraction. Also, we know Y1 and Y2's latitudes so it's simple to get the "height" of Y0 and Y3: H0 and H3. Getting back to radians from "height" is as simple as Y0 = log(tan(H0) + 1/cos(H0)). All that's left is to convert the answer from radians back into degrees. Simples.