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.

And, look: shaded relief:

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.

How bad is it?

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.

QGIS 3 and Vector map tiles

If you’ve ever tried printing a map from the Open Street Map tile service (available in QGIS through QuickMapServices plugin) you’ve experienced the annoyance of having all your labels come out tiny. In other words, you see this on the screen…

But then you get this out of your 300 dpi print layout….

The reason those labels and icons come out so small with raster tiles is that your screen display is something like 96 pixels/inch, but the print layout has to pack them down to 300 pixels/inch. The print layout also winds up asking for tiles at a higher zoom level. Text labels and icons, like everything else, are “baked” into the raster tiles at a reasonable size for display at 96 dpi. At 300 dpi, everything is much smaller.

In 2018 I described, at the end of a post, a way to work around this, but now there’s a better option. If you work instead with vector tiles, the text labels are not baked into the tile: they are rendered by QGIS. This means that if the label is supposed to be 2.5 mm high, it’ll be 2.5 mm high no matter what dpi you export at out of the print layout.

QGIS 3.14 and above can handle vector tile servers as data sources. This means the server sends tiles of vector data, which are then styled at the user’s end by QGIS.

How they work

Vector map tiles (VMTs) were not developed for GIS use. Their main use is for people viewing maps through browsers. You commonly have a map that’s been embedded in a web page that includes the javascript libraries for mapbox rendering, such as here.

But vector tiles offer the cartographer a lot. As raster tiles, OpenStreetMap data was ideal for areas where you just could not get decent data—if, say, you were mapping the streets of a city in Turkmenistan. But now, as vector tiles, this data allows you to control colours, line weight, font, and, of course the visibility of elements. Imagine your own version of OpenStreetMap, with pink highways, no railroads and emphasized parks, and you begin to get the idea.

There are multiple formats for tiles, but here we’ll just look at the Mapbox vector tiles format, since that is what QGIS can support.

Because vector map tiles are designed to be viewed over the internet, the most basic principle of VMTs is that you typically to have a tile server, which is firing off vector tiles of data, and a GL style file that understands what data layers are delivered in this server’s tiles, and tells the browser how to render that data. Often these are at different URLs.

The URLs of vector tile servers look like this:

The z, y and x in curly brackets get replaced by the zoom level, row and column, and the data comes as a PBF (Protobuf) file or a MVT (Mapbox vector tile) file.

The GL style file, on the other hand, is a text file in JavaScript Object Notation (JSON). It might begin like this, with information on the sources for tiles, sprites and glyphs (VMT lingo for markers and fonts):

