Sometimes I create hillshades using the QGIS hillshade plugin and then overlay the original DEM over it. I set the DEM to have a false colour pallette and set it to be semi-transparent to produce something like this:

That is all well and good but a bit impractical. It would be much nice to have the colour pallette permanetly assigned to the hillshade. Also I want to be able to clip and mask the resulting raster to the boundary specified in a shapefile. Fortunately, GDAL provides all the tools you need to make this happen. There are already some nice articles (like this one) that describe parts of this workflow but I am writing this because I wanted to note the additional steps that I took to make this work well for me.

Before you begin

Before you begin you should have:

1. a raster containing digital elevation data (DEM) – in my case its called ‘alt.tif’
2. a vector layer (typically a shapefile) containing the area of interest for your final product – in my case its called ‘tanzania.shp’

Create the hillshade image

The first thing we need to do is generate a hillshade. There is a nice python plugin for doing this in QGIS, you can do it in GRASS using the QGIS-GRASS plugin. But in this article I’m going for an all-GDAL approach so we will be using the handy gdaldem command line application.  I won’t be explaining the parameters I have used here as they are well explained in the gdaldem manual pages.

So to create our hillshade we do something like this:

gdaldem hillshade alt.tif shade.tif -z 5 -s 111120 -az 90

Which will produce something like this:

Our initial hillshade image

Create the terrain map

Here we want to create a pleasing colour palette for the terrain. There is a lot of content on colour theory available out there on the internet – googling for ‘colorbrewer‘ would be a good place to start if you want to learn more.  Another favourite site of mine is colourlovers.com and for this tutorial I decided to use a pallette from there to see how it would turn out.

Once you have selected a suitable colour pallette (around 5 or 6 colour classes should do it), the next thing you need to do is get some information about the range of values contained in your DEM. Once again you can easily point and click your way to this in QGIS, but here is the way to get it in gdal from the command line:

gdalinfo -mm alt.tif

The resulting output includes the computed min/max for Band 1 – which is what we are after:

Band 1 Block=1334x3 Type=UInt16, ColorInterp=Gray
 Computed Min/Max=1.000,5768.000
 NoData Value=65535

Ok so our minimum height is 1m and maximum is 5768m – Tanzania is the home of Kilimanjaro after all! So lets split that range up into 5 classes to match the ‘Landcarpet Europe‘ colourlover pallette I selected. I set nearly white as an additional colour for the highest altitude range.

65535 255 255 255
5800 254 254 254
3000 121 117 10
1500 151 106 47
800 127 166 122
500 213 213 149
1 218 179 122

The value in the first column is the base of the scale range. The subsequent values are RGB triplets for that range. I saved this as a text file called ‘ramp.txt’ in the same directory as my ‘alt.tiff’ dataset. You will notice that I made the value ranges closer together at lower altitudes and wider appart at higher altitudes. You may want to experiment a little to get pleasing results – especially if you have a relatively small number of high lying terrain pixels and the rest confined to lower lying areas.

Also note that I assigned the highest terrain ‘nearly white’ so that I could reserve white (RGB: 255 255 255) for the nodata value (65535) in this dataset. We will use the fact that white is only used for nodata to our advantage when we do a clip further on in these instructions.

Ok now we can use gdaldem to generate the terrain map:

gdaldem color-relief alt.tif ramp.txt relief.tif

This is what my relief map looked like:

Don’t worry about the fact that it does not resemble the colour pallette you have chosen – it will do in the next step!

Merging the two rasters

The next step is to combine the two products. I used Frank Warmerdam’s handy hsv_merge.py script for this purpose.

./hsv_merge.py relief.tif shade.tif colour_shade.tif

Which produced this:

You may have noticed that it is only at this point that the colours of our product resemble the original pallette we used.

One little gotcha with the hsv_merge.py script is that it throws away our nodata values, so what was sea before (and nodata in my original alt.tif dataset) is now white (RGB: 255 255 255).

Clipping and masking

You may have everything you need from the above steps. Alternatively you can also clip and mask the dataset using a shapefile.

gdalwarp -co compress=deflate -dstnodata 255 -cutline Tanzania.shp \
          colour_shade.tif colour_shade_clipped.tif

My final image is now a compressed tif with nodata outside of the country of Tanzania and looks like this:

A final note

One of the things I love about the command line is the repeatable and reusable workflows it allows for. I can distill everything in this blog post into a sequence of commands and replay them any time. If you are still stuck doing things in a GUI only way, give BASH a try – you can really do powerful geospatial data processing with it!