Sunday 1 July 2012

Terrain mapping with Mapnik

In the previous posts we’ve created three different GeoTIFF images from a single Digital Elevation Model (DEM); hillshade, slopeshade and color relief. It’s time to combine the images into a single terrain map of Jotunheimen.


Mapnik is a great open source toolkit for map rendering, and we’ll use the RasterSymbolizer to blend the GeoTIFFs. The styles and layers (images) can be defined in an XML file (jotunheimen_terrain.xml):

The map projection (UTM 32N) is defined in the first line. Mapnik uses the PROJ.4 library, and you’ll find the srs code by clicking on the Proj4 link on spatialreference.org. Mapnik can’t alter the projection of raster images, so the images have to be in the desired projection before feeding them into Mapnik.

Next, I’ve defined two styles for this map with the RasterSymbolizer. The hillshade is going to be blended with the color relief image, using the multiply blend mode. The color numbers for each pixel of the hillshade are multiplied with the corresponding pixel for the shaded relief image. I’m also setting the opacity of the hillshade to 0.6, making it less dominant. The map layers are defined below the styles. Mapnik will render the layers in order, starting from the top.

You can use a simple Python script (jotunheimen_terrain.py) to render the map image (you can also use tools like TileCache and MapProxy):

#!/usr/bin/env python
import mapnik
map = mapnik.Map(3134, 3134)
mapnik.load_map(map, 'jotunheimen_relief.xml')
map.zoom_all() 
mapnik.render_to_file(map, 'jotunheimen_relief.png')

This script will create this PNG image:

python jotunheimen_terrain.py

[ Download image ]

Let’s add the slopeshade image to our XML file (jotunheimen_terrain2.xml):

I’ve adjusted the opacity of the slope- and hillshade to create this blend:

[ Download image ]

Which terrain map you prefer is a matter of personal preference. The first map has more glow on the light facing slopes. In the last map, the steep mountains are more defined and the light facing slopes are more detailed.

On both maps, I miss two important features for Jotunheimen - lakes and glaciers. Adding land cover data is the topic of the next blog post.  

6 comments:

Jochen said...

Thanks for sharing your knowledge. This article has helped me to improve my maps on http://alpenkarte.eu and http://skitourenkarte.eu.

Have added a link to this article on my blog: http://alpenkarte.blogspot.de/

Unknown said...

Hi Bjorn,

Thanks for those precious tutorials.

I'm running into this error when running the python script :

File "jotunheimen_terrain.py", line 4, in
mapnik.load_map(map, 'jotunheimen_terrain.xml')
RuntimeError: Could not create datasource for type: 'gdal' encountered during parsing of layer 'color relief' in Layer at line 14 of 'jotunheimen_terrain.xml'

It seems that Mapnik (maybe not at the time you've written the blog post) doesn't have the gdal plugin installed by default. (https://github.com/mapnik/mapnik/wiki/PluginArchitecture)

But I installed Mapnik with brew, not from source, so I don't see how I could be able to run Scons to install plugins...

Is there something I missed? I'm a bit stuck here :(

Thanks for your help.

David McNamara said...

@Erik If you haven't figured this out yet, here is how you can install mapnik with gdal support:

$brew reinstall mapnik --with-gdal

More info:
https://github.com/mapnik/mapnik/wiki/MacInstallation_Homebrew

Muzaffer Tunç said...

Hi,
The generated image has no CRS.
Is there any way for generated tiff image to have CRS from the source geotiff image/s.

Muzaffer Tunç said...

I used gdalinfo to get source image's extensions then using GDAL_edit.py to copy it over mapnik image.
Also used gdalcopy.py to copy original image's projection.

Anonymous said...

Muzaffer Tunç,

Using gdal_edit I cannot put back the crs into the blended image, which does not have any after being created by Mapnik. Could you specify the parameters of gdal_edit that worked for you?.

So far, the only thing that works for me is by copying and paste the xml of any of the individual images and re-name it with the blended image´s name. Actually, if I perform gdal_translate over the blended image the output xml, the SRS line differs from those of the individual images, which may be the cause why is not working at all.