# Making shaded relief directly from DEMs projected in degrees?

In a presentation I gave at GeoIgnite 2021, I was explaining the process of using QGIS to make shaded relief from DEMs, more or less as detailed in this post. There wasn’t time to explain the process of re-projection, so I rather boldly asserted that you don’t need to re-project your DEMs before making shaded relief. Instead, I said, you can just use the magic scaling number 111120 to convert from degrees to metres.

Really? Is this true?

Well, yes, in many cases it works pretty well. In QGIS you pull in an unreprojected DEM (its native CRS using degrees as units) but set the Project CRS to some reasonable metric projection, such as UTM. Then, opening the GDAL hillshade tool from the Toolbox, you use the DEM as your Input layer, set scale to 111120, and enter the parameters of your choice for Z factor, Azimuth of the light and Altitude of the light.

This seems so much simpler than having to first re-project your DEM to square cells measured in metres. But there are some caveats. This technique fools with the terrain a bit, and could potentially distort how it looks. If you understand how it distorts it, you’ll be able to evaluate if you want to use it in your map.

### The magic number

The magic number 111120 comes from the GDAL documentation, specifically the documentation for the gdaldem tool:

The gdaldem generally assumes that x, y and z units are identical. If x (east-west) and y (north-south) units are identical, but z (elevation) units are different, the scale (-s) option can be used to set the ratio of vertical units to horizontal. For LatLong projections near the equator, where units of latitude and units of longitude are similar, elevation (z) units can be converted to be compatible by using scale=370400 (if elevation is in feet) or scale=111120 (if elevation is in meters). For locations not near the equator, it would be best to reproject your grid using gdalwarp before using gdaldem.

https://gdal.org/programs/gdaldem.html

Gdaldem is the tool that I recommend using for creating hillshades in QGIS, and there in the docs they say it: unless your locale is near the equator, it’s best to re-project your grid before using gdaldem.

The obvious problem with entering a scale factor of 111,120 metres per degree is that as you move away from the equator, lines of longitude become closer to one another. For example, in Canada at 54° north, one degree of latitude is close to 111,0000 metres, but a degree of longitude is only about 65,000 metres. In other words, a DEM cell is not square on the ground. So when you enter 111,120 as the scaling factor to convert degrees to metres, you will get a hillshade calculated on the basis of a square cell that has been created by stretching a rectangular cell in an east–west direction.

For example, with a 1″ DEM, each cell is 0.0002777777° (that’s 1 ÷ 3600) on a side. If we scale both axes at 111,120, the hillshading algorithm thinks it has a square cell 30 metres on a side. However the reality at 54° north is that on the ground the cell is about 30 m north–south and 18 metres east–west.

You’d expect that because the landscape is stretched east–west, the slope would be under-calculated for cells that slope east–west. Basically, these cells would appear to be on a slope less steep than it really is.

And slope is important in the hillshading algorithm: it plays a role in determining how bright a hillshade pixel is. Pixels whose aspect and slope indicate that they are orthogonal to the sun direction are the brightest white, while those facing away from the sun direction are darkest black.

So if, for example, you had a 45° slope facing west, and the sun in the west at an elevation of 45 °, you’d expect pixels on the slope to have the highest possible hillshade value, 255. But if you calculated the hillshade from a degree-DEM using a scale of 111,120, that slope would be under-estimated; maybe it would appear to be 35°. As a result it would not get assigned a 255, but something less, maybe 240. There might be a nearby area, also facing west, whose true slope was 55°, but it now gets assigned the brightest value because it appears to be 45°.

In general then, for east or west-facing surfaces, we’d expect the hillshade built from the un-re-projected degree-DEM to distribute its brightest and darkest patches differently, particularly when the sun is due east or west.

That’s just what you see.

However, when the sun is not due east or west it’s hard to see this difference in slope calculation.

You can also experiment with other, lower scale factors, such as 90,000, 80,000 or 70,000. These, in effect, cause slopes to be less under-estimated on east–west aspects, but over-estimated on north–south aspects. In other words you are distributing the error around.

This suggests that in many cases, making a hillshade directly from a DEM projected in degrees is fine.

### Other questions

A second distortion arising from creating a hillshade using a degree-DEM and a scale factor of 111,120 concerns the sun’s azimuth. If you are mapping in a local UTM projection, you know that UTM grid north is usually different from true north. This was visible in the figure above with the contour lines: the cell that was “square” in degrees turns out to be rectangular in a UTM projection, and also a bit tilted to the side.

