Difference between revisions of "Generating Wetness Indices for Watersheds in GRASS"

From CUOSGwiki
Jump to navigationJump to search
Line 2: Line 2:
 
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.
  
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 flow accumulation accuracy.
+
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 of that region, in order to ensure accuracy in contributing upslope area.
−
In the sample analysis, wetness index was calculated for the Carp watershed based on [WetnessIndex=ln(A/tan(B))] where A represents the upslope contributing area and B represents slope.
  
   
   
Line 13: Line 12:
 
[[File:DEM.png]]
  
[[File:DEM.png]]
   
−
Figure 1: 50m DEM for the Ottawa-Gatineau region.
 +
Figure 1: A 50m DEM for the Ottawa-Gatineau region that is used in the sample analysis.
   
 
=Setting Up Your Workplace=
  
=Setting Up Your Workplace=
Line 25: Line 24:
   
 
=Generating Watershed Boundaries:=
  
=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 different ways.
+
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''===
  
===''r.watershed''===
Line 38: Line 37:
 
1. Open the ''r.watershed'' module.
  
1. Open the ''r.watershed'' module.
   
−
1. Select the DEM as the input map.
 +
2. Select the DEM as the input map.
   
−
2. Select a proper input value. The input value (or threshold in command line) refers to the minimum size in image pixels a watershed can contain, for it to be delineated.
 +
3. Select a proper input value. The input value (or threshold in command line) refers to the minimum size in image pixels a watershed can contain, for it to be delineated.
   
−
3. Specify a name for 'Output Map: Unique label for each watershed basin'.
+
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.
 
   
  +
- This watershed output image will include ''all'' watersheds that exist within the DEM, given the threshold parameter that was specified.
−
- Every watershed in the image is assigned a unique pixel value.
  
   
  +
- Every watershed in the image is assigned a unique pixel value.
   
  +
