GDAL Terrain analysis I (Arctic DEM)

Introduction

From a Digital Elevation Model (DEM), a range of terrain metrics can be extracted. Perhaps the most common metrics are slope and aspect where slope is the steepness and aspect the horizontal direction of the slope. To illustrate the 3D character of a DEM, the hillshade can be computed and set as an opacity-transparency mask. Karttur’s GeoImagine Framework can use either GDAL, GRASS or python scripting to derive different terrain indices. GDAL is the most convenient and well documented alternative, but also the least flexible/powerful. This post introduces the Framework interface to gdaldem.

Overview

The processes GdalDemTiles and GdalDemRegion can calculate the following terrain indexes (as defined in gdaldem:

slope
aspect
hillshade
Topographic Position Index (TPI)
Topographic Ruggedness Index (TRI)
roughness
color-relief

This post deals with the processing of hillshade, slope and aspect. It covers how to use Karttur’s GeoImagine Framework for calling the standard DEM derivates that can be calculated with GDAL DEM utility (gdaldem). The next post introduces GDAL derived metrics that can be used for capturing elevation variations and positions.

All GDAL related DEM processing can be done from either original tiles, virtual tiles with mosaicked (expanded) coverage, or regions. The construction and use of virtual tiles is recommended and outlined in a previous post. You must prepare the mosaicked tiles beforehand, then just set the parameter mosaic to true when calling the process GdalDemTiles.

By default GdalDemTiles (and GdalDemRegions) applies a 3x3 kernel for calculating DEM derivates. You can change that for all derived indexes with the parameter radiuscsv, and a comma separated value list in the form of a text string: The default value “1” applies a (1+1+1 = 3 cells square kernel) 3x3 kernel for calculating the DEM index for the central pixel. A value of “2” aggregates (by averaging) all 2x2 pixels to a single cell and then calculates the index on the average of 3x3 cells originally derived from 6x6 cells; a value of “3” calculates the index on the average of 3x3 cells, originally derived from 9x9 cells; etc. The process NumpyDemTiles and Grass1to1Tiles (or GrassOnetoManyTiles can be set to a more fine tuned spatial aggregation. GRASS also offers much more resourceful alternatives for deriving DEM indexes.

All destination layers have default names that build on the input DEM and the spatial aggregation. You can override some parts, but it is strongly recommended to let the Framework set the names automatically and use the parameterizations as illustrated in the following examples.

Hillshade

Hillshade produces an 8-bit raster with a shaded relief effect useful for visualizing the terrain. You can optionally specify the azimuth and altitude of the light source, as well as a vertical exaggeration factor and a scaling factor to account for differences between vertical and horizontal units.

Framework process and json commands

Framework process for GDAL Hillshade (other alternatives exist): GdalDemTiles.

The json parameters relate to the parameters in the GDAL DEM utility (gdaldem).

Slope

With mode set to slope GdalDemTiles outputs an 32-bit Float raster with a slope in degrees unless the -p (percent) parameter is set.

Framework process and json commands

Framework process for GDAL Slope (other alternatives exist): GdalDemTiles.

The json parameters relate to the parameters in the GDAL DEM utility (gdaldem).

Aspect

With mode set to aspect GdalDemTiles outputs an 32-bit Float raster with azimuth in degrees; if the parameter -trigonometry is given, the output is in angle degrees with zero set to the East (rather than North).