TARDEM, A suite of programs for the Analysis of Digital Elevation Data

David G. Tarboton
June 9, 1999

Acknowledgements

These have been developed during the course of my research over the years with support from a variety of sponsors, whose support is gratefully acknowledged. Specific sponsors include:
1. Massachusetts Institute of Technology, research assistantship under Rafael Bras, for my Sc.D. research where this all got started. Some remnants of the code from this work still remain.
2. National Science Foundation grant EAR-9318977 for the development of the D¥ approach (Tarboton, D. G., 1997).
3. Forest Renewal of British Columbia, for the development of Terrain Stability Mapping methodology and Arcview Implementation, in a collaborative project involving Canadian Forest Products Ltd., Vancouver, British Columbia, Terratech Consulting Ltd., British Columbia (Bob Pack), and Craig Goodwin.
4. National Science Foundation grant INT-9724720 and NIWA New Zealand for the work on methods for mapping and identification of flow methods from digital elevation data.

These programs are distributed freely with the following conditions.
1. In any publication arising from the use for research purposes the source of the program should be properly acknowledged and a pre-print of the publication sent to David Tarboton at the address below.
2. In any use for commercial purposes or inclusion in any commercial software appropriate acknowledgment is to be included and a free copy of this software made available to David Tarboton.


The programs and what they do

These are task oriented command line programs that need to be run from a DOS prompt on a PC, or shell on UNIX.  The program can be run with all input on the command line, or if no input is given the programs will prompt for input or let you know the format of command line input required.  Where grid file input is required the programs work with BOTH the ESRI binary grid format and ASCII export grid format.  If you specify a file name in one or the other the code will detect which it is and open it appropriately.  (The programs try the ESRI binary format first and if that fails assume the ASCII format.)  TARDEM has no display capability.  For this ArcView (or some other GIS or display package) needs to be used.  ESRI binary grids can be added as theme's to a view once each program has completed.  ASCII grids need to be imported then added as themes.

Flood
This takes as input an elevation data grid and outputs a grid file ‘fel’ with pits filled, using the flooding algorithm.

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).

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.

Aread8
This takes as input a D8 flow directions file ‘p’ and outputs the contributing area. Optional command line arguments for the outlet coordinate result in only the area contributing to the designated outlet being calculated.  The result is the number of grid cells draining through each grid cell.

Areadinf
This takes as input a Dinf angle file ‘ang’ and outputs the specific catchment area. Optional command line arguments for the outlet coordinate result in only the area contributing to the designated outlet being calculated.  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.

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:

file_name_prefix -m  4 .4 .1 .05 20 -xy outlet_x Outlet_y
A 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 5 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
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].

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.  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.


Background on working with Digital Elevation Models

The data storage structures available to digitally encode topography comprise: (1) Grid Digital Elevation Models (DEMs); (2) Triangular irregular networks (TINs); and (3) contour based storage structures. Grid DEMs consist of a matrix data structure with the topographic elevation of each pixel stored in a matrix node. TINs store the X-Y location as well as elevation at irregularly spaced nodes. Contour based data structures store vector data along contour lines. Grid DEMs are readily available and simple to use and hence have seen widespread application to the analysis of hydrologic problems. Slope, flow directions, and contributing area are the primary hydrologic quantities derived from DEMs. Other useful quantities are derived from these three. Here grid DEM are used due to their availability and simplicity. The grid DEM processing routines used are based upon methods described by O'Callaghan and Mark (1984), Marks et al. (1984), Band (1986), Jenson and Domingue (1988), Tarboton (1989, 1997) and Garbrecht and Martz (1997). The steps involved are: (1) Pit filling corrections, (2) Computation of slopes and flow directions; (3) Computation of contributing area and specific catchment area and (4) Channel network extraction and computation of other quantities.

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

Contributing Area

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.

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.


Data formats

Grid data
The programs are written to directly access the ESRI binary grid format accessible through the gridio application programmers interface that is part of Spatial Analyst (version 1.0a or higher) with ArcView (version 3.0a or higher).  The programs can also access ASCII grid data files in the format used by ESRI for export of files from ArcView and Arc/Info. This comprises a few lines of header data followed by lists of cell values. The header data includes the following keywords and values: For example,
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
etc
The 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.
p.    D8 drainage directions.
sd8.  D8 slopes.
ad8.  D8 contributing area’s.
slp.  Dinf slopes.
ang.  Dinf flow directions.
sca.  Dinf contributing area’s.
plen. Longest path length to each grid point along D8 directions.
tlen. Total path length to each grid point along D8 directions.
gord. Strahler order for grid network defined from D8 flow directions.
src.  Network mask based on channel source rules.
ord.  Grid with Strahler order for mapped stream network.
w.    Subbasins mapped using subbasinsetup.
The .asc extension is used if the data is ASCII. Otherwise it is assumed to be ESRI’s proprietary grid format and accessed using the gridio library supplied with ArcView. Command line C programs are aware of this suffix convention so only need be given a command line argument of the form nnn[.asc]. The suffix sss is appended automatically.  DO NOT USE SUFFIXES ON PROGRAM COMMAND LINES.

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 NUMBER
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
In 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.

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)
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)
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.

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.


