From OpenStreetMap

In February 2008 contours were implemented on the cycle map (a Mapnik rendered Slippy Map). A similar approach was then used for OpenPisteMap, which was then copied by the cycle map in turn.

In both cases, a "toolchain" was set up, such that ongoing OSM data improvements are re-rendered along with the contour graphics at regular intervals. Rendering happens for the entire set of tile images (at several zoom levels), which makes it quite an intensive process in terms of CPU/memory and storage. This page describes the toolchain, and various other considerations.

Contents 1 Toolchain 1.1 Downloading data from NASA 1.2 The shapefiles approach 1.2.1 Preparing for gdal_contour 1.2.2 Use gdal_contour to generate shapefiles 1.2.3 Index the shapefiles 1.2.4 Generate the mapnik layers 1.3 The PostGIS approach 1.3.1 Importing the SRTM3 data 1.3.2 Setting up Mapnik 2 Other considerations 2.1 Colours, styles etc 2.2 Void filling 2.3 Other sources 2.4 Tools 2.5 Related 3 Storage requirements

Toolchain

In summary

• Download [srtm] data from NASA
• Prepare for gdal_contour
• Run gdal_contour to create contour lines
• Store the contours in either shapefiles or postgis
• Add layer definitions for mapnik

Downloading data from NASA

We're interested in srtm v2 SRTM3 from [1]. Version 1 of the srtm data had contours over the oceans - version 2 is clipped to land regions. SRTM3 covers the whole world between 60S and 60N - SRTM1 is higher resolution for the US but we're sticking with SRTM3 even in North America.

The data comes in zip files covering 1degree x 1degree sections of the planet at 3601x3601 pixels of spot-height data. So download each .hgt.zip file you are interested in. At the time of writing we had 248 zip files.

wget ftp://e0srp01u.ecs.nasa.gov/srtm/version2/SRTM3/Eurasia/N54W001.hgt.zip
etc...

The shapefiles approach

There are two approaches in use to render contours onto the map. The cycle map orignially stored the contours as shapefiles, as described here, but has since moved to storing the contours in the postgis database alongside the main osm data.

Preparing for gdal_contour

The .hgt tiffs aren't quite the right format for gdal_contour - mainly revolving around the projection. Given that the projection meta-data comes directly from the name of the zip file a script is run to create the correct ancillary files.

Grab srtm_generate_hdr.sh from http://mpa.itc.it/rs/srtm/ to run that over the .zip files.

Use gdal_contour to generate shapefiles

For the cycle map, we want to render the 10m, 50m and 100m contours separately, so for every .zip file we'll generate 3 shapefiles. The gdal_contour command is along the lines of:

gdal_contour -i 10 -snodata 32767 -a height N54W001.tif N54W001c10.shp

-i sets the contour interval. -snodata is needed to override the void detection, otherwise it thinks the voids are actually 32km high spikes (so you get very dense circular contours around the voids). -a height means that the height of each line will end up in the shapefile, otherwise it'll just be a collection of lines with no actual data.

The previous two steps can easily be scripted, as per process.sh:

#!/bin/bash

cd /home/osm/shape-resources/srtm/
for X in *.hgt.zip
do
    yes | /home/osm/andy/srtm/srtm_generate_hdr.sh $X
done

rm *.hgt #remove unzipped .hgt files, not needed

rm *.shp #remove old shapefiles so that gdal_contour doesn't barf
rm *.shx
rm *.dbf

