Terrain Analysis Using Digital

Elevation Models (TauDEM)

David Tarboton1, Dan Watson2,
Rob Wallace3
1UtahWater Research Laboratory, Utah State
University, Logan, Utah
2Computer Science, Utah State University,

Logan, Utah
3US Army Engineer Research and Development

Center, Information Technology Lab, Vicksburg,

Mississippi

http://hydrology.usu.edu/taudem dtarb@usu.edu

This research was funded by the US Army Research and

Development Center under contract number W9124Z-08-P-0420
Research

Theme
To advance the capability for hydrologic prediction
by developing models that take advantage of new
information and process understanding enabled by
new technology
Topics
Hydrologic information systems (includes GIS)
Terrain analysis using digital elevation
models(watershed delineation)
Snow hydrology
Hydrologic modeling
Streamflow trends
Hydrology and Ecology
2
 My Teaching
Physical Hydrology CEE6400 (Fall)

GIS in Water Resources CEE6440 (Fall)

A multi-university course presented on-line with shared
video in partnership with David Maidment at the
University of Texas at Austin.

Engineering Hydrology CEE3430 (Spring)

Rainfall Runoff Processes (online module)

http://www.engineering.usu.edu/dtarb/rrp.html
 Deriving hydrologically useful information from
Digital Elevation Models

Raw DEM Pit Removal (Filling)

Channels, Watersheds, Flow

Flow Field Related Terrain Information

Watersheds are the most basic

hydrologic landscape elements
flowing to downslope
neighboring
TauDEM - Channel Network and Watershed
grid cell 4 is
Proportion flowing to
neighboring grid cell 3
1/(1+2) is 2/(1+2)
Delineation
3
Software Pit removal (standard flooding approach)
4
2 Flow directions and slope
2 1 D8 (standard)
D (Tarboton, 1997, WRR 33(2):309)
1
Flat routing (Garbrecht and Martz,
5 1997, JOH 193:204)
Drainage area (D8 and D)
Network and watershed delineation
6 8 Support area threshold/channel
7
maintenance coefficient (Standard)
Flow direction measured as
counter-clockwise angle
Combined area-slope threshold
from east. (Montgomery and Dietrich, 1992,
Science, 255:826)
Local curvature based (using Peuker
and Douglas, 1975, Comput.
Graphics Image Proc. 4:375)
Threshold/drainage density selection by
stream drop analysis (Tarboton et al.,
1991, Hyd. Proc. 5(1):81)
Other Functions: Downslope Influence,
Upslope Dependence, Wetness index,
distance to streams, Transport limited
accumulation

Developed as C++ command line executable functions

MPICH2 used for parallelization (single program multiple data)
Relies on other software for visualization (ArcGIS Toolbox GUI)
The challenge of increasing Digital Elevation
Model (DEM) resolution
1980s DMA 90 m
102 cells/km2

1990s USGS DEM 30 m

103 cells/km2

2000s NED 10-30 m

104 cells/km2

2010s LIDAR ~1 m
106 cells/km2
 Website and Demo
http://hydrology.usu.edu/taudem
Model Builder Model to Delineate Watershed
using TauDEM tools
The starting point: A Grid
Digital Elevation Model
720 720
Contours
740

720

700

680

740 720 700 680

Grid Data Format Assumptions
Input and output grids are
uncompressed GeoTIFF
Maximum size 4 GB
GDAL Nodata tag preferred (if not
present, a missing value is
assumed)
Grids are square (x= y)
Grids have identical in extent, cell
size and spatial reference
Spatial reference information is not
used (no projection on the fly)
 The Pit Removal Problem

DEM creation results in artificial pits in the

landscape
A pit is a set of one or more cells which has
no downstream cells around it
Unless these pits are removed they become
sinks and isolate portions of the watershed
Pit removal is first thing done with a DEM
 Pit Filling

Increase elevation to the pour

point elevation until the pit
drains to a neighbor
 Pit Filling