With UTM projections, this difference is not more than 3°. At worst, then, you might ask for a sun azimuth of 315° and instead get one of 312°. I doubt whether you can notice this difference.

The last problem we might wonder about is whether shadows calculated by the hillshading algorithm are distorted. Are they foreshortened in an east–west direction when the hillshade is re-projected into a metric projection? The answer to this may be surprising: the GDAL hillshading algorithm does not calculate shadows. It calculates shadowed surfaces, but this is purely on the basis of the slope and aspect of the cell it is looking at (and its eight neighbours).

In other words, a piece of terrain at the bottom of a deep gorge will be rendered purely on the basis of its slope and aspect, as if the sun could shine on it right through the surrounding mountains.

Hence, no, we do not need to worry that shadows will be shortened.

### So…

All of these things considered, I come to the surprising (for me) conclusion that yes, you actually can often just hillshade a DEM that’s still projected in degrees, using the scale factor 111,120. It seems heretical, but when I can’t tell the difference in the results, I have to admit that for practical purposes it’s a good solution.

# Shaded relief with BlenderGIS (2019), part 3

## Why do this?

The shaded relief images (also often called a hillshade) that GIS software produces are pretty good. Here’s one made with QGIS.

Hillshading algorithms typically take three parameters: the elevation of the sun, the azimuth (direction) of the sun, and the vertical exaggeration of the landscape. (60°, 337° and 1X, in this case.) They do some geometrical calculations to figure out the intensity of incident light on all the pieces of the landscape, and then produce a greyscale image. (They do not, typically, however, calculate actual thrown shadows.)

You can tweak the brightness and contrast in your GIS software, or stretch the histogram to your liking. This can do a lot to lighten up the darkest shadows, or to turn the hillshade into a ghostly wash that merely suggests relief without dominating everything else on the map.

Much of the time what we can get out of GIS software meets our needs for a hillshade, which, after all, is not the final map but a merely a layer of the map. (Although occasionally it is the pièce de résistance…)

But perhaps you don’t like the glossy, shiny quality of that hillshade above. You think the mountains look like they were extruded in plastic. They remind you a little too much of Google maps Terrain layer. Maybe you’d like to see actual shadows thrown by precipitous cliffs. Maybe you’d like something that looks more like it was chisled from stone, like this…

Or perhaps you are interested in…

This might be why you are investigating Blender.

The big settings that make a difference (besides the conventional settings of azimuth, elevation and vertical exaggeration) are…

• material
• multiple lights and their colours
• amount of light bounce
• denoising the render

## Material

In animation modelling, Material is what reflects, absorbs and scatters light. Blender spends much of its time tracing light rays and deciding how the surfaces they encounter affect them.

With your plane selected, go to the Material tab and hit New. The default material that comes up has these Surface properties. (And this is what Blender used for your basic render.)

Principled BSDF is a sort of super-versatile material that allows you to have all of these properties (subsurface scattering, metallic look, sheen, transmission of light) that in earlier versions of Blender were assigned to specific surface types, like “Diffuse BSDF,” “Glossy BSDF” or “Subsurface scattering.”

(BSDF stands for “bidirectional scattering distribution function.”)

These various surfaces are actually shaders, which are pieces of software that render the appearance of things (and may actually run on your GPU, not your CPU). If you click on “Principled BSDF” next to the word “Surface” a list of all of the possible shaders comes up.

You can learn a lot more about shaders in the Blender manual, but Blender essentially wants to give you the tools to be able to simulate any material, from water to hair, and some of the effects you can get applying this to shaded relief are pretty weird.

The main shader you probably want to play with is the Mix Shader, which allows you to mix the effects of two different shaders. The Mix Shader’s  factor (from 0 to 1) determines how much the results are influenced by the second shader.

The other aspect of material that is worth experimenting with is value. Blender’s default material, the Principled BSDF, has a default colour of near-white, and its value is 0.906 (on a scale of 0 to 1).

By darkening this you create a more light-absorbing material.

You might think that the effect of a darker material is just to shift the distribution of pixel values down in the render, and that you could get the same effect by turning down the brightness of the original. But if you look at the histogram of each image you see the pixel values are distributed differently, and even when both are displayed with “histogram stretch,” they are distinct from each other.

And of course if it works for your map you can give the material real colour. (Don’t forget to save the render as RGB rather than BW). This, however, is very similar to colorizing your shaded relief in GIS software.

## Multiple lights and their colours

To add another light to your scene, go Add>Light>Sun. You can also experiment with adding other types of lights: point lights (which shine in all directions from a specific place), spot lights (which send out a specific cone of light) and area lights.

