Figure 1. Flow direction defined as steepest downward slope on planar triangular facets on a block centered grid.
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.
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_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 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.
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.
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. 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 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.
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.
| 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 referenceSettings/Fortran/Libraries
External Name Interpretation: Lowercase
String Length Argument Passing: After all argumentsUse C debug library (if C is using C debug library)Settings/Link/GeneralInclude 'dfor.lib' at the beginning of the Object Library Modules searched to avoid a clash between matherr defined in LIBCD.lib.
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.