A Gentle Introduction to GDAL Part 5: Shaded Relief
In my previous posts on GDAL (written more than five years ago!) I covered how to open and interpret maps and images with embedded geographic information; how to transform maps from one projection to another; some of the complexities introduced working with highly detailed maps; and how to read and manipulate satellite imagery. This post and the next one will cover using GDAL for visualizing other types of data: measurements like elevation, cloud cover, city lights, and vegetation.
I’ll start with elevation, showing off some of the specialized tools GDAL has for rendering shaded relief maps, like multidirectional shading, the ability to calculate the direction a slope faces, and the ability to color-code by altitude (a function that can be appropriated for other types of data, as well).
But first, a slight digression. In the first part of this series I recommended a now-outdated method for installing GDAL. These days (late summer, 2023) I think the best way to install GDAL is through Conda. If you’re not familiar with it, Conda is a tool that helps one set up and maintain a programming environment and all the myriad libraries (like GDAL) that perform specific tasks. Installing Conda is beyond the scope of this post, but it shouldn’t be too complicated if you’re familiar with the command line on Windows, MacOS, or LINUX. If you’re just getting started I’d recommend the regular installation (which comes with lots of bundled libraries) and installing GDAL from conda-forge.
With that out of the way, go ahead and download an elevation dataset. If you don’t have one of your own handy, the following examples were all created from a high resolution digital elevation model of Mount St. Helens in Washington state. (That should be a direct link, if it doesn’t work and you didn’t have anything in mind you can find & download data from the USGS National Map, the Copernicus Land Monitoring Service, the ALOS Global Digital Surface Model, or a similar archive.) A digital elevation model (DEM) is a technical name for a file containing topographic data. You might also run into the terms digital surface model (DSM), which is the elevation of the surface as seen from above, including things like tops of trees and buildings; and digital terrain model (DTM) which is the elevation of the bare earth. Any of these variations will allow you to make a hill shaded representation of the topography, which is a great foundation to build a map on.
If you grabbed the Mount St. Helens data, you’ll have a file called USGS_1M_10_x56y512_WA_FEMAHQ_2018_D18.tif which I renamed mount_st_helens_USGS_1m_dem.tif for my own sanity.
Once you have a dataset, navigate to the directory where you have it stored in a command line window, make sure you have a Conda environment with GDAL running, and use gdalinfo to take a look at the metadata:
gdalinfo mount_st_helens_USGS_1m_dem.tif
This command will display metadata for the file, including details like the map projection, resolution, number of bands, etc.
Driver: GTiff/GeoTIFF
Files: mount_st_helens_USGS_1m_dem.tif
Size is 10012, 10012
Coordinate System is:
PROJCRS["NAD83 / UTM zone 10N",
BASEGEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
CONVERSION["UTM zone 10N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-123,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["North America - between 126°W and 120°W - onshore and offshore. Canada - British Columbia; Northwest Territories; Yukon. United States (USA) - California; Oregon; Washington."],
BBOX[30.54,-126,81.8,-119.99]],
ID["EPSG",26910]]
Data axis to CRS axis mapping: 1,2
Origin = (559993.999978157342412,5120006.000014467164874)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
LAYOUT=COG
PREDICTOR=3
Corner Coordinates:
Upper Left ( 559994.000, 5120006.000) (122d13'19.06"W, 46d13'51.54"N)
Lower Left ( 559994.000, 5109994.000) (122d13'23.64"W, 46d 8'27.18"N)
Upper Right ( 570006.000, 5120006.000) (122d 5'31.69"W, 46d13'48.09"N)
Lower Right ( 570006.000, 5109994.000) (122d 5'37.03"W, 46d 8'23.74"N)
Center ( 565000.000, 5115000.000) (122d 9'27.85"W, 46d11' 7.71"N)
Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
NoData Value=-999999
Overviews: 5006x5006, 2503x2503, 1251x1251, 625x625, 312x312Yeah, I know it’s a lot. The important things here are that the Coordinate System (map projection) is UTM Zone 10N (you can look up EPSG:26910 for more info) with units in metres (note the international spelling), Pixel Size is 1 by 1 meters (with lots of significant figures) and that Band 1 (the only band in the dataset) is Type=Float32.
Why are these details important? Firstly, GDAL needs to know how the horizontal units in the projection relate to the vertical units of the elevation data. Many projections appropriate for large scale (high resolution) maps will have units of feet or meters. This is convenient because in that case the elevation data in a file is probably going to match the coordinate system (yes, you read that right, some U.S. datasets still use feet, so their could be a mismatch between imperial and metric units). For small scale (low resolution, wide area) datasets the coordinate system is likely to use angular units (degrees) which need to be scaled to generate accurate shaded relief. For data with horizontal units in degrees (like the equirectangular projection commonly used for global topographic & bathymetric data) you’ll need to include -s 111120 in your GDAL command. (Don’t worry, I’ll show an example later.)
Secondly, the data within a digital elevation model can come in a few flavors, so it’s good to know the data type. In addition to Float32 (floating point — 32 bit data that can be positive or negative), common data types for elevation are UInt16 (unsigned integers — 16 bit data with only positive values) or Int16 (signed integers — 16 bit data with negative values representing bathymetry or land elevations below sea level). Encoding the elevation with one of these data types means there doesn’t need to be any scaling to convert from data to real-world units, and there’s usually plenty of precision to avoid banding in the data (which would look like terraces in shaded relief — yuck).
Unfortunately, floating point and signed integer data probably won’t display properly in consumer image processing applications like Photoshop or GIMP, but require something like QGIS, ArcGIS, ENVI, or the Geographic Imager plugin for Photoshop. If you don’t have any of those installed you can generate a simple grayscale image of the DEM with GDAL. You can get the minimum and maximum values from the file by including the flag -mm (“Force computation of the actual min/max values for each band in the dataset.”) when you run gdalinfo, like so:
gdalinfo -mm mount_st_helens_USGS_1m_dem.tif
This command will return Computed Min/Max=550.701,2537.783 near the bottom of the info dump.
Then convert the data into an 8-bit grayscale file with gdal_translate. Use the -scale parameter with the values for max and min in the data and max and min of the output data type in the format: -scale original_minimum original_maximum scaled_minimum scaled_maximum. You’ll probably want the minimum and maximum to span the full range of the output data, which will be 0 to 255 for a typical single-channel 8-bit image. You can scale from 0 to 65535 and use -ot Uint to convert to a16 bit unsigned integer, but your display is probably only 8 bits anyways, so why bother? Here’s the command I used:
gdal_translate -scale 550.701 2537.783 0 255 -ot Byte mount_st_helens_USGS_1m_dem.tif mount_st_helens_USGS_1m_dem_8bit_gray.tif -co COMPRESS=LZW
As a reminder, gdal_translate takes the input filename first, output filename second, and the options can got pretty much anywhere. The “creation option” -co COMPRESS=LZW losslessly compresses the data to minimize the file size while preserving 100% of the information. (Creation options are specific to the export file format, so if you’re exporting as something aside from TIFF keep that in mind.)
The resulting image looks like this:
If the data is formatted properly (the value for “no data” should be specified) and you don’t want to muck about running gdalinfo and copy/pasting some numbers by hand, gdal_translate will scale the maximum and minimum values in the data to 0 to 255 automatically. Use -scale and omit all the max and min values afterwards. You just have to remember to set the output data type to byte with -ot Byte, or else you’ll end up with data that remains in the original format, but no longer in real world units.
gdal_translate -scale -ot Byte mount_st_helens_USGS_1m_dem.tif mount_st_helens_USGS_1m_dem_8bit_gray_auto.tif -co COMPRESS=LZW
That was a lot of steps to get to a kinda plain grayscale image, but datasets will come in a variety of data types and file formats and it’s good to know your way around them. Now it’s time for the fun stuff.
GDAL includes a tool — gdaldem — that creates hill shaded images from DEMs (which can also be appropriated to apply color palettes to numerical data — but I’ll get to that later). Here’s the generic command to create a grayscale hill shade with the default options:
gdaldem hillshade input_dem output_hillshade
gdaldem is the program (equivalent to gdalwarp & gdal_translate)
hillshade is the mode (I’ll describe a few other modes later)
input_dem is the source data
and output_hillshade is the resulting image (format is guessed from the extension on the output filename — I usually create GeoTIFFs by using .tif).
Here’s the command with the Mount St. Helens elevation data:
gdaldem hillshade mount_st_helens_USGS_1m_dem.tif mount_st_helens_USGS_1m_dem_hillshade.tif -co COMPRESS=LZW
Here’s what that looks like:
That’s a plain hill shaded version of the Mount St. Helens elevation data with GDAL’s defaults illumination angles — a 315˚ azimuth (the terrain is lit as if the light source is coming from the upper left (45˚ counter-clockwise from the top of the image)) and a 45˚ altitude (the light is pitched 45˚ (midway) between the horizon and directly overhead). The 315˚ azimuth takes advantage of the human brain’s preference for interpreting shape and texture when lighting is coming from the upper left, and the 45˚ elevation angle is a good compromise that highlights both steep and shallow slopes. Note that natural illumination from these angles is impossible at this latitude! Sunlight will always be coming from the south this far north (46.2˚) of the equator.
GDAL allows you to set these angles directly for different effects (or to match the sun’s position in a satellite image) with the -az (azimuth) and -alt (altitude) parameters. Azimuth is measured clockwise from the top of the image, and altitude is measured from the horizon (0˚) to the directly overhead (90˚).
Here’s what the terrain looks like with the sun coming from the upper left (315˚), upper right (45˚), lower right (135˚), and lower left (225˚):
And here’s what it looks like at four different altitudes — 15˚, 30˚, 60˚, and 75˚ above the horizon (clockwise from upper left). Azimuth is constant at 335˚. Notice how the overall image gets brighter as the illumination source gets closer to the zenith. The brightest slopes will be perpendicular to the light source, and the darkest will be angled 90˚ or more away.
This is all well and good, but what if you want to implement some of the more advanced techniques described by Tom Patterson on the Shaded Relief site? GDAL’s got you (partially) covered. There are a few alternate shading methods available for the hillshade mode that you can select by adding a flag into the gdaldem command:
-multidirectional blends light sources from several different angles, clustered around the default 315˚. This helps make sure that linear features aren’t over- or under-estimated if they are aligned perpendicular or parallel to the light source.
-combined is “a combination of slope and oblique shading” (from the gdaldem documentation) I’m not quite sure what this is doing, but to my eye looks like it’s emphasizing texture more than slope.
-igor is a more subtle type of hillshading designed to be used in combination with additional layers of data. It’s similar to -combined but lower-contrast.
In practice, I find these different shading algorithms work best blended together in image processing or GIS software, with brightness and contrast adjustments tweaked to suit the specific map I’m making. But that’s probably the topic of another tutorial.
So far I’ve only shown data with matching horizontal and vertical units. What happens when you try to make a map with units in degrees, which is typical of most global datasets, like GEBCO combined topography & bathymetry?
gdaldem hillshade GEBCO_7200px.tif GEBCO_7200px_hillshade.tif -co COMPRESS=LZW
Oh … oh no. The Earth isn’t nearly that bumpy.
Fortunately there’s a relatively easy solution, if you know the magic number. Add -s 111120 to the command, like this:
gdaldem hillshade -s 111120 GEBCO_7200px.tif GEBCO_7200px_hillshade_scaled.tif -co COMPRESS=LZW
Better! But since the scale of the Earth is a lot bigger than the highest relief on land or in the oceans, a realistic shaded relief map at global scale is pretty underwhelming. As with -s, there’s a single-command that will adjust vertical exaggeration: -z.
gdaldem hillshade -s 111120 -z 8 GEBCO_7200px.tif GEBCO_7200px_hillshade_scaled.tif -co COMPRESS=LZW
For a global map displayed on screen, a vertical exaggeration value of 8 looks pretty good to me, but the exact value will depend on what a map is trying to show and how it will be displayed. It might even be beneficial to vary the vertical exaggeration between land and oceans to balance the apparent relief. Maps are, after all, tools, and it’s best to create maps that communicate effectively rather than maps that are 100% literal but misleading or hard to interpret.
Aside from hill shading, what else can gdaldem do with topographic data? This is where the different modes come in. Slope and aspect are relatively straightforward, in my opinion.
Slope is the steepness of the terrain, by default calculated in degrees with a range of 0˚ — 90˚, but can be set to percent with the -p flag, in which case the values will be go from 0% (flat) to up to 100% (vertical).
Aspect is the direction the slope is facing. It ranges from 0˚ (directly north) clockwise — 45˚ is northeast, 90˚ east, 135˚ southeast, etc. Keep in mind that the numbers wrap, so a slope of 359˚ is oriented mostly north, angled just a tiny bit to the west.
Notice that both slope and aspect return floating point results, not 8-bit (grayscale) like the images from hillshade. Here’s what maps of these two parameters look like (at least after they’ve been scaled and converted from data to imagery):
So far the examples I’ve shown have been a little, well gray. What if I wanted to bring a little flair? A bit of pizzazz? A dash of color? Color-relief will convert a plain old DEM into a map with hypsometric tints, cartographic jargon for color-by-numbers. With color-relief you’ll need two input files. The digital elevation model and a text file that tells the algorithm which colors go with which numbers. The text file is formatted like so:
elevation red green blue alpha(optional)Put the elevation in the first column, then the value for red, then green, then blue, with an optional column for alpha (transparency). Elevation can be an integer or floating point value, and each color component is an integer from 0 to 255 (one byte’s worth). Alpha also ranges from 0–255, with 255 being fully opaque. (The default is opaque if there’s no value for alpha.) Column separators are fairly flexible — commas, spaces, tabs, or colons. Here’s an example, modified from the “cold humid” color palette in The Development and Rationale of Cross-blended Hypsometric Tints, by Patterson & Jenny:
500 112 147 141
600 120 159 152
700 130 165 159
900 145 177 171
1100 180 192 180
1300 212 201 180
1500 212 184 163
1800 212 193 179
2100 212 207 204
2400 220 220 220
3000 235 235 237
4000 245 245 245Here’s how to run the generic version of the command:
gdaldem color-relief input_dem color_text_file output_color_relief_map
Substituting the elevation data for Mount St. Helens and my custom color palette file (all on one line, Medium is wrapping the text):
gdaldem color-relief mount_st_helens_USGS_1m_dem.tif st_helens_hypsometric.txt mount_st_helens_USGS_1m_cold_humid_hypsometric.tif -co COMPRESS=LZW
Which yields a color image (below, left). By itself it’s not all that interesting — it may even be less informative than a plain grayscale image of the mountain’s elevation (like the image towards the top of this post). But combined with shaded relief, even a simple hillshade with a 335˚ azimuth, the map starts to come alive (below, right).
Note: I made the combined elevation/shaded relief image by using blend modes in QGIS, accessible from the Layer Rendering options. The overlay blend mode applied to the hillshade layer makes areas of the layer below lighter when brightness is above 50% and darker when brightness is below 50%. As a result, slopes facing towards the illumination source get lighter, while those facing away get darker.
Color-coded elevation data plus a grayscale hillshade are the fundamental components of a shaded-relief map. But there are a wealth of additional techniques cartographers use to illustrate topography — and it’s possible to replicate some of these advanced techniques with the tools provided by GDAL.
For example, Swiss style cartography and other elegant topographic maps often employ subtle shifts in hue that emulate the contrast between directly sunlit surfaces — which exhibit warm (yellow) hues — and shadowed surfaces — which are tinted with cool (blue) colors. Conveniently enough, the output from GDAL’s aspect mode has the information needed to replicate this effect, and the color-relief mode can be applied to any numerical data — not just elevation!
As I mentioned earlier, a slope’s aspect is merely the direction it faces. So if one color-codes aspect so that slopes facing the illumination source are given a slight yellow tint, and slopes facing away from the illumination are colored a bit blue, the resulting map takes on a nice touch of realism.
Just like creating colors from elevation data, both a data file (in this case floating point values generated from the digital elevation model with the aspect mode of gdaldem) and a text file (that specifies what colors go with which values) are necessary to run color-relief.
The tricky part is that aspect varies from 0˚ to 360˚, wrapping around so that both 0˚ and 360˚ should be represented by the same color, or else there will be an obvious discontinuity east and west of north. (I mapped the elevation data sequentially from low to high, so those colors are intentionally distinct.) With a typical direction of illumination of 315˚, the slope colors should look something like this:
The variation in hue is symmetric, so that colors going clockwise from 315˚ to 135˚ are the same as those going counterclockwise. One thing to note is that the lightness of all these colors should remain constant while hue and saturation vary — so I used the HCL Wizard Palette Creator to make a color ramp. (If you want to learn more about human perception, color models, and how they apply to data visualization you can read Lisa Charlotte Muth’s superb posts on color on the Data Wrapper blog, or my own series Subtleties of Color.)
Here’s what the text file associating colors with aspect angle looks like:
0 243 241 231
22.5 242 241 235
45 241 241 239
67.5 240 241 243
90 239 241 246
112.5 238 240 250
135 236 240 254
157.5 238 240 250
180 239 241 246
202.5 240 241 243
225 241 241 239
247.5 242 241 235
270 243 241 231
292.5 245 241 227
315 246 241 222
337.5 245 241 227
359.999 243 241 231We’re not all that great at distinguishing hues, especially desaturated (pale) hues, so the raw output (below, left) can be underwhelming. It looks better (and is much more interpretable) with hillshade (below, right).
GDAL’s final three modes are on the technical side, but they can add some subtle details to shaded relief maps so I think they’re worth going over. TRI (Terrain Ruggedness Index), TPI (Topographic Position Index), and roughness are all measures of the variability in elevation in the area immediately surrounding each pixel. Terrain Ruggedness Index and roughness are so similar as to be almost indistinguishable, even viewed side by side. So I’ll just show roughness (below, left). Note that I’ve inverted the image, so the roughest areas are black and the smoothest white. Topographic Position Index is at least different, but mostly gray with a few lighter & darker details where the topography is very complex.
The implementations of these algorithms in GDAL are all limited to a very small region (the 9 pixel square formed by a central pixel and the 8 pixels immediately surrounding it) so other software may be more useful if you need these parameters for analytic work. But the outputs share some similarities with the ambient occlusion and texture shading techniques that are becoming increasingly popular in cartography — so I think they’re worth experimenting with.
The most elegant and informative maps don’t typically rely on a singe technique. Rather, they are the end result of multiple elements blended together. Some obvious, some subtle. Here’s my map of Mount St. Helens, incorporating six different layers: roughness, hypsometric tints, Igor hillshade, multidirectional hillshade, simple hillshade with illumination from 315˚, and aspect colored yellow and blue. Pretty intricate coming from a single dataset!
Now that I’ve gone through all this you may be wondering — why bother with GDAL? Can’t all this be done in QGIS, or ArcGIS, or any number of other very capable applications? Well, yeah, it can (and QGIS even uses GDAL under the hood). But GDAL is wonderfully scriptable. I converted most of the web-friendly PNGs from high-res TIFFs with gdal_translate:
for file in *.tif
do
gdal_translate -outsize 1200 0 -r bilinear -of PNG $file ${file%.tif}.png -co worldfile=YES -co zlevel=8
doneAnd you can even script the input parameters, like this:
for azimuth in {1..360}
do
gdaldem hillshade -az $azimuth ../mount_st_helens_USGS_5m_dem.tif 00${azimuth}_mount_st_helens_USGS_5m_hillshade.tif -co COMPRESS=LZW
doneWhich outputs 360 files with illumination azimuths coming from 1–360˚. Perhaps a bit of a gimmick here, but potentially useful. On the other hand, maybe you’d like to create custom multidirectional hillshading, automate the production of dozens of maps, or show a client how light would vary in a park over the course of a year — all possible with a bit of creativity.
In my next post I’ll show how to use color-relief with other datasets, that have nothing at all to do with elevation or hillshading. Plus some tips on doing calculations, reading scientific formats like HDF, and integrating GDAL with Python. With any luck it’ll take fewer than 6 years!
- A Gentle Introduction to GDAL
- Map Projections & gdalwarp
- Geodesy & Local Map Projections
- Working with Satellite Data
- Shaded Relief (you are here)
- Visualizing Data
- Transforming Data
- Reading Scientific Data Formats
Thanks go to Frank Warmerdam, who got me over my fear of GDAL, and Tom Patterson, for sharing his deep knowledge of mountain cartography, and Bernice Rogowitz, for opening my eyes to the possibilities of color.