for X in /home/osm/shape-resources/srtm/*.tif
do
    echo $X
    gdal_contour -i 10 -snodata 32767 -a height $X ${X%%.tif}c10.shp
    gdal_contour -i 50 -snodata 32767 -a height $X ${X%%.tif}c50.shp
    gdal_contour -i 100 -snodata 32767 -a height $X ${X%%.tif}c100.shp
done

Index the shapefiles

The shapefiles can be given a spatial index, which significantly speeds up the rendering (since mapnik can check the index and ignore shapefiles that are outwith the bounds of the particular tile it's rendering. If you've got mapnik installed, you'll have the shapeindex utily, and you'll be able to

index.sh:

#!/bin/bash

for X in /home/osm/shape-resources/srtm/*.shp
do
   shapeindex ${X%%.shp}
done

Generate the mapnik layers

You'll need a layer definition for every shapefile. Given that we're going to have different styles for each set of contours (10, 50, 100m) we need to correlate them. So the two styles (for lines and text) are suffixed by the set number. The y iterator simply makes sure the layers have unique names.

generatexml.py:

#!/usr/bin/python
import glob

f = open('../osm-cycle/contour-layers-c100.include', 'w')
y = 0
for X in glob.glob('/home/osm/shape-resources/srtm/*c100.shp'):
  f.write("<Layer name=\"srtm100"+`y`+"\" status=\"on\" srs=\"+proj=latlong +datum=WGS84\">\n")
  f.write("  <StyleName>contours100</StyleName>\n")
  f.write("  <StyleName>contours-text100</StyleName>\n")
  f.write("  <Datasource>\n")
  f.write("    <Parameter name=\"type\">shape</Parameter>\n")
  f.write('    <Parameter name="file">' + X[0:-4] + "</Parameter>\n")
  f.write("  </Datasource>\n")
  f.write("</Layer>\n")
  y = y + 1
f.close

f = open('../osm-cycle/contour-layers-c50.include', 'w')
y = 0
for X in glob.glob('/home/osm/shape-resources/srtm/*c50.shp'):
  f.write("<Layer name=\"srtm50"+`y`+"\" status=\"on\" srs=\"+proj=latlong +datum=WGS84\">\n")
  f.write("  <StyleName>contours50</StyleName>\n")
  f.write("  <StyleName>contours-text50</StyleName>\n")
  f.write("  <Datasource>\n")
  f.write("    <Parameter name=\"type\">shape</Parameter>\n")
  f.write('    <Parameter name="file">' + X[0:-4] + "</Parameter>\n")
  f.write("  </Datasource>\n")
  f.write("</Layer>\n")
  y = y + 1
f.close

f = open('../osm-cycle/contour-layers-c10.include', 'w')
y = 0
for X in glob.glob('/home/osm/shape-resources/srtm/*c10.shp'):
  f.write("<Layer name=\"srtm10"+`y`+"\" status=\"on\" srs=\"+proj=latlong +datum=WGS84\">\n")
  f.write("  <StyleName>contours10</StyleName>\n")
  f.write("  <StyleName>contours-text10</StyleName>\n")
  f.write("  <Datasource>\n")
  f.write("    <Parameter name=\"type\">shape</Parameter>\n")
  f.write('    <Parameter name="file">' + X[0:-4] + "</Parameter>\n")
  f.write("  </Datasource>\n")
  f.write("</Layer>\n")
  y = y + 1
f.close

Take the three include files and put them into the osm.xml file at whichever point you want the contours to show up - usually just after the coastline shapefiles.

The PostGIS approach

The alternative approach, developed by OpenPisteMap, is to store the contours data in PostGIS. This only requires one set of layers to be set up for the whole data set, rather than one set per shapefile.

Importing the SRTM3 data

Put the zip files from NASA in a single directory, make sure you have srtm_generate_hdr.sh from http://mpa.itc.it/rs/srtm/ in your path and run this import script:

#!/bin/bash

PREP_TABLE="1"
for X in *.hgt.zip; do
	yes | srtm_generate_hdr.sh $X
	rm -f "${X%%.zip}"
	
	# Import 10m contours
	rm -f "${X%%.hgt.zip}.shp" "${X%%.hgt.zip}.shx" "${X%%.hgt.zip}.dbf"
	gdal_contour -i 10 -snodata 32767 -a height "${X%%.hgt.zip}.tif" "${X%%.hgt.zip}.shp"
	[ "$PREP_TABLE" ] && shp2pgsql -p -I -g way "${X%%.hgt.zip}" contours | psql -q gis
	shp2pgsql -a -g way "${X%%.hgt.zip}" contours | psql -q gis
	
	rm -f "${X%%.hgt.zip}.shp" "${X%%.hgt.zip}.shx" "${X%%.hgt.zip}.dbf"
	rm -f "${X%%.hgt.zip}.bil"
	rm -f "${X%%.hgt.zip}.hdr"
	rm -f "${X%%.hgt.zip}.prj"
	rm -f "${X%%.hgt.zip}.tif"
	unset PREP_TABLE
done

This will import all of the contours data for 10m interval contours into a PostGIS table called contours.

Setting up Mapnik

You need to add layers for the contour intervals you are interested in:

<Layer name="srtm_10" status="on" srs="+proj=latlong +datum=WGS84">
        <StyleName>contours10</StyleName>
        <StyleName>contours-text10</StyleName>
        <Datasource>
                <Parameter name="type">postgis</Parameter>
                <Parameter name="host"></Parameter>
                <Parameter name="port"></Parameter>
                <Parameter name="user"></Parameter>
                <Parameter name="password"></Parameter>
                <Parameter name="dbname">gis</Parameter>
                <Parameter name="estimate_extent">false</Parameter>
                <Parameter name="table">(select way,height from contours WHERE height::integer % 10 = 0 AND height::integer % 50 != 0 AND height::integer % 100 != 0) as "contours-10"</Parameter>
                <Parameter name="extent">-180,-90,180,89.99</Parameter>
        </Datasource>
</Layer>
<Layer name="srtm_50" status="on" srs="+proj=latlong +datum=WGS84">
        <StyleName>contours50</StyleName>
        <StyleName>contours-text50</StyleName>
        <Datasource>
                <Parameter name="type">postgis</Parameter>
                <Parameter name="host"></Parameter>
                <Parameter name="port"></Parameter>
                <Parameter name="user"></Parameter>
                <Parameter name="password"></Parameter>
                <Parameter name="dbname">gis</Parameter>
                <Parameter name="estimate_extent">false</Parameter>
                <Parameter name="table">(select way,height from contours WHERE height::integer % 50 = 0 AND height::integer % 100 != 0) as "contours-50"</Parameter>
                <Parameter name="extent">-180,-90,180,89.99</Parameter>
        </Datasource>
</Layer>
<Layer name="srtm_100" status="on" srs="+proj=latlong +datum=WGS84">
        <StyleName>contours100</StyleName>
        <StyleName>contours-text100</StyleName>
        <Datasource>
                <Parameter name="type">postgis</Parameter>
                <Parameter name="host"></Parameter>
                <Parameter name="port"></Parameter>
                <Parameter name="user"></Parameter>
                <Parameter name="password"></Parameter>
                <Parameter name="dbname">gis</Parameter>
                <Parameter name="estimate_extent">false</Parameter>
                <Parameter name="table">(select way,height from contours WHERE height::integer % 100 = 0) as "contours-100"</Parameter>
                <Parameter name="extent">-180,-90,180,89.99</Parameter>
        </Datasource>
</Layer>

And you need some styles:

<Style name="contours10">
        <Rule>
                <MaxScaleDenominator>51185</MaxScaleDenominator>
                <MinScaleDenominator>1599</MinScaleDenominator>
                <LineSymbolizer>
                        <CssParameter name="stroke">#9cb197</CssParameter>
                        <CssParameter name="stroke-width">0.5</CssParameter>
                </LineSymbolizer>
        </Rule>
</Style>
<Style name="contours50">
        <Rule>
                <MaxScaleDenominator>204741</MaxScaleDenominator>
                <MinScaleDenominator>51185</MinScaleDenominator>
                <LineSymbolizer>
                        <CssParameter name="stroke">#9cb197</CssParameter>
                        <CssParameter name="stroke-width">0.6</CssParameter>
                </LineSymbolizer>
        </Rule>
        <Rule>
                <MaxScaleDenominator>51185</MaxScaleDenominator>
                <MinScaleDenominator>1599</MinScaleDenominator>
                <LineSymbolizer>
                        <CssParameter name="stroke">#747b90</CssParameter>
                        <CssParameter name="stroke-width">0.6</CssParameter>
                </LineSymbolizer>
        </Rule>
</Style>
<Style name="contours100">
        <Rule>
                <MaxScaleDenominator>409483</MaxScaleDenominator>
                <MinScaleDenominator>204741</MinScaleDenominator>
                <LineSymbolizer>
                        <CssParameter name="stroke">#9cb197</CssParameter>
                        <CssParameter name="stroke-width">0.7</CssParameter>
                </LineSymbolizer>
        </Rule>
        <Rule>
                <MaxScaleDenominator>204741</MaxScaleDenominator>
                <MinScaleDenominator>51185</MinScaleDenominator>
                <LineSymbolizer>
                        <CssParameter name="stroke">#747b90</CssParameter>
                        <CssParameter name="stroke-width">0.7</CssParameter>
                </LineSymbolizer>
        </Rule>
        <Rule>
                <MaxScaleDenominator>51185</MaxScaleDenominator>
                <MinScaleDenominator>1599</MinScaleDenominator>
                <LineSymbolizer>
                        <CssParameter name="stroke">#855d62</CssParameter>
                        <CssParameter name="stroke-width">0.7</CssParameter>
                </LineSymbolizer>
        </Rule>
</Style>
<Style name="contours-text50">
        <Rule>
                <MaxScaleDenominator>51185</MaxScaleDenominator>
                <MinScaleDenominator>1599</MinScaleDenominator>
                <TextSymbolizer name="height" face_name="DejaVu Sans Book" size="8" fill="#747b90" halo_radius="1" placement="line" />
        </Rule>
</Style>
<Style name="contours-text100">
        <Rule>
                <MaxScaleDenominator>102370</MaxScaleDenominator>
                <MinScaleDenominator>51185</MinScaleDenominator>
                <TextSymbolizer name="height" face_name="DejaVu Sans Book" size="8" fill="#747b90" halo_radius="1" placement="line" />
        </Rule>
        <Rule>
                <MaxScaleDenominator>51185</MaxScaleDenominator>
                <MinScaleDenominator>1599</MinScaleDenominator>
                <TextSymbolizer name="height" face_name="DejaVu Sans Book" size="8" fill="#855d62" halo_radius="1" placement="line" />
        </Rule>
</Style>

Other considerations

Colours, styles etc

If you are using the shapefiles approach, one thing to consider when you are coming up with your styles is that the 100m contours will be rendered 3 times over - the 100m contours are in the 10, 50 and 100m shapefiles. So it's a good idea to make them slightly wider than the 50s, which are slightly wider than the 10s. It doesn't matter if they are all the same colour, but when they are different colours they will go muddy due to anti-aliasing. The PostGIS approach avoids this problem since the select statements in the layers will exclude overlapping contours (i.e. the 10m layer won't contain contours at 50m and 100m heights and the 50m layer won't contain contours at 100m heights).

Void filling

Although the voids aren't rendered, they still exist as gaps in the contours. Although single-pixel gaps and other small voids could be filled by extrapolation, there are larger areas that need other sources of data.

Other sources

Given the voids, and the lattitude limits on the data, future work will include using the lower-resolution SRTM30 dataset to fill in the gaps.

Tools

See Srtm2Osm to convert srtm data to osm files.

Related

HikingBikingMaps shows how to render images with hill shading in addition to contour lines.

Storage requirements

Home are some numbers from the mailinglist on disk space usage for those interested. I had a go at importing 10 metre contour lines for the whole of Eurasia into PostGIS - latitudes of 0 - 46 degrees North required about 110 gig of disk space for the Postgres table and amounted to around 105 million contour lines. (At this point before I ran out of space) - the SRTM3 data set extends up to 60 degrees North).

Between 0 and 46 degrees north across Eurasia amounts to 3244 1x1 degree tiles. So this averages you around 35MB of disk space to import a 1x1 degree tile into PostGIS (obviously dependent on the terrain the tile covers), giving rough estimated numbers of:

Location	est. storage req.	number of tiles
Africa	111 GB	(3250 tiles)
Australia	36 GB	(1060 tiles)
Eurasia	202 GB	(5902 tiles)
Islands	5 GB	(141 tiles)
N America	82 GB	(2412 tiles)
S America	62 GB	(1807 tiles)
WHOLE WORLD	498 GB	(14572 tiles)

So with half a terabyte of disk you can import the whole lot... There is also the higher resolution SRTM1 data set covering North America - I'm not clear on how using those data would affect these numbers - probably not much since you'll probably have roughly the same number of contour lines, they will just be positioned more accurately.