−
Alternatively, this can be done in command line based on: ''[r.watershed elevation=DEM threshold=# basin=Watersheds]''
  
  +
===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. In the mask, areas within the watershed are assigned a value of 1, and all other areas are assigned a value of 0.
  +
  +
A watershed mask can be produced in command line based on: [r.mapcalc WatershedMask=if (Watersheds==’unique watershed label’]
   
   
Line 77: Line 81:
 
7. Specify the easting and northing of the watershed outlet point in 'The Map E Grid Coordinates' and 'The Map N Grid Coordinates'.
  
7. Specify the easting and northing of the watershed outlet point in 'The Map E Grid Coordinates' and 'The Map N Grid Coordinates'.
   
−
8. Specify a name for "Name of Raster Map To Contain Results'.
 +
8. 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]
−
- The watershed output image will only include the watershed for which the outlet point was specified.
  
  +
[r.water.outlet drainage=DrainageDirection basin=Watersheds easting=x northing=y]
   
−
- The watershed output image is a mask of the watershed, where areas within the watershed are assigned a value of 1, and all other areas are assigned a value of 0.
 +
- 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.
   
−
Alternatively, this can be done in command line based on:
  
   
−
[r.watershed elevation=DEM drainage=DrainageDirection]
  
   
  +
[[File:Watersheds2.png]]
−
[r.water.outlet drainage=DrainageDirection basin=Watersheds easting=x northing=y]
 
   
  +
Figure 2: Left: The watershed boundaries generated in r.watersheds, using a threshold parameter of 10000. Right: The watershed boundaries generated in r.water.outlet, using x,y coordinates obtained from a topographic map.
−
[[File:TwoWatersheds.png]]
  
   
−
===Generating a Watershed Mask:===
  
−
In order to make images watershed-specific, a mask for the watershed of interest is required. In this case, the mask is used to prepare images for visual display, and to eliminate extraneous data in output images, allowing for proper color schemes to be used. In r.mapcalc, generate a mask for the watershed of interest based on
 
−
[r.mapcalc WatershedMask=if (Watersheds==24]
 
−
Where 24 is the cell value representing the Carp watershed in the watersheds map.
  
   
   
−
=Generating Slope Maps:=
 +
=Generating Wetness Indices:=
  +
==Slope==
−
To generate slope maps of the watershed that will be used to calculate wetness index, use the r.slope.aspect module. Report the input DEM and output slope-map names, and set the format by which the module reports slope (important depending on the wetness index equation that will be used).
 
  +
This section will explain the steps involved in generating and preparing a slope image for the wetness index calculation.
−
In the sample analysis, slope maps were set to report slope in degrees. Slope maps were generated based on [r.slope.aspect elevation=DEM slope=Slope format=degrees].
 
   
  +
To generate a slope map in r.slope.aspect:
−
===Eliminating Zero Slope Areas:===
  
−
Since slope appears in the denominator of the wetness index equation, regions with slope values of 0 will result in undetermined results for wetness index. To avoid null areas in the wetness index images, cell values of 0 should be assigned a value of 0.0001.
 
−
To eliminate cells with a value of 0, use r.mapcalc to multiply the mask of the watershed by 0.0001. The slope mask is generated based on
 
−
[r.mapcalc SlopeMask = WatershedMask*0.0001]
  
   
  +
1. Open the r.slope.aspect.slope module.
−
Next, in r.mapcalc, add the slope image and the slope mask together. This will replace values of 0 degrees with values of 0.0001 degrees, while producing only a negligible increase in slope in all other areas.
 
−
In the sample analysis, the final slope map was generated based on:
 
−
[r.mapcalc Slope_wNoZeros = Slope + SlopeMask]
  
   
  +
2. Select the DEM for ‘Name of Elevation Raster Map’.
−
The slope image, and the color scheme used to display the image, is illustrated below in Figure 3, and Table 1, respectively. Note that in Figure 3 the slope image has been multiplied by the watershed mask for visual display.
  
   
  +
3. Specify a name for ‘Name for Output Slope Raster Map’.
−
[[File:Slope.png]]
  
  +
  +
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.
   
−
=Generating Accumulated Flow Maps:=
  
−
To generate maps of accumulated flow, use the r.watershed module. The same threshold parameter that was originally used to delineate the watersheds should be used to generate accumulated flow.
 
−
In the sample analysis, accumulated flow maps were generated based on
 
−
[r.watershed elevation=DEM threshold=10000 accumulation=AccumFlow].
  
   
  +
===Eliminating Zero Slope Areas in Slope Maps:===
−
The accumulated flow image, and the color scheme used to display the image, is illustrated below in Figure 4, and Table 3, respectively. Similar to Figure 3, the accumulated flow image has also been multiplied by the watershed mask for visual display.
  
  +
To avoid conflict, 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:
−
[[File:AccumFlow.png]]
  
   
  +
1. Open the r.mapcalculator module.
   
  +
2. Select the output slope map as ‘A’.
−
=Generating Wetness Index Maps:=
  
−
Using the slope and flow accumulation maps from Sections 2 and 3, calculate the wetness index of the watershed using r.mapcalc. As discussed earlier, wetness index is calculated based on a region’s contributing upslope area and slope, where contributing upslope area is calculated by multiplying the flow accumulation image by the cell size.
  
−
In the sample analysis, wetness index was calculated based on
 
−
[r.mapcalc WetnessIndex = log(((AccumFlow*AREAcell)/tan(Slope)), 2.718281828)]
  
   
  +
3. Enter ‘A+0.0001’ as the ‘Formula’.
−
The wetness index image, and the color scheme used to display the image, is illustrated in Figure 5, and Table 2, respectively. Finally, the image was multiplied by the watershed mask.
  
   
  +
4. Specify a name for ‘Name for Output Raster Map’.
−
[[File:WetnessIndex.png]]
  
  +
  +
  +
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.
  +
  +
[[File:Slope1.png]]
  +
  +
Figure 3: The slope image that was used to calculate wetness index. The user-defined color rules that were used to display the image are also illustrated.
  +
  +
  +
=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 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 a single cell 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’]
  +
  +
  +
[[File:AccumFlow1.png]]
  +
  +
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)]
  +
  +
  +
