Automated river-course, ridge and basin delineation from digital elevation data

Journal of Hydrology, 86 (1986) 299-314

299

Elsevier Science Publishers B.V., Amsterdam -- Printed in The Netherlands

[1] A U T O M A T E D RIVER-COURSE, RIDGE AND BASIN D E L I N E A T I O N F R O M D I G I T A L E L E V A T I O N DATA

O.L. PALACIOS-VI~LEZand B. CUEVAS-RENAUD Centro de Hidrociencias, Colegio de Postgraduados, Chapingo, M~x. 56230 (Mexico)

(Received August 28, 1985; accepted after revision March 26, 1986)

ABSTRACT Palacios-V~lez, O.L. and Cuevas-Renaud, B., 1986. Automated river-course, ridge and basin delineation from digital elevation data. J. Hydrol., 86: 299-314. Starting from a Digital Elevation Model (DEM) formed by a triangular irregular network (TIN) of facets, an algorithm that allows for the tracking of the path of steepest slope, downstream or upstream, from any point on the topographic surface has been developed. The analysis shows that when starting such a path in a downstream (upstream) direction from point P, down (up) to a point M, and then from point M in an upstream (downstream) direction and ifthe path returns to point P, then the segment P M corresponds to a "pure" streamline. Otherwise it corresponds to a river course (ridge).This property allows an easy way to delineate the network of river courses and ridges. Likewise, the repetitive use of this algorithm allows for the definition of basin and sub-basins, also called "contributing areas", to any point on a river course. Finally, itis shown h o w a regular grid D E M can be treated as a particular case of the TIN-DEM, which gives generality to the presented procedures.

INTRODUCTION In different kinds of studies, such as hydrological, erosion and surface drainage it is important to be able to delineate river-course and ridge networks, contributing areas to any point on a river course and paths of water, starting from any point on the topographic surface. The manual handling of these problems is highly laborious, time-consuming and prone to many errors. This is the main reason t h a t has motivated the development and use of automated procedures based on the analysis of Digital Elevation Models (DEM). The purpose of this paper is to present one such procedure, but before going into detail, some topics about the DEM's, their generation and the automated channel ridge and basin delineation will be briefly reviewed.

0022-1694/86/$03.50

© 1986 Elsevier Science Publishers B.V.

300 DIGITAL ELEVATION MODELS A DEM is defined as a computer data structure which represents the shape of some topographic surface. The most common form of DEM is the square regular grid, also known as the "altitude matrix". Sample points are located at the intersections of two orthogonal sets of regularly spaced parallel lines. Only the altitude of the surface at each sample point is measured and stored in the computer. The geographical locations are determined by the grid spacing and are implicit in the sequential position of the altitude value within the computer storage array. An advantage of this type of DEM is that the neighbors of a given point are also implicit in the position of the altitude value. Several considerations have motivated the development of an alternate concept of DEM. They are discussed by Peucker et al. (1976) and Mark (1975). The main disadvantage of the regular grid is its tendency toward redundancy in the ~'smoother" sub-areas of the study area. This redundancy behaves as "noise" that makes more difficult the extraction of river course and ridge information. A more flexible approach is the Triangulated Irregular Network (TIN). The surface is modelled as a set of contiguous non-overlapping triangular facets whose vertices are located adaptively on the "surface-specific" locations, that is, locations where the slope changes abruptly (Peucker et al., 1978; Tarvydas, 1984). Within each triangular facet, planar interpolation is used to estimate the altitude. According to some estimates (Mark, 1975; Peucker et al., 1976) from 14 up to more than 250 points of a regular grid are needed instead of only one point in the TIN approach, to achieve the same accuracy when measuring some geomorphological features. It should be noted, however, that in the TIN system each recorded point requires more storage space, because the TIN has to store the x and y coordinates of the point, as well as the z coordinate. It is also necessary to predefine and store the neighbors of each point, an operation that takes some computer time, but for which a series of et~icient computer programs has already been developed (Fowler, 1977; Palacios-V~lez, 1986). There are at least two ways to define the structure of a TIN model. The first is to store the labels of all points (the coordinates are stored separately) which are linked with each other point, in a counterclockwise (or clockwise) direction. Since the average number of neighbors per point cannot exceed six, the expected number of pointers is six times the number of points. The second approach uses a storage by triangular method (Tarvydas, 1984), which requires three pointers to the component nodes and three pointers to adjacent triangles, for a total of six pointers per triangle. This relationship is illustrated in Fig. 1. Since the number of triangles is about twice the number of points, the total number of pointers is twelve times the number of points or twice the storage requirements of the first approach. This form of data organization is, however,