Original DEM Pits Filled

7 7 6 7 7 7 7 5 7 7 7 7 6 7 7 7 7 5 7 7
9 9 8 9 9 9 9 7 9 9 9 9 8 9 9 9 9 7 9 9
11 11 10 11 11 11 11 9 11 11 11 11 10 11 11 11 11 9 11 11
12 12 8 12 12 12 12 10 12 12 12 12 10 12 12 12 12 10 12 12
13 12 7 12 13 13 13 11 13 13 13 12 10 12 13 13 13 11 13 13
14 7 6 11 14 14 14 12 14 14 14 10 10 11 14 14 14 12 14 14
15 7 7 8 9 15 15 13 15 15 15 10 10 10 10 15 15 13 15 15
15 8 8 8 7 16 16 14 16 16 15 10 10 10 10 16 16 14 16 16
15 11 11 11 11 17 17 6 17 17 15 11 11 11 11 17 17 14 17 17
15 15 15 15 15 18 18 15 18 18 15 15 15 15 15 18 18 15 18 18

Pits Pour Points

 Some Algorithm Details
Pit Removal: Planchon Fill Algorithm
Initialization 1st Pass 2nd Pass

Planchon, O., and F. Darboux (2001), A fast, simple and versatile algorithm to fill
the depressions of digital elevation models, Catena(46), 159-176.
 Parallel Approach
MPI, distributed
memory paradigm
Row oriented slices
Each process
includes one buffer
row on either side
Each process does
not change buffer row
 Parallel Scheme

Communicate
Initialize( Z,F)
Do
for all grid cells i
if Z(i) > n Z denotes the original elevation.
F(i) Z(i) F denotes the pit filled elevation.
Else n denotes lowest neighboring elevation
F(i) n i denotes the cell being evaluated
i on stack for next pass
endfor
Send( topRow, rank-1 ) Iterate only over stack of changeable cells
Send( bottomRow, rank+1 )
Recv( rowBelow, rank+1 )
Recv( rowAbove, rank-1 )
Until F is not modified
D8 Flow Direction Model
- Direction of steepest descent
30

80 74 63 4 3 2

69 67 56 5 1

60 52 48 6 7 8

67 52
0.50 Slope = Drop/Distance
30
67 48 Steepest down slope direction
0.45
30 2
Grid Network
Contributing Area (Flow Accumulation)
1 1 1 1 1 1 1 1 1 1

1 3 3 3 1 1 3 3 3 1

1 1 11 1 2 1 1 2
1 11
2 1 1 15 1 2 1 1
1 15
1 5 2 20 2 1 5 2 2
20

The area draining each grid cell includes the grid cell
itself.
 Stream Definition
Flow Accumulation Stream Network for
> 10 Cell Threshold 10 cell Threshold
Drainage Area
1 1 1 1 1 1 1 1 1 1
1 3 3 3 1 1 3 3 3 1

1 1 11 1 2 1 1 1 2
11
2 1 1 15 1 2 1 1
1 15
1 5 2 20 2 1 5 2 2
20
Watershed Draining to Outlet
Watershed and Stream Grids
 DEM Delineated Catchments and
Stream Networks
For every stream
segment, there is a
corresponding
catchment
Catchments are a
tessellation of the
landscape
Based on the D8
flow direction model
 Edge contamination

Edge contamination arises due to 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. The algorithm
recognizes this and reports "no data" resulting in streaks of "no data" values extending
inwards from boundaries along flow paths that enter the domain at a boundary.
 Representation of Flow Field
Proportion Steepest direction
Steepest flowing to downslope
neighboring Proportion flowing to
single grid cell 4 is neighboring grid cell 3
direction 1/(1+2) is 2/(1+2)
3 2
4 2 1
48 52 Flow
direction.

56 67 5 1