[[File:WetnessIndex1.png]]
  +
  +
Figure 5: The output wetness index image and the user-defined color rules that were used to display the image.
   
 
===Standardizing Wetness Index Maps:===
  
===Standardizing Wetness Index Maps:===
−
A common approach when comparing wetness indices is to standardize values between 0-1. If the wetness index image is going to be compared with other indices, standardize index values using r.mapcalc and r.info. Standardized wetness images can be generated based on [r.mapcalc WetnessIndex_Standardized= ((WetnessIndex – minWetnessIndex)/(maxWetnessIndex – minWetnessIndex))].
 +
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:=
  
=Displaying Maps:=
−
Due to large flow values in rivers, streams, and creeks within the watershed, flow accumulation and wetness index maps are often ‘washed out’ due to poor color schemes. To allow proper interpretation of these maps, generate user defined color rules based on the distributions of values in the map. To acquire distributions, use r.stats –c to view the counts of each value that exist in the image, or use the histogram function provided for gdal rasters in QGIS. Apply color rules in r.colors based on
 +
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 other images required user-defined color rules for proper visual display.
  +
−
[r.colors map=FlowAccum rules=c:/CarpWatershed/Location/Mapset/FlowColors.txt]
  
  +
To generate user-defined color rules:
  +
  +
1. Obtain the distributions for values in the image using:
  +
\tA. r.stats –c FileName – returns counts for each value that exists in the image.
  +
\tB. QGIS Histogram function.
  +
  +
2. Based on the distributions, determine suitable color breaks for the color rules.
  +
  +
3. Create color rules in a text file (.txt) based on:
  +
  +
[[File:ColorRules.png]]
  +
  +
  +
To apply user-defined color rules:
  +
  +
1. Open the r.colors module.
  +
  +
2. Select the image in ‘Name of Input Raster Map’.
  +
  +
3. 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:=
  
=DEM Filtering:=
−
In order to achieve a specific level of detail in the final wetness index image, 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.
+
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.
−
By limiting filtering to the DEM only, all output maps are derivatives of a specific elevation model. For example, a slope output image from a filtered DEM represents slopes in that specific DEM. However, if the slope image was filtered after output, it would be difficult to tell what DEM the slope image is reporting slope for.
 
−
The type of filter that is applied to the DEM, and the number of times the filter(s) are run, will have large effects on both the appearance and usability of slope, flow accumulation, and wetness index images. Generally trial and error is required to find the proper series of filters to produce the desired level of detail, and ‘smooth’ the DEM enough to reduce future errors in output maps.
 
−
In the sample analysis, filters were applied to the DEM based on
 
−
[r.mfilter input=DEM output=DEM_filtered filter= ‘path to filter file’ repeat= ‘number of repeats’].
  
   
  +
Although slope, flow accumulation and wetness index can be filtered directly after output, by limiting filtering to the DEM only, all output maps are known to be derivatives of a specific elevation model. In other words, by limiting filtering to the DEM, it is very clear which DEM the maps are reporting slope, flow, or wetness for.
−
Figures 6, 7, and 8 show wetness index images that were generated from a DEM that was filtered a varying number of times with a 5x5 average filter. All images have the same scale as the wetness index image in Figure 4.
  
   
  +
In this 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 four DEMs that were filtered 0 times, 5 times, 10 times and 15 times prior to generating slope and flow accumulation images.
−
[[File:5x5AVGx5.png]]
  
   
  +
To apply a filter to a DEM:
−
[[File:5x5AVGx10.png]]
  
   
  +
1. Create a filter file which contains the series of filters that are required.
−
[[File:5x5AVGx15.png]]
  
  +
  +
- In the sample analysis, the filter file contained a single 5x5 average filter; however multiple filters can be included in one filter file.
  +
  +
- In the sample analysis, the filter file was as follows:
  +
  +
[[File:FilterFile.png]]
  +
  +
  +
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’].
  +
  +
  +
