Ex Nihilo

A Processing Sketch Blog

Archive for the ‘Maps’ Category

GPS to Mercator Projection

Thursday, November 19th, 2009

A geography issue I’ve been struggling with over the past little while is correcting my assumptions about how latitude / longitude coordinates work. When I originally built Elevation I decided to take a shortcut and ignore Earth curvature, treating my GPS points as a Cartesian grid.

And it worked, sort of. It gave me land forms that come close to representing the actual terrain. In hindsight it’s pretty obvious how distorted they were, but for the first few iterations of the program they were good enough.

I’ve now started looking into representing paths with more accuracy. At the moment Elevation’s speed scale is totally relative to the data, but I’d like to peg actual km/h values to it. I’d also like to indicate the scale of the map with sliding meter/kilometer markers that adjust based on zoom level. As I started going down this road I quickly realized the math needed a closer look.

Without a clue where to begin, I turned to Wikipedia and various open source Java/Javascript libraries to see if they could offer a clue about converting my GPS coordinates into a proper grid. Long story short, what I needed was a map projection system.

Mercator is probably the most familiar projection out there, and it’s the same one Google, Yahoo, Microsoft etc. use for their online mapping services. While it seriously distorts extreme high and low latitudes (Greenland isn’t really the same size as Africa, after all), it has the advantage of treating smaller local areas more uniformly. Mercator won’t work at the North Pole or in the Antarctic, but at a regional level like city or state/province it’s a fairly uniform distortion so the maps just look right; since Elevation is intended for those smaller areas, that’s perfect for my needs.

So, how does one go about converting GPS lat/long coordinates to actual meters? That’s where the math comes in, and Wikipedia handily lists a few formulas for it. Processing doesn’t have sinh, cosh, or sec functions, so only the first two functions for y will work; I chose the second:

x = λ - λo

y = 0.5 * ln * ((1 + sin(φ)) / (1 - sin(φ)))

The x value basically ends up being your raw longitude coordinate, though it’s handy to offset it from the central x axis of your map. The y value requires a bit more heavy lifting. In Processing, the above functions end up looking something like this:

float lambda = gpsLong - (mapWidth / 2);

float phi = radians(gpsLat);
float adjustedPhi = degrees(0.5 * log((1 + sin(phi)) / (1 - sin(phi))));

While Wikipedia doesn’t explicitly say anything about radians, I found it necessary to convert y before the math would work in Processing, then convert the result back to degrees.

The resulting values still don’t represent real units on the surface of the globe; to get that you have to multiply each by the degree length, a variable value that corresponds to the surface distance between a single degree change in latitude. At the equator there are approximately 111km between degrees, but this number shrinks the closer you get to the poles.

Solving this problem is something I’m currently working on. An unfortunate red herring was that the degree length in Vancouver appears to be almost exactly 1.6 times less than the value at the equator, which looked for all the world like a flubbed imperial/metric conversion. If I were to start with the equator distance and divide by 1.6:

x = (adjustedPhi * 111319.9) / 1.6;

y =(lambda * 111319.9) / 1.6;

I would get values that look fairly darn close to accurate in Vancouver. I thought that had solved it, but I now suspect the math at different latitudes is still very much wrong. I’ll update this post once I know more.

Update: Thanks Microsoft, for publishing this piece on MSDN about how your Bing Maps tile system works. This bit in particular seems like it contains what I’m looking for:

The ground resolution indicates the distance on the ground that’s represented by a single pixel in the map. For example, at a ground resolution of 10 meters/pixel, each pixel represents a ground distance of 10 meters. The ground resolution varies depending on the level of detail and the latitude at which it’s measured. Using an earth radius of 6378137 meters, the ground resolution (in meters per pixel) can be calculated as:

ground resolution = cos(latitude * pi/180) * earth circumference / map width

= (cos(latitude * pi/180) * 2 * pi * 6378137 meters) / (256 * 2^level pixels)

The map scale indicates the ratio between map distance and ground distance, when measured in the same units. For instance, at a map scale of 1 : 100,000, each inch on the map represents a ground distance of 100,000 inches. Like the ground resolution, the map scale varies with the level of detail and the latitude of measurement. It can be calculated from the ground resolution as follows, given the screen resolution in dots per inch, typically 96 dpi:

map scale = 1 : ground resolution * screen dpi / 0.0254 meters/inch

= 1 : (cos(latitude * pi/180) * 2 * pi * 6378137 * screen dpi)
         / (256 * 2^level * 0.0254)

Posted in Maps, Projects | 1 Comment »