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

This is an archive of an old version of TARDEM whose use is no longer recommended. For the latest version see http://hydrology.usu.edu/taudem

David G. Tarboton                                                        May 5, 2000
Utah State University
http://www.engineering.usu.edu/dtarb/


Distribution and Copyright

Copyright (C) 2000  David Tarboton, Utah State University

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


Quick Downloads

This is an archive of an old version of TARDEM whose use is no longer recommended. For the latest version see http://hydrology.usu.edu/taudem

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 latest versions have not been compiled or tested for 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.  The standalone version only works with ASCII format grids and is intended for use when you do not have ArcView.

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


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.

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.


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) (This capability is only available with the Full Version and requires you to have ArcView).  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 capability is available in both the Full and Standalone versions.  The standalone version is intended for use when you do not have ArcView). The ASCII grid data file format 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.  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. 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.

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.


Installation and Setup

Compiled Windows 95/98/NT/2000 executable programs are distributed in the following zip files.  You should choose the full version if you have ArcView, the Standalone version, otherwise: 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. (With the standalone version you do not need to include an ESRI directory on the 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/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.

  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. Add the shape themes test.shp and testw.shp to the View within ArcView to see the subwatersheds mapped in vector form.
  18. Experiment with some of the other programs and look at datafiles.
The illustrative data given above is from the Grey River watershed on the west coast of the South Island of New Zealand, based on contours supplied by Land Information New Zealand.

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.  For the latest release most of  the code was switched to C++ from C and grouped into fewer files, with most of the substantial work in the file tardemlib.cpp.  This was done to facilitate linkage into an OCX control for use with a Visual Basic program TMDL toolkit.  There are functions in tardemlib.cpp that access Numerical Recipe, Fortran DLL and Shape library routines.  The result is that (somewhat counter to my goal) compilation and linking is more complex, against requiring library functions that are not actually used for each program.  I have grouped the files needed for compilation into library sets, and indicated which library set is required for each executable.  The following table gives the library sets required for each program, followed by a table giving the files in each library set.

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
1.  Supplies as part of Spatial Analyst GRIDIO Package.  For the standalone use, the file gridio_null.cpp may be substituted.  The programs will then not work with ESRI grid files.

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.


History and Old Versions

The current version is version 4

Changes in version 4. - released May 5, 2000.

  1. Code restructures into C++.
  2. Addition of capability to follow existing streams input as shape files.
  3. Addition of capability to output shape files.
  4. Addition of options to aread8 and areadinf to use weight files and override edge contamination checks.
  5. Addition of hack program to output area and length data to determine the exponent in Hack's law.
  6. Distributed under GNU General public license.
Changes in version 3. (released June 9, 1999)
  1. Consolidation of SOURCEDEF, NETEX, NETPROP into a single program NETSETUP.
  2. Addition of the programs GRIDNET, and SUBBASINSETUP.
  3. Revision of the interfaces of all programs to ( hopefully) make easier to use and better documented.

Version 2 is still available (at your peril).

Changes in Version 2
  1. Use of ArcView/ArcInfo grid file format.  This works with ASCII files and direct access of Arc binary grid files using the ESRI application programmers interface.
  2. Edge contamination checking in contributing area computations.
  3. More channel network extraction programs.  Channel networks can be extracted using a number of algorithms, such as Peuker Douglas curvature criterion, Slope and Area function, in addition to the constant support area threshold that was used earlier.  Extracted channel networks can be saved in an Arc export format for ready import to ArcInfo or Arcview.

Version 1.

The actual code used to generate results in Tarboton (1997).


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.
5. Idaho National Engineering and Environmental Laboratory for work on the adaptation of these codes for use with the TMDL Toolkit, and integration of flow with existing channel networks.

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.