Week 6: Surface Analysis

Back To Week 5 and Raster <-> Vector Conversion | On to Week 7 and Feature-level Data Preparation

Overview:


Watershed Delineation

The definition of biologically meaningful management units has become an important topic in natural resources management. The most commonly identified management unit has become the watershed.

A single watershed is a part of the landscape, which drains through a common drainage outlet. All precipitation which falls on a watershed will drain out through this common outlet. A watershed, thus, is defined by its outlet. For example, the Snake River Watershed is defined as the area which drains through a point just upstream of the confluence of the Snake and Columbia Rivers.

Grid provides a set of tools to make watershed delineation quite simple, assuming you get all the steps in the correct order and know a few "tricks." The basic steps are:

  1. FILL. Eliminate areas of internal drainage ("sinks") in the elevation grid. If there are sinks in the elevation grid, watershed delineation will encounter an endless loop and never finish. Use the FILL command to fill sinks in a grid.
  2. FLOWDIRECTION. Create a grid whose value is the direction of flow for each cell in the elevation grid. The flow direction grid will route the flow of water, which is precisely what defines watersheds. Use the FLOWDIRECTION function for this.
  3. FLOWACCUMULATION. Create a grid whose value is the accumulation of the number of upstream cells from each cell in the elevation grid. Knowing where the highest flow cells are will aid in the creation of a hydrologically correct "stream" network, and will also aid in the selection of outlet points. Use the FLOWACCUMULATION function here.
  4. POUR POINTS. Create outlet ("pour") points. These will be used to define upstream areas for each point.
  5. WATERSHEDS. Build watersheds. Use the WATERSHED function.


FILL

The first step in watershed delineation is the filling of sinks:

Usage: FILL <in_grid> <out_grid> {SINK | PEAK} {z_limit} {out_dir_grid}

You can optionally fill sinks (default) or chop off peaks. The {z_limit} allows you to limit the depth of sinks to be filled or the height of peaks to be cut. You can also specify an optional output flowdirection grid.

Grid: fill dem demfill
Filling...
Number of sink(s): 219
Number of sink(s): 39
Number of sink(s): 5
Number of sink(s): -1


FLOWDIRECTION

Direction of flow is calculated in a 3 x 3 cell neighborhood. For the central cell, every one of the 8 adjacent cells is tested for horizontal distance and elevation difference. The direction of the neighboring cell with the greatest downward descent is chosen for the value of the center cell. The directions are assigned according to this 7-bit matrix:

32
64
128
16
 
1
8
4
2

So when a cell flows to the southwest, its flowdirection value is 8.

Usage: FLOWDIRECTION (<surface_grid>, {o_drop_grid}, 
{NORMAL | FORCE})

You can optionally create a {o_drop_grid}, which, "(f)or adjacent cells, ... is analogous to the percent slope between cells. Across a flat area, the distance becomes the distance to the nearest cell of lower elevation. The result is a map of percent rise in the path of steepest descent from each cell." The {NORMAL | FORCE} option will either cause edge cells to flow as though there were no edge (default), or will force edge cells to flow off the edge of the grid (FORCE).

Grid: flowdir = flowdirection (demfill)


FLOWACCUMULATION

Once direction of flow for each cell is known, then a flow accumulation network can be created. This is in essence like placing a marble on each cell and letting all the marbles roll downhill. Imagine a flow gauge at each cell, counting the number of marbles passing through that cell. The number of marbles passing through that cell is the flow accumulation.

Usage:  (F) FLOWACCUMULATION (<dir_grid>, {weight_grid})

FLOWACCUMULATION should always be performed on a flowdirection grid that has been created with the FLOWDIRECTION function, otherwise it can go into an infinite loop and never stop processing. An optional {weight_grid} can be applied, which uses the values from another grid rather than the default value of 1 per cell, and accumulates this flow.

Grid: flowacc = flowaccumulation (flowdir)

 

There are a few extra steps to take in defining flow accumulations. First, we need to create a hydrologically correct stream network in order to select pour points. To do this, we should try several different flow accumulation values, and see what "looks" best.

Grid: mapextent dem
Grid: gridpa dem # linear nowrap gray
Grid: $$display = con (flowacc gt 100, 4)
Grid: arcl str2 3

looks a little too dense

Grid: gridpa dem # linear nowrap gray
Grid: $$display = con (flowacc gt 500, 4)
Grid: arcl str2 3

looks a little better, but I really want the flowaccumulation-generated streams to look similar in extent to the mapped streams, so I know I'm selecting a path with a high amount of flow volume. Let's zoom in a bit and take a closer look

Grid: mape 1185458 554872 1190607 560324
Grid: $$display = con (flowacc gt 1000, 4)
Grid: arcl str2 3

let's try one more

Grid: $$display = con (flowacc gt 5000, 4)
Grid: arcl str2 3

 

OK. If this does not work, we can still do it over

Grid: flow5000 = con (flowacc gt 5000, 1)
Grid: stream5000 = streamline (flow5000, flowdir)

Now we have a grid network representing high flow paths and also a line coverage representing the same data.


POUR POINTS

