Watershed Analyses

To determine a stream network, basins, or flow characteristics (amount, timing, etc) from a DEM, a series of steps are taken in that allow the GIS user to establish the characteristics of water flow across the DEM surface. These several steps need to be taken, in order, to complete a hydrological analysis of from a DEM dataset.

  1. Eliminate pits
  2. Determine flow direction
  3. Calculate flow accumulation
  4. Determine streams
  5. delineate watersheds

Some tools in QGIS combine many of these steps into one process (e.g., GRASS r. fill.dir does 1 & 2, SAGA Channel Network and Drainage Basins does them all).

  1. Eliminating “pits” or “sinks” in the DEM
    1. Where are they? Undrained depressions prevent most watershed algorithms from working. Sinks can be natural, such as “sinkholes” or reservoirs behind dams, but they’re also due to random errors in datasets.
    2. Most GIS routines “fill” the sinks numerically such that computational “water” will drain across the sink to the next lowest point. Other GIS routine involve a combination of fill and breach, wherein the depression is filled as the cell in the steepest path downhill is abraded.
    3. QGIS filling routines are in GRASS, SAGA, and PCRaster and vary in their method. But mostly they are the “fill only” variety.
    4. ArcGIS’s approach: If you look up the process for ArcGIS is says that it does a combination of the two, but I don’t think the ArcGIS routine “remove peaks” as shown below. Filling tends to make stream profiles into stair steps. Using the tools in the Hydrology toolbox within the Spatial Analyst toolbox, you can execute a “fill pits” routine that iteratively calculates flow direction.

      Here is an example of a river profile using raw DEM data containing pits (from Harbor et al., 2005 paper on my website). This “route” or stream path follows cell to cell the lowest downstream outflow from each cell to create a stream path, but it is full of errors.
  2. Determining flow direction (from original surface and/or “filled” DEM)
    1. Most routines follow the steepest decent, rather than just the lowest cell in the nine cell neighborhood. The direction of flow is determined by finding the direction of steepest descent from each cell–change in z value / distance, where the distance is determined between cell centers. Therefore if the cell size is 1, the distance between two orthogonal cells is 1 and the distance between two diagonal cells is 1.414214.”
    2. For flow direction, most GIS routines the “D-8 method” (selecting one of the nearest 8 cells), but the way in which they’re coded in the output flow direction layer varies greatly.
      alternate flow direction numbering schemes
      Because they differ, it is important to stick with one set of tools as you work through the flow routing routines!
    3. Pits are easily found with a flow direction map, because flow is “trapped” where two or more cells point at each other, as shown below.
    4. Does this approach work everywhere? How about for this common geomorphic feature in Death Valley?

      https://goo.gl/maps/9JdjTEhe6xAuSgPCA
      There are several approaches to flow direction calculation in the literature. Some partition flow to only one cell, some divide it up based on the proportion of flow into each cell (usually 3 or less neighboring cells).
    5. It is possible, using what is known as the “D-infinity” method developed by David Tarboton at Utah State as part of TauDEM to partition the flow incrementally to many cells in the downslope direction. And TauDEM used to have a QGIS plugin, but I think it is toast. The GRASS routine r.fill.dir uses D* but r.watershed has the option to use multiple flow directions (MFD) based on D-infinity.
  3. Computing Flow Accumulation
    1. Flow accumulation is the total number of cells that drain into a given cell. The value of zero or one (ArcGIS vs Q) is given to the ridge lines, where flow originates. This determined from sequential analysis of flow direction (here created by PCRaster lddcreate tool). This routine filled the pit in the SE valley to create the outflow see in the mesh arrows.
    2. not all algorithms do this the same; the “SAGA Flow Accumulation (One Step)” tool forces edge cells off the map, and therefore ignores them in the calculation (even when they flow in to the map!)
      SAGA flow accumulation
      so this returns an area of 6400m2, or 64 cells draining to the SW vs 96 for the PCRaster version of flow accumulation.
  4. Determining Streams – Two strategies will get you to “streams” from the above inputs based on a threshold value of Flow Accumulation or Strahler order. Flow accumulation threshold has a better basis in geomorphic understanding of what defines the start of a stream.
    1.  Using Flow and/or accumulation We use a simple boolean statement on flow accumulation to define a stream as greater than a number of cells or area, depending on the output of the tool used. Based on PCRaster output of flow accumulation in cells , we can use this statement in the raster calculator
      (“accum@1” >= 10) / (“accum@1” >= 10)
      This works because where the value >=10, it returns 1 or true, so 1/1 =1
      but where the value is <10, it returns a 0 for false, and 0/0 cannot be calculated so it becomes “null.”
      stream delineation in QGIS
      Those streams can then be made into vector lines using GRASS r.to.vect (perhaps preceded by r.thin where two streams run side by side to make a cell path greater than one cell wide).
      lots of streams in QGIS
      –note: r.watershed will do most of these steps, AND MORE, in one pass. See the explanation here.
    2. Using Strahler Order The stream order method of stream “accounting” was developed from paper maps originally, but is now applied to cell-to-cell accounting in raster GIS. Any stream that starts on a hillside is a “first-order” stream, where two first orders meet, they form a 2nd order. More 1st order streams can join the 2nd and it stays a 2nd until it meets another 2nd to create a 3rd order, etc, etc. The above image has PCRaster’s stream order representation of our tiny watershed. Whereas this gives a 4th order for the stream leaving the SW corner, SAGA’s Strahler Order tool yields a 2nd order stream because 1) the area is smaller because edge cells weren’t counted, 2) only cells where two “stream lines” converge are 1st order and the rest are zero-order, and 3) I didn’t use the pit-filled DEM for this series of calculations so the pit gets some of the territory.
      saga gis stream order
  5. ARCGIS: Determining streams 
    1. Our hydrology expert needs to decide a “threshold” value of accumulation required to create a stream; how much area is needed to create enough water to have a stream. Here this calculation is accomplished in one of two ways using the “raster calculator” tool.
    2. where SetNull( conditional statement, value if false) –if the conditional statement is true, the cell is set to “NoData” or “Null” (where NoData is actually a number you don’t expect to find in the data.  See the properties box to find the NoData value)
      or
      or using the Con function: “Con(conditional statement, value if true, value if false).
      If you omit the “value if false” and the conditional statement is false, the cell value is set to NoData or “Null” .
    3. The “Basin” tool finds the watersheds draining to different end points in the DEM, the “Watershed” tool finds basins above selected points (help file)
    4. GIS operators can also define a “pour point” where the water leaves one or more watersheds and the GIS will define the watershed above each pour point.

Hydrologic modeling has important uses in

  • flood warning (real-time rainfall estimates is derived from RADAR and delivered to a GIS-based watershed model).
  • channel and navigation studies
  • water quality modeling (nutrient loading to the Chesapeake Bay, for example)
  • geomorphic study of erosion and uplift