301

No. of Triongle

NO. of Vertices

!

I

2

312

42

3

6 4

4

452

No. of Adj. TrJongles 423 156

~

I

? 8

~0 |

9

9

Fig. 1. The TIN model and the structural matrix.

somewhat easier to prepare and is more efficient for "triangle-by-triangle" processing and will be used in the following.

GENERATION OF TIN MODELS

There are two principal phases in the generation of TIN models: the selection of the data points and their connection into triangular facets. These two operations may occur either separately or in a unified process. In addition, they may be performed either manually or they may be automated. (Peucker et al., 1978). Sometimes the altitude and x, y coordinate data are the product of a transit and stadia survey in which the rodman selects ground points at which there is a change in slope (Collins, 1975). At other times these data are extracted manually from contour maps. The use of a digitizer may facilitate the determination of the x, y coordinates. In light of the increasing availability of machine readable terrain data, particularly in the form of very dense raster models, produced by automated orthophoto machines, it is becoming desirable to use automated point selection in conjunction with automated triangulation (Peucker et al., 1978). At least one system of this type has already been developed (Fowler and Little, 1979).

AUTOMATED DEM's

CHANNEL, RIDGE AND

BASIN D E L I N E A T I O N F R O M

REGULAR

GRID

Peucker and Douglas (1975) outlined automatic methods extracting ridge and channel information from a regular grid D E M . In one method, for each point with an (I,J) location in the square grid, the heights of the three cells at (I,J + 1),(I + 1,J) and (I + 1,J + 1) are examined (see Fig. 2).The lowest and

302

highest cells of the four are marked as not being ridge and channel candidates, respectively. When the scanning is completed a separate procedure is used to connect these points into ridges and channels. The following is a description of the ridge finding process, a similar procedure with the operations inverted is used to construct channels. The connection process starts at pass points and climbs to the neighboring ridge candidate at the highest altitude. Passes are found by examining all candidates and accepting as "passes" all those which are lower than all neighboring candidates. From these starting points the neighboring candidate at the highest altitude is added to the growing ridge, continuing until the ridge terminates at a peak or joins an existing ridge. To yield a simpler characterisation of the ridge and channel lines a generalization procedure is used to drop those points that are approximately located over a spacial straight line. Carroll (1983) presents a different procedure based on the analysis of height variation in four directions: two diagonals, as well as a vertical and a horizontal, so that for each point of the regular grid DEM its eight neighbors are examined. When the central height is a minimum (maximum), the location and orientation of the gully (ridge) are defined. Collins (1975) outlined a method by which all drainage basins in a DEM could be identified. His solution requires that all grid points first be sorted into ascending order of elevation. The lowest point is the outflow of a drainage basin. If the second lowest point is not a neighbor of the lowest point, then it

\

J+l

" "\\ Z

\

/

A.,c,

0Y

/KA* t

I

I+!

Fig. 2. Subdivision of each cell of the regular grid into four t r i a n g u l a r facets.

303 must be a separate drainage basin. By connecting neighbors in this sequence, all points on the grid are assigned to drainage basins. Marks et al.,(1984) present a recursive "order N" solution to this problem. In this method the eight neighbors of each point or "cell" of the regular grid D E M , that has been previously defined as part of the basin, are examined. A neighbor is considered "upstream" and therefore "in" if its slope is "nearly" horizontal or if the direction of its slope, or "exposure" that can have one of eight discrete values, passes through the cell.The algorithm starts at the basin outlet and travels recursively through the grid until it encounters either a basin boundary or the edge of the grid. It then backtracks repeating the procedure until the entire basin has been defined. Palacios-V~lez and Sailer (1985) define the contributing area to any point on a river course as the locus of cells whose paths of steepest slope pass through the given point on the river course. In order to calculate the paths of steepest slope the regular grid D E M is further subdivided into a network of triangular facets by adding a virtual point in the center of each cell,with an elevation equal to the average of the elevations at the corners and joining it with the corners, as is usually done when defining contour lines (see Fig. 2). As far as the writers k n o w all these problems have been solved only for a regular grid D E M , although Heil and Brych (1978) remark that "... with a properly constructed (TIN) terrain model, it is a relatively easy process to define the exit route of a drop of water fallingon the area at any location...", but give no details about the process or its possible applications for rivercourse, ridge and basin delineation.

