Hack 29. Plot Arbitrary Points on a World Map

Of the dozens of common cartographic projections used to depict the entire world, the Equidistant Cylindrical projectionalso known as the rectangular projectionis by far the most hackable, making it easy to plot points wherever you want.

Since our lovely world is round but our maps are not, the would-be map hacker has a choice to make when selecting a cartographic projection for a given map. The range of available choices are described in detail in [Hack #28] . Sometimes, the best projection, however, is no projection at all. The Equidistant Cylindrical projection is made by treating every degree of latitude and longitude as equal in length, resulting in a perfectly rectangular map with a 2:1 width-to-height ratio. First attributed to Eratosthenes of Cyrene in the third century BC, this projection does not preserve area, shape, direction, distance, or scale (except along meridians), but it has one supreme advantage. It's easy to hack!

Here's why: since meridians and parallels in this projection meet at equally spaced right angles, latitude and longitude can simply be treated as rectangular coordinates, shortcutting the tortuous spherical trigonometry that usually plagues the art of cartographic projection. By consequence, the Equidistant Cylindrical projection is also referred to as equirectangular, or simply rectangular. (You'll also see it referred to as the Plate Carrée projection, which is French for "square plane.") Naturally, some cartographers think it's cheating to simply set x = R and y = R, and so maps of this kind are also referred to as unprojected.

3.9.1. Finding Base Maps on the Web