A second light can bring out features in a very nice way.

By giving colour to lights, you can create differential lighting and colour.

There is one more light in your scene that often goes unnoticed, and this is what Bender called the world background colour. You can think of this as a very far away surface that nonetheless does contribute a bit of light to your scene—from all directions. You can see the world colour in action if you render your scene with the sun strength turned down to 0.

If there are places in your scene that sunlight does not reach, the world background still contributes to their illumination—and colour.

By default the word background colour is 25% grey (value = 0.25), but you can change this to a dark blue if you would like dark blue light to collect in your shadows.

## Amount of light bounce

Part of the charm of Blender is that it calculates how much light reflects off surfaces and then hits other surfaces, which is why even shadowed areas have some light. You can control, however, the number of bounces a light ray has before it expires. This is on the Render tab, in the Light Paths section.

By default, the light paths settings are as above. With the material I am using, the Max Bounces for Diffuse materials is the important part, and by default it is 4.

The icon to the right of “Light Paths” indicates that there are presets available. Clicking here reveals three key presets:

• Direct Light: light gets few bounces, and none off diffuse material
• Limited Global Illumination: light gets one bounce off diffuse material
• Full Global Illumination: light gets 128 bounces off diffuse material

Changing among these does not make a huge difference, but in scenes with deep shadows it is visible.

## Denoising the render

Blender offers the post-processing feature of denoising the render. I liken this to a blanket of snow on your landscape: it erases tiny details.

To activate denoising, go to the View Layer tab, scroll to the very last section, Denoising, and check the Denoising box.

I do not find that denoising does much if you keep the default settings. I tend to turn both Strength and Feature Strength up to 1.0.

## In conclusion

This is only been the tip of the iceberg. Some things I have not talked about are

• Applying a subdivision surface modifier to your mesh so that Blender is interpolating and smoothing it on the fly.
• BlenderGIS’s ability to read your DEM As DEM Texture, which creates a plane with subdivision and displace modifiers instead of a plane whose every vertex corresponds to a DEM cell.
• How to focus the orthographic camera down on one small part of your DEM where you can do quick test renders while you work out lighting, material, etc.
• The node editor for complex materials
• Using a perspective camera to shoot a scene that is not straight down
• Cutting DEMs into tiles for Blender to work with, when the whole DEM is simply too much for the RAM on your computer.

Have fun exploring Blender’s many features!

# Shaded Relief with BlenderGIS (2019), part 2

In the overall process we are now at step 3, but things will go faster now.

2. Read the DEM into Blender as DEM raw data build
3. Adjust Z scaling (vertical exaggeration)
4. Create and adjust a georef camera
5. Correct the final pixel dimensions of the output image to match the DEM
6. Set final image type to be TIFF
7. Turn the light into a Sun, and adjust its properties
8. Do a test render
9. Do a full render

At the end I’ll cover some of the variations on this process and extra tweaks you can do.

## Adjust Z scaling (vertical exaggeration)

The first thing to adjust on your plane is what we usually call the vertical exaggeration, but Blender thinks of as the Z scaling.

With the plane selected, go to the Object tab. Here you should see a number of sections, the first of which is the Transform section.

Under Scale increase the Z value if you want to create vertical exaggeration. You should see your plane change shape.

It’s important to increase Z scaling before setting up the georef camera, because if you increase Z scaling afterwards, you may accidentally have set the camera below the tops of the highest features.

In my case the terrain has plenty of relief, so I’ll leave Scale Z at 1.000.

## Create and adjust georef camera

Make sure the plane is still selected (in the Outliner) and go GIS>Camera>Georef.

A new camera (“Georef cam”) will appear in the Outliner. You do not need to delete the old camera.

When the georef cam is selected, then in the 3D Viewport you’ll see something like this (you may have to pull back).

This is an orthographic camera (meaning, in effect, that all points of its lens look directly down: there is no perspective distortion) placed over the plane.

To check that your camera is properly pointed at your landscape, and sees all of it, hover the mouse over the 3D Viewport, and hit 0 on the numeric keypad (or go View>Viewpoint>Camera). You should see what the camera sees (“camera orthographic view”)

Hit NumPad-0 again (or go View>Viewpoint>Camera again) to go back to regular (“User perspective”) view.

## Output tab: Correct the final pixel dimensions of the output image to match the DEM, and set final image type to be TIFF

Go to the Output tab, and under Dimensions check the Resolution X and Y that were set up when you created the georeferenced camera. Also note the percentage (%), which is 100% at this point.