DEFINITION OF PATHS OF STEEPEST SLOPE As will be seen later on, the automatic delineation of channels, ridges and the contributing area to any point on a channel, including the definition of whole watersheds, can be based on the analysis of paths of steepest slope, that theoretically would follow a drop of water falling on any location of the topographic surface. It will also be shown that solving these problems for the more general case of TIN models, the solution can be easily extended to the regular grid models, that are treated as a special case of TIN models. As mentioned above, for processing TIN models the coordinates of its points as well as their structural matrix should be provided. The matrix has 2 N + 2 rows, one for each triangular facet (where N is the number of points) and six columns, three for the number of vertices and three for the number of neighboring facets (see Fig. 1). In regular grid D E M ' s the determining the direction of m a x i m u m slope is a straight operation, since each cellis divided into two (by means of a diagonal) or four (by introducing a virtual point into the middle of the cell and then joining it with the vertices of the cell) right-angled triangles. In this case it su~ces to determine the slopes along the legs of the righ~angled triangle and then sum them up vectorially.

3o4 In T I N - D E M ' s w h e r e the t r i a n g l e s are n o t in g e n e r a l right-angled, a n o t h e r a p p r o a c h is proposed. Suppose we w a n t to join t h e point P ( u , v, w ) with the s t r a i g h t line O - M defined by: y

=

ax + b

(1)

z

=

cx + d

(2)

w h e r e t h e p a r a m e t e r s a, b, c and d are c a l c u l a t e d from the c o o r d i n a t e s of the two vertices t h a t form a n edge of a given t r i a n g u l a r facet, as it is i l l u s t r a t e d in Fig. 3. T h e problem is to define a p o i n t like M ( x , y , z), o v e r t h e line given by eqns. (1) and (2), in such a way t h a t the g r a d i e n t is m a x i m u m or minimum. T h e g r a d i e n t of line P - M is defined as: g

z - w D

=

(3)

where D

= x/(x-

u) 2 + ( y -

v) 2

By s u b s t i t u t i n g eqns. (1) and (2) in (3), c a l c u l a t i n g t h e d e r i v a t i v e of g with respect to x and e q u a t i n g it to zero, one obtains a f t e r some simplification: x* =

~ + ( y ° - v)2 x o -- x 1

xlx0

\ \\~,,

\3

