Difference between revisions of "Generating Wetness Indices for Watersheds in GRASS"
| (173 intermediate revisions by 40 users not shown) | |||
| Line 1: | Line 1: | ||
| + | =Introduction:= |
||
| − | GENERATING WETNESS INDICES FOR WATERSHEDS IN GRASS |
||
| + | 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. In other words, watersheds 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. Thus, as the outlet of a 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 slope and contributing upslope area, 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 the 'contributing upslope area' parameter. |
||
| − | This tutorial will illustrate the processes involved in generating a topological wetness index from a digital elevation model using GRASS modules through Quantum GIS. Each step involved in generating a wetness index for a watershed will de explained, and example maps from a sample analysis will be provided. In the sample analysis, a wetness index was generated for the Carp River watershed in Ottawa, Ontario using a DEM of the Ottawa-Gatineau region (Fig.1). |
||
| + | =Procedure:= |
||
| − | [[File:DEM.pngx150px]] |
||
| + | This exercise was originally written for use by Carleton University for the purpose of providing a learning tool in the Department of Geography and Geomatics, under the course code GEOM 4008. The first version was written in 2010, and it was checked and updated in 2017. The updated screenshots seen below are from a Windows 7 Operating System. If you do not currently have these programs on your desktop please see the following links: For Quantum GIS download: [http://www.qgis.org/en/site/forusers/download.html] For GRASS GIS download: [https://grass.osgeo.org/download/] |
||
| + | =Outline:= |
||
| − | Once the DEM is loaded in QGIS, create a GRASS mapset with an appropriate projection and workplace extent. Using r.in.gdal.qgis, import the DEM into the new mapset. |
||
| + | 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 index images. |
||
| + | Images from a sample analysis are also provided for each step 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 DEM of the Ottawa-Gatineau region, with 50m cell size, was used for the sample analysis. |
||
| − | 1.0 BACKGROUND INFORMATION: |
||
| − | 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. |
||
| + | [[File:DEM.png|500px]] |
||
| − | 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. |
||
| − | 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. |
||
| + | ''Figure 1: A 50m DEM for the Ottawa-Gatineau region that was used in sample analysis.'' |
||
| − | 2.0 GENERATING WATERSHED BOUNDARIES: |
||
| − | In order to identify the spatial extent of the watershed of interest, use r.watershed to generate an image of watershed boundaries across the DEM. The r.watershed module uses the DEM and a threshold parameter as input. In this case, the threshold parameter refers to the minimum size in cells a watershed must be for it to be delineated. It may take some experimentation to acquire a proper threshold parameter that will produce watershed boundaries that correspond with the watershed outlet that was chosen in your research. |
||
| − | Once the watershed boundaries are output, use them to focus the GRASS extent around the watershed of study (Fig.2). This will reduce processing times and eliminate extraneous data from calculations and visual representation. |
||
| + | =Setting Up Your Workplace= |
||
| − | In the sample analysis, watershed boundaries were delineated based on |
||
| + | 1. In QGIS, open the DEM using the 'Add Raster Layer' button in the 'Manage Layers' toolbar. |
||
| − | [r.watershed elevation=DEM threshold=10000 basin=Watersheds] |
||
| − | [[File: |
+ | [[File:Addraster.jpg]] |
| + | 2. Once the DEM is loaded in QGIS, create a new GRASS mapset with an appropriate projection and workplace extent. |
||
| − | 2.1 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. |
||
| + | 3. Using r.in.gdal.qgis, import the DEM into the new mapset. |
||
| + | [[File:Uploadingtool.png|500px]] |
||
| − | 3.0 GENERATING SLOPE MAPS: |
||
| − | 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). |
||
| − | 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]. |
||
| + | =Generating Watershed Boundaries:= |
||
| − | 3.1 ELIMINATING ZERO SLOPE AREAS: |
||
| + | 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. |
||
| − | 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] |
||
| + | ==''r.watershed''== |
||
| − | 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. |
||
| + | * This method is designed for analysis in which the delineation of all watersheds over a specified size is required. |
||
| − | In the sample analysis, the final slope map was generated based on: |
||
| − | [r.mapcalc Slope_wNoZeros = Slope + SlopeMask] |
||
| + | * If a watershed with a ''specific outlet point'' is required, use the r.water.outlet module. |
||
| − | 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. |
||
| − | [[File:Slope.png]] |
||
| + | To delineate watersheds in r.watershed: |
||
| + | 1. Open the ''r.watershed'' module. |
||
| − | 4.0 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]. |
||
| + | [[File:Rwatershed.png|400px]] |
||
| − | 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. |
||
| + | 2. Select the DEM as the input map. |
||
| − | [[File:AccumFlow.png]] |
||
| + | 3. Select a proper input value. The input value (or threshold in the command line) refers to the minimum size in image pixels a watershed must contain to be delineated. |
||
| + | [[File:Watershed1.PNG|400px]] |
||
| − | 5.0 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)] |
||
| + | 4. Specify a name for 'Output Map: Unique label for each watershed basin'. |
||
| − | 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. |
||
| − | [[File:WetnessIndex.png]] |
||
| + | Alternatively, this can be done in the command line: |
||
| − | 5.1 STANDARDIZING WETNESS INDEX MAPS: |
||
| + | :r.watershed elevation=DEM threshold=# basin=Watersheds |
||
| − | 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))]. |
||
| + | * This watershed output image will delineate ''all'' watersheds that exist within the DEM, given the threshold parameter that was specified. |
||
| − | 6.0 USER DEFINED COLOR RULES: |
||
| + | * Each watershed in the image is assigned a unique pixel value. |
||
| − | 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 |
||
| − | [r.colors map=FlowAccum rules=c:/CarpWatershed/Location/Mapset/FlowColors.txt] |
||
| + | '''Watershed Masks:''' |
||
| − | 7.0 FILTERING THE DEM: |
||
| − | 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. |
||
| − | 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’]. |
||
| + | Since this method produces a watershed map consisting of numerous watersheds, if analysis on a particular watershed is required, a watershed mask will also be required. |
||
| − | 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. |
||
| − | [[File: |
+ | [[File:Findingmapcalc.png|400px]] |
| − | [[File: |
+ | [[File:MapcalcAF.PNG|400px]] |
| + | A watershed mask can be produced in command line based on: |
||
| − | [[File:5x5AVGx15.png]] |
||
| + | :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''== |
||
| − | 8.0 VERIFYING OUTPUTS: |
||
| + | * Unlike r.watershed, this method will delineate a watershed based on the coordinates of its outlet point. |
||
| − | 8.1 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. |
||
| − | 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. |
||
| + | * This allows practitioners to customize the scale of output watersheds based on their research needs, without having to tamper with a threshold parameter. |
||
| − | [[File:WatershedVerification.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. |
||
| + | To delineate watersheds in r.water.outlet: |
||
| + | 1. Open the r.watershed module. |
||
| − | 8.2 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. |
||
| − | [[File: |
+ | [[File:Rwatershed.png|400px]] |
| + | 2. Select the DEM as the input map. |
||
| − | [[File:TopoVarification2.png]] |
||
| + | 3. Specify a name for 'Output Map: Drainage Direction'. |
||
| − | Based on Figures 9 and 10, the wetness index map seems to be 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. |
||
| + | |||
| − | 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. |
||
| + | * The output image quantifies the drainage direction of cells throughout the DEM. |
||
| + | |||
| + | 4. Open the r.water.outlet module. |
||
| + | |||
| + | [[File:Rwateroutlet.png|400px]] |
||
| + | |||
| + | 5. Select the Drainage Direction map as 'Name of Input Raster Map'. |
||
| + | |||
| + | [[File:Watershed2.PNG]] |
||
| + | |||
| + | 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. |
||
| + | |||
| + | |||
| + | |||
| + | [[File:WatershedsComp.png|800px]] |
||
| + | |||
| + | ''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. |
||
| + | |||
| + | [[File:Slope aspect.png|400px]] |
||
| + | |||
| + | 2. Select the DEM for ‘Name of Elevation Raster Map’. |
||
| + | |||
| + | [[File:Slopeopen.PNG|400px]] |
||
| + | |||
| + | 3. Specify a name for ‘Name for Output Slope Raster Map’. |
||
| + | |||
| + | [[File:Slopeopen2.PNG|400px]] |
||
| + | |||
| + | |||
| + | 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 or percent slope, according to the user's choice. |
||
| + | |||
| + | * 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 wetness index images, zero slope areas in 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. |
||
| + | |||
| + | [[File:Findingmapcalc.png]] |
||
| + | |||
| + | 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. |
||
| + | |||
| + | ==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. |
||
| + | |||
| + | [[File:Rwatershed.png|400px]] |
||
| + | |||
| + | 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. |
||
| + | |||
| + | [[File:Findingmapcalc.png|400px]] |
||
| + | |||
| + | 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|600px]] |
||
| + | |||
| + | ''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 is used to calculate wetness index, accumulated flow is illustrated in Figure 4 due to a time constraint. However, due to the nature of the upslope contributing area, 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. |
||
| + | |||
| + | [[File:MapcalcAF.PNG|400px]] |
||
| + | |||
| + | 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|600px]] |
||
| + | |||
| + | ''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. |
||
| + | |||
| + | [[File:Findingrinfo.png|400px]] |
||
| + | |||
| + | 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 and upslope contributing area images are often hard to interpret without user-defined color rules. 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 '' 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 [http://grass.fbk.eu/grass62/manuals/html62_user/r.colors.html r.colors] for syntax of color rules. |
||
| + | |||
| + | |||
| + | 4. Open the r.colors module. |
||
| + | |||
| + | [[File:Findingrcolours.png|400px]] |
||
| + | |||
| + | 5. Select the image in ‘Name of Input Raster Map’. |
||
| + | |||
| + | [[File:Workingrcolours.PNG|400px]] |
||
| + | |||
| + | 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 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 [http://grass.osgeo.org/gdp/html_grass64/r.mfilter.html r.mfilter] for the syntax of filter files. |
||
| + | |||
| + | |||
| + | 2. Open the r.mfilter module. |
||
| + | |||
| + | [[File:Findingrmfilter.png|400px]] |
||
| + | |||
| + | 3. Select the DEM in ‘Name of Input Raster Map’. |
||
| + | |||
| + | [[File:Mfilter1.PNG|400px]] |
||
| + | |||
| + | 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. |
||
| + | |||
| + | [[File:Mfilter2.PNG|400px]] |
||
| + | |||
| + | 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. All indices were generated using the methodologies described in this tutorial. Each was generated from a DEM that was filtered by a 5x5 average filter a varying number of times. The DEMs were filtered 0 times, 5 times, 10 times and 15 times prior to generating wetness index images. Screenshots of each wetness index are provided to allow an investigation into changes in wetness index that result from changes in DEM filtering processes. |
||
| + | |||
| + | Although these screenshots will convey the general effects of DEM filtering, trial and error will likely be required to find the proper series and types of filters for your research needs. |
||
| + | |||
| + | |||
| + | [[File:WetnessIndexComparison.png|600px]] |
||
| + | |||
| + | ''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]'' |
||
| + | |||
| + | |||
| + | [[File:NorthernWetnessIndices.png|850px]] |
||
| + | |||
| + | ''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]'' |
||
| + | |||
| + | |||
| + | [[File:SouthernWetnessIndices.png|850px]] |
||
| + | |||
| + | ''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. Note that, in Figure 9, the watershed delineated by r.water.color is more accurate that that delineated by r.watershed. This is likely due to an improper threshold parameter in r.watershed. |
||
| + | |||
| + | |||
| + | [[File:WatershedVerification2.png|850px]] |
||
| + | |||
| + | ''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. Superimpose the wetness index image (with the transparency of ~60%) over the topographic map. Based on the 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. |
||
| + | |||
| + | |||
| + | [[File:WetnessIndexVerification.png|850px]] |
||
| + | |||
| + | ''Figure 10: A screenshot during wetness index (from 5x5 average filtered DEM x5times) verification, in the northern region of the Carp watershed.'' |
||
| + | |||
| + | |||
| + | [[File:WetnessindexVerification2.png|850px]] |
||
| + | |||
| + | ''Figure 11: A screenshot during wetness index (from 5x5 average filtered DEM x5times) verification, in the southern region of the Carp watershed.'' |
||
| + | |||
| + | =Additional:= |
||
| + | |||
| + | For a practical and fully visual experience of working with QGIS and GRASS GIS in relation to the delineation of watersheds please see this two-part video from YouTube to further knowledge: |
||
| + | |||
| + | https://www.youtube.com/watch?v=oXLPmsvcVLY |
||
| + | |||
| + | https://www.youtube.com/watch?v=DZRhFCMwrgo |
||
| + | |||
| + | |||
| + | |||
| + | [[Category:Tutorials]] |
||
Latest revision as of 15:14, 28 November 2017
Contents
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. In other words, watersheds 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. Thus, as the outlet of a 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 slope and contributing upslope area, 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 the 'contributing upslope area' parameter.
Procedure:
This exercise was originally written for use by Carleton University for the purpose of providing a learning tool in the Department of Geography and Geomatics, under the course code GEOM 4008. The first version was written in 2010, and it was checked and updated in 2017. The updated screenshots seen below are from a Windows 7 Operating System. If you do not currently have these programs on your desktop please see the following links: For Quantum GIS download: [1] For GRASS GIS download: [2]
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 index images.
Images from a sample analysis are also provided for each step 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 DEM 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.
- If a watershed with a specific outlet point is required, use the r.water.outlet module.
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 the command line) refers to the minimum size in image pixels a watershed must contain to be delineated.
4. Specify a name for 'Output Map: Unique label for each watershed basin'.
Alternatively, this can be done in the command line:
- r.watershed elevation=DEM threshold=# basin=Watersheds
- This watershed output image will delineate 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 will also be 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, without having to tamper with a threshold parameter.
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 or percent slope, according to the user's choice.
- 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 wetness index images, zero slope areas in 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.
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 is used to calculate wetness index, accumulated flow is illustrated in Figure 4 due to a time constraint. However, due to the nature of the upslope contributing area, 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 and upslope contributing area images are often hard to interpret without user-defined color rules. 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 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 r.colors 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 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 r.mfilter for the 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. All indices were generated using the methodologies described in this tutorial. Each was generated from a DEM that was filtered by a 5x5 average filter a varying number of times. The DEMs were filtered 0 times, 5 times, 10 times and 15 times prior to generating wetness index images. Screenshots of each wetness index are provided to allow an investigation into changes in wetness index that result from changes in DEM filtering processes.
Although these screenshots will convey the general effects of DEM filtering, trial and error will likely be required to find the proper series and types 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. Note that, in Figure 9, the watershed delineated by r.water.color is more accurate that that delineated by r.watershed. This is likely due to an improper threshold parameter in r.watershed.
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. Superimpose the wetness index image (with the transparency of ~60%) over the topographic map. Based on the 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 (from 5x5 average filtered DEM x5times) verification, in the northern region of the Carp watershed.
Figure 11: A screenshot during wetness index (from 5x5 average filtered DEM x5times) verification, in the southern region of the Carp watershed.
Additional:
For a practical and fully visual experience of working with QGIS and GRASS GIS in relation to the delineation of watersheds please see this two-part video from YouTube to further knowledge:
