For Information that Supplements the Contents of this Unit:

USGS Digital Elevation Models (Delaney) -- Data accuracy, acquisition, applications and availability.

Environmental Modeling and Visualization with GRASS (US Military) -- Interacting fields; surface modeling; multidimensional modeling; 3D scattered data interpolation; terrain analysis; etc.

Terrain Analysis and Erosion Modeling (US Military) -- (Graphics) Elevation surface; slope angle; aspect angle; profile and tangential curvature; topographic potential for net erosion/deposition; etc.

B. ESTIMATING ELEVATION

C. ESTIMATING SLOPE AND ASPECT

D. DETERMINING DRAINAGE NETWORKS

surfaces, such as the surface of the earth, are continuous phenomena (fields) rather than discrete objects

to fully model the surface, would need an infinite amount of points

there are various ways of representing continuous surfaces in digital form using a finite amount of storage

this lecture will look at digital elevation models as one way of representing surfaces and will examine some important algorithms based on DEMs

the term digital elevation model or DEM is frequently used to refer to any digital representation of a topographic surface

however, most often it is used to refer specifically to a raster or regular grid of spot heights

this is the definition that is used here

digital terrain model or DTM may actually be a more generic term for any digital representation of a topographic surface, but it is not so widely used

the DEM is the simplest form of digital representation of topography and the most common

a variety of DEMs are available, including coverage of much of the US from the US Geological Survey

the resolution, or the distance between adjacent grid points, is a critical parameter

the best resolution commonly available is 30 m, with a vertical resolution of 1 m

coverages of the entire globe, including the ocean floor, can be obtained at various resolutions

at 30m resolution there are approximately 10¹² points required to cover the globe

several different methods have been used to create DEM series like those from the USGS

see USGS (1987) for more details on the following

conversion of printed contour lines

existing plates used for printing maps are scanned

the resulting raster is vectorized and edited

contours are "tagged" with elevations

additional elevation data are created from the hydrography layer

i.e. shorelines provide additional contours

finally, an algorithm is used to interpolate elevations at every grid point from the contour data

draw four lines through the point, interpolate profiles along each line, and average

by photogrammetry

this can be done manually or automatically (soft photogrammetry):

manually, an operator looks at a pair of stereophotos through a stereoplotter and must move two dots together until they appear to be one lying just at the surface of the ground

automatically, an instrument calculates the parallax displacement of a large number of points

e.g. for USGS 7.5 minute quadrangles, the Gestalt Photo Mapper II correlates 500,000 points

extraction of elevation from photographs is confused by flat areas, especially lakes, and wherever the ground surface is obscured (buildings, trees)

there are two techniques for choosing sample points when using manual photogrammetry:

1. profiling

the photo is scanned in rows, alternately left to right and right to left, to create profiles

a regular grid is formed by resampling the points created in this process

because the process tends to underestimate elevations on uphill parts of each profile and overestimate on downhill parts, the resulting DEMs show a characteristic "herringbone" or "firth" effect when contoured

2. contour following

contour lines are extracted directly from stereopairs during compilation of standard USGS maps

contour data are processed into profile lines and a regular grid is interpolated using the same algorithms used for manual profiling data

DEMs from each source display characteristic error artifacts

e.g. effects of mis-tagged contours in the products of scanned contour lines

determining attributes of terrain, such as elevation at any point, slope and aspect

finding features on the terrain, such as drainage basins and watersheds, drainage networks and channels, peaks and pits and other landforms

modeling of hydrologic functions, energy flux and forest fires

the principal components of a drainage basin are its topographic form and the topologic structure of its drainage network

the quantification of these components is tedious and time consuming when accomplished manually

the automated determination of these components is an ideal application of GIS technology

watersheds comprise one method of completely partitioning space and many environmental phenomena can be related to watersheds

determination of the drainage network and the associated drainage divides provides an important first step in the creation of a hydrologic information system

registration and segmentation of digital imagery can be enhanced if use is made of the drainage basin information

knowledge of the drainage divides and of the drainage network can be used to provide better estimates of slopes and aspects (e.g., slopes should break at divides and at channels)

in this unit we look at a number of simple algorithms for DEMs

to estimate the elevation of some point, we need to know first whether the point of interest is exactly at a point in the raster, or in between