(Z - W ~ \\\\

2

/

Fig. 3. Calculation of gradient between point P and line OM.

(4)

305

where: u + a(v Xo

=

1 +

-

b)

a2

(5)

Yo = axo + b

By the way, it is noted that point Q(xo, yo), on the line defined by eqn. (1) in the X - Y plan, is the closest to W(u, v). xl =

w-d c

(7)

y* and z* are calculated by eqns. (1) and (2). The sign of the slope can be calculated by substituting x, y and z by x*, y* and z* in eqn. (3). It should be noted that when x0 = xl the direction of optimum slope is parallel to eqn. (1) and therefore another edge of the triangular facet should be examined. When x* falls outside the examined edge, it should take the value of the abscissa of the vertex closest to x*. In this case point (x*, y*, z*) coincides with a vertex. When a is very large (even if it is only greater than 1) a similar set of formulas is used, but in which y, not x, is the variable. Based on this formulation, the path of steepest slope, in the downstream or upstream direction, can be determined as follows: (1) Start from a point given by its coordinates (u, v, w) and the number of the triangular facet. (2) Select an edge and calculate the (x*, y*, z*) point. If the direction of the optimum slope is parallel, or the sign of the gradient is not the one desired, select the next edge. (3) If there is no slope or no other facet, stop. (4) If the new point coincides with one vertex, go to step 6. (5) Update the coordinates and number of triangular facet and repeat step 2. In this step the use of the structural matrix is indispensable. (6) Analyse the opposite edges of all triangular facets which have the given vertex in common. Select the (x*, y*, z*) point with the biggest or smallest gradient. (7) If there is no slope, stop. (8) If the new point coincides with another vertex repeat step 6. (9) Update the coordinates and number of triangular facet where the new point is located and repeat step 2.

RIVER-COURSE A N D R I D G E D E L I N E A T I O N F O R TIN-DEM's

Let us examine a segment W M of a channel in a downstream direction (see Fig. 4). Considering that channels are lines of convergent slopes, it is evident

306

..........

.......

...........................

Fig. 4. Definition of river-course and ridge segments.

that for each point on that segment, for instance for point P, the lowest of all surrounding points at a arbitrary small distance, such as M, R, W and Q will be precisely on that segment and will be point M, otherwise it would not be a channel through which the water flows. Observe, however, that the highest of all surrounding points at an arbitrary small distance, will never be on the segment, but somewhere outside it (could be Q or R), which explains the convergence of the slope in a transverse direction A completely similar analysis can be made with a segment of a ridge in an upstream direction. In other words, following the path of steepest slope from any point P in a downstream (upstream) direction, down (up) to a point M, and then from this point in the upstream (downstream) direction, if the path returns to the initial point P, then the segment PM belongs to a "pure" streamline, with no convergence or divergence of slope in the transverse direction (segments of "pure" streamlines in Fig. 4 are QP and RP). Otherwise the segment ~ appertains to a channel (ridge). This property provides an easy, direct and general way of defining the channel (ridge) network, by analysing the paths of steepest slope in a downstream (upstream) direction. In a topographic surface made up of planar triangular facets, the edge between two adjoining facets, such as QPM and R I M in Fig. 5, will be a segment of a channel if, when analysing facet QPM, the point of minimum slope in relation to vertex P results outside the facet and below the other facet, as point T. It will be also necessary that when analysing facet RPM the point of minimum slope in relation to P will be outside RPM and below QPM, as point S. A similar analysis is made when searching the ridge segments.

307

I

t

i/ S

I

,\ T ~

I

Fig. 5. The common edge between two facets as a channel segment.

An algorithm to obtain the segments of the channel (ridge) network would consider the following steps: (1) For each vertex of the DEM, analyse in pairs the adjoining facets that have vertex in common. (2) For each pair of facets, calculate the point of optimum slope downstream (upstream) from the vertex. If for both cases this point falls outside the analysed facet and below (above) the other facet, the common edge will be a segment of the channel (ridge) network. A separate procedure can then join all of these segments to form the channel or ridge network.

DEFINITION OF THE CONTRIBUTING AREA TO ANY POINT ON A CHANNEL

The contributing area to a given point in a channel or river course, is the locus of all points in the watershed, whose paths of steepest slope pass through the given point. In the simplest case the contributing area is bounded by: (1) a segment of a "pure" streamline that climbs from the given point up a ridge; (2) a ridge line, sometimes alternating with segments of "pure" streamlines; and (3) a second segment of a "pure streamline" that descends from the ridge line and arrives at the given point in the channel. If there are points and closed ridge lines within the assumed contributed area, they should be substracted.

308 Clearly in order to apply this approach for the definition of contributing areas, a good routine, capable of tracking the ridge lines in ascending as well as descending direction should be developed first. However, the existence of segments of "pure" streamlines alternating with ridge segments would further complicate the procedure. A much simpler approach follows the given definition of a contributing area. The main steps to be considered, are: (1) Subdivide the whole watershed or studied area in "small" unit areas, for instance small squares. (2) Follow the path of steepest slope from the center of each unit area. (3) If the path passes through the given point on the channel, the unit area is assumed to be part of the contributed area. Otherwise it is not. The efficiency of this simple procedure can be significantly improved by: (4) Discarding from the beginning all unit areas whose maximum altitude is inferior to the given point on the channel. (5) Suspending the tracking procedure and discarding the unit areas whose paths descend to an elevation equal or less than the altitude of the given point, without reaching it. (6) Marking all vertices by which the path passes the first time, with - 1, if the path is finally discarded and with + 1 if the path reaches the given point. (7) Suspending the tracking procedure when a marked vertex is reached, and counting the unit area as part of the contributing area if the mark is + 1. Observe that the existence of internal pits and ridges is considered automatically.

REGULAR GRID M O D E L S AS A SPECIAL CASE OF TIN M O D E L S Usually different algorithms that lead to different computer codes are used w h e n dealing with T I N and regular grid DEM's. In our case, however, all of the FORTRAN codes, developed to solve the problems under discussion can be used for the T I N as well as for regular grid models. Since this generalization m a y be of independent interest, the changes in the computer programs that allow to treat a regular grid as a special case of T I N models are presented in the following. (1) Instead of three coordinate arrays and the structural matrix, only an altitude matrix ZA(I,J); I = 1, 2 .... NC; J - 1, 2 .... N R (where N C is the number of columns and N R is the number of rows) is required. Observe that while the number of points in the altitude matrix is N P = N C * N R , the number of cells is only N L = (NC - 1)*(NR - 1). There are, however, N L additional points in the center of each cell. Therefore the total number of points is N T P = N P + NL. The label of each point is K, K = 1, 2 .... NTP, the relationships between K and the I-column and J-row are as follows: K = (J - I)*NC + 1 or, using FORTRAN integer operations:

(8)

309 J

=

(K - 1)]NC + 1

