Mapnik/Hillshading using Mapnik, GDAL and SRMT data
_-_The_Noun_Project.svg){kind=small}
The given reason is: This guide is out of date and relies on tools whose sources are gone (domain grabbed). Adding raster layers like hill shade is documented in the Mapnik documentation, their generation (the hill shade itself) is covered elsewhere on the web in a better way. --Nakaner (talk) 14:50, 13 August 2021 (UTC).
Since I have been doing a decent amount of piste mapping this winter, I wanted to see my results rendered in a visually appealing way. For skiers, it is important to know whether you are going uphill or downhill, hence a piste map needs elevation data. Personally, I find hillshading more intuitive than contours - at the cost of absolute height information and, to an extent, precision, but I can live with that. Here's how I did it.
Prerequisites
My original setup included the following:
- Ubuntu 10.10
- Mapnik - it is available through the Ubuntu package manager, and even the Mapnik viewer is included, which is very convenient for quick map debugging. I use an OSM file as the data source, which saves me the effort of maintaining a PostGIS database.
I assume you are already familiar with Mapnik and have managed to render a simple map. Also, I assume you are using Ubuntu - things should look pretty much the same on other Linux distributions and even on other Unices, but if you decide to try this on a different OS, be aware that you'll be a pioneer...
Setup
We will need a few extra tools:
- GDAL (available through the Ubuntu package manager, you need gdal-bin and libgdal1-dev)
- srtm_generate_header.sh, available at [1]
- PerryGeo (GDAL-based DEM utilies, available at http://www.perrygeo.net/wordpress/?p=7 (now dead link)
This goes for Ubuntu 10.10, which includes GDAL 1.6.0. GDAL 1.7.0 includes all of the above functionality, so the other two tools will not be needed any more.
srtm_generate_header.sh
srtm_generate_header.sh just needs to be downloaded to a convenient location (~/bin is a good idea). If you are going to use SRTM3 data (or anything similar with 3" resolution), you can use that script as it is.
If you are planning to use SRTM1 data, with a resolution of 1", you will need to modify some hard-coded values in the script. I recommend creating a copy named srtm1_generate_header.sh and editing that - you can thus process both SRTM1 and SRTM3 data, depending on which you can get.
At line 88, you will find the following code:
echo "BYTEORDER M LAYOUT BIL NROWS 1201 NCOLS 1201 NBANDS 1 NBITS 16 BANDROWBYTES 2402 TOTALROWBYTES 2402 BANDGAPBYTES 0 NODATA -32768 ULXMAP $ULXMAP ULYMAP $ULYMAP XDIM 0.000833333333333 YDIM 0.000833333333333"> $TILE.hdr
Change that to
echo "BYTEORDER M LAYOUT BIL NROWS 3601 NCOLS 3601 NBANDS 1 NBITS 16 BANDROWBYTES 7202 TOTALROWBYTES 7202 BANDGAPBYTES 0 NODATA -32768 ULXMAP $ULXMAP ULYMAP $ULYMAP XDIM 0.000277777777778 YDIM 0.000277777777778"> $TILE.hdr
This sets the number of rows and columns to three times the value (3600 instead of 1200) and changes the number of bytes per row accordingly. Conversely, the X/Y resolution is divided by three.
PerryGeo
The PerryGeo utilities need to be built by hand. We just need one file, hillshading.cpp, from the sources - download it from SVN or extract it from the archive. The build instructions on the home page are outdated and do not work on Ubuntu 10.10 - we will need to adapt some parameters. Simply go to the directory into which you downloaded the source file and run
g++ hillshade.cpp `gdal-config --cflags` `gdal-config --libs` -o hillshade
Again, move the binary to a convenient location such as ~/bin.
Preparing Data
Getting Raw DEM Data
With the tools set up, it's time to get some data. First, decide on the area you want to try out.
SRTM data is arranged in 1°×1° tiles. At least north of the equator and east of Greenwich (e.g. most of Europe), the tile name is determined by its southwestern corner - in other words, round the coordinates off to integer values. I decided to run my tests on the 3 Vallées ski area at 45.3211N/6.5934E, so the tile I want is N45E006. Each tile comes as a ZIP file, which contains a single .hgt file.
I used Viewfinder's DEM for the Alps. See SRTM#Obtaining for sources of DEM data.
From Raw DEM to Hillshading
Now comes the moment of truth: we're converting the raw elevation data to a nice image file. First we need to convert the zipped .hgt file to GeoTIFF. Remember to use srtm1_generate_header.sh (the modified version) if you are working with 1" resolution data - else you'll get little more than garbage out of the process.
srtm_generate_hdr.sh N45E006.hgt.zip
This will create a bunch of new files in the directory. The .tif file is what we're interested in. It's basically a georeferenced image file with elevation data represented as different brightness values. You can try opening it in GIMP and normalize the colors to see it (don't save the results back to the original file, though).
First we need to make sure the SRS is specified in the file (which is not always the case):
gdal_translate -of GTiff -co "TILED=YES" -a_srs "+proj=latlong" N45E006.tif N45E006_adapted.tif
Now warp the image into whatever projection you intend to use for your map - in most cases this is spherical Mercator.
gdalwarp -of GTiff -co "TILED=YES" -srcnodata 32767 -t_srs "+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m" -rcs -order 3 -tr 30 30 -multi N45E006_adapted.tif N45E006_warped.tif
Now, create the hillshading:
hillshade N45E006_warped.tif N45E006_hillshade.tif
You can open the hillshade file in your image editing program - you should able to clearly make out mountain features.
GDAL 1.7
GDAL 1.7.0+, which will be included in Ubuntu 11.10, includes the hillshade command as gdaldem hillshade, along with the other utilities from PerryGeo. From GDAL 1.5.0, .hgt files are supported without using srtm_generate_header.sh, you only need to unzip them.
The whole procedure for GDAL 1.7+ is:
unzip N45E006.hgt.zip gdal_translate -of GTiff -co "TILED=YES" -a_srs "+proj=latlong" N45E006.hgt N45E006_adapted.tif gdalwarp -of GTiff -co "TILED=YES" -srcnodata 32767 -t_srs "+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m" -rcs -order 3 -tr 30 30 -multi N45E006_adapted.tif N45E006_warped.tif gdaldem hillshade N45E006_warped.tif N45E006_hillshade.tif
Tweaking the Results
If the bounding box you need is smaller than the original DEM tile, you can use the -projwin option for the gdal_translate call.
To get a better resolution, try the following parameters in the gdalwarp call:
-tr 15 15 -wt Float32 -ot Float32 -wo SAMPLE_STEPS=100
The -tr option sets the resolution in meters per pixel. The original setting of 30m/px is already higher than the resolution of the original data, but interpolation will still make it look smoother than a bitmap rendered at lower resolution and scaled to size. The -wt and -ot options can be useful for removing artefacts. Note that Mapnik also has an option to interpolate the rendered relief - leaving the interpolation to Mapnik allows you to keep the raster files for input small.
If you find the relief a little "flat", try
-z 2
in the call to hillshade. The -z parameter (known as ZFactor) will multiply heights with the factor specified, thus setting it to a value greater than 1 will give the expression of steeper slopes and/or a clearer sky. (A ZFactor between 0 and 1 will flatten the relief accordingly.) You may want this if you are going to overlay parts of the relief with semitransparent map features.
By default, hillshade will simulate light coming from the top left (i.e. northwest). This is mostly for psychological reasons - areas whose top/left border is brighter than the bottom/right border are perceived as bumps, areas whose top/left border is darker are perceived as pits. (Just look at the GUI controls of your operating system.) However, on the northern hemisphere this is an unnatural position for the sun - slopes which in reality are always in the shadow will now appear in bright light and vice versa. Some extra command line options allow you to control the position of the "sun":
-az controls the azimuth of the sun - 0° means light coming from the north, 90° from the east, 180° from the south and so on. The default is 315°, which is north-west.
-alt sets the altitude of the sun, as an angle. Default is 45°.
Rendering Hillshading
Now it's time to get the relief into Mapnik.
A First Run
Since the hillshading file is opaque, it has to be the bottom layer, i.e. the first layer you define in your map file.
Add this to your map file, before the first layer definition:
<Style name="raster"> <Rule> <RasterSymbolizer> <CssParameter name="scaling">bilinear8</CssParameter> </RasterSymbolizer> </Rule> </Style> <Layer name="dem-N45E006" status="on"> <StyleName>raster</StyleName> <Datasource> <Parameter name="type">gdal</Parameter> <Parameter name="file">N45E006_hillshading.tif</Parameter> <Parameter name="format">tiff</Parameter> </Datasource> </Layer>
Refinements
You may notice that areas cover up the relief entirely. There are a few approaches to this:
- Render areas semitransparently to keep the relief underneath visible.
- Render areas underneath the relief and make the raster layer semitransparent by including the following line in your RasterSymbolizer element:
<CssParameter name="opacity">0.7</CssParameter>
- Another approach is described in Hillshading using the Alpha Channel of an Image - here a semitransparent image is overlaid over the map. Note that in this case, the relief will shade not just areas, but also icons and text labels.
Another area for improvement has to do with the fact that each SRTM tile requires its own data source in the map file. You may want to keep these in a separate XML file, which you generate on the fly with a script, and include a reference to that in your main map file.
Finally, you may also want to create a batch script if you intend to create hillshading images for large regions.