in the first case, the elevation can be taken directly from the database

in the second case, we need to develop some method of interpolation, or estimation of elevation

can use the elevation of the nearest point, but this leads to sharp changes of elevation halfway between points

instead, the normal approach is to fit a plane to the nearby raster points, and use it to estimate elevation at any point

the plane passing through these points is represented as:

z = a + bx + cy

since a plane will generally not pass exactly through all the points

 
the plane which minimizes the sum of squared elevation differences between the plane and the data at each of the nearby points is often used

can determine the equation of the plane as follows:

use the four nearest grid points (known as the "neighborhood" of the point or the "2x2 window" around the point)

define an origin in the middle of the 2x2 window, and give the four neighboring points the coordinates (-1,-1), (-1,1), (1,-1) and (1,1)

since the four points are evenly spaced, the coefficients in the equation can be calculated from the following:
 

a = [z(-1,-1) + z(-1,1) + z(1,-1)+ z(1,1)]/4

b = [-z(-1,-1) - z(-1,1) + z(1,-1) + z(1,1)]/4

c = [-z(-1,-1) + z(-1,1) - z(1,-1) + z(1,1)]/4


note: the coefficients can be solved using larger neighborhoods, e.g. the 9 points in a 3x3 window

 
the points are at (-1,-1), (-1,0), (-1,1), (0,-1), (0,0), (0,1), (1,-1), (1,0), (1,1)

a = [z(-1,-1) + z(-1,0) + z(-1,1) + z(0,-1) + z(0,0) + z(0,1) + z(1,-1) + z(1,0) + z(1,1)]/9

b = [-z(-1,-1) - 2z(-1,0) - z(-1,1) + z(1,-1) + 2z(1,0) + z(1,1)]/8

c = [-z(-1,-1) - 2z(0,-1) - z(1,-1) + z(-1,1) + 2z(0,1) + z(1,1)]/8


having determined the coefficients, the elevation (z) can be determined from:
 

z = a + bx + cy


a major disadvantage of this approach is that the planes do not pass through all the data points

 
the planes simply come as close as possible to all 4 or all 9


bilinear interpolation avoids this problem

 
see diagram

slope and aspect can be calculated from the fitted plane

to estimate these at a raster point, a 3x3 window centered on the point is usually used

slope is calculated from:

(b² + c²)^1/2

aspect is calculated from:
 
tan^-1 c/b


normally a "slope map" or "aspect map" will display the attribute values generalized over areas (regions) instead of at points, such that within each area, all slopes fall into a certain range (e.g. 10-15%) or all aspects fall into a certain quadrant (e.g. NW)

 
to generate such a map, slope or aspect is determined at each raster point, and then these values are aggregated into polygons based on a set of predefined ranges

this way of representing slope or aspect is not as accurate as the original raster form


since both slope and aspect are derivable from elevation by a simple process, is there any need to store them as separate layers?

a raster DEM contains sufficient information to determine general patterns of drainage and watersheds

think of each raster point as the center of a square cell

the direction of flow of water out of this cell will be determined by the elevations of surrounding cells

algorithms to determine the flow direction generally use one of the following cases:

