Shaded Relief Maps using GDAL and Open Data
A Beginner’s Guide
Recently I embarked on a journey to learn more about GIS. The goal was to learn how to process rasterized GeoTIFF files in order to build a shaded relief map that can be used to contribute to OpenStreetMap. I decided to compile my notes into a single article in the hopes that it can be of use to someone else. Please note that I am not a GIS specialist.
✉️ If you spot a mistake or have a suggestion, feel free to contact me! My Twitter handle is @dmoraisferreira.
Why bother with relief maps instead of just using orthoimagery?
Mapping trails in wooded areas can be a challenging task, particularly when there are no GPX tracks and no reliable datasets that can be referenced. Orthoimages, or orthophotos, are aerial pictures of a terrain that are georeferenced. These images are two dimensional and only captures the features with the highest elevation. Therefore, anything covered by trees or other obstacle will not be visible. In contrast to the orthorectified imagery, LiDAR uses laser light to illuminate the terrain and measures the reflections to build a three-dimensional point dataset. This laser light is able to pass through tree canopy and reflect on the ground surface back into the plane-mounted sensors. The resulting points are classified and can be converted into a Digital Elevation Model (DEM), from which a shaded relief map, or hillshade, can be computed. The higher the resolution of the point cloud is, the better the resolution of the DEM will be. As the LiDAR data points are classified, it is possible to remove classes such as buildings or trees before generating the DEM. Depending on which classes are chosen, the resulting DEM will be referred to as Digital Terrain Model (DTM) or as a Digital Surface Model (DSM).
For our purposes, DTMs capture the bare-earth elevation, that is the elevation of the ground without trees or buildings, while DSMs capture the elevation of all features, including buildings and pylons. Generally, trees are never included, as this would defeat the purpose of modelling the elevation of large wooded areas.
Furthermore, hillshades can reveal information that is not easily visible to the naked eye. Guillaume Rischard tweeted a great example: the old dismantled fortress walls in Luxembourg City can still be detected in a hillshade. Mark Walters tweeted a couple of interesting visualisations of Gallo-Roman villas. Rouven Meidlinger tweeted a very interesting video that shows the Celtic oppidum on the Titelberg.
📚 Recommended reading: DEM, DSM & DTM Differences — A Look at Elevation Models in GIS, GISGeography.
Prerequisites
You will need to install the GDAL library on your computer. There are a couple of ways to achieve this:
- install GDAL
- install QGIS, which comes with GDAL installed
- docker pull osgeo/gdal
- If you cannot make persistent changes, use OSGeoLive, a self-contained bootable Linux distribution
For this tutorial, I will be using QGIS 3.12.1 and GDAL 3.0.4.
📚 Recommended reading: The GDAL documentation pages for the different programs.
Step 0: Data
In 2017, the Administration of Air Navigation (ANA) in Luxembourg released a high-resolution DTM that was computed from the first LiDAR survey flight in the same year. As of January 2020, the raw LiDAR dataset from the 2019 flight has been published on the OpenData Portal. This LiDAR dataset is still in its raw form and will have to be converted into a DEM first. This project is currently being worked on by Guillaume Rischard.
Download the data
wget https://data.public.lu/fr/datasets/r/b21141d7-afc3-42a1-9c65-4587169ef3a7 -O ANA_LUREF_NGL_DTM.zip
unzip ANA_LUREF_NGL_DTM.zip -d ANA_LUREF_NGL_DTM
cd ANA_LUREF_NGL_DTM
find Zone* -type f -iname ‘*.tif’ > input-files-epsg2169.txt
What are the individual files?
- tif: File format for storing raster graphic images
- prj: Contains the Well-known text (WKT) for the GeoTIFF file projection
- tfw: TIFF World File, contains translation of image coordinates to real-world coordinates
What do the filenames stand for?
The filenames are labelled using a grid in EPSG:2169 (LUREF) coordinates. You can find a visualization of this grid by clicking here.
What are we working with?
david@desktop$ gdalinfo Zone1_Zone2/MNT_Luref_NGL_1m_100000_101000.tif
Driver: GTiff/GeoTIFF
ERROR 1: PROJ: proj_create_from_database: datum not found
ERROR 1: PROJ: proj_create_from_database: prime meridian not found
Files: Zone1_Zone2/MNT_Luref_NGL_1m_100000_101000.tif
Size is 1000, 1000
Coordinate System is:
PROJCRS["Luxembourg 1930 / Gauss",
BASEGEOGCRS["unknown",
DATUM["unnamed",
ELLIPSOID["International 1924",6378388,297.000000028404,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",49.8333333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",6.16666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",80000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",100000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",2169]]
Data axis to CRS axis mapping: 1,2
Origin = (100000.000000000320142,100999.999999999155989)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
AREA_OR_POINT=Point
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 100000.000, 101000.000) ( 6d26'40.92"E, 49d50'31.17"N)
Lower Left ( 100000.000, 100000.000) ( 6d26'40.73"E, 49d49'58.80"N)
Upper Right ( 101000.000, 101000.000) ( 6d27'30.97"E, 49d50'31.04"N)
Lower Right ( 101000.000, 100000.000) ( 6d27'30.77"E, 49d49'58.68"N)
Center ( 100500.000, 100500.000) ( 6d27' 5.85"E, 49d50'14.92"N)
Band 1 Block=1000x1 Type=Float32, ColorInterp=Gray
NoData Value=-32767
Unit Type: metreThis tells us that our projection is EPSG:2169, that we have one grayscale raster band and that our data is compressed using LZW, a loss-less data compression algorithm. The image structure metadata also describes what data format is used to store pixel and band information. You can read more about multi-band representation here.
Furthermore, gdalinfo also tells us what the pixel size is. This is useful when deciding on what resolution parameter to choose when processing the dataset.
The two PROJ errors are raised by the GTiff driver, read more here. This is no big deal, as we will be warping the data later anyways.
If you want to get rid of the error messages, you can execute the following command:
gdal_translate -a_srs “EPSG:2169” -co “TFW=YES” input.tif output.tif
Step 1: Build a Virtual Dataset (VRT)
Since our dataset consists of many smaller TIF files, we should first combine then into a mosaic using gdalmerge or gdalbuildvrt. While the former generates a single large TIF file, the latter generates a VRT file that only contains references to our input files.
What resolution to use?
As the output of the following command shows, the pixel size is identical for all our TIF files. Therefore, we can use the default “average” resolution.
david@desktop$ < input-files-all.txt xargs -I{} -d'\n' gdalinfo {} | grep "Pixel Size" | sort | uniq -c
1 Pixel Size = (0.999999999999992,-1.000000000000000)
3840 Pixel Size = (1.000000000000000,-1.000000000000000)
78 Pixel Size = (1.000000000000007,-1.000000000000000)
1 Pixel Size = (1.000000000000012,-1.000000000000000)What resampling algorithm to use?
We want to preserve thin lines as much as possible, therefore it would be best to select resampling algorithms that don’t distort the pixels too much. For our purposes, I choose nearest. The lanczos algorithm would also have been a good choice, as this algorithm preserves local contrast.
gdalbuildvrt -r nearest -a_srs “EPSG:2169” ANA_LUREF_NGL_DTM_epsg2169.vrt -input_file_list input-files-epsg2169.txt
📚 Recommended reading: Resampling methods useful in a spatial context?
Step 2: Generate the Hillshade
There are a handful of parameters that can influence the final appearance of the hillshade. I ended up using the following parameters, as this seems to result in a hillshade that is useful for our purposes.
gdaldem hillshade ANA_LUREF_NGL_DTM_epsg2169.vrt ANA_LUREF_NGL_DTM_epsg2169_hillshade.tif -co BIGTIFF=YES -co TILED=YES -co COMPRESS=DEFLATE -of GTiff -z 1.0 -s 0.5 -multidirectional
Vertical Exaggeration Factor z
The elevation values are multiplied with the z-Factor. A higher value will result in hillier terrain.
📚 Recommended reading: Don’t be shy — exaggerate your maps!
Example A: Flat Terrain
As can be seen for z=4, the features really stand out. Unfortunately, this leads to very dark slopes in hilly terrain (Example B).
Example B: Hilly Forest
By exaggerating the scales, the shades become much darker. This can obscure small trails. As can be seen below, z=4 produces a very dark slope.
📚 Recommended reading: Setting the Z Factor parameter correctly.
Scale Parameter s
This parameter defines the ratio of vertical to horizontal units. High values result in a flatter hillshade.
By changing the default s=1 to s=0.5, we gain a lot of detail in both examples.
Example A: Flat Terrain
Example B: Hilly Forest
📚 Recommended reading: The gdaldem documentation page.
Unidirectional and Multidirectional Light Sources
The light azimuth and altitude can be defined by using -az and -alt. When using a single light source, most programs default to an azimuth of 315 degrees and an altitude of 45 degrees. According to this thread, this is due to how the human eye perceives shadows.
Multidirectional lighting is also supported by GDAL. In the following example, four light sources were added from different angles (225, 270, 315 and 360 degrees).
Step 3: Slope Map and Shading
gdaldem slope ANA_LUREF_NGL_DTM_epsg2169.vrt ANA_LUREF_NGL_DTM_epsg2169_slope.tif
touch slope-ramp.txt && printf ‘%s\n%s\n’ ‘0 255 255 255’ ’90 0 0 0' >> slope-ramp.txt
david@desktop$ cat slope-ramp.txt
0 255 255 255
90 0 0 0The slope ramp maps terrain angles to colors. The first column represents the angle and the second, third and fourth column represent red, green and blue respectively. Flat terrain (0 degrees) will be white, while 90-degree slopes will be black. Values in-between those boundaries are interpolated.
gdaldem color-relief ANA_LUREF_NGL_DTM_epsg2169_slope.tif slope-ramp.txt ANA_LUREF_NGL_DTM_epsg2169_slope-shade.tif
This is what the resulting slope shade looks like in normal blending mode:
Adding the slope shade to the hillshade, we can see that the much brighter slopes regain some information.
📚 Recommended reading: Blend modes.
Step 4: Color Relief
Adding color to a grayscale relief makes hills and valleys stand out even further.
Choosing an appropriate color palette
When working on the hillshade, I tried out a couple of different color palettes that were all based on the rainbow palette 🌈. However, this video convinced me to use the viridis color palette. These palettes are easier to read by people who are color blind. They are also perceptually more uniform.
While I created my color ramp using the QGIS wizard, it is always useful to know the highest and lowest point of a country. For Luxembourg, the lowest point is 129.9 meters and the highest point is 559.9 meters.
touch viridis-gdal.txt && printf ‘%s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n’ ‘133.863,68,1,84,255,133.86’ ‘242.697,57,84,139,255,242.69’ ‘297.115,42,117,142,255,297.11’ ‘351.532,31,149,140,255,351.53’ ‘405.95,47,180,124,255,405.94’ ‘460.367,110,206,88,255,460.36’ ‘510.599,188,223,39,255,510.59’ ‘552.458,253,231,37,255,552.45’ >> viridis-gdal.txt
david@desktop$ cat viridis-gdal.txt
133.863,68,1,84,255,133.86
242.697,57,84,139,255,242.69
297.115,42,117,142,255,297.11
351.532,31,149,140,255,351.53
405.95,47,180,124,255,405.94
460.367,110,206,88,255,460.36
510.599,188,223,39,255,510.59
552.458,253,231,37,255,552.45The format is similar to that of the slope ramp, the first column represents the elevation, the second, third and fourth column represent red, green and blue respectively, the fifth column represents the alpha (opacity) and the sixth column is a text label.
gdaldem color-relief ANA_LUREF_NGL_DTM_epsg2169.vrt viridis-gdal.txt ANA_LUREF_NGL_DTM_epsg2169_color-relief.tif
Note that you can also generate a VRT file instead of a TIF, this has the advantage of being a lot faster to compute and doesn’t generate a large TIF file.
gdaldem color-relief ANA_LUREF_NGL_DTM_epsg2169.vrt viridis-gdal.txt ANA_LUREF_NGL_DTM_epsg2169_color-relief.vrt -of vrt
If you are feeling patriotic, you could always use the Flag of Luxembourg ;)
I also experimented with a more traditional color palette. While it does look like the “official” map I grew up with, I still find the viridis palette more useful for OpenStreetMap.
Step 5: Visualization
Color Relief
In QGIS, we add the color relief and the hillshade layers. The color relief blending mode is set to multiply and the opacity is set to 85%.
Color Relief and Slope Shading
By further multiplying the slope shade onto out color relief and hillshade, the depth perception further increases. For OpenStreetMap mapping, it seems that adding a slope shade can add detail to trails on the slope.
Step 6: Warping and Compression
Serving EPSG:2169 can be suboptimal, as most, if not all, users of this relief map will be using EPSG:3857. By using gdalwarp, we can warp, i.e. reproject, and compress the hillshade in one step.
gdalwarp -t_srs epsg:3857 -r lanczos -multi -wo NUM_THREADS=ALL_CPUS -co COMPRESS=DEFLATE -co TILED=YES -co ZLEVEL=9 -co NUM_THREADS=ALL_CPUS -co BIGTIFF=YES ANA_LUREF_NGL_DTM_epsg2169_hillshade.tif ANA_LUREF_NGL_DTM_epsg3857_hillshade-compressed.tif
Ideally, you only want to reproject the data once. If you have to use lossy compression algorithms, try to use them as late as possible in the processing pipeline.
Step 7: Building Overview Images
You should build overview images for each file that you intend to serve locally or over the internet. For example, if you want to publish the VRT file we build earlier in MapServer, you could execute the following command:
gdaladdo -ro — config COMPRESS DEFLATE — config COMPRESS_OVERVIEW DEFLATE — config ZLEVEL 9 — config BIGTIFF_OVERVIEW IF_SAFER — config GDAL_TIFF_OVR_BLOCKSIZE 512 -r nearest ANA_LUREF_NGL_DTM.vrt 4 16 64 256 1024 4096
This generates overview tiles in TIF format using the lossless DEFLATE algorithm (ZLEVEL=9 sets the compression to the highest value possible). Furthermore, we generate overview tiles that are level-n times smaller than the input image (a level of 4 generates an overview that is a quarter of the input image). GDAL_TIFF_OVR_BLOCKSIZE defines the tile dimensions of 512x512 pixels.
Alternatively, you can build JPEG overview files using the following command:
gdaladdo -ro — config COMPRESS_OVERVIEW JPEG — config INTERLEAVE_OVERVIEW PIXEL — config BIGTIFF_OVERVIEW IF_SAFER — config GDAL_TIFF_OVR_BLOCKSIZE 512 -r average ANA_LUREF_NGL_DTM.vrt 4 16 64 256 1024 4096
Since we did not specify JPEG_QUALITY_OVERVIEW, the default value 75 is used. By using IF_SAFER for the BIGTIFF_OVERVIEW, we tell GDAL to use the BigTIFF format if the filesize exceeds 4GB.
Step 7: Making the Imagery available
Luckily the Operations Working Group of OpenStreetMap provided me with a MapServer instance to host the data. To host the map locally, you could use TileMill or QGIS. Alternatively, you could use gdal_calc.py to compute a single TIF file.
The final hillshade can be viewed by clicking here. The TMS URL is:
Note that the final hillshade might change over time and does not necessarily reflect the contents of this article. I’ll try and keep this article updated if any changes are made.
Acknowledgments
- Firefishy for helping with GDAL and MapServer
- Guillaume Rischard for the back and forth of information, stay tuned for his LiDAR hillshade using 2019 LiDAR data 😉
- The OpenStreetMap Operations Working Group for hosting the hillshade and making it available to the OSM Luxembourg community
- Chris Henrick for his awesome GDAL tutorial
