Figure 1. Flow direction defined as steepest downward slope on planar triangular facets on a block centered grid.
David G. Tarboton
May 5, 2000
Utah State University
http://www.engineering.usu.edu/dtarb/
This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License version 2, 1991 as published by the Free Software Foundation.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
A copy of the full GNU General Public License is included in file gpl.html.
This is also available at:
http://www.gnu.org/copyleft/gpl.html
or from:
The Free Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA.
If you wish to use or incorporate this program (or parts of it) into
other software that does not meet the GNU General Public License conditions
contact the author to request permission.
David G. Tarboton
Utah State University
8200 Old Main Hill
Logan, UT 84322-8200
USA
http://www.engineering.usu.edu/dtarb/
email: david.tarboton<at symbol>usu.edu
All programs take command line arguments. When run without a command line argument, they tell you what arguments they are expecting, or prompt for input.
Streamtogrid
This takes as input an EPA style river reach file in 'shape file' format
and a DEM grid, and outputs a grid of flow directions along the rivers
to enforce drainage along delineated river channels.
Tdprepro
A batch file that runs Flood, D8, Aread8 and Gridnet. This takes
as input an elevation data grid and optionally a river flow enforcement
direction grid.
Flood
This takes as input an elevation data grid and outputs a grid file
‘fel’ with pits filled, using the flooding algorithm. The river flow
enforcement direction grid is an optional input, and if input, pits will
be filled consistent with drainage along existing streams.
D8
This takes as input pit filled elevation data file 'fel' and outputs
D8 flow directions 'p' and slopes 'sd8' for each grid cell. The following
coding is used for direction of flow from a pixel to ONE of its eight neighbors:
4 3 2
5 1
6 7 8
In flat areas flow directions are assigned away from higher ground
and towards lower ground using the method of Garbrecht and Martz (1997).
The river flow enforcement direction grid is an optional input, and if
input, directions will be set consistent with drainage along existing streams.
Dinf
The Dinf approach assigns a flow direction based on steepest slope
on a triangular facet (Tarboton, 1997). This is saved as an angle ‘ang’
in radians anti-clockwise from east. In flat areas the D8 flow directions
are converted to angles and used. The river flow enforcement direction
grid is an optional input, and if input, directions will be set consistent
with drainage along existing streams.
Aread8
This takes as input a D8 flow directions file ‘p’ and outputs the contributing
area. The result is the number of grid cells draining through each grid
cell. Optional command line arguments for the outlet coordinate result
in only the area contributing to the designated outlet being calculated.
By default the program checks for edge contamination.
The edge contamination checking may be overridden with the optional command
line argument -nc.
Areadinf
This takes as input a Dinf angle file ‘ang’ and outputs the specific
catchment area. Specific catchment area is defined as contributing area
per unit contour length. Here the contour length is taken as the
grid cell size. The result has length units the same as grid cell
size. Optional command line arguments for the outlet coordinate result
in only the area contributing to the designated outlet being calculated.
By default the program checks for edge contamination.
The edge contamination checking may be overridden with the optional command
line argument -nc.
Gridnet
This takes as input a D8 flow directions file 'p' and outputs three
grid files:
'plen' Each grid cell contains the path length from the furthest
cell that drains to each cell.
'tlen' Each grid cell contains the total length of all paths
draining to each cell.
'gord' Each grid cell contains the Strahler order associated
with that cell for a flow network defined using the D8 flow directions
and including each grid cell. Strahler order is defined as follows.
Cells that don't have any other grid cells draining in to them are order
1. When two (or more) flow paths of different order join the order
of the downstream flow path is the order of the highest incoming flow path.
When two (or more) flow paths of equal order join the downstream flow path
is increased by 1. Algorithmically this is implemented as:
Order = Max(Highest incoming flow path order, Second highest incoming
flow path order + 1)
These outputs are used by netsetup in some of the options for mapping
channel networks.
Netsetup
This takes input from D8 flow directions and contributing area grids,
the pit filled elevation grid and needed grids from 'plen', 'gord', 'sd8'
depending upon the method used to map channel networks. Typical command
line input is:
netsetup file_name_prefix -m 4 .4 .1 .05 20 -xy outlet_x Outlet_yA method as well as outlet coordinates must be specified. I suggest identifing outlet coordinates using the identify tool within a View in ArcView. The methods available range from a specified contributing area threshold to the use of upward curved grid cells (e.g. Tarboton and Shankar, 1998) and grid network pruning by order (e.g. Peckham, 1995 as used in RiverTools). My current favorite method uses upward curved grid cells with the parameters above. If a constant support area threshold is used I recommend selecting it objectively using the methods of Tarboton et al. (1991, 1992). Output consists of two grids and two ascii files defining the channel network. The grid files are the 'src' file (essentially an intermediate file) containing the contributing area computed using only pixels that meet the channel network extraction threshold for the method being used, and the 'ord' file containing the Strahler stream order of grid cells on the channel network. The ASCII files are the 'tree.dat' file containing a list of links (channel segments) comprising the channel network, with connectivity and other information, and the 'coord.dat' containing the vector points along each channel segment.
The methods for defining channel networks 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, Auc >= 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]
Arclinks
This takes input from the 'tree.dat' and 'coord.dat' files that define
a channel network and imports them into ArcView as a Vector feature file.
The vector feature file has a single item for each link (stream segment
between junctions). An intermediate *.e00 file is generated and the 'import71'
program run. Attribute information is included with these files allowing
stream properties to be identified. The program also allows a subnetwork,
or network with different support area threshold to be imported through
the
specification of a link number (from the links field in the attribute table)
using option -l and alternative support area threshold, using option -am.
Arcstreams
This is similar to arclinks in that it takes input from the 'tree.dat'
and 'coord.dat' files that define a channel network and imports them into
ArcView as a Vector feature file. The vector feature file has a single
item for each Strahler stream (a sequence of stream segments with the same
Strahler order). An intermediate *.e00 file is generated and the
'import71' program run. Attribute information is included with these
files allowing stream properties to be identified. The program also
allows a subnetwork, or network with different support area threshold to
be imported through the specification of a link number (from the links
field in the attribute table) using option -l and alternative support area
threshold, using option -am.
Linkan
Program to compute link properties associated with channel network.
Input is 'coord.dat' and 'tree.dat' files Output is an ASCII file
'li.dat' listing properties associated with each link. These results
are used in the analysis to objectively determine a support area threshold
to use to map channel networks following the procedure of Tarboton et al.
(1991, 1992).
Streaman
Program to compute Strahler stream properties. Input is 'coord.dat'
and 'tree.dat' files. Output is an ASCII file 'st.dat' listing properties
associated with each Strahler stream. These results are also used in the
analysis to objectively determine a support area threshold to use to map
channel networks following the procedure of Tarboton et al. (1991, 1992).
Asfgrid
Program to tabulate the area and slope data for an area-slope analysis
to objectively determine a support area threshold to use to map channel
networks following the procedure of Tarboton et al. (1991, 1992).
Input are the files 'sd8' and 'ad8' if D8 is being used and 'sca' and 'slp'
if D¥ is being used. Output is two
ASCII text files suffixed by 'saraw.dat' containing each slope and
area value, and 'sa.dat' containing the data binned into a specified bin
size (default 100). The 'saraw.dat' file can be very large (15 MB
for an 800 x 700 DEM) so be sure you have sufficient disk space before
execution of this program.
Subbasinsetup
Program to define subbasins draining to each link (stream segment)
in a channel network controlled by the Strahler order threshold parameter
given. Input is the flow directions 'p' grid and 'tree.dat' and 'coord.dat'
files. Output is a grid (suffix 'w') that has a separate value for
each subbasin over which parameter values can be averaged for distributed
hydrologic modeling using these subbasins as model elements.
Streams of order lower than the threshold given are stripped off the network
before mapping resulting in a coarser watershed delineation.
The output also includes two vector shape files, suffix '.shp' a stream
network shape file and suffix 'w.shp' a shape file of subwatershed boundaries.
In these shape files the 'wsno' (watershed number) identifier in '.shp'
corresponds with the 'id' in 'w.shp' and is the same as the values in the
'w' grid allowing database linking and cross referencing.
Hack
Program to tabulate the contributing area and longest upstream path
at the downstream end of each link in a channel network. Empirically
the relationship between Area and mainstream length has been reported to
be L= A^H where H is Hack's exponent around 0.55 to 0.6. This program
allows this relationship to be tested. Input is the files defining
the channel network, with suffixes 'tree.dat' and 'coord.dat' as well as
grid file 'plen' that gives the longest flow path to each grid cell.
Output is a table of contributing area and length.
Pit Filling Corrections
Pits in digital elevation data are defined as grid elements or sets of grid elements surrounded by higher terrain that, in terms of the DEM, do not drain. These are rare in natural topography and generally assumed to be artifacts arising due to the discrete nature and data errors in the preparation of the DEM. They are eliminated here using a ‘flooding’ approach. This raises the elevation of each pit grid cell within the DEM to the elevation of the lowest pour point on the perimeter of the pit .
Slopes and Flow Directions
Working with grid DEMs slope may be computed as the difference in elevation between two adjacent cells divided by the distance between them. In dealing with flow this is usually done in a forward downwards direction. The slope associated with a cell is the slope from the cell to a downslope neighbor. This makes sense because it is where water will go. Radiation computations sometimes use slope based upon central finite difference methods. The earliest and simplest method for specifying flow directions is to assign flow from each grid cell to one of its eight neighbors, either adjacent or diagonally, in the direction with steepest downward slope. This method, designated D8 (8 flow directions), was introduced by O'Callaghan and Mark (1984) and has been widely used. The D8 approach has disadvantages arising from the discretization of flow into only one of eight possible directions, separated by 45° . These have motivated the development of other methods comprising multiple flow direction methods , random direction methods and grid flow tube methods . Tarboton (1997) discusses the relative merits of these.
In the D¥ method, the flow direction angle measured counter clockwise from east is represented as a continuous quantity between 0 and 2p. This angle is determined as the direction of the steepest downward slope on the eight triangular facets formed in a 3 x 3 grid cell window centered on the grid cell of interest as illustrated in figure 1. A block-centered representation is used with each elevation value taken to represent the elevation of the center of the corresponding grid cell. Eight planar triangular facets are formed between each grid cell and its eight neighbors. Each of these has a downslope vector which when drawn outwards from the center may be at an angle that lies within or outside the 45o (p/4 radian) angle range of the facet at the center point. If the slope vector angle is within the facet angle, it represents the steepest flow direction on that facet. If the slope vector angle is outside a facet, the steepest flow direction associated with that facet is taken along the steepest edge. The slope and flow direction associated with the grid cell is taken as the magnitude and direction of the steepest downslope vector from all eight facets. This is implemented using equations given in Tarboton (1997).
Figure 1. Flow direction defined as steepest downward slope on planar triangular facets on a block centered grid.
In the case where no slope vectors are positive (downslope), the flow direction is set using the method of Garbrecht and Martz (1997) for the determination of flow across flat areas. This makes flat areas drain away from high ground and towards low ground. The D¥method is preferred for the computation of flow directions on hillslopes where D8 grid bias is significant the the calculation of specific catchment area. D8 is still used for the definition of channel networks because we can not (have not yet learned to) work with channel networks that bifurcate in a downwards direction
Upslope area (counted in terms of the number of grid cells) is calculated for both single and multiple flow directions using a recursive procedure that is an extension of the very efficient recursive algorithm for single directions (Mark, 1988). The upslope area of each grid cell is taken as its own area (one) plus the area from upslope neighbors that have some fraction draining to it. The flow from each cell either all drains to one neighbor, if the angle falls along a cardinal (0, p/2, p, 3p/2) or diagonal (p/4, 3p/4, 5p/4, 7p/4) direction, or is on an angle falling between the direct angle to two adjacent neighbors. In the latter case the flow is proportioned between these two neighbor pixels according to how close the flow direction angle is to the direct angle to those pixels, as illustrated in Figure 1. Specific catchment area, a, is then upslope area per unit contour length, taken here as the number of cells times grid cell size (cell area divided by cell size). This assumes that grid cell size is the effective contour length, b, in the definition of specific catchment area and does not distinguish any difference in contour length dependent upon the flow direction.
The contributing area programs check for edge contamination. This is defined as 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 or areas with no data values for elevation. The algorithm recognizes this and reports no data for the contributing area. It is common to see streaks of no data values extending inwards from boundaries along flow paths that enter the domain at a boundary. This is the desired effect and indicates that contributing area for these grid cells is unknown due to it being dependent on terrain outside of the domain of data available. Edge contamination checking may be overridden by the option -nc in aread8 or areadinf in cases where you know this is not an issue or want to ignore these problems, if for example the DEM has been clipped along a watershed outline.
Channel Networks
When a map of contributing area is viewed using a threshold, the channel
networks stand out as those cells with contributing area greater than a
threshold of contributing area. It is an issue to decide the most appropriate
threshold, or whether some other quantity such as slope should be part
of the threshold. This is discussed at length in my research papers, Tarboton
et al. (1991, 1992). One approach that has some theoretical justification
is to look for a break in the plot of slope versus contributing area .
Once a threshold has been established the channel network can be defined
(mapped) as all those grid cells with contributing area greater than the
threshold.
ncols 480 nrows 450 xllcorner 378923 yllcorner 4072345 cellsize 30 nodata_value -32768 43 3 45 7 3 56 2 5 23 65 34 6 32 etc 35 45 65 34 2 6 78 4 38 44 89 3 2 7 etc etcThe first row of data is at the top of the data set, moving from left to right. Cell values should be delimited by spaces. No carriage returns are necessary at the end of each row in the data set. The number of columns in the header is used to determine when a new row begins. The number of cell values must be equal to the number of rows times the number of columns.
Grid naming convention.
The following default naming convention is suggested and used by the programs with command line input. Any file names may be used with interactive input, but I suggest sticking to this convention to avoid confusion. File names are:
nnnnsss[.asc]
nnnn comprises the name of the dataset. Maximum length is operating
system dependent.
sss comprises the suffix used to designate the data type as follows:
| 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 |
Vector Data
The following files are used to represent channel networks.
Network connectivity file, nnnntree.dat.
This is essentially a list of links comprising a channel network. It
is text with 7 columns as follows:
1 LINK NUMBERIn fortran programs this is generally stored in the array CNET. This file is produced by the program ‘netsetup’. The second and third columns refer to point coordinates, vectors along each link, from upstream to downstream, stored on the network coordinate, or ‘coord.dat’, file.
2 START POINT NUMBER IN COORD
3 END POINT NUMBER IN COORD
4 NEXT (DOWNSTREAM) LINK NUMBER IN CNET
5&6 PREVIOUS (UPSTREAM) LINK NUMBERS IN CNET
7 STRAHLER ORDER OF LINK IN CNET
Network coordinate file, nnnncoord.dat.
This is a list of coordinates defining the points along each channel
link. It is text with 5 coulmns of data as follows:
1 X COORDINATE (metres)In fortran programs this is generally stored in the array COORD. This file is produced by the program ‘netsetup’. The coordinates are based on the coordinate system (and projection) implicit in the header bounding box information in the raster grid file. Coordinates are the centers of grid elements (pixels) corresponding to each channel network link. This file is only useful in conjunction with the ‘tree’ file which gives the start and end position (line or record) in this file of each channel network link.
2 Y "
3 DISTANCE ALONG CHANNELS TO GAUGE (metres or whateve units grid size is in)
4 ELEVATION (metres or whateve units the DEM is in)
5 CONTRIBUTING AREA (meter2 or whatever units grid size is in)
Arc Export/Import ‘.e00’ file.
This is an ASCII file that is my interpretation of the format used
by Arc/Info and Arcview for vector export files. It is written by the programs
‘arclinks’ and ‘arcstreams’ and can be read into Arcview using the IMPORT71
utility. Then it is available as a ‘feature data source’ for viewing as
a ‘theme’. There is an ‘arc’ for each link or Strahler stream that has
various attributes.
Shape Files
This is an open ESRI data format that stores vector data in DBF files.
It is described in a white
paper. TARDEM reads EPA reach files in Shape file format to enforce
flow directions to follow existing streams where desired. TARDEM
also outputs the delineated channel network and reach subwatersheds in
Shape file format. The header information in these files is designed
to be self explanatory.
Windows 95/98
SET PATH=C:\WINDOWS;C:\WINDOWS\COMMAND;C:\DOS;C:\TARDEM\EXEC;C:\ESRI\AV_GIS30\ARCVIEW\BIN32
Windows NT/2000
path=c:\tardem\exec;c:\esri\av_gis30\arcview\bin32
There may be other 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. It is common in Win2000 for systems not to have an autoexec.bat file at all. If this is the case then you need to create one (c:\autoexec.bat) that has the single line giving path entry similar to above.
Test data.
The enclosed file test.asc includes a
small grid dataset that can be used to unsure you have things installed
right. The following sequence of commands illustrates the use of
the suite.
Link library sets required for each executable program.
| Name | Main | Library Sets or files required for linking |
| Aread8 | aread8.cpp | Tardem |
| Areadinf | areadinf.cpp | Tardem |
| Arclinks | arclinksmn.c | arcls.f, gridio |
| Arcstreams | arcstreamsmn.c | arcls.f, gridio |
| Asfgrid | aread8.cpp | asfgrid.f, gridio |
| D8 | d8.cpp | Tardem |
| Dinf | dinf.cpp | Tardem |
| Flood | floodmn.cpp | Tardem |
| Gridnet | gridnet.cpp | Tardem |
| Hack | hack.cpp | Tardem |
| Linkan | linkanmn.cpp | linkan.f, gridio |
| Netsetup | netsetupmn.cpp | Tardem |
| Streaman | streamanmn.cpp | streaman.f, gridio |
| Streamtogrid | streamtogrid.cpp | shapetogrid, gridiocpp |
| Subbasinsetup | Subbasinsetupmn.cpp | Tardem, wshedtopoly |
Link library set files
| Library set | Files |
| Tardem | tardemlib.cpp, nrutil.cpp, feedback_print.cpp, avcalls_null.cpp, gridio.cpp, avgridio.lib1,cell.cpp, dbf.cpp, exception.cpp, field.cpp, point.cpp, shape.cpp, shapefile.cpp, shapemain.cpp, shp_point.cpp, shp_line.cpp, shp_polygon.cpp, asfgrid.f, streaman.f, tardemflib.f |
| gridio | gridio.c, avgridio.lib1 |
| gridiocpp | gridiocpp, avgridio.lib1 |
| wshedtopoly | areatree.cpp, boundbox.cpp, Fileheader.cpp, linklist.cpp, myshp.cpp, point1.cpp, polygon.cpp, polygonmap.cpp, polyline.cpp, shpopen.cpp, shputils.cpp, touchpoint.cpp |
The tardemsource.zip file also
contains the necessary directory structure and Visual Studio workspace
to compile and link these programs on a Windows PC, with Microsoft Visual
Studio and Digital/Compaq Visual Fortran.
Changes in version 4. - released May 5, 2000.
Garbrecht, J. and L. W. Martz, (1997), "The Assignment of Drainage Direction Over Flat Surfaces in Raster Digital Elevation Models," Journal of Hydrology, 193: 204-213.
Jenson, S. K. and J. O. Domingue, (1988), "Extracting Topographic Structure from Digital Elevation Data for Geographic Information System Analysis," Photogrammetric Engineering and Remote Sensing, 54(11): 1593-1600.
Mark, D. M., (1988), "Network models in geomorphology," Chapter 4 in Modelling in Geomorphological Systems, Edited by M. G. Anderson, John Wiley., p.73-97.
Marks, D., J. Dozier and J. Frew, (1984), "Automated Basin Delineation From Digital Elevation Data," Geo. Processing, 2: 299-311.
O'Callaghan, J. F. and D. M. Mark, (1984), "The Extraction of Drainage Networks From Digital Elevation Data," Computer Vision, Graphics and Image Processing, 28: 328-344.
Peckham, S. D., (1995), "Self-Similarity in the Three-Dimensional Geometry and Dynamics of Large River Basins," PhD Thesis, Program in Geophysics, University of Colorado.
Peuker, T. K. and D. H. Douglas, (1975), "Detection of surface-specific points by local parallel processing of discrete terrain elevation data," Comput. Graphics Image Process., 4: 375-387.
Tarboton, D. G., (1989), "The analysis of river basins and channel networks using digital terrain data," Sc.D. Thesis, M.I.T., Cambridge, MA, (Also available as Tarboton D. G., R. L. Bras and I. Rodriguez-Iturbe, (Same title), Technical report no 326, Ralph M. Parsons Laboratory for Water resources and Hydrodynamics, Department of Civil Engineering, M.I.T., September 1989).
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.
Tarboton, D. G., R. L. Bras and I. Rodriguez-Iturbe, (1992), "A Physical Basis for Drainage Density," Geomorphology, 5(1/2): 59-76.