Watershed analysis

At this point of the tutorial it may be useful to experiment with some more advanced GRASS functionality for hydrologic analysis.

First of all reset the region settings with g.region making them match the DTM map elevation.dem (from the window g.region set elevation.dem in the field Set region to match this raster map, making sure that the resolution is 30 meters).
For this exercise the map elevation.dem will be used as DTM, aspect as aspect and slope as slope.

Sub basins extraction

Display the map elevation.dem (DTM). The area can be divided into a set of basins and sub basins. This operation can be performed automatically using the r.watershed (Raster -> Hydrologic modeling -> Watershed analysis) module to identify all the basis with a minimum area of 1 square km. The r.watershed module can perform further analysis which will be skipped now: some of the parameters can be left to their default values and some of the fields can be left blank (the corresponding maps will not be created).
Conversely, the following parameters must be provided:

r.watershed module interface

r.watershed can be also run in background; you can check the r.watershed module options using g.manual (note that in some cases r.watershed performs better on DEMs which have been smoothed using a neighborhood filter using the mean operator. This can be done using the r.neighbors module; this operation is not necessary in our case.).

basin raster map

The output map from r.watershed identifies each basin with a different integer number, assigned to all the cells within the basin; this value can be detected using the already described procedure for querying raster maps.

r.watershed gives to the map a "random" color table. Using the r.colors module (Raster -> Manage map colors -> Set colors to predefined color tables), it is possible to change the color table, for example using the "rainbow" color table like in the figure below (remember to refresh the display by clicking the icon on the Map Display), the color sequence hints at how the basins are numerated starting from the gauge point of each main drainage.

r.colors interface

basins raster map and its sequential sorting of the watershed basins

Sometimes it is useful to better understand the results to "merge" the output map with the aspect map while displaying. To do this, in the GIS Manager use the raster aspect as Base map and, in this case, basins as Drape map.

GIS.m: shadowing effect

Merge of the aspect map and basins

River network extraction

The accumulation map evaluates the number of cells draining in each cell. Negative values indicate that a cell receives flux from an area outside the current region.

The accumulation map

It is possible to detect the main channels, but not much more. Using r.mapcalc it is possible to compute the logarithm of the absolute value of the accumulation, an useful parameter for hydrologic calculations. On the console, write:

> r.mapcalc 'log_accumulation=log(abs(accumulation)+1)'

Display, then, the new map (log_accumulation):

Accumulation logarithm

Use again r.mapcalc to extract from this map the river network: a good threshold value is 6. The value 1 is assigned to the cells with a value whose logarithm is larger than 6, 0 to the other ones.

> r.mapcalc 'inf_rivers=if(log_accumulation>6)'

The inf_rivers raster map

Finally, to create a vector map of the river network, use r.thin from the menu, selecting Raster -> Transform features -> Thin linear features to thin out to a single pixel width the lines on the binary raster map (inf_rivers): the new raster map thin_rivers will be created.

Using the module r.to.vect (selecting File -> Map type conversions -> Raster to vector map on the menu) it is possible to create a new vector map of the river network, named here rivers.

r.thin module interface

Display the vector map rivers using the icon in the Display Manager (as seen before when transforming contour lines from raster to vector).

rivers vector map

The rivers map contains, at the upper and lower borders, two "horizontal" lines which are not part of the river network: they are due to the fact that the two outermost basins are not close in the current region. To avoid this pathologic situation, only the basins falling completely inside the region should be evaluated. Use the v.digit module: after selecting the rivers vector map in the GIS Manager click on the icon in the Map Display window. The following new window appears:


Select the icon and, clicking twice on the horizontal lines, delete them. At the end of the operation, click with the right button in the digitizing window and exit from the digitizing module by clicking on the icon and closing all the windows which have been created by the v.digit module.

Overlay this new vector map of the river network to the basins (raster) map: it is possible now to detect the watersheds. The result is shown below:

Basins and river network

It may be interesting to create a vector map showing the basins boundaries.
Select on the menu File -> Map type conversions -> Raster to vector map and set the parameters as shown in the next figure (note the line feature type, using the area feature the basins areas and not the boundaries would be created in the output vector map):

r.to.vect interface

Display the new vector map basins: try to hide the "centroids" and "boundaries" elements and to assign random colors to the areas (in the next figure the map rivers have been overlayed to give a more complete view).

Display of the basins and rivers vector maps


The suggested resolution for this exercise is 30 meters: it may be interesting to try the same analysis with different resolutions, for example 10 or 50 meters.