Installation and Setup

Compiled Windows 95/98/NT executable programs are distributed in a zip file.  These should be extracted and placed in any convenient directory, e.g. c:\tardem\exec.  To effectively 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 easily met by including the following line in the autoexec.bat file.

Windows 95/98
SET PATH=C:\WINDOWS;C:\WINDOWS\COMMAND;C:\DOS;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 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.

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.

  1. Start ArcView.  Load Spatial Analyst extension.
  2. Import the data source (file 'test.asc') as an ASCII grid.  Name the grid 'test'.
  3. Open a DOS prompt and navigate to the folder where you saved the imported ESRI grid.
  4. flood test
  5. d8 test
  6. dinf test
  7. aread8 test
  8. areadinf test
  9. gridnet test
  10. View the results by loading some of the layers created, particularly 'testsca', 'testad8' and 'testgord'.
  11. Use the ArcView identify button to select the outlet of the subbasin to work with.  With the test data a good one is x: 2439530  y: 5857710
  12. netsetup test -m 4 .4 .1 .05 20 -xy 2439530  5857710
  13. arclinks test -i
  14. Add the feature theme testli to the View within ArcView to see the channel network mapped.
  15. subbasinsetup test 1
  16. Add the grid theme testw to the View within ArcView to see the subwatersheds mapped.
  17. Experiment with some of the other programs and look at datafiles.

Windows 95/98/NT executable programs, tardem.zip.


Source Code, tardemsource.zip.

These programs have all been compiled and tested under Microsoft Windows using the Microsoft Visual C++ (version 6.0) compiler and Digital Visual Fortran (version 6.0).  Many (though not the latest) have also been compiled and tested under UNIX.  The following table gives the components necessary to link together for each program, in case you wish to change anything or implement on a different system.
Name Main Files to Link
Flood floodmn.c flood.c
avcalls_null.c
callocate.c1
gridio.c
lsmcom.c
avgridio.lib2
D8 d8.c flood.c
avcalls_null.c
callocate.c1
gridio.c
lsmcom.c
setdir.c
avgridio.lib2
Dinf dinf.c flood.c
avcalls_null.c
callocate.c1
gridio.c
lsmcom.c
setdir.c
avgridio.lib2
Aread8 aread8mn.c aread8.c
avcalls_null.c
callocate.c1
gridio.c
lsmcom.c
avgridio.lib2
Areadinf areadinf.c area.c
avcalls_null.c
callocate.c1
gridio.c
lsmcom.c
avgridio.lib2
Gridnet gridnet.c callocate.c1
gridio.c
lsmcom.c
avgridio.lib2
Netsetup Netsetup.c sourcedef.c
netex.f
netprop.f
callocate.c1
gridio.c
gridio.f
lsmcom.c
avgridio.lib2
Arclinks arclinksmn.c arclinks.f
callocate.c1
gridio.c
lsmcom.c
avgridio.lib2
Arcstreams arcstreamsmn.c arcstreams.f
callocate.c1
gridio.c
lsmcom.c
avgridio.lib2
Linkan linkanmn.c linkan.f
callocate.c1
gridio.c
lsmcom.c
avgridio.lib2
Streaman streamanmn.c streaman.f
callocate.c1
gridio.c
avgridio.lib2
Asfgrid asfgridmn.c asfgrid.f
callocate.c1
gridio.c
avgridio.lib2
Subbasinsetup subbasinsetup.c linksource.f
nrutil.f
callocate.c1
gridio.c
avgridio.lib2
All Include files avexec32.h1
gioapi.h1
gridio.h
lsm.h
Notes.
1. Supplied as part of Spatial Analyst Gridio package.
2. In the directory c:\esri\av_gis30\arcview\gridio. To compile programs for use in standalone mode without ArcView substitute the file gridio_null.c for this library. This disables the capability to directly access ESRI binary grid format files, but use with ASCII grid files is possible.

To link FORTRAN and C together on a PC with the compilers mentioned I needed to use the following settings.

Settings/Fortran/External Procedures
Argument Passing Conventions:  C by reference
External Name Interpretation:  Lowercase
String Length Argument Passing:  After all arguments
Settings/Fortran/Libraries
Use C debug library (if C is using C debug library)
Settings/Link/General
Include 'dfor.lib' at the beginning of the Object Library Modules searched to avoid a clash between matherr defined in LIBCD.lib.

References

Band, L. E., (1986), "Topographic partition of watersheds with digital elevation models," Water Resources Research, 22(1): l5-24.

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.

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.

Tarboton, D. G. and U. Shankar, (1998), "The Identification and Mapping of Flow Networks from Digital Elevation Data," Invited Presentation at AGU Fall Meeting, San Francisco, December 6 to 10.