{
  "version": 8,
  "name": "Positron",
  "metadata": {...  },
  "sources": {
    "openmaptiles": {
      "type": "vector",
      "url": "https://api.maptiler.com/tiles/v3/tiles.json?key={key}"
    }
  },
  "sprite": "https://openmaptiles.github.io/positron-gl-style/sprite",
  "glyphs": "https://api.maptiler.com/fonts/{fontstack}/{range}.pbf?key={key}",

and then it will continue with a list of the layers in the tiles, and how they should be styled.

  "layers": [
    {
      "id": "background",
      "type": "background",
      "paint": {"background-color": "rgb(242,243,240)"}
    },
    {
      "id": "park",
      "type": "fill",
      "source": "openmaptiles",
      "source-layer": "park",
      "filter": ["==", "$type", "Polygon"],
      "layout": {"visibility": "visible"},
      "paint": {"fill-color": "rgb(230, 233, 229)"}
    },
    {
      "id": "water",
      "type": "fill",
      "source": "openmaptiles",
      "source-layer": "water",
      "filter": [
        "all",
        ["==", "$type", "Polygon"],
        ["!=", "brunnel", "tunnel"]
      ],
      "layout": {"visibility": "visible"},
      "paint": {"fill-antialias": true, "fill-color": "rgb(194, 200, 202)"}
    },     ...

The styling information here is written in something called Mapbox GL style, and here are some GL Style file URLs:

Vector tiles are relatively new, and almost no one gives them away for free, but we do have a few heroic companies offering this. Most require API-keys or access-tokens, but these are frequently free.

The first question you might have is what you get if you ask for vector tiles from the tile server, but don’t combine them with a GL Style file. VMTs without a style file aren’t useless—but what you get is data where all points share a single style, as do all lines and polygons.

Data from the Nextgen vector tile server without styling

A simple VMT service brought in using QGIS3’s native abilities

Our best candidate at the moment for free OpenStreetMap (OSM) data delivered through vector tiles is the ESRI OSM vector tile layer. (Thank you ESRI!). The tile server is here:

https://basemaps.arcgis.com/arcgis/rest/services/OpenStreetMap_v2/VectorTileServer/tile/z/y/x.pbf

and the style file is here:

https://www.arcgis.com/sharing/rest/content/items/3e1a00aeae81496587988075fe529f71/resources/styles/root.json?f=pjson

We can add this tile server to QGIS by

  1. opening the Browser panel
  2. Right-clicking the new Vector Tiles category and choosing New ARCGIS vector tile service connection
  3. Filling it out like this:

Notice that the Service URL ends with “VectorTileServer.” There’s no “tile/{z}/{y}/{x}.pbf.” QGIS will append that automatically before making a tile request.

When we add this layer to our map canvas, and zoom to Baku, Azerbaijan at 1:8,028 scale, we get this pleasing map that looks remarkably like OpenStreetMap:

However, we also get a pile of errors, under the banner “Vector tiles: Style could not be completely converted.” If you click on “Details,” some will be like this:

Referenced font Arial Unicode MS Regular is not available on system

You can repair these font errors by installing the missing font – or simply accepting the local font that has been substituted.

But there are also a pile of errors along the lines of

amenity area/bowling alley: Could not retrieve sprite 'amenity area/bowling alley'

I will explore how to solve these sprite problems below. First let’s talk about what you’ve got so far.

Styling the vector tile layers

If I open the styling panel on the right of the map I get something truly new.

Wow, that’s the longest list of rule-based styling I’ve ever seen!

Can I literally uncheck or change the style for every type of layer? Are these layers within a layer?

Yes and yes. Vector map tiles contain many layers. (The set of layers is called the schema of a tile service.) Out of each layer (e.g., “landcover”), features are selected using a filter, which is quite similar to a rule in rule-based styling. For example, in the image above there is a specific fill style for features in the landcover layer for which _symbol = 48,and the zoom is 9 or higher. It’s called landcover/park.

I can uncheck landcover/park and get a map without parks…

Or I can double-click on the sample fill swatch and change the parks to a pink/white gradient fill.

In fact I can apply any QGIS style to these polygons: shapeburst fill, line fill, marker fill… I don’t have full control—I can’t change the min and max zooms that define this style, or the filter rule—but within the predefined filters, the possibilities are mind-boggling.

Once I like the way in which I have tweaked the styling, I can save the style for this layer to a .qml file. (Open the layer’s Properties, find the Style button, and choose Save Style…) This will allow me to apply this custom styling to the same vector tile layer, but in another QGIS project.

My powers over the styling of the VMTs are clumsy—a GL Style file typically defines scores of styles, and if I want to change, say, all my water styles from blue to grey, I will have to do that to each one individually. But it’s still very powerful.

There are a few other things to point out. Note at the bottom of the styling panel you have a checkbox for Visible rules only. Visible rules only cuts out styles that are not displayed at the current zoom. (It does not remove styling for features not visible in the current map canvas.)

The current zoom is also displayed there. This is always shown as an integer, but with vector map tiles you are not restricted to integer zoom levels. There’s no “Zoom to Native Resolution (100%)” like there is when you right-click a raster tile layer. The zoom level shows as 14, but the map canvas is at 1:10,000, which is somewhere between 1:17,061 (the break point for vector tiles changing from zoom 13 to 14) and 1:8,530 (the break point for vector tiles changing from zoom 14 to 15). If you do the math, you can conclude we’re at something like zoom level 14.17. This might make you queasy! And, welcome to the world of vector tiles.

[Note that zoom levels for vector tiles are one more than the zoom levels for typical raster map tiles (e.g., Google Maps, OSM ). This is a consequence of MapBox vector tiles being 512 pixels on a side. As well, I have a suspicion that, in QGIS, the detail shown on ESRI OSM vector tiles at zoom level 14 is the same as shown on OSM raster tiles at zoom level 13. But I’m not quite sure why this is.]

Styles are implemented in order, from the top of the list down, and you can drag them to reorder.

Finally, it’s important to understand that the styling rules you see here in QGIS are not the GL Style document itself. They are a translation of that document. You cannot (at, least, as far as I know) translate back from QGIS styling to a GL Style file.

Labels

If you move to the labelling tab in the style panel, you see this:

So, you also have extensive control over the labelling. Double click on a label’s style swatch, and change the font, size, colour… all of the usual label parameters are available.

For example, if I would like street names to be blue (bold-italic) and subway stations to be labelled in red, I can just do it by altering three styles.

Finally, note that while you can’t open an attribute table on your vector map tiles layer, you can use the Information tool. You will get results for every single layer under the mouse click.

I can now print at 300 dpi from this, and I get just what I see on the screen!

Sprites

All that might have been enough for you, and now you’re off to try customizing OpenStreetMap yourself. However, if you’re interested in fixing those sprite errors we encountered above, continue reading.

When you add a vector tile service to QGIS you frequently get this:

When you check out the Details, many of the errors refer to “could not retrieve sprite.”

At the same time, in the QGIS’s Network Log window, you’ll see something like this:

2020-12-31T20:27:24     WARNING    Error transferring https://cdn.arcgis.com/sharing/rest/content/items/684018c92e9b424ca8e900da8dad8b23/resources/styles/../sprites/sprite@2x.json - server replied: Not Found
2020-12-31T20:27:24     WARNING    Error transferring https://cdn.arcgis.com/sharing/rest/content/items/684018c92e9b424ca8e900da8dad8b23/resources/styles/../sprites/sprite.json - server replied: Not Found

So what’s going on?

If I paste the URL for the GL Style file into my browser (I’m using Firefox), the result is that whole code for the styling gets delivered to me in an unformatted blob.

It’s a bit dense reading. But near the top I can pick out this piece:

"sprite": "https://cdn.arcgis.com/sharing/rest/content/items/3e1a00aeae81496587988075fe529f71/resources/styles/../sprites/sprite"

Sprites, in vector map tile parlance, are markers. QGIS needs to use this URL look up the sprite image (a PNG file with all the sprites packed into it) and its partner sprite JSON file (which explains which piece of the image a given sprite is). But, for QGIS, the “styles/..” part of this particular URL is confusing. We might understand it as meaning descend into the styles sub-directory and then pop back up out of it, but QGIS doesn’t know this.

I can make a better address for the sprites by saving the GL Style file locally, editing it so the URL of the sprites is simpler, and then telling QGIS to now use my local version of the GL Style file.

I save the GL Style file locally by selecting all the text out of my browser, and copying and pasting into my text editor. Then I change the sprite URL by removing the “styles/..” part, so it now reads

"sprite": "https://cdn.arcgis.com/sharing/rest/content/items/3e1a00aeae81496587988075fe529f71/resources/sprites/sprite"

Once I’ve saved the file as ESRI_OSM_Style.json, I return to QGIS, where I edit the Vector Tiles connection I made for this server, and specify this new file as the “Style URL.” I could post it on a server somewhere and give a http:// URL to it, but, with Ubuntu, I can also just use the file:// protocol:

Sprites make quite a difference at high zoom levels.

No sprites, zoom level 17
With sprites, zoom level 17

A few final points

There are three different types of support for vector tiles in QGIS. The original one was the Vector Tiles Reader plugin, which received the vector tiles, converted them to geoJSON data, and made a different layer for each kind of features in the tiles. There is the new native support, which I have used here: it receives the tiles and applies filters, or rules, to style the many different kinds of features in them. Finally, there is the MapTiler plugin, which has a great deal of potential but requires a free account with MapTiler/OpenMapTiles, and this has an unfortunate data cap of 100,000 tile requests/month. (It seems like a lot, but while learning to use and style vector map tiles, I went through my 100,000 tiles in eleven days.)

The QGIS implementation of vector tiles permits you to add “New generic connection” if you have URLs for the server and the GL Style file. Yet outside of MapTiler and ESRI I have found precious few servers that also provide GL Style files. Here’s a quick survey:

  • The city of Wien (Vienna)—the data is limited to Austria. [Server URL] [GL Style URL]
  • Ordnance Survey UK—the data is limited to the U. K. and requires a free API key. [Server URL to which you have to append your free API key] [GL Style URL to which you have to append your free API key] Because of the API-key, QGIS does not yet seem to be able to make the request for the sprites correctly.
  • MapTiler/OpenMapTiles (OSM data)—beware the 100,000 tile limit. (Literally every time you zoom, pan or resize the screen you can trigger a tile request.) [Server URL to which you have to append your free access key] [GL Style URL for the “Basic” style]
  • Nextzen (OSM data)—lovely data, but I haven’t found a GL Style file. [Server URL to which you append your API key.]
  • ESRI generously offers another vector tile server based on more data sources than just OpenStreetMap. [Server URL] It comes with many different styles which you can find here. To obtain the GL Style file for one of these, you typically Open the style in Map Viewer and then in the upper right is a button called View style. Copy the link out if it and use it, as we did above, in the Style URL field. (Generally speaking, you’ll have to also do the same hack with the sprites as I showed above.)
  • Mapbox, if you can decode their mysterious URLs (mapbox://mapbox.mapbox-streets-v11), serves up vector tiles, and the public access key is free (without data limits), but I have not been able to find a GL Style file for its schema that QGIS can fully translate. [Server URL for Mapbox-streets-v8] [GL Style URL for MapBox-streets-v11]
  • Thunderforest has nice vector tile servers, such as their Transport layer, and their Outdoors layer, that you can have access to with a free API key—but I have not been able to find GL Style files for the schemas of these servers. [Server URL to which you have to append your free API key]

I understand that converting data to vector tiles, and keeping it up to date, costs money, so I don’t begrudge these companies charging fees or imposing limits. On the other hand, I really appreciate those which don’t!

As time passes there will no doubt be more vector map tile servers, and the way schemas and styling work will change. Nonetheless working with vector map tiles is addictive, and it will be exciting to watch this evolve.

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.

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

Resolution (degrees)Resolution (metres)SourceLimitationsRecommended scale (screen/print)Notes
0.75″ or 0.0002083333333°20CDEM or CDSM at https://maps.canada.ca/czs/index-en.htmlCanada only1:60,000/1:200,000Alternately, can be downloaded in a metric projection
1″ or 0.0002777777°30SRTM1 at https://dwtkns.com/srtm30m/60° South to 60° North1:90,000/1:300,000Comes in 1° x 1° tiles
3″ or 0.000833333°90SRTM3 at https://dwtkns.com/srtm/60° South to 60° North1:270,000/1:900,000Comes in 5° x 5° tiles
15″ or 0.00416667°450SRTM15+ at https://topex.ucsd.edu/WWW_html/srtm15_plus.html1:1,300,000/1:4,500,000The world is divided into 33 tiles.

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.

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.

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:

Make a static hillshade

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

Displaying hillshades

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.

Two common solutions are to adjust the opacity of your hillshade, or to increase its brightness.

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.

What I like about this is that I then know that no part of my hillshade is darker than 50% grey.

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.

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.