Making shaded relief from digital elevation models (DEMs) in QGIS: a British Columbia perspective

Who would not agree with Tom Patterson, creator of the fabulous website, when he said, “There is no more important component of a map than the shaded relief.”

But the topic of creating your own shaded relief from a DEM is rather complex, so I’ve made a few assumptions in this how-to. I’m assuming:

  • you know your way around QGIS (I am using QGIS 2.18 for this post)
  • you’re familiar with the idea of projections, and re-projecting raster data
  • you know that each projection/datum combination has an EPSG number, which is a convenient way to refer to it. For example
    • Lat/Long/WGS84 = 4326
    • UTM Zone 9N/WGS84 = 32609
    • BC Albers = 3005

If that’s the case, there are really only three things you need to know in order to made your own hillshades:  where to get the data, how to transform it, and what pitfalls to avoid.

Why even make your own hillshade?

Usually when I want to include shaded relief in a map that is in British Columbia, the first place I will turn is the WMS service at

001_WMS hillshade example

The 315 black-and-white hillshade from The black-and-white hillshade with the sun at an azimuth of 315° is typically the most useful of the layers on this WMS server.

It’s got some limitations:

  • there are only two sun azimuths to select from: 315° (northwest) and 225° (southwest). More flexibility would be good, because each landscape seems to have a different ideal azimuth to bring out the landforms that you want to bring out. More about this down below.
  • there’s no height exaggeration possible
  • occasionally there are odd artifacts, like straight lines running across the hillshade

On the other hand, you can adjust its brightness and contrast, and use the Multiply blend mode, which means you can do some nice things with it.

002_WMS hillshade example with overlay

When you set the blend mode of a hillshade to MULTIPLY, it shows through upper layers.

But, if you want to adjust the azimuth or exaggerate height, you’ll need to find a DEM and make your own hillshade.  It’s well-established (though not well-explained) that the human eye needs to see light from above, preferably from above and to the left. If your map is, say,  south-up, you will need a hillshade where the light comes from the southeast (which will be in the upper left).

004_WMS hillshade example rotated 180

Now south is at the top (map is rotated 180°) and the shaded relief appears flat or inverted. For this map you would need a hillshade where the sun’s azimuth is southeast, or 135°.

005_DEM hillshade rotated 180

Here the map is still rotated 180°, but the hillshade has been manufactured with an azimuth of 135°. The mountains look like mountains again.

006_DEM hillshade rotated 180 3x

A 3x vertical exaggeration of the terrain emphasizes the detail in the valley bottoms.

The overall process

007a_flow chart

Let’s talk about the three principles.

  1. The hillshading algorithms require a DEM in a metric projection. That means that DEMs projected in degrees won’t work: you have to re-project them first. Unfortunately, just about all DEMS come projected in degrees.
  2. The scale of your final map determines what sort of cell size you want in your re-projected DEM. A DEM with 10m cells is far too detailed for a map at 1:500,000, and the file would be enormous. On the other hand, a cell size of 500m would make a very coarse hillshade at 1:500,000. As a general guideline, divide the denominator of your scale by 5000 to get roughly what cell size you want. So if your map is 1:100,000, you’d be looking for a DEM with (roughly) 20m cells.
  3. It is an option to style any DEM as “Hillshade” (other options are singleband grey, multiband colour, paletted, etc.) but the hillshading algorithm in the toolbox (Processing>Toolbox) produces a better result.

Acquiring DEM data

The first thing you need to do is pick your resolution. Because DEMs usually come projected in 4326, DEM resolution is typically expressed in arc-seconds, or seconds of latitude. Because degrees of latitude in BC are bigger than degrees of longitude, these cells are not square. They are upright rectangles.

What you want to know is how these non-square cells measured in arc-seconds will convert to square cells measured in metres. Here’s a handy chart.

resolution degrees resolution metres pixel size in degrees source limitations recommended scale Notes
1/3 “ 6 9.26E-05 USA only ~ 1:30,000 It comes as 1°x 1° tiles. Files have names like “,” which would be the 1°x 1° tile northeast of 44°N, 110°W.
1” 17 0.000208333 SRTM 1 for all of North America at
Europe at
North America and Europe ~ 1:100,000 From recordlist these come as 1°x1° tiles in HGT format, each about 25MB. The file N55W128.hgt would be north and east of 55°N, 128°W
3” 45 0.000833333 This coverage, based on SRTM data, is available for the world in two different versions.

worldwide~ 1:300,000v4 data is better, but it comes as 5°x5° tiles. The SRTM 3 comes as 1°x1°tiles in HGT format. The old CDED dems were this resolution.