(9)

I --- K - (J - 1)*NC

(10)

and:

(2) T h e c o o r d i n a t e s of t h e K-point are c a l c u l a t e d as follows: (2.1) if K = NP: X(K)

=

(I-

Y(K)

=

(J-

Z(K) =

1)*DX

(11)

1)*BY

(12)

ZA(I,J)

(13)

w h e r e DX a n d DY are t h e spacing b e t w e e n c o l u m n s and rows in t h e a l t i t u d e matrix. (2.2) if K > N P (the p o i n t is in t h e c e n t e r of t h e cell): (a) it s h o u l d be n o t e d t h a t I K = K - N P is t h e n u m b e r of t h e cell; a n d (b) t h e n u m b e r of t h e lower left v e r t e x of t h e cell, is: KA

=

[(IK - 1)*NC]/(NC - 1) + 1

(14)

T h e r e f o r e t h e c o r r e s p o n d i n g row and c o l u m n are: JA

=

(KA - 1)/NC + 1

(15)

IA

=

K A - (JA - 1)*NC

(16)

T h e c o o r d i n a t e s are: X(K)

=

DX]2. + (IA - 1)*DX

(17)

Y(K)

=

DY/2. + (JA - 1)*DY

(18)

Z(K)

-- [ZA(IA,JA) + ZA(IA + 1, JA) + ZA(IA + 1,JA + 1) + ZA(IA,JA + 1)]]4

(19)

(3) I n s t e a d of t h e s t r u c t u r a l m a t r i x MTRI(IT,L); IT = 1, 2 . . . . NF; L = 1, 2 . . . . 6, a c o r r e s p o n d i n g f u n c t i o n s h o u l d be definedl Observe t h a t t h e t o t a l n u m b e r o f facets is f o u r times t h e n u m b e r of cells: N F = 4* NL, since e a c h cell is subdiVided into four triangles. G i v e n t h e n u m b e r of t r i a n g u l a r facet IT, t h e n u m b e r o f cell will be: IK

=

(IT

-

1)]4 + 1

(20)

and t h e n u m b e r of t h e l o w e r left vertex, c o r r e s p o n d i n g row and c o l u m n m a y be c a l c u l a t e d by eqns. (14)--(16). W h e n t h e n u m b e r of t h e lower left vertex, KA, is k n o w n , t h e n t h e cell n u m b e r is: IK --

[KA*(NC - z)]INC + 1

(21)

T h e n u m b e r o f facets, up to t h e p r e v i o u s cell is: NFA(KA)

= 4 * [ K A t ( N C - 1)]/NC