Change Resolution X and Resolution Y to match the pixel dimensions of your DEM. In my case, the initial DEM was 839 x 702 pixels, so I enter these two numbers for X and Y. With the dimensions of the hillshade matching those of the DEM I can, at the end, apply the DEM’s world file to the Blender output, and georeference it.

It’s possible at this stage in the process to set percentage to something small, like 20%. This way you can do some quick test renders: Blender will produce a rendering whose pixel dimensions are only 20% of the overall Resolution numbers you set. (Before the final render you’ll be back here and set this to 100% again.)

Lower down on the same tab, under Output, set the parameters for the type of image you want.

I typically set these to File Format = TIFF, Color=BW and Color Depth = 8, but you can get JPG, PNG, etc. If you choose later to assign colour (other than grey or white) to the plane’s material, the world background or the sun, you will probably want to set Color to RGB or RGBA.

## Turn the light into a Sun, and adjust its properties

Your plane is built, your camera is set: now all that is left to arrange is the Sunlight.

Select the Light in the outliner, and go to the Object Data tab.

Initially you should see something like this under Light:

Click on Sun, and then change Strength to 1.

Click on the Use Nodes button below, and set its Strength to 2. It should now look like this.

You will want to play around with the Strength setting (in the Nodes section) when you do test renders. The strength of the light affects how bright your hillshade is.

The Angle setting for the sun is important. The Angle is how many degrees across the sun’s face looks from the ground. A sun with a smaller angle produces shadows with sharper edges.

A 1° sun is essentially a point source, and it casts sharp shadows. A 10° sun is bigger and the edges of shadows are diffuse. Take your pick.

Note that if you decide to use a coloured sun later (see “Part 3“) this is where you will change its colour: Light>Object Data>Nodes>Color.

Now go to the Object tab, and consider the location, rotation and scale of the light.

The location of the light does not matter: a “Sun” type light is considered to shine from an infinite distance regardless of its location. Scale should be left as all 1’s.

The rotation of the light however is all-important to us. It is here that we set the sun elevation and azimuth. These are both counter-intuitive in Blender, so here’s how it works.

Rotation X is always 0

Rotation Y is the complement of the familiar elevation angle (that is, 90 minus elevation). E.g., if you want a sun 60° above the horizon, set Rotation Y to 30°. Think of a light in a theatre clamped to a lighting pipe suspended over the stage. The pipe is the Y axis. The light is initially pointing straight down (0°) and you are swinging it up by 30°.

Rotation Z is the azimuth of the light, but instead of beginning at 0° for North and increasing clockwise (as we are familiar with from compass bearings), it begins at 0° for East and increases counterclockwise. (Think of angles measured on a coordinate plane.) It’s generally quickest to figure out by sketching on the back of an envelope, but the actual formula is that the Rotation Z = 450 – Azimuth, and subtract 360 from your answer if it is greater than 360.

A few common angles: If I want a sun 50° above the horizon, coming from azimuth 337° (North-northwest), my light’s Rotation numbers will look like this:

## Do a test render (F12)

There are two good ways to do this.

The first is to switch into camera view (Numpad-0) and change the viewport shading from Solid to Render. In the 3D Viewport, locate this group of icons in the upper right:

and note the group of four circles on the right end. At present, the second one, Solid viewport shading, is selected. Click on the fourth one to change to Render viewport shading. You will see Blender begin to render your scene with the camera and lights you have configured.

The other way to do a test render is to set percentage on the Output tab to something small (e.g., 10% or 20%) and and hit F12 to get a test render at a fraction of the size of the final render. A little Render window opens and after a pause you (hopefully!) get something like this.

Either way, doing a test render often catches common mistakes, like forgetting to set up your light or camera. It also reveals how the overall image will look, and so you can adjust your light strength, azimuth and elevation to increase or decrease shadows and the overall brightness of the image.

When finished, be sure to either switch back to Solid viewport shading or set the output percentage to 100%.

## Render

On the Output tab, adjust Sampling. To the right of the word “Sampling,” note this menu icon:

This icon indicates that presets exist for this setting. Clicking on it, select Final. This increases the sampling numbers for the render, which will improve image quality.

Hit F12 to render.

Note that you can interrupt a render by hitting ESC, or by closing the render window.

I get this result:

Image>Save As allows you to save the image. If I saved this as “Blender shaded relief 1x 50 337.tif” I would then make a copy of the TFW world file I created back at the beginning and name it “Blender shaded relief 1x 50 337.tfw” This makes something georeferenced that I can pull into my GIS.