15”4500.00416667 1:1,500,000  250m and 500m (Password: ThanksCSI!) worldwide~ 1:1,000,000 and 1:2,000,000re-samplings of the finer resolution data30”1 km0.00833333 Or 1:3,000,000at this point you should be considering the shaded relief at’10 km 1:30,000,000 

You’ll notice that I favour the DEMs created from the Shuttle Radar Topography Mission (SRTM). These have a few advantages over the DEMs produced by Natural Resources Canada or the provincial mapping agencies:

  • they are usually more recent (dating to about 2000)
  • they represent actual measurements, as opposed to a grid generated from contour lines
  • they go right across provincial and international boundaries

From here, I’ll demonstrate how this works with an actual example. In this case I want shaded relief for a series of maps that I’m making. All the maps are in the same area, and they range in scale from 1:45,000 to 1:130,000. Looking at the handy chart above, this scale suggests I’m going to want to use SRTM 1 DEM.

007c determining the tiles you want2

Displaying an area in a print composer with grid of lines in CRS 4326 is one way to figure out which tiles you’ll need.

So I go to Recordlist, enter the bounding box and select SRTM1.

008_recordlist selecting DEMs

Downloading DEM tiles from Recordlist

Once they are downloaded and unzipped, I’ll read them into QGIS to confirm that they cover the right area.

Merge and Clip

You want to merge all of the individual DEMs into one, using Raster>Miscellaneous>Merge. (Incidentally, they do not have to be read into QGIS to do this.)

016_merging DEMs

QGIS’s raster merge dialogue

I tend to name the resulting DEM with “_4326” on the end so that later I will know what projection it’s in.

It’s tempting to display as hillshade now, but don’t. The hillshade styling is not meant for DEMs projected in degrees.

040_displaying as hillshade when projected in degrees

Styling the DEM as “hillshade” produces ugly results when the DEM is projected in degrees

If you want to clip the merged DEM, now is the time to do it. Remember that with Raster>Extraction>Clipper you will need to change the QGIS projection to the DEM’s native projection (4326) before you draw the clipping box. Be sure to check, once you come back to your map’s projection, that the clip you made covers your whole print composer.

Re-project and re-sample

Reprojection is the process of giving raster cells coordinates in a new projection. Resampling, on the other hand, is the process of creating new cells based on old cells at a coarser or finer resolution. The two are essentially inseparable, since as you reproject from 4326 to, say, 32609, you will also want to go from the rectangular cells of the degree-projected DEM to nice square cells in the UTM projection.

The first thing to do is to right-click the DEM and choose Save As… so you can see the dimensions of the original DEM cells.

050_Save As dialogue

Note that once you set the CRS for the saved copy to be 32609, you get a suggested resolution of (roughly) 17 x 31m cells. That’s the native cell size of this SRTM DEM at this latitude. Note down the “17” (the smaller dimension)  somewhere, and close this dialogue.

You can reproject and resample using QGIS’s Save As… feature, but you don’t get control over the resampling algorithm used, and the results, once you get to the final hillshade, are ugly.

060_hillshade produced with QGIS Save As

When you reproject using Save As…, the final hillshade sometimes has these weird cross-hatching artifacts in it

Instead, you want to re-project and re-sample with the Warp tool in the toolbox. (Go Processing>Toolbox and search on “warp.”)

084_reprojecting using gdalwarp through QGIS toolbox

The toolbox’s Warp tool

In this dialogue…

  • Source SRS should be 4326
  • Destination SRS should be 32609 (or whatever metric projection you are making your final map in)
  • Output file resolution can be whatever you want, but I get good results with the smaller of the two cell dimensions I saw in the Save As dialogue — in this case, 17.
  • Trial and error with making hillshades has convinced me that the best resampling method to choose here is “bilinear.”
  • I name the resulting DEM with a “_32609” on the end so I will know its projection in the future.

The new 32609 DEM should look the same as the 4326 DEM when read into QGIS, but if you go to the Metadata tab in Layer Properties you’ll see it has a cell size of 17 metres, and quite different pixel dimensions.

For your final hillshade, the one you use in your map, you will want something better that what you get if you just style this DEM as “hillshade.”  But for now, go ahead and style it as “hillshade.” This enables you to play with the sun azimuth and elevation, and the vertical exaggeration, to see what is going to work best for your terrain. The human eye probably wants an azimuth around north-northwest (337.5) but there are a lot of azimuths on either side that will work.

Notice how changing azimuth changes what’s brought out in this piece of terrain:

Make a static hillshade