67 52 D8 D
0.50 6 8
30 7

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.)
Proportion Steepest direction
flowing to downslope D-Infinity Slope, Flow Direction
neighboring
grid cell 4 is
Proportion flowing to
neighboring grid cell 3
and Contributing Area
1/(1+2) is 2/(1+2)
3 2
4 2 1
Flow
direction.

5 1

6 8
7
 Pseudocode for Recursive Flow
Accumulation
Global P, w, A,
FlowAccumulation(i) Pki
for all k neighbors of i
if Pki>0
FlowAccumulation(k)
next k
Ai wi PkiA k
{k:Pki 0}
return
General Pseudocode Upstream Flow
Algebra Evaluation
Global P, ,
FlowAlgebra(i) Pki
for all k neighbors of i
if Pki>0
FlowAlgebra(k)
next k
i = FA(i, Pki, k, k)
return
 Example: Retention limited runoff
generation with run-on
Global P, (r,c), q
FlowAlgebra(i)
qk r
for all k neighbors of i
if Pki>0 c qi
FlowAlgebra(k)
next k
q i max( P q
{k:Pki 0}
ki k ri ci ,0)

return
 Retention Capacity
Retention limited
runoff with run-on
A C r=5
r=7 c=6
c=4 0.6 qin=1.8
q=3 q=0.8

1
r=4 0.4 Runoff from uniform input of 0.25
r=4
c=6 c=5
1
q=0 qin=2
q=1
B D

q i max( P q
{k:Pki 0}
ki k ri ci ,0)
Decaying Accumulation
A decayed accumulation operator DA[.]
takes as input a mass loading field m(x)
expressed at each grid location as m(i, j)
that is assumed to move with the flow field
but is subject to first order decay in
moving from cell to cell. The output is the
accumulated mass at each location DA(i,j).
The accumulation of m at each grid cell
can be numerically evaluated

DA[m(x)] = m(i, j)2

+ pkd(ik , jk )DA(ik , jk )
k contributing neighbors

Here d(i ,j) is a decay multiplier giving the

fractional (first order) reduction in mass in
moving from grid cell (i,j) to the next
downslope cell. If travel (or residence)
times t(i,j) associated with flow between
cells are available d(i,j) may be evaluated Useful for a tracking
as d(i, j) exp( t (i, j)) where is a first contaminant or compound
order decay parameter.
subject to decay or attenuation
 Transport limited accumulation
Supply Capacity Transport Deposition

S Tcap ca 2 tan(b ) 2 Tout min{S Tin , Tcap} D S Tin Tout

Useful for modeling erosion and sediment delivery, the spatial

dependence of sediment delivery ratio and contaminant that adheres to
sediment
 Parallelization of Contributing Area/Flow Algebra
1. Dependency grid
Queues
andempty
resulting
Decrease in onsountil
socross
new exchange
D=0 border
completion info.
cellsdependency
partition on queue
Executed by every process with grid flow field
P, grid dependencies D initialized to 0 and an
empty queue Q.
FindDependencies(P,Q,D) A=1
D=0 A=1
D=0 A=1
D=0
for all i
for all k neighbors of i
if Pki>0 D(i)=D(i)+1
if D(i)=0 add i to Q
next A=1.5
D=1
D=0 D=0
D=3
D=2
D=1
A=3 D=1
D=0
A=1.5
2. Flow algebra function
Executed by every process with D and Q B=-1
B=-2 B=-1
initialized from FindDependencies.
FlowAlgebra(P,Q,D,,) A=1
D=0 D=2
A=5.5
D=0 D=1
A=2.5
D=0
while Q isnt empty
get i from Q
i = FA(i, Pki, k, k)
for each downslope neighbor n of i
if Pin>0
D(n)=D(n)-1 D=0
A=1 A=6
D=3
D=2
D=1 A=3.5
D=1
if D(n)=0
add n to Q
next n
end while
swap process buffers and repeat
 Capabilities Summary
11 GB

Capability to run larger problems

Grid size 6 GB
Processors
Theoretcal Largest 4 GB
used
limit run
2008 TauDEM 4 1 0.22 GB 0.22 GB 1.6 GB