Images are provided from each output wetness index, to allow investigation of the differences in wetness index output resulting from filtering the DEM with a 5x5 average filter. Generally trial and error is required to find the proper series of filters for your research needs. However, the following images should give an idea on the basic effects of average filtering.
  +
  +
[[File:WetnessIndexComparison.png]]
  +
  +
Figure 6: Wetness indices that were generated from differently filtered DEMs.
  +
  +
  +
[[File:NorthernWetnessIndices.png]]
  +
  +
Figure 7: Wetness indices, generated from differently filtered DEMs, in the northern region of the Carp watershed.
  +
  +
  +
[[File:SouthernWetnessIndices.png]]
  +
  +
Figure 8: Wetness indices, generated from differently filtered DEMs, in the southern region of the Carp watershed.
   
   
 
=Verifying Outputs:=
  
=Verifying Outputs:=
 
==Verifying Watershed Boundaries:==
 
==Verifying Watershed Boundaries:==
−
To assess the accuracy of the watershed boundaries output, obtain digital maps (either raster or vector) of verified watershed boundaries from local conservation authorities. This will allow visual and statistical comparison of the output watershed and the verified watershed boundaries.
+
To assess the accuracy 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. Furthermore, statistical comparison of the verified and output watersheds is possible, if need be.
−
In the sample analysis, watershed boundaries were compared to a watershed boundary vector file that was published by Ministry of Natural Resources. The two sets of watershed boundaries are illustrated in Figure 9.
 
   
−
[[File:WatershedVerification.png]]
  
   
  +
[[File:WatershedVerification2.png]]
−
Based on the watershed comparison, r.watershed delineated the Carp watershed fairly accurately in the southern and mid portions of the region. However, in the northern portion of the watershed delineation was inaccurate due to a variety of potential errors. One potential source of error is the use of an inappropriate threshold parameter in r.watersheds while delineating watersheds. Another potential source of error is the possibility that Ministry of Natural Resources used a different DEM, different DEM preprocessing methods, or different software when generating the watershed boundaries.
  
  +
  +
Figure 9: Two watershed outputs (r.watershed, r.water.outlet) with a superimposed watershed verification vector. In the sample analysis, the watershed verification vector was published by Ministry of Natural Resources.
  +
   
   
 
==Verifying Wetness Index:==
  
==Verifying Wetness Index:==
−
If digital topographic maps are available, they provide a good means to evaluate the potential accuracy of output wetness index images. With wetness index and topographic maps on screen, explore relatively wet areas on the wetness index image. Based on the landscapes that are reported by the topographic map, evaluate the accuracy of the wetness index image.
  
−
In the sample analysis, the wetness index display in QGIS was set to be ~50% transparency. A topographic map covering the extent of the watershed was placed behind the wetness index image. Figures 10 and 11 illustrate screenshots during verification. The wetness index image in Figures 10 and 11 use the same color scheme as those in Figures 5, 6, 7, and 8, however since the image is set to be partially transparent, colors appear washed out.
 
   
  +
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 relative accuracy of each wetness index can be evaluated.
−
[[File:TopoVarification1.png]]
  
  +
To evaluate wetness index accuracy, the types of land cover (based on topographic maps) that fall in relatively wet or dry areas (based on wetness index) should be investigated. 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 that an accurate wetness index image was generated. 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.
  +
  +
  +
[[File:WetnessIndexVerification.png]]
  +
  +
Figure 10: Screenshot during wetness index verification, in the northern region of the Carp watershed.
   
−
[[File:TopoVarification2.png]]
  
   
  +
[[File:WetnessindexVerification2.png]]
−
Based on Figures 9 and 10, the wetness index map seems to be fairly accurate in identifying spatial surface moisture patterns. The topographic map verified that most relatively wet features in the wetness index image represent creeks, streams or swampy areas within the watershed. Furthermore, the majority of relatively dry areas seem to be located in topographically high areas. Although most areas that are shown to be wet in the topographic map are indicated as such in the wetness index image, the wetness index appears to have misidentified some wet areas in the landscape (ex: wetlands near Corkery in Fig.10). However, these inaccuracies are likely due to the color scheme that was used to report wetness index, and are not due to errors in the wetness index itself.
 
   
  +
