The flow direction map can be visualized with d.rast.arrow.
The method adopted to filter the elevation map and rectify it is based on the paper titled "Software Tools to Extract Structure from Digital Elevation Data for Geographic Information System Analysis" by S.K. Jenson and J.O. Domingue (1988).
The procedure takes an elevation layer as input and initially fills all the depressions with one pass across the layer. Next, the flow direction algorithm tries to find a unique direction for each cell. If the watershed program detects areas with pothholes, it delineates this area from the rest of the area and once again the depressions are filled using the neighborhood technique used by the flow direction routine. The final output will be a depressionless elevation layer and a unique flow direction layer.
This (D8) flow algorithm performs as follows: At each raster cell the code determines the slope to each of the 8 surrounding cells and assigns the flow direction to the highest slope out of the cell. If there is more than one equal, non-zero slope then the code picks one direction based on preferences that are hard-coded into the program. If the highest slope is flat and in more than one direction then the code first tries to select an alternative based on flow directions in the adjacent cells. r.fill.dir iterates that process, effectively propagating flow directions from areas where the directions are known into the area where the flow direction cannot otherwise be resolved.
The flow direction map can be encoded in either ANSWERS (Beasley et.al, 1982) or AGNPS (Young et.al, 1985) form, so that it can be readily used as input to these hydrologic models. The resulting depressionless elevation layer can further be manipulated for deriving slopes and other attributes required by the hydrologic models.
In case of local problems, those unfilled areas can be stored optionally. Each unfilled area in this maps is numbered. The -f flag instructs the program to fill single-cell pits but otherwise to just find the undrained areas and exit. With the -f flag set the program writes an elevation map with just single-cell pits filled, a direction map with unresolved problems and a map of the undrained areas that were found but not filled. This option was included because filling DEMs was often not the best way to solve a drainage problem. These options let the user get a partially-fixed elevation map, identify the remaining problems and fix the problems appropriately.
r.fill.dir is sensitive to the computational region settings. Thus the module can be used to generate a flow direction map for any sub-area within the full map layer. Also, r.fill.dir is sensitive to any raster MASK in effect.
In some cases it may be necessary to run r.fill.dir repeatedly (using output from one run as input to the next run) before all of problem areas are filled.
r.fill.dir input=ansi.elev output=ansi.fill.elev direction=ansi.asp
North Carolina sample dataset example: The LiDAR derived 1m elevation map is sink-filled. The outcome are a depressionless elevation map, the flow direction map and an error map.
# set computational region to elevation map g.region raster=elev_lid792_1m -p # generate depressionless DEM and related maps r.fill.dir input=elev_lid792_1m output=elev_lid792_1m_filled \ direction=elev_lid792_1m_dir areas=elev_lid792_1m_error # generate elevation map of pixelwise differences to see obtained terrain alterations r.mapcalc "elev_lid792_1m_diff = elev_lid792_1m_filled - elev_lid792_1m" r.colors elev_lid792_1m_diff color=differences # assess univariate statistics of differences r.univar -e elev_lid792_1m_diff # vectorize filled areas (here all fills are of positive value, see r.univar output) r.mapcalc "elev_lid792_1m_fill_area = if(elev_lid792_1m_diff > 0.0, 1, null() )" r.to.vect input=elev_lid792_1m_fill_area output=elev_lid792_1m_fill_area type=area # generate shaded terrain for better visibility of results r.relief input=elev_lid792_1m_filled output=elev_lid792_1m_filled_shade d.mon wx0 d.shade shade=elev_lid792_1m_filled_shade color=elev_lid792_1m_filled d.vect elev_lid792_1m_fill_area type=boundary color=red
Last changed: $Date: 2017-11-21 19:33:36 -0800 (Tue, 21 Nov 2017) $
© 2003-2018 GRASS Development Team, GRASS GIS 7.2.3svn Reference Manual