Partial
Sept
implement- 8 4 GB 1.6 GB
2009 0.22 GB
ation
June
TauDEM 5 8 4 GB 4 GB
2010
Multifile on
Sept Hardware
48 GB 4 6 GB
2010 limits
RAM PC
Multifile on
Sept cluster with Hardware
128 11 GB
2010 128 GB limits
RAM

Single file size limit 4GB

At 10 m grid cell size
 Improved runtime efficiency
Parallel Pit Remove timing for NEDB test dataset (14849 x 27174 cells 1.6 GB).

2000
ArcGIS T ~ n 0.03
500 1000

Total
Compute

Seconds
Seconds

500
T ~ n 0.44 Total
C ~ n 0.56

200
200

C ~ n 0.69 Compute

1 2 3 4 5 7 1 2 5 10 20 50
Processors Processors

8 processor PC 128 processor cluster

Dual quad-core Xeon E5405 2.0GHz PC with 16GB 16 diskless Dell SC1435 compute nodes, each with 2.0GHz dual
RAM quad-core AMD Opteron 2350 processors with 8GB RAM
 Improved runtime efficiency

Parallel D-Infinity Contributing Area Timing for Boise River dataset (24856 x 24000 cells ~ 2.4 GB)

Total

500
500

Compute T ~ n 0.18
T ~ n 0.63 to 48 proc. Total
Seconds

Seconds
100 200
Compute
C ~ n 0.95
200

C ~ n 0.93
100

to 48 proc.
50
1 2 3 4 5 7 10 20 50 100
Processors Processors
8 processor PC 128 processor cluster
Dual quad-core Xeon E5405 2.0GHz PC with 16 diskless Dell SC1435 compute nodes, each with 2.0GHz
16GB RAM dual quad-core AMD Opteron 2350 processors with 8GB RAM
 Scaling of run times to large grids

Number of PitRemove D8FlowDir

Dataset Size Hardware Processors (run time seconds) (run time seconds)
(GB) Compute Total Compute Total
GSL100 0.12 Owl (PC) 8 10 12 356 358
GSL100 0.12 Rex (Cluster) 8 28 360 1075 1323
GSL100 0.12 Rex (Cluster) 64 10 256 198 430
GSL100 0.12 Mac 8 20 20 803 806
YellowStone 2.14 Owl (PC) 8 529 681 4363 4571
YellowStone 2.14 Rex (Cluster) 64 140 3759 2855 11385
Boise River 4 Owl (PC) 8 4818 6225 10558 11599
Boise River 4 Virtual (PC) 4 1502 2120 10658 11191
Bear/Jordan/Weber 6 Virtual (PC) 4 4780 5695 36569 37098
Chesapeake 11.3 Rex (Cluster) 64 702 24045

1. Owl is an 8 core PC (Dual quad-core Xeon E5405 2.0GHz) with 16GB RAM
2. Rex is a 128 core cluster of 16 diskless Dell SC1435 compute nodes, each with 2.0GHz dual quad-core
AMD Opteron 2350 processors with 8GB RAM
3. Virtual is a virtual PC resourced with 48 GB RAM and 4 Intel Xeon E5450 3 GHz processors
4. Mac is an 8 core (Dual quad-core Intel Xeon E5620 2.26 GHz) with 16GB RAM
 Scaling of run times to large grids
100000 100000

PitRemove run times D8FlowDir run times

10000
Time (Seconds)

Time (Seconds)
10000
1000

100 Compute (OWL 8) Compute (OWL 8)

