Tutorial: Calculate the Topographic Wetness Index from a DTM

1. Introduction

2. Interpolate Voids in the DTM

The provided DTM is derived from LIDAR data from the Netherlands (AHN3). The data in the Dutch Amersfoort projection (EPSG: 28992). The DTM, however, has voids. Voids are pixels with nodata. These voids need to be interpolated before we can continue.

1. Start QGIS.

2. Add the provided DTM (dtm.tif) to a blank project.

You can clearly see the voids.

3. Go to the Processing Toolbox: in the main menu choose Processing | Toolbox.

4. In the Processing Toolbox select GDAL | Fill nodata.

5. In the Fill nodata dialogue choose dtm as Input layer. Keep other settings as default (you can play with these if the voids are bigger). Save the result as dtm_voidfilled.tif.

6. Click Run. Close the dialogue after processing.

7. Remove the dtm layer from the Layers panel.

In the result you can clearly see that the voids have been interpolated.

8. Style the DTM with Singleband Pseudocolor. You can use a Topography ramp from the cpt-city catalogue. Blend the DTM with the hillshade (see other tutorials on how to do this).

Now the DTM is ready we can proceed with calculating the slope.

Note that in QGIS 3.14 you need to set the projection of the raster layer after each calculation. This is not the case in the 3.10 LTR versions.

3. Calculate Slope in Degrees and in Radians

In order to calculate TWI we need the slope in radians. With the raster analysis tools we can calculate the slope in degrees and percentages. So we're first going to calculate the slope in degrees and then convert the values to radians.

1. In the main menu go to Raster | Analysis | Slope...

2. In the Slope dialogue choose dtm_voidfilled as Input layer. Check the box to use the ZevenbergenThorne formula instead of Horn's one. The ZevenbergenThorne formula works better with less relief, while we use Horn's formula for steep areas. Choose slopedegrees.tif as output slope layer.

3. Click Run. Close the dialogue after processing.

Because later we're going to calculate the natural logarithm (ln) we need to check if there are pixels with a slope of 0 degrees.

4. Click right on slopedegrees and choose Properties.

5. In the Layer Properties dialogue go to the Information tab and check the minumum value (STATISTICS_MINIMUM).

In this case the minimum slope is larger than 0, so we can proceed. However, if the minumum slope is 0 you need to do the following steps, that you can also apply to this data, because it will not have any effect.

6. Click Cancel to close the Layer Properties dialogue.

7. In the main menu go to Raster | Raster Calculator.

8. In the Raster Calculator create the following equation:

( "slopedegrees@1" <= 0 ) * 1 + ( "slopedegrees@1" > 0 ) * "slopedegrees@1"

This means: if cells have a slope less than or equal to 0, then give value 1 (True) and multiply by 1. For the cells that are larger than 0, assing value 1 (True) and multiply with its slope value. The result will be 1 for cells <= 0 and the original slope for cells > 0.

9. Save the result to slopedegrees_modified.tif and click OK.

To convert from degrees to radians we need to multiply the pixel values by 0.01745

10. Go back to the Raster Calculator and calculate the following equation:

"slopedegrees_modified@1" * 0.01745

11. Save the result to sloperadians.tif and click OK.

12. Style the result from yellow to red.

13. Remove the other slope layers.

The next step is to calculate the contributing upslope area for each pixel of the DTM.

4. Calculate the Contributing Upslope Area

In this section we're going to calculate the contributing upslope area for each pixel in the DTM.

1. Go to the Processign Toolbox and choose SAGA | Simulation | Flow accumulation (qm of esp).

2. In the dialogue choose dtm_voidfilled as DEM, keep the default Preprocessing to fill sinks temporarily and save the Contributing Area to contributing_upslope_area.sdat.

3. Click Run. Close the dialogue after processing.

The result is for each pixel the amount of pixels that are upstream.

Now we're ready to calcute TWI.

5. Calculate TWI

Now we have the slope in radians and the upslope contributing area (in pixels) we can calculate the Topographic Wetness Index in the Raster Calculator.

1. Go to the Raster Calculator.

2. Write the following equation:

ln ( ( "contributing_upslope_area@1" * 5 * 5 ) / tan ( "sloperadians@1" ) )

Note that we multiply the contributing upslope area (amount of pixels) with the area of the pixel (5 m x 5 m).

3. Save the result to twi.tif.

Note that if the Output CRS is unknown, you need to set it for the input layer.

4. Click OK.

5. Style the result with a red to yellow to blue colour ramp and use blending with the hillshade layer.

We can also visualise the result in 3D.

6. Visualise TWI in the 3D view