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:

Typical usage of a hillshade with false colour overlay

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

The terrain colour map I produced

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:

The result of merging my hillshade and my terrain colour map

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:

Final result: A false coloured elevation map for Tanzania

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!

by Tim Sutton at 8:39 PM under fossgis , linux , qgis , ubuntu