Figure 11: Screenshot during wetness index verification, in the southern region of the Carp watershed.
−
In the sample analysis, verification method was used to evaluate the relative accuracies of the four wetness index images (Figures 5, 6, 7, 8). Generally, the wetness index image that was generated from the DEM, filtered by a 5x5 average x5 times, was found to be the most accurate in predicting spatial moisture patterns throughout the Carp watershed.
  

Revision as of 13:59, 5 December 2010

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 of that region, 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.

DEM.png

Figure 1: A 50m DEM for the Ottawa-Gatineau region that is used in the 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 can 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. 

- Every 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. In the mask, areas within the watershed are assigned a value of 1, and all other areas are assigned a value of 0.

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


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. Identify the x,y coordinates of the watershed outlet point.

5. Open the r.water.outlet module.

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

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

8. 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.


File:Watersheds2.png

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


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 conflict, 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.

Slope1.png

Figure 3: The slope image that was used to calculate wetness index. The user-defined color rules that were used to display the image are also illustrated.


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 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 a single cell 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’]


AccumFlow1.png

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)]


WetnessIndex1.png

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 other images required user-defined color rules for proper visual display.

To generate user-defined color rules:

1. Obtain the distributions for values in the image using: \tA. r.stats –c FileName – returns counts for each value that exists in the image. \tB. QGIS Histogram function.

2. Based on the distributions, determine suitable color breaks for the color rules.

3. Create color rules in a text file (.txt) based on:

ColorRules.png


To apply user-defined color rules:

1. Open the r.colors module.

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

3. 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:

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.

Although slope, flow accumulation and wetness index can be filtered directly after output, by limiting filtering to the DEM only, all output maps are known to be derivatives of a specific elevation model. In other words, by limiting filtering to the DEM, it is very clear which DEM the maps are reporting slope, flow, or wetness for.

In this 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 four DEMs that were filtered 0 times, 5 times, 10 times and 15 times prior to generating slope and flow accumulation images.

To apply a filter to a DEM:

1. Create a filter file which contains the series of filters that are required. 

- In the sample analysis, the filter file contained a single 5x5 average filter; however multiple filters can be included in one filter file.

- In the sample analysis, the filter file was as follows:

File:FilterFile.png


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’].


Images are provided from each output wetness index, to allow investigation of the differences in wetness index output resulting from filtering the DEM with a 5x5 average filter. Generally trial and error is required to find the proper series of filters for your research needs. However, the following images should give an idea on the basic effects of average filtering.

WetnessIndexComparison.png

Figure 6: Wetness indices that were generated from differently filtered DEMs.


NorthernWetnessIndices.png

Figure 7: Wetness indices, generated from differently filtered DEMs, in the northern region of the Carp watershed.


SouthernWetnessIndices.png

Figure 8: Wetness indices, generated from differently filtered DEMs, in the southern region of the Carp watershed.


Verifying Outputs:

Verifying Watershed Boundaries:

To assess the accuracy 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. Furthermore, statistical comparison of the verified and output watersheds is possible, if need be.


WatershedVerification2.png

Figure 9: Two watershed outputs (r.watershed, r.water.outlet) with a superimposed watershed verification vector. In the sample analysis, the watershed verification vector was 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 relative accuracy of each wetness index can be evaluated. To evaluate wetness index accuracy, the types of land cover (based on topographic maps) that fall in relatively wet or dry areas (based on wetness index) should be investigated. 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 that an accurate wetness index image was generated. 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.


WetnessIndexVerification.png

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


WetnessindexVerification2.png

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

Retrieved from "https://dges.carleton.ca/CUOSGwiki/index.php?title=Generating_Wetness_Indices_for_Watersheds_in_GRASS&oldid=330"