-- 4*(IK - 1)

(22)

310 The labels of the four triangular facets of the IK-cell will be: N F A ( K A ) + 1, N F A ( K A ) + 2, N F A ( K A ) + 3 and N F A ( K A ) + 4, as shown in Fig. 2. The number of the central point is K C = N P ÷ IK, the vertices and numbers of neighboring facets are shown in Fig. 2. These values should be returned w h e n the function is called: the number of vertices, w h e n L = I, 2, 3 and the number of neighboring facets w h e n L = 4, 5, 6. While these are the only indispensable changes other changes are convenient, because they improve the efficiency of the computer program. For instance, the location of a point inside the network of triangular facets is a straight operation for the regular grid model, while a m u c h more complicated search procedure must be used in the TIN model. In this case the method of "triangular coordinates" (Gold et al., 1977) is used.

APPLICATIONS A 651 point TIN model, as well as a regular grid 62 × 21 altitude matrix were manually extracted from a contour map of the Walnut Gulch Experimental Watershed, Arizona. For the regular grid an average of about one minute was necessary to interpolate the height of each point. For the TIN model it took about five minutes to locate one point, interpolate its height and determinate its coordinates. Figures 6 and 7 show the "three-dimensional" block diagrams for the TIN as well as the regular grid model.

Fig. 6. 3-D block diagram of the Walnut Gulch Experimental Watershed TIN-DEM.

311

Fig. 7. 3-D block diagram of the Walnut Gulch Experimental Watershed regular grid DEM.

Figures 8 and 9 show the river-course networks for both the TIN and the regular grid models. From these figures it draws the attention to the lack of continuity of some segments. Actually this lack of continuity, much greater in the regular grid model, is due to an insufficient number of data points in the models, that not always permits the definition of convergent transversal slope, indispensable for the formation and identification of channels. In these somewhat artificial situations, the water flows, however, from the end of one segment of a channel to the beginning of the next separated one through a "pure" streamline, as defined previously. There are also some apparent loops. They are not truly loops because they include saddle points in which the water flows in opposite

Fig. 8. River-course network and contributing area to a point (TIN-DEM).

312

Fig. 9. River-course network and contributing area to a point (regular grid DEM).

direction. All these circumstances, unexpected at the beginning of the research, somehow reduce the utility of the river-course network. A more clear and complete idea of the water flow may be obtained by drawing lines of steepest slope from different locations, as shown in Fig. 10. Figures 8 and 9 also show the contributing areas or sub-basins to arbitrary points on the main river course. Once again, when a channel apparently crosses the boundary of the contributing area, what really happens is that two different river courses are formed at the ridge. It should be noted that the less the unit areas are, the greater will be the accuracy or '~resolution" of the procedure, but the computer time will also increase. There may be two kind of errors in this process: (1) it is possible that portions of a unit area, defined as part of the sub-basin, actually fall outside the sub-basin; and (2) it is also possible that parts of a unit area, defined as not being part of the sub-basin actually are part of it, as could be demonstrated, simply by defining a smaller unit area and repeating the whole procedure. During a test with the same TIN model the whole watershed was subdivided in unit areas of different size, in order to study the influence of this size in the calculation of the contributing area and the required computer time (for an HP-1000 with RTE-VIB operating system). The results are presented in Table 1. From this table it can be seen t h a t when the contributing area was defined only by a few tens of relative "big" unit areas, a fairly good estimate

Fig. 10. Lines of steepest slope t h a t define the water flow.

313 TABLE

1

Changes in the calculation of contributing area as a function of the unit area size Contributing area

Number of unit areas

C P U (sec)

Surface (ha)

No. of full paths

Total

Per unit area

6.24 12.50 20.31 42.36 72.10 158.86 288.29

O.1058 0.0735 0.0582 0.0545 0.0501 0.0483 0.0491

No. ofpa~s 59 170 349 777 1439 3289 5878

44.11 43.98 45.29 44.82 44.80 44.90 44.83

O.1515 0.1579 0.1332 0.1410 0.1380 0.1374 0.1359