Now if this hillshade has no strange artifacts, you’re done. But if you want a smoother results, it’s time to take the azimuth, sun elevation and vertical exaggeration you’ve chosen, and go over to the Hillshade tool in the toolbox (go Processing>Toolbox and search on “hillshade”). This will produce a raster hillshade that you just display as singleband grey. It also tend to be about half the file size of the DEM itself.

086 running toolbox hillshade2

The toolbox’s Hillshade tool

And here’s the result

090_hillshade produced with gdalwarp gdaldem hillshade

Look, no artifacts!

Displaying hillshades

Typically you put the hillshade as the bottom layer in your map, and adjust brightness and contrast — because hillshades tend to be a bit darker than you want. Usually it’s brightness up, contrast down. Sometimes you will also make it semi-transparent.

More important, in terms of getting the rest of you map layers to appear to drape over your hillshade, is to set the blend mode of the hillshade to Multiply.

091_adjusting brightness etc

Multiply causes the values (lightness, darkness) in the hillshade to be adopted by layers above it.

So you can get something like this.

095_final result with multiply

Overlaying a hillshade whose blend mode is set to Multiply

That’s it. Using these techniques, you should be able to manufacture hillshades with any azimuth, pretty much all over the world. It gets the most challenging when you are mapping at large scales like 1:20,000.

Finally, if you master this, and you are really in love with shaded relief, you will want to experiment with making shaded relief in Blender.


Buffering (or haloing) text over complex backgrounds using the Screen blend mode

This is a common problem when you are working over a complicated, multi-colour background. Such backgrounds tend to swallow text!

Here’s an example, where text lies over a background with a lot of variability in value (darkness).text over complex background

Here’s the same image in greyscale, just to show the pixel values. Out of a total range of 0 to 255, the pixels here range from 6 (almost black) to 240 (almost white). Black text over light portions reads easily, but it’s hard to read the letters, M, N and T because of the way they overlap dark areas.value map of background

The classic, or perhaps I should say default solution, is to buffer the text in white. This works fine, in that you can read the text more easily, but it’s a crude solution and it’s a bit loud. It lacks subtlety. white 1mm buffer

A somewhat more elegant solution is to have the buffer match the background colour. In this example the text stretches across many background colours, so no matter what colour you select it’ll be a compromise. Still,  it looks better than white. The problem is that across the map as a whole each piece of text will probably get a different buffer colour, depending on what is dominant in its background.1mm colour buffer

What we really want is a way to pick up the background colour for each pixel, and colour the buffer that way. We can actually approximate this by using the blend mode of Screen.

Screen is the opposite of Multiply. Where Multiply transfers dark values onto another layer, Screen uses dark values to lighten another layer.

In this case, you set your buffer to something fairly dark, like 80% grey. Then set the blend mode of the buffer to Screen. You get this.screen 80 percent grey

Here the buffer colour echoes the underlying colour. The buffer around the T is orange-y, while the buffer around the O is grey. The buffer for the M uses both colours.

How would you do this in QGIS and Inkscape?

For text labels displayed in QGIS, this impressive piece of software can do this automatically. On the Buffering tab of the Labels dialogue, check Draw text buffer and specify an 80% grey for the buffer fill. (Note that 80% grey is RGB 51/51/51, HSV 0/0/20, or HTML #333333.)

Then set the blend mode to Screen. There are a couple things to note. One is that the sample shown will look ghastly; you can ignore that. The other is that the blend mode of the label is different from the blend mode of the buffer. The label itself probably has the default blend mode of Normal, and this can be checked on the Text tab.QGIS_buffer_screen

Here’s the result: Tiltusha example

In Inkscape, where you might be doing more complex text, text set on curved paths with adjusted kerning, you can also do this. Blend modes in Inkscape are applied to an entire layer, so you can’t give one piece of text a different blend mode than another. But by introducing a special layer that has a blend mode of Screen just beneath your text, you can put the buffers there.

The procedure goes something like this:

  1. Create a new layer below your text layer and set its blend mode to Screen.
  2. On the text layer, select your text and duplicate it with Ctrl-D. duplicate
  3.  Turn on Stroke for the duplicated text. (Text normally has fill but no stroke) and set the stroke colour to 80% grey. (You can do both of these steps simultaneously with a shift-click on the 80% grey swatch in the Palette at the bottom of the screen). Set its width as you see fit (1mm in this example). It should look pretty ghastly.buffer stroke on
  4. Press Shift-PgDown to send that duplicate piece of text to the layer below, the layer with blend mode of Screen. It should now look pretty nice. buffer moved to screen

I should point out one disadvantage of this two-layer approach: you can’t group the text with its buffer, a practice I normally do because it makes it easier to move labels around later.