Creating a Terrain Wetness Index from a DEM

High-value areas of a Terrain Wetness Index in mainland British Columbia, east of the northern tip of Vancouver Island

Areas identified as high value within the Terrain Wetness Index. Located in mainland British Columbia, east of the northern tip of Vancouver Island.

From the archive

This post was originally published in April 2011 on my old "Mostly a GIS Guy" blog. It has been preserved here, lightly edited, as part of the blog migration to Little Earth.

A Terrain Wetness Index (TWI) is a means for measuring drainage in a landscape — in particular, areas with poor drainage. The purpose is to isolate dips and hollows that maintain water for longer periods, exclusive of permanent water bodies like lakes and swamps.

A Digital Elevation Model (DEM) was acquired at 25 metre resolution. It was in an older ESRI grid format and separated into tiles.

I'm too cheap to buy ArcInfo… so I use GRASS. First job is to import the grids:

Import the DEM tiles

r.in.gdal input=tdem092k.asc output=tdem092k
r.in.gdal input=tdem092l.asc output=tdem092l
...

Next is to merge them together:

Merge tiles into a seamless DEM

r.series input=tdem092k,tdem092l,tdem092m,tdem092n,tdem093c,tdem093d,tdem093e,tdem093f,tdem102i,tdem102p,tdem103a,tdem103h output=dem method=average

Now we have a seamless dataset to work with called dem. Time to do the real work.

Elevation Variance

Areas of high negative elevation variance.

The simplest way to calculate TWI is to run an elevation variance algorithm on the DEM. Simply put, this is a neighbour function that looks at the average elevation of all surrounding cells, then calculates variance based on itself and the average.

To do this, two separate commands were required. The first calculated mean elevation for each cell:

Mean elevation in a 3×3 window

r.neighbors -c input=dem output=mean3 size=3

This produced a dataset (mean3) holding the mean elevation for a window of 3×3. Considering the cells are 25 metres, this resulted in a 75 by 75 metre window. The r.neighbors command defaults to mean as the function, so it isn't necessary to specify it on the command line.

The second step is to calculate variance by comparing each elevation value in the original DEM with the new mean values. GRASS's extremely powerful and versatile r.mapcalc command is used for this.

Difference from the local mean

r.mapcalc mean3diff = dem - mean3

We are looking for the dips only, not the peaks. Peaks have positive values whereas dips have negative values. So we need to select for the negative variance values:

Keep only negative variance (dips)

r.mapcalc mean3dip = "if(mean3diff < 0, mean3diff, 0)"

These steps were repeated for different sized windows — 7 and 15. Then the three were added up:

Sum the multi-scale dips

r.series input=mean3dip,mean7dip,mean15dip output=dips method=sum

The result is illustrated above.

You may have noticed the image clearly identifies steep river valleys. This isn't exactly what we are after. We want to find low-lying areas of poor drainage. We got low-lying areas of high drainage.

Calculating mean elevation variance alone results in something like the following:

Open valley shape isolated by mean elevation variance, not yet a closed basin

Terrain shape isolated by Mean Elevation Variance.

We are after a closed shape. Kind of like this:

Closed bowl-shaped dip produced by combining mean elevation variance and mean aspect variance

A closed dip resulting from the combination of Mean Elevation Variance and Mean Aspect Variance.

To achieve this we need to add another variable: aspect variance.

Aspect Variance

Aspect is a term for terrain direction, or the direction in which a slope is facing. If we apply the process above to aspect, we can expect to isolate areas of high aspect variance.

Generate an aspect layer from the DEM with the r.slope.aspect command:

Derive aspect from the DEM

r.slope.aspect elevation=dem aspect=aspect

This creates a layer called aspect which contains continuous, floating point values of aspect in degrees. It is hard to correctly measure variance with continuous data. So the next step is to cluster the values into four classes — North, South, East, and West (identified numerically as 1–4). Mapcalc takes care of that:

Reclassify aspect into four directions

r.mapcalc aspect_rc = "if(aspect > 45 && aspect <= 135, 1) + if(aspect > 135 && aspect <= 225, 2) + if(aspect > 225 && aspect <= 315, 3) + if(aspect <= 45 && aspect > 0, 4) + if(aspect > 315 && aspect <= 360, 4) + if(aspect == 0, null())"

Now we can run the variance procedure:

Multi-scale aspect variance, then sum

r.neighbors -c input=aspect_rc output=variance3aspect size=3 \
method=variance
r.neighbors -c input=aspect_rc output=variance7aspect size=7 \
method=variance
r.neighbors -c input=aspect_rc output=variance15aspect size=15 \
method=variance
r.series input=variance3aspect,variance7aspect,variance15aspect \
output=aspects_3_15 method=sum

Creating the Final Terrain Wetness Index

After reviewing both datasets, aspects_3_15 and dips, it was determined to select areas based on the following:

Aspect variance greater than or equal to 10
Mean elevation variance less than -10

Another mapcalc statement takes care of this:

Combine the two criteria into the final TWI

r.mapcalc twi = "if(aspects_3_15 >= 10 && dips < -10, dips, null())"

When the criteria are met, the resulting cells are populated with elevation variance. This allows for a gradient of depth for dips and hollows. All cells not selected are populated with null values.

The result is our final Topographic Wetness Index and looks like this:

Final Topographic Wetness Index model.

This has proven to be a valuable input to a model I am currently working on. If you want to see it in action, take a look at the Yellow Cedar Decline Mashup.

— Jamie

← All Field Notes