Now, to consider the many other things we can fiddle with in Blender—material, number of lights, denoising, etc.—go to Part 3.

# Shaded relief with BlenderGIS (2020), part 1

This tutorial replaces my tutorial from five years ago,  “Shaded Relief with BlenderGIS” tutorial. At this point (March 2020) I am using Blender 2.82 and the most recent BlenderGIS addon.  Daniel Huffman’s tutorial, which uses a different technique, has also been updated.

It’s rather an understatement to say that Blender is complex. To keep things from getting out of hand, I’m going to take a straight path, right through the centre of the software, with a focus on getting out the other side with a completed render of shaded relief. But at the end, I will return and explore some of the interesting side trails that lead to the features that make the Blender hillshades so interesting.

Here’s the basic procedure I will follow to create a hillshade, once Blender is installed with the BlenderGIS addon. This can function as a checklist once you’ve become familiar with the process:

2. Read the DEM into Blender as DEM raw data build
3. Adjust Z scaling (vertical exaggeration)
4. Create and adjust a georef camera
5. Correct the final pixel dimensions of the output image to match the DEM
6. Set final image type to be TIFF
7. Turn the light into a Sun, and adjust its properties
8. Do a test render
9. Do a full render

Before you experiment with using Blender to make shaded relief, you will want to try the relatively simpler step of producing your own shaded relief in GIS software like QGIS. Consequently, this tutorial won’t explain a host of things that you would learn in that process: what digital elevation models (DEMs) are, the importance of cell sizes, nodata values and projections, or the art of combining shaded relief with other layers, stretching histograms and adjusting brightness and contrast. I will assume that you already know that, at its most basic level, creating shaded relief involves specifying a sun elevation and azimuth, and selecting a vertical exaggeration for your terrain.

I use the free GDAL command-line tools gdal_translate and gdalwarp to do re-projection and re-sampling, as well as produce world files. If the command line makes you queasy, QGIS offers a graphical front-end to these tools as well. (Processing>Toolbox, and search on “translate” or “warp.”)

Your first step will be to go to https://www.blender.org/ and obtained the free animation software Blender 2.8.

## Installation of BlenderGIS

You only have to do this once, and then Blender is prepared to accept geographic data.

The BlenderGIS addon has installation instructions in its wiki at https://github.com/domlysz/BlenderGIS/wiki/Install-and-usage. Basically, they go as follows.

Go to the BlenderGIS site on github and hit the Clone or Download button, and then Download ZIP. You will receive a file called BlenderGIS-master.zip which you can store pretty much anywhere. Once you’ve installed BlenderGIS, you won’t need this file any more.

Within Blender,  go Edit>Preferences, and select the Addons tab. Click Install…, select the BlenderGIS-master.zip file, and click Install Add-on from File.

Once it is installed, be sure to check the box next to 3D view: BlenderGIS to enable it.

Two more tweaks will make Blender easier to use.

1. Note the default cube that is always there in a new workspace. To delete it, hover the mouse over the cube and hit the Delete key.
2. Find the Render tab on the right, and set Render Engine to Cycles. (If you have a good graphics card, you might also want to set Device to GPU compute)

Now go File>Defaults>Save Startup File. This means that in the future, Blender will open with no cube, and with Cycles as the render engine.

1. Note down the dimensions of your DEM
2. Convert (re-project) to a metric projection
3. Create a world file

First I will note down the pixel dimensions of my DEM, whcih in this case are 839 x 702. I will use these numbers later to tell Blender what size image to produce.

Now, the most important principle of making hillshades is that your vertical and horizontal units should be the same. If the DEM data came projected in degrees, you need to convert that to a projection measured in metres.

You’ll notice later that BlenderGIS claims that it can read a DEM projected in “WGS84 latlon” —in other words, EPSG 4326. I have found this works only when you are reading that DEM into a Blender scene already georeferenced in a metric projection.  You have to do a kind of complicated head-stand to make this occur, so we’ll take the simpler route and feed to BlenderGIS a DEM that is projected in metres. Typically for me (in northern British Columbia) this would be UTM Zone 9 North/WGS84 (EPSG 32609) or the BC Albers Equal-Area projection/WGS84 (EPSG 3005).

I will use gdalwarp to re-project, and (optionally) re-sample. In the following example, I am re-projecting out of  lat/long/WGS84 (EPSG code 4326) into UTM Zone 9 north/WGS84 (EPSG code 32609), and re-sampling to 16 metre square cells. I will assume your DEM is a geotiff file.