Because the Equidistant Cylindrical projection is so hackable, lots of interesting mapping apps like xplanet [Hack 46] and worldKit [Hack #39] rely on it, making good base maps easy to find on the Web. For example, NASA's beautiful Blue Marble imagery is provided in this form. You can find free images on NASA's Earth Observatory at http://earthobservatory.nasa.gov/Newsroom/BlueMarble/BlueMarble.html. Another great source of world maps in this projection is the flatplanet maps catalog, which can be found at http://flatplanet.sourceforge.net/maps/. We couldn't find a rectangular projection base map with political borders that we liked, so we made one of our own, using the Blue Marble imagery, GMT, and the GIMP [Hack #32] . You can download our Blue Marble, political-borders base map at http://mappinghacks.com/maps/.

3.9.2. The Code

Here's a sample Perl script called plot_points.pl that uses Arnar Hrafnkelsson's excellent Imager module to plot locations on a rectangular-projection world map. You can get Imager from the CPAN at http://search.cpan.org/~addi/Imager/:

#!/usr/bin/perl
 
use Imager;
use strict;
 
my ($in, $out) = @ARGV;
die "Usage: $0  
" unless $in and $out;
 
my ($top, $left) = ( 90, -180 );
my ($bottom, $right) = ( -90, 180 );
 
my $map = Imager->new;
$map->open( file => $in ) or die $map->errstr;
 
my $x_scale = $map->getwidth / ($right - $left);
my $y_scale = $map->getheight / ($bottom - $top);
 
while () {
 my ($lat, $lon) = /(-?d+.?d*)/gos;
 next if $lon > $right or $lon < $left or $lat < $bottom or $lat > 
$top;
 
 my $x = ($lon - $left) * $x_scale;
 my $y = ($lat - $top) * $y_scale;
 
 $map->circle( color => "white", r => 5, x => $x, y => $y );
 $map->circle( color => "red", r => 4, x => $x, y => $y );
}
 
$map->write( file => $out ) or die $map->errstr;

plot_points.pl takes two arguments on the command line: a filename to read the world map from, and another filename to write the map plus the points into. The script hinges on the fact that, in an Equidistant Cylindrical projection, each pixel will take up a fixed number of degrees of both latitude and longitude, no matter where it lies on the map. (This is not the case with other projections, of course.) So we start by getting the image width and height, and then calculate the scale in pixels per degree for both the x- and y-axes.

Although, by definition, the x- and y-scales should be the same, calculating them separately is only a little extra work, and it guarantees that our script will still work on rectangular maps that have been squashed or stretched consistently in one dimension or the other.

plot_points.pl then tries to read coordinates from its standard input, one pair per line, latitude first, in decimal degree notation, and does a little bounds checking. Then, for each axis, it calculates the distance to the edge of the map in degrees and converts this value to pixels by multiplying by the scale values calculated earlier. And that's all the heavy lifting there is! A red circle is drawn inside a slightly larger white circle to obtain a nice border effect (though there are actually better ways of doing this in Imager), and then, after all the points are read, the whole thing is written to disk.

Note that there's nothing particularly Perl-ish about our script. The same technique should be easy to implement in any language. In fact, in [Hack #41], we'll show off the exact same technique in JavaScript, with a few added wrinkles.

3.9.3. Running the Code

Here's a sample run of plot_points.pl, run on the Pathfinder image from the flatplanet gallery and fed points for London, New York, and Tokyo. You can hit Ctrl-D after you're done entering points by hand. Figure 3-29 shows the result:

$ perl plot_points.pl PathfinderMap.jpg three_cities.jpg
51 0
40 -74
35 140

Figure 3-29. London, New York, and Tokyo, plotted on an Equidistant Cylindrical world map

Of course, if you have a file full of points, you can just use shell redirection to dump the list right into the script. If the file contains other things, such as labels for each point, that's okay, so long as the first two decimal numbers on each line are the latitude and longitude, respectively. For example, suppose you have a tab-separated file of the world's 7,000 largest cities, sorted by population, in a file called cities.txt:

place name country current population latitude longitude
Mumbai India 12383146 18.96
 72.82

Buenos Aires Argentina 12116379 -34.61
 -58.37

"Karāchi" Pakistan 10537226 24.86
 67.01

Manila Philippines 10232924 14.62
 120.97

Dilli India 10203687 28.67
 77.21

...

If you wanted to plot just the 50 largest cities, you couldn't feed it directly into plot_points.pl, because (a) you'd get all 7,000 that way and (b) the population figure appears before the latitude and longitude in each line. But you could feed it through the Unix head(1) and cut(1) commands, and pipe that into plot_points.pl:

$ head -50 cities.txt | cut -f4,5 | perl plot_points.pl PathfinderMap.jpg 
biggest.jpg

If you don't have such a file of your own, you can get ours from http://mappinghacks.com/data/cities.txt. The map thus generated looks something like Figure 3-30. Isn't the Unix shell neat?

Figure 3-30. The world's 50 largest cities, plotted on an Equidistant Cylindrical map

 

3.9.4. Why Use Any Other Projection, if This One Is So Simple?

Well, sadly, even though the Equidistant Cylindrical projection looks fine at a global scale, and would even look all right for places near the equator, where lines of latitude and longitude really are about the same length, all that stuff near the Poles is pretty badly distorted. On the maps we've been looking at, that doesn't really matter much, but the Equidistant Cylindrical projection definitely won't make for such great maps of smaller regions farther away from the equator.

The sinusoidal projection offers one alternative. Instead of treating latitude and longitude as equal, distances between two lines of longitude are scaled to the cosine of the latitude, as would be depicted on a globe. Instead of using the following code to determine the x-coordinate:

 my $x = ($lon - $left) * $x_scale;

we could use this calculation to plot a sinusoidal projection of our data:

 my $x = ($lon - $left) * $x_scale * cos( $lat );

Although this distorts the shapes of large areas, it's easy to calculate, and it preserves the relationship of latitude and longitude better at temperate latitudes. For this reason, the TIGER Map Server makes use of the sinusoidal projection, as discussed in [Hack #14] .

We examine still other cartographic projections in [Hack #28] . Still, no other projection beats the Equidistant Cylindrical projection for ease of use, and there's absolutely something to be said for that.

3.9.5. See Also

• [Hack #39]
• [Hack #28]
• [Hack #14]
• [Hack #41]
• http://flatplanet.sourceforge.net/maps/
• http://mappinghacks.com/maps/