Pour points will be digitized using the ArcEdit module. In essence, ArcEdit is the most common module used to bring coordinate data into ARC/INFO. ArcEdit is complex, but we will be using some very basic functionality, so we do not need any major introductions.

Grid: arc ae	/* starts ArcEdit
Arcedit: create pour_pts stream500	/* create a new coverage, based on existing tics & boundary
Arcedit: image flow5000 transparent 4	/* draw the grid in the background
Arcedit: bc str2 6	/* draw the "real" streams in the background
Arcedit: bc stream5000 3	/* draw the derived stream net in the background
Arcedit: be arc	/* set what background features to draw
Arcedit: de point	/* set the foreground features to draw
Arcedit: sds 2	/* selected features in red
Arcedit: draw

This is what you will see when you draw. Zoom in closer with a zoom window (zoom to the box shown in the above image), and you will see the actual data sets (flow5000 in blue, stream5000 in cyan, and str2 in purple):

Arcedit: ef point

Before you add any points, zoom in even closer. I've added points along high-flow paths, but far enough away from each other so as not to affect the SNAPPOUR function. The points are added where indicated by the arrows.

Arcedit: add

Once you have added the points you want, quit and save.

Arcedit: q
Keep all edit changes (Y/N) : y
This will replace all changes back into the original coverage(s) and/or INFO file(s)
Do you really want to do this (Y/N) : y
Saving changes for E:\ARC_LMS\POUR_PTS
Saving labels...
4 label(s) written to E:\ARC_LMS\POUR_PTS
  from the original 0, 2 added and 0 deleted
   BND replaced into E:\ARC_LMS\POUR_PTS
Leaving ARCEDIT...

Now that we have pour points, we can create watersheds.


WATERSHEDS

To create watersheds, we need a grid of pour points and a flowdirection grid. We already have the flowdirection grid, and we can create the pour point grid on the fly by embedding the POINTGRID function within the WATERSHED function. We will also use the SNAPPOUR function to make sure that the pour points are actually on high-flow cells.

Grid: setcell dem
Grid: sheds = watershed (flowdir, snappour (pointgrid (pour_pts), flow5000, 10))
Grid: gridsh sheds
Grid: arclines stream5000 3

(note: this image's colors have been altered)

Once you have created watersheds, take a look at them. If they do not look correct, it is most likely because of poor pour point placement. Kill the watershed grid, start ArcEdit again, and edit the pour points. Then rebuild the watersheds grid. Once you have watersheds delineated, you can perform any of the other typical analyses we've already seen, such as road and stream densities, summary of forest stands, etc. on a watershed-by-watershed basis.


Slope and Aspect Calculation

Slope and aspect are easily calculated in Grid.

The build-in SLOPE function calculate slope in either degrees (default) or percent.

Usage:  (F) SLOPE (<grid>, {DEGREE | PERCENTRISE})
        (F) SLOPE (<grid>, <z_factor> {DEGREE | PERCENTRISE})

The output grid will consist of values which represent the slope of a plane tangent to the input grid at each lattice mesh point. You can specify a <z_factor> which will multiply a constant to the cell values in the output grid.

Once you have a slope grid, you may classify it as you wish to create another grid. If you want to create a classified slope grid, first create a slope grid and then classify the output:

Grid: slp = slope (dem, percentrise)
Grid: slp50 = con (slp ge 0 and slp lt 50, 1, 2)

rather than placing the SLOPE function within the conditional:

Grid: slp50 = con(slope (dem, percentrise) ge 0 and slope (dem, percentrise) lt 50, 1, 2)

Although the latter will work, it causes the Grid processor to calculate slope twice rather than once. This may not make a difference in a 2-class output, but it will certainly make a difference in a multi-class output. The same rule holds especially true for any of the Grid functions that are processor-intensive.


Creating a grid of aspect is just as easy.

Usage:  (F) ASPECT (<grid>)

The aspect grid will have values equal to the compass direction of the steepest line on a plane tangent to each lattice mesh point.

Grid: asp = aspect (dem)

Although aspect is easy for the user, it is processor-intensive, so the above rule about efficient use of processor time also holds true.

The only tricky part of using aspect is in generating zonal means. Because aspect values are not absolute, the means cannot be calculated using an arithmetic mean. We will perform this task next week.


It is possible to create a series of surface grids from a single function: CURVATURE.

Usage:  (F) CURVATURE (<grid>,  {o_pro_curv}, {o_plan_curv},
                       {o_slope}, {o_aspect})

The default output from the CURVATURE function is a grid representing the shape of the grid at each mesh point. "A positive curvature indicates that the surface is upwardly convex at that cell. A negative curvature indicates that the surface is upwardly concave at that cell. A value of zero indicates that the surface is flat." In addition to the slope curvature output, you can also create optional outputs of profile curvature, planform curvature, slope (in degrees), and aspect.


Exercises

Based on the digital elevation model of Pack Forest, create watersheds at the confluences of streams which have contributing areas of at least 100,000,000 m^2 (flowaccumulation = 100,000).

What is the area of each watershed?
What is the road and stream density within each of these watersheds?
What is the timber age distribution within each watershed?


Back To Week 5 and Raster <-> Vector Conversion | On to Week 7 and Feature-level Data Preparation