`gdalwarp  -s_srs EPSG:4326 -t_srs EPSG:32609 -r bilinear -tr 16 16 -of GTiff -co "TFW=YES" myDEM.tif myDEM_32609_16x16.tif`

Where

-s_srs
the Spatial Reference System (projection & datum) of the original (the “source”)
-t_srs
the SRS of the result (the “target”)
-r
resampling method (bilinear, cubic, etc.)
-tr
resolution, in metres (give two values, one for X and one for Y)
-of
output format
-co
additional creation options (TFW=YES means create a world file)

The terms “projection/datum,” SRS (Spatial Reference System), and CRS (Coordinate Reference System) are interchangeable for our purposes.

[I’m not going to suggest converting your DEM to Float32 data type. This makes a difference only with large scale mapping (e.g., 1:25.000 or larger), where you might actually see “steps” in your hillshade where the elevations change by whole metres. If you are doing such mapping, consider floating your DEM first.]

Finally, you will want a world file for the DEM, so you can georeference the hillshade Blender produces. World files used to be pretty common, but now that raster data so often comes with its georeferencing built in, I will explain what they are.

The world file was a clever invention: a six-line text file that contained the cell dimensions and the coordinates of the image origin (upper left corner). It was paired with the image file by having the exact same filename, but with a TFW extension (for TIFF images). For example, you could have myDEM.tif and its associated world file myDEM.tfw. This allowed the georeferencing information to be held externally to the image file.

When an image file has no georeferencing (i.e., it is an ordinary TIFF image, not a GeoTiff), the world file is all your GIS software needs to correctly place and scale the image. The only additional piece that it will have to ask you for is the projection of the image, in order to make sense of the world file.

For JPG and PNG images, the world file extensions are JGW and PNW, respectively.

To make a world file I use gdal_translate as follows:

`gdal_translate -co "TFW=YES" myDEM.tif deleteme.tif`

This copies myDEM.tif to a dummy file called deleteme.tif (which I will immediately delete), and in the process makes a world file called deleteme.tfw. For now I’ll rename this to myDEM.tfw. Later I’ll change its name again to match the image created by Blender.

So at the end of this preparation step I have three things:

• a DEM in a metric projection (e.g., myDEM_32609_16x16.tif)
• a corresponding world file (e.g., myDEM_32609_16x16.tfw)
• I know the EPSG projection code for the projection the DEM is in (e.g., 32609)

## Read the DEM into Blender as RawDEM

On the GIS menu, go Import>Georeferenced raster, and navigate to your DEM file.

In the right margin, set Mode to DEM raw data build (slow). (If you do not check the Build faces box, you will get a point cloud.)

You also have the option here of selecting the correct CRS (coordinate reference system) for your DEM. You need to do this only if you plan on bringing other georeferenced data into Blender to lay atop your DEM. If not, you can leave this at the default value (WGS84 latlon) even though that is incorrect.

[If the CRS you want to use is not on the dropdown menu, you can add it by clicking the “+” button, checking the Search box, typing in the EPSG code for your CRS into the Query box, and hitting Enter. It should appear in the Results box, and then you just check Save to addon preferences and click OK.]

After a pause, during which Blender builds a plane mesh out of your DEM, you will see a grey, 3D rendering of your terrain, floating in a coordinate space.

## A quick tour of Blender

Let’s take a brief detour here and learn some of the the language of Blender, as well as some of its controls and unexpected behaviours.

Blender, being 3D animation software, is a world of vertices (points in space), edges (which connect vertices) and faces (which fill in between edges). In this case, the cells of the DEM have been translated into a square array of vertices, each with the appropriate height for the DEM cell it represents. These vertices have been connected with edges that form a square lattice like the cells of the DEM. Each square of four edges has a face.

The whole set of vertices, etc. is a mesh, specifically in this case a plane mesh.

In the lower right  corner, you will notice a status bar:

This tells us that we have 588,978 vertices and 587,438 faces. My DEM is 839 columns by 702 rows, and 839 x 702 = 588,978, so I could have predicted the number of vertices. This is useful to know so you can estimate before reading in a DEM, how much RAM it will require. On average (and this varies as the size of the DEM increases), every 3,500,000 vertices requires 1GB of RAM.

The number of faces is slightly less because no faces are generated by the final row and final column (839 + 702 = 1541, minus 1 face which would be generated by both the final row and final column, hence 1540 fewer faces than vertices.)

The status bar also tells us we are presently using 291MB of RAM. This will go up during the final render. Near as I can tell, available RAM is the only limit to the size of DEMs that Blender can handle.