1. assume only 4 possible directions of flow (up, down, left, right - the Rook's move directions from chess); or

2. assume 8 possible directions (the Queen's move directions)

in both cases, number the move directions clockwise from up

water is assumed to flow from each cell to the lowest of its neighbors

 
if no neighbor is lower, the cell is a "pit" and gets code "0"


combinations which would be hydrologically impossible, such as a 4 to the left of a 6 in the 8 move case, are logically impossible in this scheme

Flow directions (4 moves)
 

3	3	3	3
3	3	3	3
2	3	3	4
2	0	4	4

Flow directions (8 moves)
 
 

4	4	5	6
4	4	5	6
4	5	6	6
3	0	7	7

Determining the watershed

a watershed is defined here as an attribute of each point on the network which identifies the region upstream of that point

to find a watershed

begin at the specified cell and label all cells which drain to it, then all which drain to those, etc. until the upstream limits of the basin are defined

the watershed is then the polygon formed by the labeled cells

to draw the drainage network, connect the moves with arrows

a zero on the edge of the array is interpreted as a channel which flows off the area

since in natural systems, small quantities of water generally flow overland, not in channels, we may want to accumulate water as it flows downstream through the cells so that channels begin only when a threshold volume is reached

accumulation of volume proceeds as follows:

start by setting each cell to zero

then beginning at each cell, add one to it and all cells downstream of it, following the directions indicated in the network

to simulate actual stream channels, assume a channel begins only when the accumulated water passing through a cell reaches some critical value

this means that small tributaries in the examples above will be deleted

in the example, channels start only when the flow reaches a volume of 2

Accumulated flow
 

1	1	1	1
1	2	4	1
1	2	8	1
1	16	3	1

the networks found by this process can be thought of as estimates of real channel networks

real networks consist of junctions or forks, links, and sources, and all of these can be identified on the simulated networks

how do networks obtained from DEMs differ from real ones?

real streams sometimes branch downstream

but this is impossible using this method, the simulated networks cannot bifurcate

DEM data contains large numbers of ties of elevation, because the vertical resolution is not very high

using this method, water cannot "flow" from one cell to an adjacent cell with the same elevation

as a result, ties can lead to large numbers of unwanted pits

e.g. in this example, using Rook's case (4 directions) central cell has no outflow direction

to avoid the problem, allow water to flow between neighbors at the same elevation, determining the direction of flow by evaluating local slope (i.e. over a larger window)

e.g. here the local slope is to the south

alternatively, deal with the problem by regarding the cell as a very small lake, and simulating its overflow (see next point)

pits occur frequently on DEMs, largely as a result of data errors

if a cell has no lower neighbors, it is a pit

the pit can be "flooded" to form a "lake" by the following process:

initiate a lake at the elevation of the cell, with a "shoreline" defined by the cell's perimeter

find the lowest cell adjacent to the shoreline, raise the lake to that level and expand the shoreline to include it

if one of the neighbors is now lower than the lake, it is the outlet: terminate the process

if the lowest neighbor is part of another lake, merge the lakes and continue

the number of streams joining at a junction, known as the valency of the junction, is almost always 3 in reality, but may be as high as 4 with the 4-move case, and 8 with the 8-move case

junction angles are determined by the cell geometry in the simulation, but in reality are a function of the terrain and the erosion process

in areas of uniform slope the technique generates large numbers of parallel streams

in reality streams tend to wander because of unevenness, and the resulting junctions reduce the density of streams in areas of approximately uniform slope

drainage density is very high in the simulations

in many types of terrain, channels are incised, and often the width of the incised channel is too small to show on the DEM

this problem can be dealt with by searching the DEM for possible channels - see Band (1986) for example

these methods do well on highly dissected landscapes of high drainage density

they do better at modeling watershed boundaries than drainage channels

therefore, ideally, a spatial database for modeling runoff and other processes related to hydrology should include both the DEM and the stream channels themselves (the "blue lines" of a topographic map)

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

Burrough, P.A., 1986. Principles of Geographical Information Systems for Land Resources Assessment, Clarendon, Oxford. Chapter 3 reviews alternative methods of terrain representation.

Evans, I.S., 1980. "An integrated system for terrain analysis and slope mapping," Zeitschrift fur Geomorphologie 36:274-95.

Marks, D., J. Dozier and J. Frew, 1984. "Automated basin delineation from digital elevation data," Geoprocessing 2:299-311.

O'Callaghan, J.F. and D.M. Mark 1984. "The extraction of drainage networks with lakes," Water Resources Research, 18(2):275-280.

Pfaltz, J.L., 1976. "Surface networks," Geographical Analysis 8:77-93. Discussion of surface-specific points and their relationship to ridge and channel lines.

USGS, 1987. Digital Elevation Models, Data Users Guide 5, US Department of the Interior, USGS, Reston, VA. Describes the creation and data structures of USGS DEMs in detail.

1. Discuss some of the problems encountered with algorithms which extract drainage networks from digital elevation models, and present some possible solutions to those problems.

2. How would the incorporation of hydrologic information-- such as drainage divides and stream networks--into a GIS assist a resource manager?

3. Discuss the problems of obtaining maps of slope and aspect from DEMs.

4. What possible ways are there for displaying a DEM on a computer screen? Discuss the advantages and disadvantages of each from the point of view of a) the users and b) the programmers.