Generating Wetness Indices for Watersheds in GRASS

Introduction:

In hydrology, a watershed refers to an area of land in which all the precipitation runoff drains into a common waterway, such as a lake, river, or ocean. Therefore, drainage basins act as natural funnels, collecting the precipitation that lands within their borders, and directing it into the common waterway outlet. Watersheds are arranged in hierarchical order, where the size of a given watershed is determined by the location of its outlet. Therefore, as the outlet of the watershed moves downstream the watershed becomes larger, and consists of more sub-watersheds. In hydrology, watersheds are extremely important since they delineate regions that contribute to flow in a given waterway.

Wetness indices are used to describe spatial moisture patterns within a given region. They are calculated based on a region’s contributing upslope area and slope, where areas with low slope and large contributing area are expected to be relatively wet and areas with high slope and small contributing area are expected to be relatively dry. When a wetness index is calculated for a given region, the index should be calculated for the entire watershed in which that region is located, in order to ensure accuracy in contributing upslope area.

Outline:

This tutorial illustrates the processes involved in delineating watershed boundaries, as well as generating topological wetness indices inside those boundaries. These processes will be performed using GRASS tools through Quantum GIS. The tutorial also discusses the effects of filtering DEMs on wetness indices, and introduces a means of verifying both watershed boundaries and wetness indices.

Images from a sample analysis are also provided for each step that is described in the tutorial. In the sample analysis, the Carp River watershed in Ottawa, Ontario is delineated, and a wetness index is generated for the watershed region. A digital elevation model of the Ottawa-Gatineau region, with 50m cell size, was used for the sample analysis.

Figure 1: A 50m DEM for the Ottawa-Gatineau region that was used in sample analysis.

Setting Up Your Workplace

1. In QGIS, open the DEM using the 'Add Raster Layer' button in the 'Manage Layers' toolbar.

2. Once the DEM is loaded in QGIS, create a new GRASS mapset with an appropriate projection and workplace extent.

3. Using r.in.gdal.qgis, import the DEM into the new mapset.

Generating Watershed Boundaries:

There are two available modules in GRASS that allow watershed delineation: r.watershed and r.water.outlet. Although these modules produce similar output, they delineate watersheds in two different ways.

r.watershed

- This method is designed for analysis in which the delineation of all watersheds over a specified size is required.

- However, if a watershed with a specific outlet point is required, the r.watershed module is not recommended.

To delineate watersheds in r.watershed:

1. Open the r.watershed module.

2. Select the DEM as the input map.

3. Select a proper input value. The input value (or threshold in command line) refers to the minimum size in image pixels a watershed must contain, for it to be delineated.

4. Specify a name for 'Output Map: Unique label for each watershed basin'.