You can see the vertices, edges and faces of your plane mesh if you zoom in and press TAB to toggle into Edit mode.

Don’t worry if you can’t figure out how to do this just yet. If you do though, be sure to press TAB again to leave Edit mode before continuing on.

The upper right corner of Blender has a panel called an outliner that shows a tree of objects in the Blender world. At present it looks like this.

Notice that the plane (called “Thomlinson 32609 16m clip” in my case) is selected. (In the 3D view this is reflected by the fact that the plane is outlined in orange.) In addition there is a camera and a light.

Below the outliner is a Properties panel. On its left margin you will see a series of tabs identified by these icons. Hovering over them reveals their names.

We won’t have anything to do with some of these (Particles, Physics, etc.) The important ones for our purposes are Render, Output, Object and Object Data.

Note: sometimes the collection of tabs changes, depending on what object in the Blender world is selected. At this moment the plane is selected, but if the light were selected some of these tabs would disappear, and the Object Data tab would display a green light bulb icon.

The rest of the screen is taken up by the large 3D Viewport panel. Navigating in this is similar to using Google Earth, but a little different.

• Rolling the mouse wheel zooms you in and out.
• Holding down the mouse wheel button while moving the mouse in this panel allows you to rotate around objects.
• To pan, hold down Shift and the mouse wheel button while moving the mouse to either side.

The other thing to know about Blender is that some keys work only when the cursor is over the 3D Viewport.

• n toggles the 3D View Properties Panel
• Numpad -. (period) centres the view on the selected object.
• Numpad-0 toggles between this view and camera view.
• TAB toggles the selected mesh in and out of Edit mode (and we won’t be using this).

Remember that Blender is a studio for animators. It is designed to allow you to sculpt new objects, and the generate many frames of animation in which various objects are lit by lights and shot from a specific camera. We are not going to create any new objects, and are going to shoot only a single frame; so there are many controls we will not use.

However, setting up our light and camera are crucial, and that’s what we do in Part 2.

# Making shaded relief from digital elevation models (DEMs) in QGIS

Updated, April, 2021

Tom Patterson, creator of the fabulous shadedrelief.com website, once said, “There is no more important component of a map than the shaded relief.” Who would not agree?

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 3.16 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.

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 openmaps.gov.bc.ca: http://openmaps.gov.bc.ca/imagex/ecw_wms.dll?

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.

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).

### The overall process

Let’s talk about the three principles.

1. The hillshading algorithm requires a DEM in a metric projection. That means that DEMs projected in degrees won’t work: you have to re-project them first. [Although maybe sometimes not.] Unfortunately, just about all DEMs come projected in degrees, so re-projection is a must. (The Geospatial data extraction site at canada.ca, however, is one that does offer DEMs projected in metres. Choose to download your data in the “Canada Atlas Lambert (EPSG: 3979)” projection.)
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 3,000 (for screen display) or 10,000 (for printing at 300dpi) to get roughly what cell size you want. So if your map is 1:100,000, and it will be printed at 300dpi, you’d be looking for a DEM with (roughly) 10m cells.
3. In QGIS it is an option to style any DEM as “Hillshade” (other options are singleband grey, multiband colour, paletted, singleband pseudocolour, etc.) but the GDAL 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 Canada are bigger than degrees of longitude, these cells are not square on the ground. They are rectangles that are taller (north–south) than they are wide (east–west).

What you want to know is how these non-square cells measured in arc-seconds will convert to square cells measured in metres. Here are some rough ideas.

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 are all in the same area, and range in scale from 1:45,000 to 1:130,000. Looking at the chart above, this scale range suggests I can get away with using SRTM1.

So I go to Derek Watkins’s 30-Meter SRTM Elevation Data Downloader, and select each of those SRTM1 tiles.

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

### Merge and Clip

You will 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.) If you made shaded relief out of each individually, there would be some artifacts along the lines where tiles meet.

I tend to name the resulting merged 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.

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 changing the resolution of the DEM. 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.

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.

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

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:

Now if this DEM-displayed-as-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 GDAL Hillshade tool in the toolbox. (Go Processing>Toolbox and search on “hillshade.” Notice that two different hillshade tools are available: a QGIS one, and a GDAL one. The GDAL one has an interface I prefer.)

Enter your vertical exaggeration under Z factor, your sun azimuth under Azimuth of the light, and sun elevation under Altitude of the light. Generally speaking you leave scale as 1.0000, because your have re-projected your DEM and the vertical units and the horizontal units are both metres.

And here’s the result

