Print OpenStreetMap with Gnuplot

From OpenStreetMap

Software described on this page or in this section is unlikely to be compatible with API version 0.5 deployed on 8 October, 2007. If you have fixed the software, or concluded that this notice does not apply, remove it.

Background: Close to me (LA2), there are no more completely white areas that need to be mapped in OSM. In some areas, however, only the major roads have been covered and I need to go back there and fill in the smaller roads. I want to avoid the roads that are already done. It would be extremely handy to download all of OSM to a handheld unit (PDA, GPS receiver, StreetPilot), but lacking that I could use a way to generate a PDF (or PostScript) file for printing on an A4/legal size paper. That is roughly 2000x3000 pixels in 300 dpi or 4000x6000 pixels in 600 dpi, compared to the 700x500 slippy map on the OSM website, so printing a screenshot is not a useful solution.

The correct way could be to use the Osmarender program. Here we present an alternative solution.

The following is a quick-and-dirty hack, based on planet.osm, the monthly static database dump. It assumes you run Linux and are comfortable with the command line and shell scripts. It has been tested on a standard SuSE 9.1 installation on a 2 GHz CPU with 512 MB RAM, and Ubuntu 6.06 with gnuplot. It also assumes the planet.osm dump has a line-oriented formatting, which was true in April and May 2006.

1. Download the planet.osm dump and uncompress it, as specified on the page planet.osm. Make sure the resulting file is named planet.osm in the current directory.

2. Run this script. It might take some minutes, but you only have to do this once each month, after you download the new planet.osm dump. The resulting file OpenStreetMap is almost as big as the planet.osm file (for May 2006, planet.osm is 93 megabytes, and OpenStreetMap is 59 megabytes).

 #!/bin/sh
 
 <planet.osm tr \' \" |
 sed '/<node/!d;s/.*id="\([0-9]*\)" lat="\([-.0-9]*\)" lon="\([-.0-9]*\)".*/\1 \2 \3/' |
 sort -k1,1 >nodes
 
 <planet.osm tr \' \" |
 sed '/<segment/!d;s/.*id="\([0-9]*\)" from="\([0-9]*\)" to="\([0-9]*\)".*/\1 \2 \3/' |
 sort -k3,3 | join -1 3 -2 1 - nodes |
 sort -k3,3 | join -1 3 -2 1 - nodes |
 sort -k3,3 -n |
 awk '{printf "\n# segment %d from %d to %d\n%s %s\n%s %s\n", $3, $1, $2, $7, $6, $5, $4; }' >OpenStreetMap

The resulting file can be read by gnuplot (step 3 and 4). It is named OpenStreetMap because it will look nice in the default gnuplot legend. The file looks something like this:

 # segment 14 from 39 to 36
 -0.146549582481384 51.5254936218262
 -0.146119650312045 51.5254364014072
 
 # segment 19 from 42 to 40
 -0.147806495428085 51.5235748291016
 -0.145178407205123 51.5240177646579
 
 # segment 27 from 44 to 47
 -0.148549228906631 51.5251197814941
 -0.147991970210569 51.5252208698207

Gnuplot uses X and Y coordinates, which is why the longitude (0.14 degrees west of Greenwich) appears before the latitude (51 degrees north) on each line. If you specify more coordinates, each on a line of its own, gnuplot will draw lines from each point to the next. The blank lines cause gnuplot to lift the pen and not draw a line to the next point. Since this is a quick hack and not highly optimized, we instruct gnuplot to draw each line segment as a separate line. This solution ignores the new concept of "ways", and also all tags, names, and classes of roads.

3. For each area you want to plot, write a gnuplot command file like the one below. This example shows an area in Sweden, from 14.9 to 15.7 degrees east of Greenwich, and 58.3 to 58.6 degrees north of the equator. Use negative numbers if you are west of Greenwich or south of the equator. Next month, when you have a new database dump, you can use this plot file again. The last part sets the linewidth to half (lw 0.5) of the default. Depending on your printer's quality, you might want to try other values.

 set terminal postscript "Helvetica" 10
 set encoding iso_8859_1
 set grid
 set output "myplot.ps"
 set label "Fornåsa"   at 15.23, 58.48
 set label "Klockrike" at 15.34, 58.49
 plot [14.9:15.7] [58.3:58.6] "OpenStreetMap" with lines lw 0.5

The default font is Helvetica at 14 points, which is somewhat large for a map. You need to specify the encoding only if you set labels (placenames) on the map that contain non-ASCII characters (e.g. Fornåsa). Make sure the file you save uses the encoding you specify. Observe that X (longitude) comes before Y (latitude) when you specify coordinates for the labels.

4. Run these commands, supposing we called the file above myplot.gnuplot. It should only take a few seconds.

 gnuplot myplot.gnuplot
 ps2pdf myplot.ps myplot.pdf
 acroread myplot.pdf

5. You should be able to print the resulting map from acroread (or your favorite PDF viewer), as shown in the screenshot below. Or you can use the PDF file in e-mails or on the web.

Map projection

1. using proj

This uses proj, see Converting_to_WGS84 for pointers.

To use proj you need to filter your data first, that means instead of using the whole planet.osm file you need to limit your area. Filter out all lat/lon pairs that falls outside you map area. If you leave these ones out you will get very funky results.

To make this work I would have to figure out how to not scale the plot that gnuplot outputs. It doesn't matter how much accuracy I can reproject that data, it will still be squezed when printed on paper. ;-)

project to Swedish Grid RT90

mv OpenStreetMap OpenStreetMap.wgs84
cat OpenStreetMap.wgs84 |
proj  +ellps=WGS84 +proj=tmerc +lat_0=0 +lon_0=015d48.377m +k=1.0000056 +x_0=1500064.1 +y_0=-668.0|
sed 's/-.*//;'>OpenStreetMap.rt90;

2. without proj

This simple solution doesn't care about map projections, but assumes that latitudes and longitudes can be used as a rectangular grid on paper. This approximation gives negligible errors for small areas (such as a city), provided that the two dimensions are correctly scaled. But scaling can be hard enough. In central Scandinavia, since the cosine of 60 degrees is 0.5, each longitude degree is only half as long as a latitude degree. A short awk script can help to compute the coordinate ranges for the plot command, given a center coordinate and map scale. When printing on A4 paper, gnuplot uses an area of roughly 200 x 150 mm for the diagram. Since latitude degrees are equally long everywhere (1/90 of the distance between pole and equator, which is 10,000 km according to the original definition of the metre, i.e. 111.1 km), this is used as the basis for the scale. First the vertical distance (yd) from the center to the edge of the map is computed, then the horizontal distance (xd).

 # Enter center latitude, longitude and demoninator of the map scale on standard input,
 # e.g. for 58.45 N and 15.30 E in scale 1:250000
 echo 58.45 15.30 250000 |
 awk '{yd = $3 * 0.075 / 111111.1; xd = yd * 4.0 / 3.0 / cos($1 * 3.1415926 / 180.0);
      printf "plot [%.5f:%.5f] [%.5f:%.5f]\n", $2 - xd, $2 + xd, $1 - yd, $1 + yd;}'

For this example, the awk script prints "plot [14.86999:15.73001] [58.28125:58.61875]", which is quite close to the example above.