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).
This is a normal part of iterative development work - I was trying to create a new workflow and making a few wrong turns is inevitable. So here was the complete revised gdal_rasterize command which I ran:
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
The first interesting thing I noticed when generating the raster was that the output size was no different from the Int16 version I created earlier. Also, you may have been worried reading the above about how gdal would map data values ranging from 0-2000 onto a Byte data set (which should only support pixel values from 0-255). Well take a look at this gdalinfo snippet:
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:
Now if we wanted to we could assign a colour pallette to it in QGIS by setting theLayer -> 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.
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:
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):
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:
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 :-)
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:
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.
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).
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.
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.
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