Hillshaded made by GDAL tend to be a bit darker than we want. Their raster values run from 1 to 255, and by default this gets displayed as a singleband greyscale from pure black to pure white. But you probably don’t want pure black areas in your map, or at least not very many. Areas of solid black can’t have much information in them.

A solution that I like is to build a colour ramp that runs from 50% grey (#808080) to white (#ffffff),

and then display the hillshade as singleband pseudocolour, using this colour ramp.

Typically you combine the hillshade with the other layers of your map by assigning it a Blending mode of Multiply.

So you can get something like this.

Be aware that if you have other layers in your stack that have partial opacity, it does make a difference whether your hillshade is above or below them. In general you get the most predictable results by using a Blending mode of Multiply and placing the hillshade at the top of your stack.

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.

# A new map of the Bulkley Valley

When I attended the ICA Mountain Cartography Workshop in Banff last spring, I was encouraged by a number of people there to consider printing and selling maps of my local area. Most of this last winter has been taken up by the production of this map, which was printed in April and is now available at Interior Stationary in Smithers, BC.There are a number of goals coming together here. From a cartography perspective, the idea was a map that would feature rich shaded relief generated in Blender. From a community perspective, the idea was to produce a map that represents the amazing topography of this area, that draws you in and makes you want to go exploring.

Every map is a story, and maps, being deliberately authored, contain bias. In this case, I wanted the story to be about the terrain — the mountains, rives and lakes — and my bias was to emphasize the romance of landscape at the expense of the technical ways we divide it up. So missing from this map are all those lines that chop up a landscape: town boundaries, ski area boundaries, recreation site boundaries, private land designations, First Nations reserves and Regional District subdivisions. Even provincial park boundaries got the axe, because while on the one hand a park indicates a piece of the landscape to be preserved from development, on the other hand is the implication that outside the boundary anything goes. Park boundaries on maps suggest that only within the line is the higher quality landscape. I wanted to avoid that suggestion.

As a concession to practicality, most roads are named and you can use it as a road map. But it’s not one of those maps where every road in the provincial database is present. I left out the vast majority of forestry roads, including only those mainline roads which are commonly driven. And here too bias plays a part, because although I’ve lived in the valley for 20 years, my notion about what roads are “commonly driven” will differ from someone else’s.

As another concession I put in some point markers for those small provincial parks, or recreation sites, where a person can camp or picnic. I did this because, like roads, these are important landmarks for people.

Trails are not shown on the map, partly because at this scale, 1:150,000, you can not use it as a trail map, but also because I didn’t want to suggest that trails are an important interpreter of the landscape. When you put trails on a map they stand out as a travel network, perceived in the mind not unlike a road network. We think that where trails go must be better than where trails don’t go. I didn’t want to bias the viewer in this way.

The colouring of the map is derived from land cover data, which is essentially a satellite or aerial photo that someone has looked at and classified into zones of different stuff. Coniferous and deciduous and mixed forest are all identified, as are ice and snow, bare rock, built-up areas and shrubland. I’m a big proponent of the idea that a map should faithfully give you the colours of a landscape before you ever arrive in the area, so the colours I chose were keyed to the summer landscape here — and consequently there’s a whole second possible project of a winter map of the valley. I think the colours work pretty well, although it’s something of an idealized summer landscape. For example, there are no lingering snowpatches such as one would actually see. It is quite clear how coniferous forests dominate upper slopes while the valley bottoms are full of cottonwood and aspen. You see how south facing slopes are different from north-facing. Given that it can’t account for how light changes throughout the day, it’s a pretty good representation of what we see.

Naming is an interesting area in map-making. If you publish an article about how such and such peak should be named after your uncle Fred, people will agree with you or not, but you’ll be known as having proposed that name. However if a map-maker puts a name on a mountain or creek, no one really questions it. Maps are seen as authoritative.

But names can be seen as impositions on the landscape like other lines. Does it improve this river that we call it the Bulkley River? Would it be a better river if it was labelled the Wedzenkwe? Or is it best not labelled at all? I made a decision to go with the naming of features that reflects the British Columbia geographic names database — in other words, the names that appear on most maps and with which most map-readers are familiar — but I’m not entirely comfortable with this. For example, I’d like to see a version of this map with Wet’suwet’en and Gitxsan names on as many features as possible.

We have a lot of unnamed features in this area — or, I should say, features without official names. I hope that this map will act as a lightning rod for name information, that people will give me local names and I can add them on to future editions, and thus substantiate them.

# 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).

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.

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.

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.

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.

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.

Here’s the result:

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

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.