of the sub-basin area was already obtained. Obviously, the delineation of the boundary is much better when the unit areas are smaller. The table also shows t h a t on the average only about 1/7 of the paths of steepest slope are fully analysed up to the boundary of the topography or up to the outlet, while 6/7 of the paths are only partially analysed up to a marked vertex, which shows the efficiency of the algorithm.

CONCLU~ONS A new procedure for defining the path of steepest slope, in downstream or upstream direction, for TIN DEM's has been developed. The analysis shows t h a t when starting such path in a downstream (upstream) direction from point P, down (up) to a point M, and then from this point in the upstream (downstream) direction, if the path returns to point P the segment PM belongs to a "pure" streamline with no convergence or divergence of transversal slope. Otherwise the segment appertains to a channel (ridge). This property provides an easy and general way to define the channel and ridge networks. Starting at the definition of "contributing area" to any point on a river course, as the locus of points whose paths of steepest slope pass through the given point on the river course, an algorithm to delineate these areas is presented. Significant ways to improve the efficiency of the algorithm are also discussed. The changes necessary to expand the computer programs for regular grid models, t h a t may be treated as special case of TIN models, are presented. All these developments are illustrated with examples for the Walnut Gulch Experimental Watershed regular grid and TIN models.

314 ACKNOWLEDGEMENT

The writers are gratefulto Dr. D.A. Woolhiser, from the A R S / S W Rangeland Watershed Research Center, Tucson, Ariz., for kindly providing the topography for the examples. They are also grateful to Mr. Juan Avila, computer operator at the Centro de Hidrociencias, for his aid in testing the computer programs.

REFERENCES Carrol, R., 1983. Automated gully delineation using digital elevation data, A C S M - A S P 49th Annu. Meet., Washington, D.C., Tech. Pap., pp. 144-151. Collins, S.H., 1975. Terrain parameters directly from D T M , Can. Surv., 29 (5):507-518. Fowler, R.J., 1977, DELTRI: An efficientprogram for producing Delaunay triangulations. Tech. Rep. 18, ONR, Contrast No. N00014-75-c.0886, Dept. Geogr., Simon Frazer Univ., Burnaby, B.C., Canada. Fowler, R.J. and Little, J.J., 1979. Automatic extraction of irregular network digital terrain models. Comput. Gr., 13 (2): 199-207. Gold, C.M., Charters, T.D. and RRmsdel, J., 1977. Automated contour mapping using triangular element data structures and an interpolant over each triangular domain, Comput. Gr., (11): 170-175. Heil, R.J. and Brych, S.M., 1978. An approach for consistent topographic representation of varying terrain. Proc. A S P - A C S M Syrup. on DTM's, St. Louis, Mo., pp. 397-411. Mark, D.M., 1975. Computer analysis of topography: A comparison of terrain storage methods. Geogr. Ann., 3-4 Ser. A.: 179-188. Marks, D., Dozier, J. and Frew, J. 1984. Automated basin delineation from digital elevation data, Geo-processing, 2: 299-311. Palacios-V~lez, O., 1986. A program for producing Thiessen Polygons and Delaunay Triangulation. In prep. Palacios-V~lez, O. and Sailer R., J.-S.,1985. An~lisis topogr~dico computarizado del escurrimiento superficial. Remitido a Agrociencia, Chapingo, M~x. Peucker, T.K. and Douglas, D.H., 1975. Detection of surface specific points by local parallel processing of discrete terrain elevation data. Comput. Gr. Image Process., 4: 375-387. Peucker, T.K., Fowler, R.J., Little, J.J. and Mark, D.M., 1976. Digital representation of threedimensional surfaces by triangulated irregular networks (TIN), Tech. Rept. 10, ONR, Contract No. N00014-75-C-0886, Dept. Geogr., Simon Frazer Univ., Burnaby, B.C., Canada. Peucker, T.K., Fowler, R.J., Little, J.J. and Mark, D.M. 1978. The Triangulated Irregular Network. Proc. ASP-ACSM Syrup. on DTM's, St. Louis, Mo. Tarvydas, A., 1984. Terrain Approximation by Triangular facets, ASP-ACSM 44th Annu. Meet., Washington, D.C., Teeh. Pap. pp. 524-532.