Creating color relief and slope shading with gdaldem

In my previous blog post we created a grayscale shaded relief (hillshade) from a Digital Elevation Model (DEM). Today, we’ll use the DEM to create color relief and slope shading maps.

Color relief or hypsometric tints depict elevation as bands of color, to enhance elevation zones so map readers can better see differences in relief. The colors selected for the tints are assumed to relate to the ground cover typically found at various elevations in the area being mapped. A typical color scheme progresses from dark greens for lower elevations up through yellows/browns, and on to grays and white at the highest elevations. I selected this scheme for my area:

Elevation	Color
0
900
1300
1900
2500

To create a color relief with gdaldem, we first have to create a text-based color configuration file (color_relief.txt), containing the association between elevation values and colors. The file generally contains 4 columns per line: the elevation value and the corresponding red, green, blue value between 0 and 255 (RGB):

0 110 220 110
900 240 250 160
1300 230 220 170
1900 220 220 220
2500 250 250 250

The above numbers define a gradient that will blend 5 colors over 2500 meters of elevation. We use this command to apply this color ramp to the Jotunheimen DEM:

gdaldem color-relief jotunheimen.tif color_relief.txt jotunheimen_colour_relief.tif

[ Download image ]

A Digital Elevation Model (DEM) can also be used to measure slope, the steepness of the terrain. Slope is estimated by comparing each pixel of the DEM to the elevations of the surrounding pixels. Creating slope shading with gdaldem is a two step process.

First, we use the DEM to create a raster where each pixel contains an angle, varying from 0 to 90 degrees as the terrain ranges from horizontal to vertical:

gdaldem slope jotunheimen.tif jontunheimen_slope.tif

We then use this raster to create slope shading by assigning a color to each pixel, using the same technique as for color relief.

Create a color configuration file (color_slope.txt) containing these two lines:

0 255 255 255
90 0 0 0

Flat terrain will be displayed in white (slope of 0°) and cliffs in black (slope of 90°), with angles in-between being various shades of grey. The command is:

gdaldem color-relief jotunheimen_slope.tif color_slope.txt jotunheimen_slopeshade.tif

[ Download image ]

In the next blog post we’ll combine slopeshade, hillshade and color relief into a single map image.

Creating hillshades with gdaldem

In the first part in this blog series we created a Digital Elevation Model (DEM) for Jotunheimen, a mountainous area in Norway. We’ll use this DEM to create hillshade or shaded relief, a popular cartographic technique to visualise terrain by modulating light and shadows on a map. 

GDAL, my favorite GIS Swiss Army Knife, allows you to create hillshade with gdaldem:

gdaldem hillshade -of PNG jotunheimen.tif jotunheimen_hillshade.png

This command will create this PNG image:

[ Download image ]

A virtual light source is placed above the DEM to calculate which areas are lightened up and which fall in the shadow. You can clearly see the numerous mountains and valleys. By default, the light source is placed in a top left position (azimuth = 315 degrees). Let’s try to move the light source to a bottom right position (azimuth = 135 degrees):

gdaldem hillshade -of PNG -az 135 jotunheimen.tif jotunheimen_hillshade_az135.png

[ Download image ]

You get the completely opposite effect, the valleys are percepted as mountain ridges and mountains appears as valleys (multistable perception illusion). I’ve computed a video with a 360 degree change of light source:


It’s easier to avoid relief inversion when the light source is changed gradually.

You can also change the altitude of light source (-alt) from 0 (ranking light) to 90 degrees (directly above the terrain). This video shows the effect:


Lastly, you can change the exaggeration (-z) used to pre-multiply the DEM elevations. This video shows a gradual change from zero to 3:1 terrain exaggeration, where mountains are three times higher:


Our hillshade is still a boring black and white image. In the next blog post we’ll create a color relief - hypsometric tints - for our terrain map.

Digital Terrain Modelling and Mapping

This is the first post in a new blog series focusing on all the fun you can do with digital terrain data. It’s going to be more related to topographic than thematic mapping, although the two map types can be used in combination. I’ll also try to restrict myself to data and software in the public domain.