Total (OWL 8) 1000 Total (OWL 8)
Compute (VPC 4) Compute (VPC 4)
10 Total (VPC 4) Total (VPC 4)
Compute (Rex 64) Compute (Rex 64)
Total (Rex 64) Total (Rex 64)
1 100
0.1 1 10 0.1 1 10
Grid Size (GB) Grid Size (GB)
1. Owl is an 8 core PC (Dual quad-core Xeon E5405 2.0GHz) with 16GB RAM
2. Rex is a 128 core cluster of 16 diskless Dell SC1435 compute nodes, each with 2.0GHz dual quad-core AMD
Opteron 2350 processors with 8GB RAM
3. Virtual is a virtual PC resourced with 48 GB RAM and 4 Intel Xeon E5450 3 GHz processors
Programming
C++ Command Line Executables
that use MPICH2
ArcGIS Python Script Tools
Python validation code to provide
file name defaults
Shared as ArcGIS Toolbox
Q based block of code to evaluate any flow algebra
expression
while(!que.empty())
{
//Takes next node with no contributing neighbors
temp = que.front(); que.pop();
i = temp.x; j = temp.y;
// FLOW ALGEBRA EXPRESSION EVALUATION
if(flowData->isInPartition(i,j)){
float areares=0.; // initialize the result
for(k=1; k<=8; k++) { // For each neighbor
in = i+d1[k]; jn = j+d2[k];
flowData->getData(in,jn, angle);
p = prop(angle, (k+4)%8);
if(p>0.){
if(areadinf->isNodata(in,jn))con=true;
else{
areares=areares+p*areadinf->getData(in,jn,tempFloat);
}
}
}
}
// Local inputs
areares=areares+dx;
if(con && contcheck==1)
areadinf->setToNodata(i,j);
else
areadinf->setData(i,j,areares);
// END FLOW ALGEBRA EXPRESSION EVALUATION
}
 Maintaining to do Q and partition sharing

while(!finished) { //Loop within partition

while(!que.empty())
{ .... // FLOW ALGEBRA EXPRESSION EVALUATION
}
// Decrement neighbor dependence of downslope cell
flowData->getData(i, j, angle);
for(k=1; k<=8; k++) {
p = prop(angle, k);
if(p>0.0) {
in = i+d1[k]; jn = j+d2[k];
//Decrement the number of contributing neighbors in neighbor
neighbor->addToData(in,jn,(short)-1);
//Check if neighbor needs to be added to que
if(flowData->isInPartition(in,jn) && neighbor->getData(in, jn, tempShort) == 0 ){
temp.x=in; temp.y=jn;
que.push(temp);
}
}
}
}
//Pass information across partitions
areadinf->share();
neighbor->addBorders();
 Python Script to Call Command Line

mpiexec n 8 pitremove z Logan.tif fel Loganfel.tif

PitRemove
Validation code to add default file names
 Multi-File approach
To overcome 4 GB file size
limit
To avoid bottleneck of
parallel reads to network
files
What was a file input to
TauDEM is now a folder
input
All files in the folder tiled
together to form large
logical grid
 Processor Specific Multi-File Strategy
Input Output
Shared file Node 1 Node 1 Shared file
store local disk local disk store

Core 1

Core 2

Node 2 Node 2
local disk local disk

Core 1
Scatter all Gather partial
input files to output from
all nodes Core 2 each node to
form complete
output on
shared store
 Summary and Conclusions
Parallelization speeds up processing and partitioned processing
reduces size limitations
Parallel logic developed for general recursive flow accumulation
methodology (flow algebra)
Documented ArcGIS Toolbox Graphical User Interface
32 and 64 bit versions (but 32 bit version limited by inherent 32 bit
operating system memory limitations)
PC, Mac and Linux/Unix capability
Capability to process large grids efficiently increased from 0.22 GB
upper limit pre-project to where < 4GB grids can be processed in
the ArcGIS Toolbox version on a PC within a day and up to 11 GB
has been processed on a distributed cluster (a 50 fold size
increase)
Limitations and Dependencies
Uses MPICH2 library from Argonne National
Laboratory
http://www.mcs.anl.gov/research/projects/mpich2/

TIFF (GeoTIFF) 4 GB file size (for single file

version)

Run multifile version from command line for > 4

GB datasets

Processor memory