Map Projections (with GNU Octave)

You have analyzed your GPS data to death and all that's left is plotting it on a world map. Thank god it's easy with GNU Octave. Since most of the Webservices (Google Maps, Bing Maps, Open Street Map) use a Mercator projection for their maps, we'll use it aswell. So according to xkcd I am not really into maps. If you'd like to know more about map projections and especially map tiling, read Troy Brant's Post on map projections with Objective-C.

Imagine you have to plot the following GPS positions:

on the map (taken from Wikipedia:File:Mercator-projection.jpg):

So first define the Mercator projection as:

function y = radians(x)
  y = x * pi / 180;
endfunction

function [x, y] = mercator(lat, lng, width, height)
  x = mod(width * (lng+180)/360, width) ;
  y = log(tan((radians(lat)/2) + (pi/4)));
  y = (height/2) - (width * y/(2*pi));
endfunction

… then load and plot the data. Note that the filesize of the resulting image is much larger than the original. If you want to serve this image from a webserver you should postprocess it.

I = imread("Mercator-projection.jpg");
[height, width, channels] = size(I);

data = [ 52.517, 13.4; # berlin
         51.473611 -0.159444; # london
         40.783333, -73.966667 # new york
         ];

fig = figure('visible','off');
imshow(I);
set(gca, 'units','normalized', 'Position', [0 0 1 1]);
set(gcf, 'units', 'pixels','position',[0 0 width height]);
hold on;

[x,y] = mercator(data(1,1), data(1,2), width, height);
plot(x,y, "go", "MarkerSize", 10, "linewidth", 3);

[x,y] = mercator(data(2,1), data(2,2), width, height);
plot(x,y, "ro", "MarkerSize", 10, "linewidth", 3);

[x,y] = mercator(data(3,1), data(3,2), width, height);
plot(x,y, "bo", "MarkerSize", 10, "linewidth", 3);

print("-dpng",sprintf("-S%d,%d",width,height),"mercator-proj.png")
replot;

and you are done: