One of the most popular posts on my blog ever is my article about creating GDAL hillshades: ‘A workflow for creating beautiful relief shaded DEMS using GDAL‘. In the techniques outlined in the aforementioned article, colours are assigned across a continuum to and superimposed on top of a hillshaded model. The GDAL tools in QGIS have come along nicely thanks to the work of Guiseppe Sucameli – to the point that you can pretty much do the workflow from the aforementioned article within the QGIS gui environment now (barring some simple text file editing activities).
Sometimes one wants to simply assign arbitrary colours to pixel values of a raster to create a colour map. You can do this with QGIS, but the results are not persistent in the dataset, which I needed because I wanted to create my own colour map and then merge that raster with a shaded relief model. There was one extra requirement I had – I wanted to create my colour map based on polygonal areas at the beginning of my workflow. In this article I describe how I went about creating a hillshaded relief model using a polygonal overlay to define my height classes. I’m also going to show the mis-steps I took to give you a more realistic view of the process of doing such work. Read on for details…
Rasterizing the vector layer
My polygon vector file represents five height classes stored in an attribute called ‘minhoyde’:
0 500 1000 1500 2000
So to make a raster from this polygon layer, I used the Raster -> Conversion -> Rasterize command in QGIS to create and run this command:
gdal_rasterize -a minhoyde -ts 59620 74579 -l hoyde /tmp/hoyde.shp /tmp/hoyde.tif
This process took rather a long time – as you can see I was generating a fairly large raster and the dataset had a lot of complex polygons to process. If I was smart at that point I would have added these options to the rasterize command before running it:
-co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 -ot Int16
The first three options would have made for a more optimally compressed file (see Rudi’s previous article on this). The last option asks GDAL to output a 16bit integer file – as it turned out, the default in the absence of this is to generate a Float64 file which is massive. Since I missed that step the first go around (and consequently landed up with a 34G file), I used the gdal_translate tool (Raster -> Conversion -> Translate in QGIS) to do this as a separate step:
gdal_translate -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 -ot Int16 hoyde.tif hoyde2.tif
After doing this my output file shrunk to a much more managable 18mb.
Refining the rasterization process
After loading up the layer in QGIS, I realised I made two important mistakes:
- I could have saved my output file as type Byte, possibly resulting in an even smaller dataset.
- I should have specified a nodata value (255 would have been a good choice given point 1 above).
gdal_rasterize -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 \ -ot Byte -a_nodata 255 -a minhoyde -ts 59620 74579 -l hoyde /tmp/hoyde.shp /tmp/hoyde.tif
Band 1 Block=59620x1 Type=Byte, ColorInterp=Gray
Min=0.000 Max=2000.000
Minimum=0.000, Maximum=2000.000, Mean=59.377, StdDev=208.406
NoData Value=255
Metadata:
STATISTICS_MAXIMUM=2000
STATISTICS_MEAN=59.377337101373
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=208.40576721446
you can see that gdal still shows the correct range of values event though it is a byte dataset. If you open the file we have produced in QGIS, it will look something like this:
Our rasterized polygon layer (click for larger image)
Now if we wanted to we could assign a colour pallette to it in QGIS by setting the Layer -> Properties -> Symbology Tab -> Colour Map option to ‘ColourMap’ and then defining a colour for each unique value on the ColourMap tab. To illustrate, I have done this to replicate the look I am aiming for when I generate my colours layer in the next step.
Assigning a colour map to our rasterized polygons in QGIS (click for larger image)
When you assign a colour map in QGIS, other applications are not aware of this colour mapping, so what we are aiming to do in this tutorial is to embed the colour map into the image using GDAL, which I show in the next step.
Generating a colour relief map
The next step in the process was to assign colours to the rasterized polygon map. To do this, I basically treat the rasterized polygon layer as if it was a DEM and use the gdaldem colour-relief tool to assign colours to it. Before I could do this, I needed to understand the values that had been assigned for the different height classes during the rasterization process above. Doing this was a simple matter of opening my ‘hoyde’ dataset (the rasterized polygon layer created above) in QGIS and clicking around until I had identified the various unique values the dataset contains. In the process I generated a little table like this in text editor:
#HoydeRaster HoydeVector
#0 0
#244 500
#232 1000
#220 1500
#208 2000
Now that I know which pixel values represent which heights, I can assign colour values (which were supplied by my client) to each height class by creating a colour lookup file that looks like this (saved as colours.txt):
0 255 242 175
244 235 210 134
232 211 187 131
220 168 163 135
208 255 255 255
The first column is the pixel value and the subsequent columns represent RGB values. Good now everything is in place and I can generate my coloured raster using Raster -> Analysis -> Dem from withing QGIS:
Generating our colours layer using gdaldem color-relief (click for larger image)
Or run from the command line:
gdaldem color-relief /home/web/utno/qgis-web-client/data/hoyde.tif \ /home/web/utno/qgis-web-client/data/colours.txt \ /home/web/utno/qgis-web-client/data/colours.tif \ -alpha -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9
After this command ran, I had a new tif image called ‘colours.tif’ – which we can then use as a basis for hsv_merge.py.
Merging the shade model and the colour model
The last step in the process is to merge the colour model with our hillshade. This process is the same as described in my previous article - basically you need to download the hsv_merge.py script and run it from the command line) like this:
hsv_merge.py colours.tif norge_hillshade_shpclip.tif colour_relief.tif
Afterwards you should have a generated file called colour_relief.tif which is a composite of the hillshade and the colour map. Here is what it looked like for me when I was done. I put a few vector layers over the top for extra credits 
Our rasterized vector layer merged with the hillshade (click for larger image)
Managing the output file size
Although we are pretty much done with our workflow, it would be worth your while to try to optimise the output file size a little (in my case the output above was 17GB). Firstly the image we created is an RGB image – but it is a good condidate for conversion to a palletted image using the gdal ‘rgb2pct.py’ application. Also we could try adding some compression options to the created pct file to make it smaller still. In QGIS you can do this using Raster -> Conversion -> RGB to PCT and filling in the dialog as shown in my example below:
Converting an RGB to Palletted Image (click for larger image)
Or on the command line:
rgb2pct.py colour_relief.tif colour_relief_pct.tif
That reduced my file down from 17GB to 4.2GB. As a final step I used gdal_translate to compress the tif (once again, see Rudi’s previous article on what these options mean).
gdal_translate -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 colour_relief_pct.tif colour_relief_pct_compress.tif
This resulted in a file size of 324MB – much more acceptible. Unfortunately that wasn’t quite the end of the story for this job. As I mentioned, doing this kind of work is a lot about trial and error. In this case converting the RBG image to a PCT caused the resulting images to have quite a grainy appearance.
Artifacts in palletted image mean that despite the smaller file size I wont use it (click for larger image)
I actually quite liked it as it looks a bit ‘arty’ but on the sum of it, the smooth tones of the original colour chaded DEM looked better to me so I used the translate step above directly on the original to compress it (forgoing the reduced colour pallette).
Compressing the output image (click for larger image)
Or from the command line (I’ve added a few more options to the command line version below, dropping the alpha channel):
gdal_translate -of GTiff -b 1 -b 2 -b 3 -co COMPRESS=DEFLATE \ -co PREDICTOR=2 -co ZLEVEL=9 colour_relief.tif colour_relief_compress.tif
This resulted in a file size of 1024MB with a much better fidelity to the original. For performance reasons, you will probably want to forgo disk space a little and build pyramids / overviews for the final raster too. A final step you may want to do is to apply a Mask to the output layer – which is left as an exercise to the reader – I’ve described the process also quite well in my original article. Note that you can now do the masking procedure from within QGIS using the raster tools using gdalwarp with the cutline option.
Conclusion
There is no excuse for having an ugly looking map on your FOSSGIS desktop GDAL and QGIS give you all the tools you need to make something beautiful and to be flexible about how you produce it. I’m not quite done with my workflow yet as I actually intend to merge other thematic layers (landuse in particular) into the colour shaded map so that I can have a single raster backdrop which depicts both landuse and height classes. But hopefully the above procedure will give you some ideas about alternative ways in which to generate shading classes for your hillshade / dem. For example you could use flood risk polygons, demographic data or other polygonal datasets you have kicking around to colour your hillshade with.
Producing a raster backdrop in this way can offer some nice rendering performance gains and help you to simplify your projects. Ultimately this work is going to be published on the web so I am trying to speed things up a little and also work around a banding issue in QGIS’s print composer when you work with semi-transparent raster data.
Addendum
One extra step that I performed that may be worth a mention was to use gdal sieve to remove ‘salt and pepper’ effect that was caused by small tranparent pixels being introduced between the polygons in the rasterization process. Example usage:
gdal_sieve.py -st 2 -of GTiff /home/web/utno/qgis-web-client/data/mask2.tif \ /home/web/utno/qgis-web-client/data/mask3.tif
I do a lot of mapping for field sales forces. Have to take into account all types of terrain which may impede travel. I use a program that is not the same or as complicated. But the concepts are similar. Your article will help me immensely when producing a viewable map with the output of territories and regions.
[...] finally, the Linfiniti Geo Blog has a great post on Creating coloured rasters with GDAL. If you haven’t checked out Linfiniti’s blog, you should – I’ve found it to [...]
[...] recently posted an article on creating coloured rasters with gdal which outlined a procedure for applying colour from polygonal areas to a raster. I used the [...]