Mapnik/Hillshading using Mapnik, GDAL and SRMT data

From OpenStreetMap Wiki
Jump to navigation Jump to search

trash bin

It has been proposed that this page be deleted or replaced by a redirect. See the discussion page for further information.
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).
Logo.png
This page describes a historic artifact in the history of OpenStreetMap. It does not reflect the current situation, but instead documents the historical concepts, issues, or ideas.


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 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.