Contents

• Reference material [This is a list for extended reading.  You are not expected to read this before doing the homework, or even at all in its entirety.  You should skim this material and be ready to search in this material for help or clarification when necessary.]
• Installation [This may require system administrator help (or at least permission) for networked configurations.]
• Exercise
• Summary of material to turn in.
• Data file, ex1.zip, with datasets used.  This is the same data as used for exercise 1.

Reference Material

• Tarboton, D. G., (1997), "A New Method for the Determination of Flow Directions and Contributing Areas in Grid Digital Elevation Models," Water Resources Research, 33(2): 309-319. (available online at http://hydrology.usu.edu/dtarb)
• TARDEM documentation.
• Tarboton, D. G., R. L. Bras and I. Rodriguez-Iturbe, (1992), "A physical basis for drainage density," Geomorphology, 5(1/2): 59-76.
• Tarboton, D. G., R. L. Bras and I. Rodriguez-Iturbe, (1991), "On the Extraction of Channel Networks from Digital Elevation Data," Hydrologic Processes, 5(1): 81-100. (Also reprinted in Terrain Analysis and Distributed Modeling in Hydrology, Edited by K. J. Beven and I. D. Moore, John Wiley & Sons, 1993, p.85-104.)

Installation.

1. SINMAP extension.  First check if SINMAP is already installed.  From the File/Extensions menu in Arcview see if SINMAP is there, and if it is there select it.  If SINMAP is not there do the following.  Download sinmaplight.zip.  This zip file contains 'sinmap.avx' and 'sinmap.dll'. sinmap.avx is the ArcView extension file used by SINMAP.  This file must be placed in the folder C:\ESRI\AV_GIS30\ARCVIEW\EXT32 (or equivalent if ArcView is installed elsewhere). sinmap.dll is the binary dynamic link library file containing SINMAP routines.  It must be placed in the folder C:\ESRI\AV_GIS30\ARCVIEW\BIN32. You should now be able to select the SINMAP extension from the file/extensions menu in ArcView.
2. TARDEM.  First check if TARDEM is already installed.  TARDEM programs are DOS style programs run from a command line, so open a MS-DOS Prompt  (Start/Programs/MS-DOS Prompt).  In this DOS window type 'flood'.  If you get the message "bad command or file name" then the programs are not (properly) installed.  If you get the message "Usage: ..." then the programs are installed.  To install download and save to disk tardem.zip.  This file contains several *.exe files that should be placed in any convenient folder, e.g. c:\tardem\myexec.  To use these programs this directory needs to be on the system search path.  Since these programs use the ESRI gridio library this library also needs to be on the system search path.  Both these requirements can be met by including the following line in the autoexec.bat file.

Windows 95/98
SET PATH=C:\TARDEM\EXEC;C:\ESRI\AV_GIS30\ARCVIEW\BIN32

Windows NT
path=c:\tardem\exec;c:\esri\av_gis30\arcview\bin32

There may be other path entries on your system depending on the software you have.  You should leave them alone and only add the entries involving TARDEM and ARCVIEW.  Once these changes have been made the system should be rebooted.  Make sure there are no blanks or trailing blanks in these paths or it won't work. The autoexec.bat file can be edited from the Windows Start Menu Start/run/sysedit.

Exercise

This exercise will use the Reynolds Creek DEM, reydem.asc from exercise 1.  This should already exist on your computer as an ESRI grid named reydem, from exercise 1.

To identify and delineate a watershed we will use TARDEM.  Open a MS-DOS prompt (Start/Programs/MS Dos Prompt) and navigate (using cd) to your working directory (e.g. cd c:\giscourse\reydata\).  Now at the prompt enter the command

tdprepro reydem
This runs a batch file that runs four programs (flood, d8, aread8 and gridnet) that fill pits, compute flow directions and slopes (by the D8 method), compute contributing area and compute grid network order and length attributes.  The files produced are identified by a suffix appended at the end of the name reydem that was given on the command line.  (You may ignore the VAT (value attribute table) warning message - we do not need a VAT, it occurs because the contributing area has a large number of different integer values - too many for ArcView to treat as attribute table categories.).    The programs could also have been run on the ASCII grid file using
tdprepro reydem.asc    (you do not need to do this)
If you do this ignore the 'invalid grid item syntax' message.  This is saying the program could not open the grid as an ESRI grid at which point it reverts to reading it as an ASCII file.  This will be a bit slower because ASCII input/output is used.  The programs may also be run individually using the commands
flood reydem   (you do not need to do these - they were done by tdprepro)
d8 reydem
aread8 reydem
gridnet reydem
This is useful if there was a problem part of the way through and you do not want to rerun work already done, or if you know you will be using a channel network delineation approach that does not use some of the output, e.g. the constant support area method does not use gridnet output.

Following are the suffixes for the grid files produces and what they contain:

no suffix.	Elevation data.
fel	Pit filled elevation data.	produced by flood	input to D8, Dinf, netsetup
p	D8 drainage directions.	produced by D8	input to aread8
sd8	D8 slopes.	produced by D8	input to netsetup, some methods
ad8	D8 contributing area’s, units are number of grid cells.	produced by aread8	input to netsetup
slp	Dinf slopes.	produced by Dinf
ang	Dinf flow directions.	produced by Dinf	input to areadinf
sca	Dinf contributing area’s, units are specific catchment area, i.e. number of grid cells times cell size.	produced by areadinf
plen	Longest path length to each grid point along D8 directions.	produced by gridnet
tlen	Total path length to each grid point along D8 directions.	produced by gridnet
gord	Strahler order for grid network defined from D8 flow directions.	produced by gridnet
src	Network mask based on channel source rules.	produced by netsetup
ord	Grid with Strahler order for mapped stream network.	produced by netsetup
w	Subbasins mapped using subbasinsetup.	produced by subbasinsetup
fdr	Flow directions enforced to follow the existing stream network	produced by streamtogrid	optional input to flood
fdrn	Flow directions enforced to follow the existing stream network after cleaning to remove any loops	produced by flood	optional input to d8, dinf

The .asc extension is used if the data is ASCII.

If you look in the directory where you are working you should see grid files (which are folders in Windows) with name suffixes fel, p, ad8, sd8, plen, tlen, and gord.  Start ArcView and load the Spatial Analyst Extension (choose Extension from File menu and check Spatial Analyst from the available extensions).  Open a view and load the grid theme reydemad8. ( add theme button). This is the contributing area grid computed using the D8 method.  Fiddle with the color scheme legend so that you can see the channel network.  Drag the legend bars for the stream reach shape file and watershed outline shape file rcout.shp to above the reydemad8 theme so that stream reaches and watershed outline are visible above the contributing area grid theme.  Note the correspondence (or lack of it) between stream reaches and large contributing area values.  It is this correspondence that is the basis for the contributing area threshold method for delineating channel networks.  Note also the black (no data) jagged edges to the contributing area theme.  This is because the contributing area computation checks for edge contamination.  This is the possibility that a contributing area value may be underestimated due to grid cells outside of the domain not being counted.  This occurs when drainage is inwards from the boundaries.  The algorithm recognizes this and reports no data.  This may be overridden by the option -nc in aread8 in cases where you know this is not an issue, if for example the DEM has been clipped along a watershed outline.

Zoom in on the outlet and use Identify  to select a grid cell on the outlet with large contributing area

Write down the outlet coordinates, in this case x=520170 y= 4789800.

Now enter the command

netsetup reydem -m 1 500 -xy 520170 4789800
This maps channel networks using method 1 (constant contributing area) with a threshold area of 500 grid cells using outlet coordinates specified.  This results in two new grid files being created, reydemord and reydemsrc.  The first contains channel Strahler stream order, and the second is an intermediate mask file used to indicate the channels mapped.  This command also results in two new text files reydemtree.dat and reydemcoord.dat.  These define the reach linkages and attributed for the mapped channels.  Details of these file formats are given in the TARDEM documentation.  Add the reydemord theme.  Experiment with some other contributingarea thresholds.  Report the highest stream order for contributing area thresholds of 100, 500, 1000.

The complete set of methods for defining channel networks via netsetup are

1. Catchment area threshold A >= p[0].
2. Area-Slope threshold A S^p[1] >= p[0].
3. Length-Area threshold A >= p[0] L^p[1]. Here L is the maximum drainage length to each cell in the plen file.
4. Accumulation area of upward curved grid cells.  The DEM is first smoothed by a kernel with value p[0] at its center, p[1] on its edges, and p[2] on diagonals.  The Peuker and Douglas (1975) method is then used to identify upwards curved grid cells and contributing area computed using only these cells.  A threshold, A_uc >= p[3] on these cells is used to map the channel network.
5. Grid order threshold  O >= p[0].
6. Use existing channel networks specified in a *fdrn file [fdrn file created from fdr file by flood.  fdr file created from shape file by streamtogrid]
Experiment with a few of these.  You may want to rename the files reydemtree.dat and reydemcoord.dat between runs to save them, as each new run overwrites the old ones.  My favourite is method 4, with parameters  .4 .1 .05 20.

With a channel network defined enter the command.

subbasinsetup reydem 1
This maps subwatersheds draining to each stream order reach controlled by the order threshold parameter given.  Streams of order lower than the threshold given are stripped off the network before mapping resulting in a coarser watershed delineation.   The output is a grid reydemw, an integer grid giving subwatersheds, reydem.shp a stream network shape file and reydemw.shp a shape file of subwatershed boundaries.  Add these themes into arcview to examine them.

Now enter the command

netsetup reydem -m 4 0.4 0.1 0.05 20 -xy 520170 4789800  (If you are using different data the outlet coordinates specified will need to be different)
This delineates a channel network according to the upwards curvature method with threshold 20.  You may run subbasinsetup to display this if you like.  Now enter the command
streaman reydem
This produces an output file reydemst.dat that contains statistical properties of the "Strahler Streams".  These are stream segments of the same order.  The file looks like:
Streams with amin =  0.00000E+00
 7
 first link
 last link
 order
 length
 drop
 g. length
 AREA
    530     0     5  0.14018E+04  0.00000E+00  0.12987E+04  0.23912E+09
     74    74     2  0.62699E+03  0.16000E+02  0.57940E+03  0.76140E+06
     58    58     1  0.70456E+03  0.63000E+02  0.56365E+03  0.45360E+06
...
The seven columns are defined in the header as first link number, last link number, order, length, drop, geometric length and area, respectively, for each Strahler stream.  We are going to perform a constant drop analysis, so focus on columns 3 and 5.  Plot the stream drop versus order, for each stream and for the mean of all streams of the same order.  To do this you will need to import this data into software such as Excel, Matlab, Splus or Mathematica (whatever you are comfortable with for plotting).  The result should look something like the figure below.

In this figure the mean drop for each order has been offset from the individual points for display purposes.  The first order stream drops seem to have a mean less than the higher orders (at least visually and discounting the single 5th order stream which is not a representative sample).  This can be tested using the t test for the difference between means.

Here and  are the means of the first and higher order streams respectively.  n_x and  n_y are the sample sizes and s_x and s_y the sample standard deviations.  With the above data this evaluates to -5.3.  Roughly speaking the difference is statistically significant when |t| (the absolute value of t) is greater than 2.  Look at a statistics book to be more precise.

Repeat this analysis for channel networks generated with thresholds of 30 and 50, i.e. with the commands

netsetup reydem -m 4 0.4 0.1 0.05 30 -xy 520170 4789800
netsetup reydem -m 4 0.4 0.1 0.05 50 -xy 520170 4789800
Find a threshold for which the t statistic measuring the significance of the difference between first order and higher order streams is not significant (i.e. |t| less than 2).  Report this threshold.  (If you are using different data the you may also need to try a few different thresholds beyond the set 20, 30 and 50, say 5, 10, 15, 80, 100, 150.  The right one depends on the drainage density or texture of the topography and resolution of the DEM).  This threshold defines is a channel network that is consistent with Horton's laws.  Other channel networks inconsistent with Horton's laws may be mapping as 'channel' parallel flow down hillslopes. Print a layout of the channel network and watershed mapped in this way.  Determine the drainage density (total channel length/watershed area) and compare the drainage density to the drainage density of the EPA reach files mapped channel network.  Report the drainage density of your mapped channel network and the EPA reach files channel network.

Summary of answers to turn in

• The highest stream order for contributing area thresholds of 100, 500, 1000.
• The stream drop versus order plots, for thresholds 20, 30 and 50.
• The threshold contributing area of upwards curved pixels that defines a channel network for which the constant drop property holds.
• A layout of the channel network mapped and watersheds delimited using the threshold on upwards curved pixels.
• The drainage density of the network mapped using the threshold on upwards curved pixels.
• The drainage density of the EPA reach files channel network.