My focus area is going to be Jotunheimen in Norway, a mountainous area popular among hikers and climbers. The 29 highest mountains in Norway are all in Jotunheimen, including the very highest - Galdhøpiggen (2469 m). Here's a few pictures I took on a hike back in July 2006 which gives you a glimpse of the beautiful landscape:

Skogadalen valley and Hurrungane mountain range

Fannaråkhytta

Bøverbreen glacier

It’s many different data formats for storing and working with digital terrain data. I’ll start looking at Digital Elevation Models (DEM) represented as rasters. A raster DEM is a grid of squares or pixels where each square represents an elevation at a geographic location. The size of each square (in meters or arc-seconds) tells us about the level of detail or data resolution.

In Norway, it’s possible to get digital terrain data with 10 meters resolution, but unfortunately it’s not yet available to the public. Instead I found a suitable dataset on Viewfinder Panoramas, created from topographic maps. You’ll find  DEMs with global coverage at various resolutions here, here and here.

I downloaded two files from Viewfinder Panoramas (N61E007.hgt and N61E008.hgt), each covering an area of 1 degree latitude and longitude marked in yellow:


Bigger map | KML

You can use the great GDAL library to combine and clip the DEMs to the area of interest. I used gdalbuildvrt to create a combined virtual dataset of the DEM files:

gdalbuildvrt jotunheimen.vrt N61E007.hgt N61E008.hgt

Then you can use gdalwarp to project and clip the DEM to the area of interest:

gdalwarp -t_srs EPSG:32632 -te 432000 6790000 492000 6850000 -r bilinear jotunheimen.vrt jotunheimen.tif

The original DEM data is unprojected and I’m projecting the data using UTM zone 32N (EPSG:32632). UTM coordinates (meters) are used to clip the dataset to an area covering 3,600 km² (60 x 60 km). The area is shown with a red square on the map above. The -r bilinear option is important when projecting elevation data because other resampling methods tend to produce odd stripes in the image.

We finally have a DEM for Jotunheimen (download). You won’t see much if you try to open this image in Windows. The image is almost all black. Use gdalinfo to see the properties of this image:

gdalinfo -mm jotunheimen.tif

Driver: GTiff/GeoTIFF
Files: jotunheimen.tif
Size is 3134, 3134
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32632"]]
Origin = (432000.000000000000000,6850000.000000000000000)
Pixel Size = (19.144862795149969,-19.144862795149969)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (432000.000, 6850000.000) (7d42'39.90"E, 61d46'36.81"N)
Lower Left  (432000.000, 6790000.000) (7d43'59.45"E, 61d14'18.21"N)
Upper Right (492000.000, 6850000.000) (8d50'54.02"E, 61d46'58.29"N)
Lower Right (492000.000, 6790000.000) (8d51'3.39"E, 61d14'39.21"N)
Center      (462000.000, 6820000.000) (8d17'9.19"E, 61d30'42.34"N)
Band 1 Block=3134x1 Type=Int16, ColorInterp=Gray
    Computed Min/Max=15.000,2466.000

What does it tell you?

• The image format is GeoTIFF, which allows georeferencing information to be embedded within a TIFF file.
• The image or grid size is 3134 x 3134 pixels or cells.
• The coordinate system is UTM zone 32N using the WGS 84 spheroid.
• The coordinate unit is meters.
• Each pixel or grid cell covers an area of approximately 19 x 19 meters.
• The extent of the area is defined by four coordinate pairs.
• The image is single band and the elevation data is stored as 16-bit integers. The minimum elevation is 15 meters and the maximum is 2,466 meters.

Let’s use gdal_translate convert the GeoTiff image to PNG showing elevation with 256 shades of gray (-scale).

gdal_translate -of PNG -ot Byte -scale 15 2466 0 256 jotunheimen.tif jotunheimen.png

This image gives you a glimpse of the elevation in Jotunheimen. Dark areas are valleys with low altitude. Light areas are mountain ridges with high altitude. Areas with uniform shades are lakes or areas with less variation in elevation.

Hillshade is a more effective technique to visualise terrain data. That’s the topic for the next blog post.