- Email: [email protected]
Automated channel ordering and node indexing for raster channel networks
PII:
Computers & Geosciences Vol. 23, No. 9, pp. 961±966, 1997 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0098-3004/97 $17.00 + 0.00 S0098-3004(97)00055-1
AUTOMATED CHANNEL ORDERING AND NODE INDEXING FOR RASTER CHANNEL NETWORKS JURGEN GARBRECHT1 and LAWRENCE W. MARTZ2 USDA, Agricultural Research Service, Grazinglands Research Laboratory, 7207 W. Cheyenne St., El Reno, OK 73036, and 2Department of Geography, 9 Campus Drive, University of Saskatchewan, Saskatoon, Saskatchewan, Canada S7N 0W0 (e-mail: [email protected])
1
(Received 1 December 1996; revised 6 April 1997) AbstractÐA numerical algorithm is proposed for automated interpretation of channel networks from raster images, indexing of network nodes, and ordering of channels by the Strahler method. Channel ordering and node indexing is fundamental to the automation of ¯ow-routing management in distributed hydrologic models and morphometric evaluation of channel-network structure. The node index numbers can also serve to link network nodes to corresponding tabulated attributes of network channels. The proposed algorithm uses a two-step approach: ®rst, the raster image of the network is interpreted and a channel attribute table is created; second, the network nodes are indexed based only on the network connectivity information contained in the channel attribute table. The use of attribute information in the second step eliminates a repeat of the time consuming cell-by-cell interpretations of the network raster image. In addition to the general case of nodes with two upstream in¯ows, the algorithm also handles complex nodes (nodes with more than two upstream in¯ows) which are rare under natural conditions, but more frequent in raster networks. The ®nal product of the presented algorithm is a table of channel orders and node indices derived from raster images of channel networks. # 1998 Elsevier Science Ltd. All rights reserved Key Words: Drainage network, Network topology, Digital elevation model, Distributed modeling, Image interpretation.
INTRODUCTION
Regional soil and water conservation studies of forest, range and agricultural lands often involve the application of distributed-parameter hydrologic models. These models require watershed segmentation and channel network delineation to account for the eects of spatial variability (Beasley, Huggins, and Monke, 1980; Knisel, 1982; Young and others, 1987; Plate, Ihringer, and Lutz, 1988; Feldman, 1995; Smith and others, 1995; Singh, 1995). Watershed segmentation has traditionally been accomplished by manual evaluation of topographic maps (Ponce, Osmolski, and Smutzer, 1985; Goodrich and others, 1994), but is increasingly replaced by automated delineation of channel network and subcatchments using Geographic Information Systems (GIS) and raster Digital Elevation Models (DEM) [Band, 1986; Jenson and Domingue, 1988; Vieux, 1991]. Hydrologic models are then applied to the delineated subcatchments, and subcatchment output is fed into the channel network for modeling downstream propagation and watershed response. For distributed modeling to be organized and automated, the channel network nodes (channel junction points) must be indexed in a hydrologically 961
meaningful way that allows the determination of the sequence in which the channel ¯ow routing in the network is to be evaluated. The automation of this task requires an algorithm that can interpret a raster image of a network and index the network nodes. Such an algorithm provides a direct link between GIS images and hydrologic models, and leads to automated processing of segmented watersheds by distributed hydrologic models. The resulting indexed network nodes can also be applied in studies of channel network morphometry and network composition. The concepts and methods underlying the algorithm are extendable to other dendritic structures such as the branching pattern of trees or certain reticulate crystal formations. In this paper an algorithm to analyze images of raster channel networks, index network nodes, and order the channels by the Strahler method is presented. The computational steps are presented in detail to permit the implementation of the algorithm by the reader. The initial network de®nition and other related raster data that are needed by the algorithm are produced by the computer program TOPAZ, an automated digital landscape analysis tool for topographic evaluation, watershed segmentation, and channel identi®cation (Garbrecht and Martz, 1995). Even though the algorithm is
962
J. Garbrecht and L. W. Martz
bering system was selected because it allows the automated determination of ¯ow routing sequence in networks (Garbrecht, 1988). Referring to Figure 1, the watershed outlet node is assigned index 1 and a trace upstream from the outlet is initiated along the channels. At the next node encountered in this trace the index number is incremented by one and assigned to that node. This is repeated until all nodes have been assigned an index number. The channel trace follows the leftmost branch at all junction nodes. At source nodes the tracing direction is reversed and channels are followed downstream until a node with at least one untraced channel is encountered. The left-most channel is then traced upstream as before, until another source node is reached. This continues until the tracing returns to the watershed outlet. The resulting node indexes are shown in Figure 1.
Figure 1. Example channel network and de®nitions.
illustrated with initial data from TOPAZ, the numerical methods can be incorporated into other GIS systems that contain similar raster data. BACKGROUND
A channel network consists of a set of channel links connected by network nodes (Horton, 1945; Shreve, 1966; Smart, 1972), hereafter simply referred to as nodes (Fig. 1). Three types of nodes are encountered in a channel network: the watershed outlet node, upstream tips of the channel network where channel links originate (source nodes), and points at which two or more channel links join (junction nodes), see Figure 1. The channels can be ordered according to Strahler (1957). The automated Strahler ordering of the channels of a raster network is the ®rst objective of the proposed algorithm. In the Strahler ordering system the channel links that originate at source nodes are ®rst-order channels. When two ®rst-order channel links join, a second-order channel is formed; when two secondorder channel links join, a third-order channel is formed; and so forth. When a lower-order channel joins a higher-order channel, the order of the higher-order channel does not change. The trunk channel through which all discharge and sediment passes from upstream sources is therefore the channel link of highest order. For sake of simplicity, the name Strahler in conjunction with the word order is omitted in the remainder of this paper. Automated node indexing of a raster network is the second objective of the proposed algorithm. Node index numbers are assigned according to the node numbering system described by Croley (1980). The origins of this numbering system can be traced to Shreve (1967) and Smart (1978). Croley's num-
INPUT DATA REQUIREMENTS
A network and ¯ow-direction raster are needed as initial input for the proposed algorithm. The channel network is identi®ed in the network raster as strings of connected raster cells having a value of 100 on a background of cells having a value of 0. The cells at the tips of the network (source nodes) are given a value of 1. In the example network raster shown in Figure 2, the cells with an arrow represent channel cells having a value of 100. This network was generated by software TOPAZ
Figure 2. Channel network and corresponding ¯ow directions. Cells with bold outline represent network nodes; channel cells with arrow have value of 100.
Channel ordering and node indexing
963
(Garbrecht and Martz, 1995) using the Critical Source Area (CSA) concept which classi®es cells with upstream channel area greater than the CSA value as being part of the network (Jenson and Domingue, 1988; Morris and Heerdegen, 1988; Martz and Garbrecht, 1992). The ¯ow-direction raster identi®es, by code numbers, the ¯ow direction at each cell as the direction of the steepest downward slope to one of the eight adjacent cells (Fair®eld and Leymarie, 1991). The ¯ow directions of the channel cells are shown in Figure 2 by arrows.
ALGORITHM OVERVIEW
The algorithm consists of two computational steps. First, the raster image of the network is processed, a Strahler order is assigned to each channel, and a table of channel attributes is created. This table contains the coordinates of the beginning and ending node of each channel, and thereby captures the topology of the network. The numerical interpretation of the raster network begins at the source nodes and proceeds in a downstream direction along the channels to the watershed outlet node. This upstream-to-downstream processing direction is required because the order of downstream channels depends on the order of upstream channels. Second, node index numbers are assigned. The numerical processing begins at the watershed outlet node, and proceeds in an upstream direction along the channels to the source nodes. This downstreamto-upstream processing direction is required by the node indexing convention previously described. This could be accomplished by a cell-by-cell tracing of the raster network from the watershed outlet node to the source nodes. However, the cell-by-cell tracing is time consuming and is done in the ®rst step of the algorithm. A faster method that relies on the previously identi®ed channel topology data in the attribute table is implemented. Using the attribute data allows tracing of the network by jumping directly from node to node without tracing cell-by-cell the channel links. Step 1: De®nition of channel order De®nition of channel order begins by identifying all ®rst-order channels (source channels), then all second-order channels, and so forth. The algorithm explained in the following is for the ®rst two channel orders only. Higher-order channels are treated the same way as second-order channels. The process ¯ow chart for the algorithm is shown in Figure 3. First order channels The network raster is scanned to ®nd a source node (cell with value 1 in Fig. 2). On ®nding a source node, the row/column coordinates of this node are stored for later use, and the channel is traced downstream following the direction indicated by the ¯ow direction. All channel
Figure 3. Process ¯ow chart for de®nition of channel orders.
cells encountered during this trace, including the initial source node cell, are given a value of ÿ1 (negative value of the order). Negative values are used to identify those channel cells that have been processed. Every cell is also checked to determine if it is a junction node. When a junction node is encountered, the trace is terminated and a value of 2 is assigned to the junction node to de®ne it as a second-order channel cell. The row/column coordinates of the junction cell and of the previous to last cell before the junction are retained. The network raster scan then resumes at the previous source cell to identify the next source node and initiate another round of processing. Once all source nodes and ®rst-order channels have been processed, the network raster consists of background values of zero, ®rst-order channels of value ÿ1, junction cells at the downstream end of ®rst-order channels of value 2, and all remaining channel cells of value 100. This situation is shown in Figure 4A. Second order channels The network raster is scanned to detect second-order nodes (cells with value 2) at which time a similar processing as for ®rst-order channels is applied, but with three major dierences: (1) during the initial scan, cells with value 2 must be checked to determine if they are at the upstream end of a second-order channel before the trace in the downstream direction is initiated; (2) at a junction node at the end of a trace, the order of the next downstream channel may remain unchanged as indicated by a value of 2 at the junction node; under these conditions the downstream trace is re-initiated from this junction node to the next, and so forth, until a junction node with value 3 or 100 is encountered; and (3) a value of 3 is assigned to the junction node at the end of a full
964
J. Garbrecht and L. W. Martz
Figure 4. Raster network: (A) after processing ®rst-order channels, and (B) after processing secondorder channels.
trace of a second-order channel segment to indicate the beginning of a third-order channel. At the end of the evaluation of all second-order channels, the network is de®ned as shown in Figure 4B. The remaining nodes and cells are identi®ed in a similar fashion. In the previous example, the third-order channels would be indicated with a value of ÿ3. Final network raster All higher-order channels are treated in the same manner as second-order channels. When all orders have been evaluated, the entire network is de®ned by negative values. The absolute value is taken and the network is de®ned by cells having positive values on a background of cells having value 0. The value of the channel cells represents the order of the channel. In addition to this raster of the network, an attribute table was created containing (for each channel link) the order and coordinates of the beginning and ending points of each channel. This attribute table is used in the following node indexing algorithm. Step 2: Network node indexing Node indices are assigned by searching and evaluating the attribute table rather than by tracing the channels cell-by-cell in the raster network. The concept of this approach is presented ®rst, followed by the algorithm for the network search and node indexing. The process ¯ow chart is shown in Figure 5. The algorithm starts at the last data record of the channel-link attribute table which represents the watershed outlet node. The watershed outlet node is assigned index number 1. Then, the channel ending
coordinates of all channels in the attribute table are searched to identify those channels that terminate at the current node. From the identi®ed in¯ow channels, the appropriate left-hand channel is selected, the node at the upstream end of this channel is given the next successive index number, and the search of the channel attribute table for upstream in¯ow is re-initiated for this new node. When a source node is encountered the search moves in a downstream direction until a partially or fully un-evaluated node is encountered, at which time the node-by-node search and indexing resumes again in an upstream direction. This proceeds until all nodes are evaluated and the watershed outlet node is reached again. Each node processed during this search is treated dierently according to the nodes type: source node, simple junction node, complex junction node, or watershed outlet node. The decision logic and indexing at each node type is detailed in the following. The watershed outlet node is encountered twice during the node search. First, at the initiation of the indexing procedure when the index number 1 is assigned to the outlet. After indexing the watershed outlet node, the upstream node of the in¯owing channel becomes the current node. Second, at the end of the indexing procedure and completion of the algorithm. If the outlet node happens to be a junction node, then the node is treated as a junction node. Source nodes are encountered only once during the network search. When a source node is encountered, the next successive index number is assigned,
Channel ordering and node indexing
965
at the current time. Non-consecutive indices indicate that a trace along the left in¯ow channel is undertaken before the index for the next in¯ow at the complex junction is assigned. This dierentiation is important to ensure consistency in channel Strahler order when decomposing complex junctions into a sequence of simple junctions. Once the channel network nodes are indexed, the morphometric properties of the network can be evaluated and the node indices can be processed by software SEQDRV (Garbrecht, 1992) to optimize the sequence of channel link evaluation for cascadetype ¯ow routing.
CONCLUSIONS
Figure 5. Process ¯ow chart for network node indexing.
and the next downstream node is selected as the current node in the search. A simple junction node is encountered three times during a network search. Once when moving up the left-most in¯ow; once when returning from the left-most in¯ow and moving up the right in¯ow; and once when returning from the right in¯ow. The ®rst time a simple junction node is encountered, the next successive index number is assigned and the left-most in¯ow channel link is followed upstream to the next node. The second time a junction node is encountered, no index number is assigned but the second in¯owing channel link is followed upstream. The third time a simple junction node is encountered, no index number is assigned and the next downstream node becomes the current node. Complex junction nodes are the result of the discrete representation of the network. During the network search they are encountered as many times as the junction has in¯ow channels plus one. The indexing is similar to that of simple junctions, except that multiple node indices are assigned to account for the extra in¯ow channels. The additional indices are in anticipation of the decomposition of the complex junction node into a sequence of simple junction nodes as required by subsequent node processing software (Garbrecht, 1988) and for a more accurate de®nition of network composition parameters (Garbrecht and Martz, 1995). The indices assigned to a complex junction node may be consecutive or non-consecutive depending on the sequence and Strahler order of the in¯owing channel links. Consecutive indices indicate that lower order channel links enter from the left and no upstream trace is undertaken along these channels
An algorithm that automatically determines channel Strahler order and node indices for raster channel networks is presented. Strahler orders are determined by a cell-by-cell trace of the raster network in a downstream direction along the channels beginning at the source nodes to the watershed outlet. During this trace, beginning and ending coordinates of each channel link are retained in an attribute table. The subsequent node indexing uses the channel topology information in the attribute table to search the network from downstream to upstream starting at the watershed outlet node and following the left hand pattern described in Croley (1980) and Garbrecht (1988). During this search, nodes are indexed consecutively. The algorithm accounts for complex junction nodes that are particular to raster channel networks and assigns indices that allow a subsequent decomposition of the complex node into a sequence of simple junction nodes. The Strahler orders and the decomposition of complex junction nodes are necessary for a more accurate determination of network composition parameters and for subsequent optimization of the sequence in which channels must be evaluated for cascade type ¯ow routing. The algorithm presented herein makes possible the automated quanti®cation of network structure from raster network images and greatly enhances the direct linkage of GIS generated channel networks and hydrologic and hydraulic models. Code is available by anonymous FTP from FTP.IAMG.ORG.
REFERENCES Band, L. E. (1986) Topographic partitioning of watersheds with digital elevation models. Water Resources Research 22(1), 15±24. Beasley, D. B., Huggins, L. F. and Monke, E. J. (1980) ANSWERS: A model for watershed planning. Transactions of American Society of Agricultural Engineers 23(4), 938±944. Croley, T. E. (1980) A micro-hydrology computation ordering algorithm. Journal of Hydrology 48, 221±236.
966
J. Garbrecht and L. W. Martz
Fair®eld, J. and Leymarie, P. (1991) Drainage networks from grid digital elevation models. Water Resources Research 27(4), 29±61. Feldman, A. D. (1995) HEC-1 ¯ood routing package. In Computer Models of Watershed Hydrology, ed. V. J. Singh, pp. 119±150. Water Resources Publication, Highlands Ranch, Colorado. Garbrecht, J. (1988) Determination of the execution sequence of channel ¯ow for cascade routing in a drainage network. Hydrosoft 1(3), 129±138. Garbrecht, J. (1992) Documentation of model and computer program SEQDRV: topologic evaluation of drainage networks for cascade routing. HYDROWAR Reports Division, Hydrology Days Publications, 1005 Country Club Road, Fort Collins, CO 80524, Report No. 92.3, 81 pp. Garbrecht, J. and Martz, L. W. (1995) TOPAZ: An automated digital landscape analysis tool for topographic evaluation, drainage identi®cation, watershed segmentation, and subcatchment parameterization. Overview. ARS Publication NAWQL 95-1, U. S. Department of Agriculture, Agricultural Research Service, Oklahoma 74702, 16 pp. Goodrich, D. C., Schmugge, T. J., Jackson, T. J., Unkrich, C. L., Keefer, T. O., Parry, R., Bach, L. B. and Amer, S. A. (1994) Runo simulation sensitivity to remotely sensed initial soil water content. Water Resources Research 30(5), 1393±1405. Horton, R. E. (1945) Erosional development of streams and their drainage basins, hydrophysical approach to quantitative morphology. Bulletin of the Geological Society of America 56, 275±370. Jenson, S. K. and Domingue, J. O. (1988) Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering and Remote Sensing 54(11), 1593±1600. Knisel, W. G. (ed). (1982) CREAMS: A ®eld-scale model for chemical runo, and erosion form agricultural management systems. USDA-Science and Education Administration, Conservation Report, 26 pp.
Martz, L. W. and Garbrecht, J. (1992) Numerical de®nition of drainage networks and subcatchment areas from digital elevation models. Computers & Geosciences 18(6), 747±761. Morris, D. G. and Heerdegen, R. G. (1988) Automatically derived catchment boundary and channel networks and their hydrological applications. Geomorphology 1(2), 131±141. Plate, E. J., Ihringer, J. and Lutz, W. (1988) Operational models for ¯ood calculations. Journal of Hydrology 100, 489±506. Ponce, V. M., Osmolski, Z. and Smutzer, D. (1985) Large basin deterministic hydrology: a case study. Journal of Hydraulic Engineering 111(9), 1227±1245. Shreve, R. L. (1966) Statistical law of stream numbers. Journal of Geology 74, 17±37. Shreve, R. L. (1967) In®nite topologically random channel networks. Journal of Geology 75, 178±186. Singh, V. P. (ed.) (1995) Computer Models of Watershed Hydrology. Water Resources Publications, Highland Ranch, Colorado, 1130 pp. Smart, J. S. (1972) Channel networks. Advances in Hydrosciences 8, 305±346. Smart, J. S. (1978) The analysis of drainage network composition. Earth Surface Processes 3, 129±170. Smith, R. E., Goodrich, D. C., Woolhiser, D. A. and Unkrich, C. L. (1995) KINEROSÐA kinematic runo and erosion model. In Computer Models of Watershed Hydrology, ed. V. J. Singh, pp. 697±732. Water Resources Publication, Highland Ranch, Colorado. Strahler, A. N. (1957) Quantitative analysis of watershed geomorphology. Transactions of the American Geophysical Union 38(6), 913±920. Vieux, B. E. (1991) Geographic information systems and non-point source water quality and quantity modeling. Hydrological Processes 5(1), 101±113. Young, R. A., Onstad, C. A., Bosch, D. D. and Anderson, W. P. (1987) AgNPS: agricultural non-point-source pollution model: a watershed analysis tool. USDAARS, Conservation Research Report, 35 pp.
Computers & Geosciences Vol. 23, No. 9, pp. 961±966, 1997 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0098-3004/97 $17.00 + 0.00 S0098-3004(97)00055-1
AUTOMATED CHANNEL ORDERING AND NODE INDEXING FOR RASTER CHANNEL NETWORKS JURGEN GARBRECHT1 and LAWRENCE W. MARTZ2 USDA, Agricultural Research Service, Grazinglands Research Laboratory, 7207 W. Cheyenne St., El Reno, OK 73036, and 2Department of Geography, 9 Campus Drive, University of Saskatchewan, Saskatoon, Saskatchewan, Canada S7N 0W0 (e-mail: [email protected])
1
(Received 1 December 1996; revised 6 April 1997) AbstractÐA numerical algorithm is proposed for automated interpretation of channel networks from raster images, indexing of network nodes, and ordering of channels by the Strahler method. Channel ordering and node indexing is fundamental to the automation of ¯ow-routing management in distributed hydrologic models and morphometric evaluation of channel-network structure. The node index numbers can also serve to link network nodes to corresponding tabulated attributes of network channels. The proposed algorithm uses a two-step approach: ®rst, the raster image of the network is interpreted and a channel attribute table is created; second, the network nodes are indexed based only on the network connectivity information contained in the channel attribute table. The use of attribute information in the second step eliminates a repeat of the time consuming cell-by-cell interpretations of the network raster image. In addition to the general case of nodes with two upstream in¯ows, the algorithm also handles complex nodes (nodes with more than two upstream in¯ows) which are rare under natural conditions, but more frequent in raster networks. The ®nal product of the presented algorithm is a table of channel orders and node indices derived from raster images of channel networks. # 1998 Elsevier Science Ltd. All rights reserved Key Words: Drainage network, Network topology, Digital elevation model, Distributed modeling, Image interpretation.
INTRODUCTION
Regional soil and water conservation studies of forest, range and agricultural lands often involve the application of distributed-parameter hydrologic models. These models require watershed segmentation and channel network delineation to account for the eects of spatial variability (Beasley, Huggins, and Monke, 1980; Knisel, 1982; Young and others, 1987; Plate, Ihringer, and Lutz, 1988; Feldman, 1995; Smith and others, 1995; Singh, 1995). Watershed segmentation has traditionally been accomplished by manual evaluation of topographic maps (Ponce, Osmolski, and Smutzer, 1985; Goodrich and others, 1994), but is increasingly replaced by automated delineation of channel network and subcatchments using Geographic Information Systems (GIS) and raster Digital Elevation Models (DEM) [Band, 1986; Jenson and Domingue, 1988; Vieux, 1991]. Hydrologic models are then applied to the delineated subcatchments, and subcatchment output is fed into the channel network for modeling downstream propagation and watershed response. For distributed modeling to be organized and automated, the channel network nodes (channel junction points) must be indexed in a hydrologically 961
meaningful way that allows the determination of the sequence in which the channel ¯ow routing in the network is to be evaluated. The automation of this task requires an algorithm that can interpret a raster image of a network and index the network nodes. Such an algorithm provides a direct link between GIS images and hydrologic models, and leads to automated processing of segmented watersheds by distributed hydrologic models. The resulting indexed network nodes can also be applied in studies of channel network morphometry and network composition. The concepts and methods underlying the algorithm are extendable to other dendritic structures such as the branching pattern of trees or certain reticulate crystal formations. In this paper an algorithm to analyze images of raster channel networks, index network nodes, and order the channels by the Strahler method is presented. The computational steps are presented in detail to permit the implementation of the algorithm by the reader. The initial network de®nition and other related raster data that are needed by the algorithm are produced by the computer program TOPAZ, an automated digital landscape analysis tool for topographic evaluation, watershed segmentation, and channel identi®cation (Garbrecht and Martz, 1995). Even though the algorithm is
962
J. Garbrecht and L. W. Martz
bering system was selected because it allows the automated determination of ¯ow routing sequence in networks (Garbrecht, 1988). Referring to Figure 1, the watershed outlet node is assigned index 1 and a trace upstream from the outlet is initiated along the channels. At the next node encountered in this trace the index number is incremented by one and assigned to that node. This is repeated until all nodes have been assigned an index number. The channel trace follows the leftmost branch at all junction nodes. At source nodes the tracing direction is reversed and channels are followed downstream until a node with at least one untraced channel is encountered. The left-most channel is then traced upstream as before, until another source node is reached. This continues until the tracing returns to the watershed outlet. The resulting node indexes are shown in Figure 1.
Figure 1. Example channel network and de®nitions.
illustrated with initial data from TOPAZ, the numerical methods can be incorporated into other GIS systems that contain similar raster data. BACKGROUND
A channel network consists of a set of channel links connected by network nodes (Horton, 1945; Shreve, 1966; Smart, 1972), hereafter simply referred to as nodes (Fig. 1). Three types of nodes are encountered in a channel network: the watershed outlet node, upstream tips of the channel network where channel links originate (source nodes), and points at which two or more channel links join (junction nodes), see Figure 1. The channels can be ordered according to Strahler (1957). The automated Strahler ordering of the channels of a raster network is the ®rst objective of the proposed algorithm. In the Strahler ordering system the channel links that originate at source nodes are ®rst-order channels. When two ®rst-order channel links join, a second-order channel is formed; when two secondorder channel links join, a third-order channel is formed; and so forth. When a lower-order channel joins a higher-order channel, the order of the higher-order channel does not change. The trunk channel through which all discharge and sediment passes from upstream sources is therefore the channel link of highest order. For sake of simplicity, the name Strahler in conjunction with the word order is omitted in the remainder of this paper. Automated node indexing of a raster network is the second objective of the proposed algorithm. Node index numbers are assigned according to the node numbering system described by Croley (1980). The origins of this numbering system can be traced to Shreve (1967) and Smart (1978). Croley's num-
INPUT DATA REQUIREMENTS
A network and ¯ow-direction raster are needed as initial input for the proposed algorithm. The channel network is identi®ed in the network raster as strings of connected raster cells having a value of 100 on a background of cells having a value of 0. The cells at the tips of the network (source nodes) are given a value of 1. In the example network raster shown in Figure 2, the cells with an arrow represent channel cells having a value of 100. This network was generated by software TOPAZ
Figure 2. Channel network and corresponding ¯ow directions. Cells with bold outline represent network nodes; channel cells with arrow have value of 100.
Channel ordering and node indexing
963
(Garbrecht and Martz, 1995) using the Critical Source Area (CSA) concept which classi®es cells with upstream channel area greater than the CSA value as being part of the network (Jenson and Domingue, 1988; Morris and Heerdegen, 1988; Martz and Garbrecht, 1992). The ¯ow-direction raster identi®es, by code numbers, the ¯ow direction at each cell as the direction of the steepest downward slope to one of the eight adjacent cells (Fair®eld and Leymarie, 1991). The ¯ow directions of the channel cells are shown in Figure 2 by arrows.
ALGORITHM OVERVIEW
The algorithm consists of two computational steps. First, the raster image of the network is processed, a Strahler order is assigned to each channel, and a table of channel attributes is created. This table contains the coordinates of the beginning and ending node of each channel, and thereby captures the topology of the network. The numerical interpretation of the raster network begins at the source nodes and proceeds in a downstream direction along the channels to the watershed outlet node. This upstream-to-downstream processing direction is required because the order of downstream channels depends on the order of upstream channels. Second, node index numbers are assigned. The numerical processing begins at the watershed outlet node, and proceeds in an upstream direction along the channels to the source nodes. This downstreamto-upstream processing direction is required by the node indexing convention previously described. This could be accomplished by a cell-by-cell tracing of the raster network from the watershed outlet node to the source nodes. However, the cell-by-cell tracing is time consuming and is done in the ®rst step of the algorithm. A faster method that relies on the previously identi®ed channel topology data in the attribute table is implemented. Using the attribute data allows tracing of the network by jumping directly from node to node without tracing cell-by-cell the channel links. Step 1: De®nition of channel order De®nition of channel order begins by identifying all ®rst-order channels (source channels), then all second-order channels, and so forth. The algorithm explained in the following is for the ®rst two channel orders only. Higher-order channels are treated the same way as second-order channels. The process ¯ow chart for the algorithm is shown in Figure 3. First order channels The network raster is scanned to ®nd a source node (cell with value 1 in Fig. 2). On ®nding a source node, the row/column coordinates of this node are stored for later use, and the channel is traced downstream following the direction indicated by the ¯ow direction. All channel
Figure 3. Process ¯ow chart for de®nition of channel orders.
cells encountered during this trace, including the initial source node cell, are given a value of ÿ1 (negative value of the order). Negative values are used to identify those channel cells that have been processed. Every cell is also checked to determine if it is a junction node. When a junction node is encountered, the trace is terminated and a value of 2 is assigned to the junction node to de®ne it as a second-order channel cell. The row/column coordinates of the junction cell and of the previous to last cell before the junction are retained. The network raster scan then resumes at the previous source cell to identify the next source node and initiate another round of processing. Once all source nodes and ®rst-order channels have been processed, the network raster consists of background values of zero, ®rst-order channels of value ÿ1, junction cells at the downstream end of ®rst-order channels of value 2, and all remaining channel cells of value 100. This situation is shown in Figure 4A. Second order channels The network raster is scanned to detect second-order nodes (cells with value 2) at which time a similar processing as for ®rst-order channels is applied, but with three major dierences: (1) during the initial scan, cells with value 2 must be checked to determine if they are at the upstream end of a second-order channel before the trace in the downstream direction is initiated; (2) at a junction node at the end of a trace, the order of the next downstream channel may remain unchanged as indicated by a value of 2 at the junction node; under these conditions the downstream trace is re-initiated from this junction node to the next, and so forth, until a junction node with value 3 or 100 is encountered; and (3) a value of 3 is assigned to the junction node at the end of a full
964
J. Garbrecht and L. W. Martz
Figure 4. Raster network: (A) after processing ®rst-order channels, and (B) after processing secondorder channels.
trace of a second-order channel segment to indicate the beginning of a third-order channel. At the end of the evaluation of all second-order channels, the network is de®ned as shown in Figure 4B. The remaining nodes and cells are identi®ed in a similar fashion. In the previous example, the third-order channels would be indicated with a value of ÿ3. Final network raster All higher-order channels are treated in the same manner as second-order channels. When all orders have been evaluated, the entire network is de®ned by negative values. The absolute value is taken and the network is de®ned by cells having positive values on a background of cells having value 0. The value of the channel cells represents the order of the channel. In addition to this raster of the network, an attribute table was created containing (for each channel link) the order and coordinates of the beginning and ending points of each channel. This attribute table is used in the following node indexing algorithm. Step 2: Network node indexing Node indices are assigned by searching and evaluating the attribute table rather than by tracing the channels cell-by-cell in the raster network. The concept of this approach is presented ®rst, followed by the algorithm for the network search and node indexing. The process ¯ow chart is shown in Figure 5. The algorithm starts at the last data record of the channel-link attribute table which represents the watershed outlet node. The watershed outlet node is assigned index number 1. Then, the channel ending
coordinates of all channels in the attribute table are searched to identify those channels that terminate at the current node. From the identi®ed in¯ow channels, the appropriate left-hand channel is selected, the node at the upstream end of this channel is given the next successive index number, and the search of the channel attribute table for upstream in¯ow is re-initiated for this new node. When a source node is encountered the search moves in a downstream direction until a partially or fully un-evaluated node is encountered, at which time the node-by-node search and indexing resumes again in an upstream direction. This proceeds until all nodes are evaluated and the watershed outlet node is reached again. Each node processed during this search is treated dierently according to the nodes type: source node, simple junction node, complex junction node, or watershed outlet node. The decision logic and indexing at each node type is detailed in the following. The watershed outlet node is encountered twice during the node search. First, at the initiation of the indexing procedure when the index number 1 is assigned to the outlet. After indexing the watershed outlet node, the upstream node of the in¯owing channel becomes the current node. Second, at the end of the indexing procedure and completion of the algorithm. If the outlet node happens to be a junction node, then the node is treated as a junction node. Source nodes are encountered only once during the network search. When a source node is encountered, the next successive index number is assigned,
Channel ordering and node indexing
965
at the current time. Non-consecutive indices indicate that a trace along the left in¯ow channel is undertaken before the index for the next in¯ow at the complex junction is assigned. This dierentiation is important to ensure consistency in channel Strahler order when decomposing complex junctions into a sequence of simple junctions. Once the channel network nodes are indexed, the morphometric properties of the network can be evaluated and the node indices can be processed by software SEQDRV (Garbrecht, 1992) to optimize the sequence of channel link evaluation for cascadetype ¯ow routing.
CONCLUSIONS
Figure 5. Process ¯ow chart for network node indexing.
and the next downstream node is selected as the current node in the search. A simple junction node is encountered three times during a network search. Once when moving up the left-most in¯ow; once when returning from the left-most in¯ow and moving up the right in¯ow; and once when returning from the right in¯ow. The ®rst time a simple junction node is encountered, the next successive index number is assigned and the left-most in¯ow channel link is followed upstream to the next node. The second time a junction node is encountered, no index number is assigned but the second in¯owing channel link is followed upstream. The third time a simple junction node is encountered, no index number is assigned and the next downstream node becomes the current node. Complex junction nodes are the result of the discrete representation of the network. During the network search they are encountered as many times as the junction has in¯ow channels plus one. The indexing is similar to that of simple junctions, except that multiple node indices are assigned to account for the extra in¯ow channels. The additional indices are in anticipation of the decomposition of the complex junction node into a sequence of simple junction nodes as required by subsequent node processing software (Garbrecht, 1988) and for a more accurate de®nition of network composition parameters (Garbrecht and Martz, 1995). The indices assigned to a complex junction node may be consecutive or non-consecutive depending on the sequence and Strahler order of the in¯owing channel links. Consecutive indices indicate that lower order channel links enter from the left and no upstream trace is undertaken along these channels
An algorithm that automatically determines channel Strahler order and node indices for raster channel networks is presented. Strahler orders are determined by a cell-by-cell trace of the raster network in a downstream direction along the channels beginning at the source nodes to the watershed outlet. During this trace, beginning and ending coordinates of each channel link are retained in an attribute table. The subsequent node indexing uses the channel topology information in the attribute table to search the network from downstream to upstream starting at the watershed outlet node and following the left hand pattern described in Croley (1980) and Garbrecht (1988). During this search, nodes are indexed consecutively. The algorithm accounts for complex junction nodes that are particular to raster channel networks and assigns indices that allow a subsequent decomposition of the complex node into a sequence of simple junction nodes. The Strahler orders and the decomposition of complex junction nodes are necessary for a more accurate determination of network composition parameters and for subsequent optimization of the sequence in which channels must be evaluated for cascade type ¯ow routing. The algorithm presented herein makes possible the automated quanti®cation of network structure from raster network images and greatly enhances the direct linkage of GIS generated channel networks and hydrologic and hydraulic models. Code is available by anonymous FTP from FTP.IAMG.ORG.
REFERENCES Band, L. E. (1986) Topographic partitioning of watersheds with digital elevation models. Water Resources Research 22(1), 15±24. Beasley, D. B., Huggins, L. F. and Monke, E. J. (1980) ANSWERS: A model for watershed planning. Transactions of American Society of Agricultural Engineers 23(4), 938±944. Croley, T. E. (1980) A micro-hydrology computation ordering algorithm. Journal of Hydrology 48, 221±236.
966
J. Garbrecht and L. W. Martz
Fair®eld, J. and Leymarie, P. (1991) Drainage networks from grid digital elevation models. Water Resources Research 27(4), 29±61. Feldman, A. D. (1995) HEC-1 ¯ood routing package. In Computer Models of Watershed Hydrology, ed. V. J. Singh, pp. 119±150. Water Resources Publication, Highlands Ranch, Colorado. Garbrecht, J. (1988) Determination of the execution sequence of channel ¯ow for cascade routing in a drainage network. Hydrosoft 1(3), 129±138. Garbrecht, J. (1992) Documentation of model and computer program SEQDRV: topologic evaluation of drainage networks for cascade routing. HYDROWAR Reports Division, Hydrology Days Publications, 1005 Country Club Road, Fort Collins, CO 80524, Report No. 92.3, 81 pp. Garbrecht, J. and Martz, L. W. (1995) TOPAZ: An automated digital landscape analysis tool for topographic evaluation, drainage identi®cation, watershed segmentation, and subcatchment parameterization. Overview. ARS Publication NAWQL 95-1, U. S. Department of Agriculture, Agricultural Research Service, Oklahoma 74702, 16 pp. Goodrich, D. C., Schmugge, T. J., Jackson, T. J., Unkrich, C. L., Keefer, T. O., Parry, R., Bach, L. B. and Amer, S. A. (1994) Runo simulation sensitivity to remotely sensed initial soil water content. Water Resources Research 30(5), 1393±1405. Horton, R. E. (1945) Erosional development of streams and their drainage basins, hydrophysical approach to quantitative morphology. Bulletin of the Geological Society of America 56, 275±370. Jenson, S. K. and Domingue, J. O. (1988) Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering and Remote Sensing 54(11), 1593±1600. Knisel, W. G. (ed). (1982) CREAMS: A ®eld-scale model for chemical runo, and erosion form agricultural management systems. USDA-Science and Education Administration, Conservation Report, 26 pp.
Martz, L. W. and Garbrecht, J. (1992) Numerical de®nition of drainage networks and subcatchment areas from digital elevation models. Computers & Geosciences 18(6), 747±761. Morris, D. G. and Heerdegen, R. G. (1988) Automatically derived catchment boundary and channel networks and their hydrological applications. Geomorphology 1(2), 131±141. Plate, E. J., Ihringer, J. and Lutz, W. (1988) Operational models for ¯ood calculations. Journal of Hydrology 100, 489±506. Ponce, V. M., Osmolski, Z. and Smutzer, D. (1985) Large basin deterministic hydrology: a case study. Journal of Hydraulic Engineering 111(9), 1227±1245. Shreve, R. L. (1966) Statistical law of stream numbers. Journal of Geology 74, 17±37. Shreve, R. L. (1967) In®nite topologically random channel networks. Journal of Geology 75, 178±186. Singh, V. P. (ed.) (1995) Computer Models of Watershed Hydrology. Water Resources Publications, Highland Ranch, Colorado, 1130 pp. Smart, J. S. (1972) Channel networks. Advances in Hydrosciences 8, 305±346. Smart, J. S. (1978) The analysis of drainage network composition. Earth Surface Processes 3, 129±170. Smith, R. E., Goodrich, D. C., Woolhiser, D. A. and Unkrich, C. L. (1995) KINEROSÐA kinematic runo and erosion model. In Computer Models of Watershed Hydrology, ed. V. J. Singh, pp. 697±732. Water Resources Publication, Highland Ranch, Colorado. Strahler, A. N. (1957) Quantitative analysis of watershed geomorphology. Transactions of the American Geophysical Union 38(6), 913±920. Vieux, B. E. (1991) Geographic information systems and non-point source water quality and quantity modeling. Hydrological Processes 5(1), 101±113. Young, R. A., Onstad, C. A., Bosch, D. D. and Anderson, W. P. (1987) AgNPS: agricultural non-point-source pollution model: a watershed analysis tool. USDAARS, Conservation Research Report, 35 pp.