Alternatively, this can be done in command line based on: [r.watershed elevation=DEM threshold=# basin=Watersheds]

- This watershed output image will include all watersheds that exist within the DEM, given the threshold parameter that was specified. 
- Each watershed in the image is assigned a unique pixel value.

Watershed Masks:

Since this method produces a watershed map consisting of numerous watersheds, if analysis on a particular watershed is required, a watershed mask is often required.

A watershed mask can be produced in command line based on: [r.mapcalc WatershedMask=if (Watersheds==’unique watershed label’]

In the mask, areas within the watershed are assigned a value of 1, and all other areas are assigned a value of 0.

r.water.outlet

- Unlike r.watershed, this method will delineate a watershed based on the coordinates of its outlet point.

- This allows practitioners to customize the scale of output watersheds, based on their research needs.

To delineate watersheds in r.water.outlet:

1. Open the r.watershed module.

2. Select the DEM as the input map.

3. Specify a name for 'Output Map: Drainage Direction'.

- The output image quantifies the drainage direction of cells throughout the DEM.

4. Open the r.water.outlet module.

5. Select the Drainage Direction map as 'Name of Input Raster Map'.

6. Specify the easting and northing of the watershed outlet point in 'The Map E Grid Coordinates' and 'The Map N Grid Coordinates'.

7. Specify a name for "Name of Raster Map to Contain Results'.

Alternatively, this can be done in command line based on: [r.watershed elevation=DEM drainage=DrainageDirection], [r.water.outlet drainage=DrainageDirection basin=Watersheds easting=x northing=y]

- The watershed output image will only include the watershed for which the outlet point was specified.
- The watershed output image is essentially a mask of the watershed.

Figure 2: Left: The watershed boundaries generated in r.water.outlet, using x,y coordinates obtained from a topographic map. Right: The watershed boundaries generated in r.watersheds, using a threshold parameter of 10000.

Generating Wetness Indices:

Slope

This section will explain the steps involved in generating and preparing a slope image for the wetness index calculation.

To generate a slope map in r.slope.aspect:

1. Open the r.slope.aspect.slope module.

2. Select the DEM for ‘Name of Elevation Raster Map’.

3. Specify a name for ‘Name for Output Slope Raster Map’.

Alternatively, this can be done in command line based on [r.slope.aspect elevation=DEM slope=Slope format=degrees]

- The slope output image reports the slopes of the DEM landscape in degrees.
- The slope image will commonly contain regions where slope equals zero. 
- Since slope appears in the denominator of the wetness index equation, zero slope areas will result in undetermined (null) regions in the output wetness index maps. 

Eliminating Zero Slope Areas in Slope Maps:

To avoid null values in the wetness index image, zero slope areas in the output slope images should be assigned a value of 0.0001 (or any other very small value).

To eliminate zero slope areas from slope maps:

1. Open the r.mapcalculator module.

2. Select the output slope map as ‘A’.

3. Enter ‘A+0.0001’ as the ‘Formula’.

4. Specify a name for ‘Name for Output Raster Map’.

Alternatively, this can be done in command line based on [r.mapcalc Slope_NoZeros = Slope + 0.0001]

- This will replace zero slope values with slopes of 0.0001.
- The increase in slope in all other areas is negligible.

Figure 3: The output slope image and the user-defined color rules that were used to display the image.

Upslope Contributing Area:

This section will explain the steps involved in generating and preparing an upslope contributing area image for the wetness index calculation.

To generate a flow accumulation map in r.watershed:

1. Open the r.watershed module.

2. Select the DEM as the input map.

3. Specify a name for 'Output Map: Number of Cells that Drain through Each Cell'.

Alternatively, this can be done in command line based on [r.watershed elevation=DEM accumulation=AccumFlow]

- The output image quantifies the number of cells that drain through each cell or, in other words, the number of cells that are upslope from a given cell.

Since flow accumulation represents the number of upslope cells relative to a given cell, to obtain upslope contributing area, multiply flow accumulation by the area of cells in the image.

To generate an upslope contributing area map:

1. Open the r.mapcalculator module.

2. Select the flow accumulation image as ‘A’.

3. Enter ‘A * cell area’ as the ‘Formula’.

4. Specify a name for ‘Name for Output Raster Map’.

Alternatively, this can be done in command line based on [r.mapcalc UpslopeContributingArea = FlowAccumulation * ‘cell area’]

Figure 4: The output accumulated flow image and the user-defined color rules that were used to display the image. Although the upslope contributing area image was used to calculate wetness index, accumulated flow is illustrated since the two images appear exactly the same.

Wetness Index:

Wetness index is calculated based on [WetnessIndex=ln(A/tan(B))] where A represents the upslope contributing area and B represents slope.

To generate a wetness index map:

1. Open the r.mapcalculator module.

2. Select the upslope contributing area image as ‘A’.

3. Select the slope image as ‘B’.

4. Enter ‘log(A/tan(B), 2.718282)’ as the ‘Formula’.

4. Specify a name for ‘Name for Output Raster Map’.

Alternatively, this can be done in command line based on [r.mapcalc WetnessIndex = log(UpslopeContArea/tan(Slope), 2.718282)]

Figure 5: The output wetness index image and the user-defined color rules that were used to display the image.

Standardizing Wetness Index Maps:

A common approach when comparing wetness indices is to standardize values between 0-1. If output wetness index images are going to be compared with other indices, standardize index values using r.mapcalc and r.info.

Standardizing images can be done in command line based on [r.mapcalc StandardizedImage= ((Image – ImageMin)/(ImageMax – ImageMin))].

Displaying Maps:

Given the extremely large flow values in creeks, streams, and tributaries within the watershed, flow accumulation images are often hard to interpret without user-defined color rules. All slope and wetness index images also required user-defined color rules for proper visual display.

To display images based on user-defined color rules:

1. Generate distributions for an image using the [r.stats –c FileName] module in GRASS or the Histogram function in QGIS.

2. Based on the distribution and your research needs, determine suitable color breaks for displaying the image.

3. Create a text file (.txt) containing the color rules for the image.

Refer to [1] for syntax of color rules.

4. Open the r.colors module.

5. Select the image in ‘Name of Input Raster Map’.

6. Specify the full path to user-defined color rules text file in the ‘Path to Rules File’ browser.

Alternatively, this can be done in command line based on [r.colors map=FlowAccum rules=‘FullPathToColorRules’]

DEM Filtering:

So far, all sample analysis images that are illustrated were generated from an unfiltered DEM. However, in order to achieve a specific level of detail in final wetness index images, filtering the DEM prior to generating the slope and flow accumulation maps may be required. Filtering can also reduce errors in future output maps that are caused by the DEM, or how the DEM was created. The DEM should be the only image that is filtered throughout the analysis. This ensures that all output images all derivative products are derivatives of a specific elevation model. However, filtering final output maps for visual display is acceptable.

To apply a filter to a DEM:

1. Create a filter file which contains the filter you have chosen.

- In the sample analysis, the filter file contained a single 5x5 average filter; however multiple filters can be included in one filter file.
- Refer to [2] for syntax of filter files.

2. Open the r.mfilter module.

3. Select the DEM in ‘Name of Input Raster Map’.

4. Specify a name for ‘Name of Output Raster Map’.

5. Specify the full path and name of filter file in ‘Name of Filter File’.

6. Specify the number of times to repeat the filter.

Alternatively, this can be done in command line based on [r.mfilter input=DEM output=DEM_filtered filter= ‘path to filter file’ repeat= ‘number of repeats’].

In the sample analysis, a total of four wetness indices were generated for the Carp watershed, each generated from a DEM filtered by a 5x5 average filter a varying number of times. All indices were generated using the methodologies described in this tutorial. The indices were generated based on 4 DEMs that were filtered 0 times, 5 times, 10 times and 15 times prior to generating slope and flow accumulation images. Images of each output wetness index are provided to allow investigation into the wetness index differences that result from filtering the DEM.

The following images should give an idea of the basic effects of average filtering.
However, Generally trial and error is required to find the proper series of filters for your research needs.

Figure 6: Wetness indices that were generated from a differently filtered DEM.
[Top Left: Unfiltered, Top Right: 5x5 filter x5 times, Bottom Left: 5x5 filter x10 times, Bottom Right: 5x5 filter x15 times]

Figure 7: Wetness indices generated from differently filtered DEMs, superimposed on topographic maps during verification.
[A: Unfiltered, B: 5x5 filter x5 times, C: 5x5 filter x10 times, D: 5x5 filter x15 times]

Figure 8: Wetness indices generated from differently filtered DEMs, superimposed on topographic maps during verification.
[A: Unfiltered, B: 5x5 filter x5 times, C: 5x5 filter x10 times, D: 5x5 filter x15 times]

Verifying Outputs:

Verifying Watershed Boundaries:

To assess the accuracies of the watershed boundary outputs, obtain digital maps (either raster or vector) of verified watershed boundaries from local conservation authorities. By superimposing the verified watersheds over the watershed outputs, the accuracy of watershed outputs can be visually evaluated.

Figure 9: The two watershed outputs (Left: r.watershed, Right: r.water.outlet) and the watershed vector file, published by Ministry of Natural Resources.

Verifying Wetness Index:

If digital topographic maps are available in your area, they can provide a good means of evaluating the potential accuracies of your output wetness index images. By superimposing the wetness index images (with transparency of ~60%) over the topographic maps, the types of land cover (based on topographic maps) that fall in relatively wet or dry areas (based on wetness index) can be investigated. Based on visual comparison between the two images, the accuracy of the wetness index can be evaluated. For example, if the majority of creeks, streams, and swampy areas on the topographic map are shown to be ‘wet’ areas on the wetness index, this would suggest a certain level of accuracy. On the other hand, if the majority of creeks, streams, and swampy areas on the topographic map are shown to be ‘dry’ areas on the wetness index, the index is likely inaccurate.

Figure 10: A screenshot during wetness index verification, in the northern region of the Carp watershed.

Figure 11: A screenshot during wetness index verification, in the southern region of the Carp watershed.