A GRASS Program for Determining the Topology of Stream Networks

John F. Stamm

Department of Geological Sciences, Case Western Reserve University,

Cleveland, OH 44106-7216, email: jfs5@po.cwru.edu

Introduction

This paper presents a GRASS (geographical resource analysis support system; Shapiro et al., 1992; United States Army Corps of Engineers, 1993) program for determining the topology of stream networks. The program inputs raster files generated by the GRASS program r.watershed. Because it determines the relationships of tributary streams, it is called r.tribs. The input files required are:

stream: a raster map of stream segments; stream segments are labelled using integer values greater than or equal to 2,

accumulation: a raster map of the number of cells that drain through each cell; absolute values is the amount of overland flow that is routed through a pixel,

drainage: a raster map of drainage directions; if there is no flow direction, a value of -1 is assigned. Otherwise, the integers 1 through 8 are assigned to the compass directions shown below:

    			5	 6	7
    			4	-1	8
    			3	 2	1

For example, a value of 2 means that the pixel drains south, a value of 5 means the pixel drains to the northwest.

Two tables are output by r.tribs. The first table list the tributaries associated with each stream segment. The second table list streams and tributaries in an order such that a stream and its tributaries are printed only if the tributaries have been previously listed. That is, first order streams are listed first, streams with tributaries that are first order streams are listed second, and so on. The program can also be run in debugging mode, in which more detailed information is printed.

The program r.tribs was written so that GRASS can be used to generate input for a runoff and erosion model called KINEROS (Smith et al., in press; Woolhiser et al., 1990). KINEROS represents a watershed as a set of related elements. Elements may be hillslopes, channels or ponds. A computational order must be specified so that boundary conditions for an element, such as the amount of water contributed by lateral hillslopes and upstream tributaries, are available. The program r.tribs provides KINEROS with that computational order.

Input Files

We assume that the program r.watershed has been run, and that the input files stream, accumulation and drainage are available. Assume these files are called stream.0, accum.0 and drain.0, respectively. The command-line version of r.tribs is:

GRASS 4.1 > r.tribs st=stream.0 ac=accum.0 dr=drain.0

Alternately, the command r.tribs can be entered on the command line:

GRASS 4.1 > r.tribs

and the user will be prompted for these files.

Methodology and Example of Output

This is followed by a listing of the tributaries associated with each stream. The tributaries of a stream are determined by using the following methodology.

a. The stream and accumulation files are scanned to determine the locations of pixels that are stream segments and have the lowest accumulation of runoff for that segment. These points mark the locations of the heads of the stream segments. Locations of the head of a stream are stored for the next step.

b. Points within one pixel of the head are scanned to determine if they are part of a different stream segment and if they drain into the head of this stream. The drainage raster file is used to determine the drainage direction. Any point meeting this criteria are listed as tributary streams.

c. All points with no tributaries are first order. Tributary values are set to zero.

d. The label of the stream and its tributaries are stored and printed to the screen.

We have used r.tribs to generate output for a small watershed (17.5 km2) in Idaho (Horse Creek) in which we used r.watershed to define a total of 50 stream segments. Note that stream segments and tributaries are labelled using even integers ranging from 2 to 102. Stream segment 2 is located at the mouth of the watershed.

The r.tribs program first outputs the number of rows (nrows) and columns (ncols) in the raster files. Then the maximum and minimum labels of stream segments are printed. This is followed by two tables. The first table lists streams and their tributaries. An example of this table is given below.

``

`nrows: 138`

`ncols: 273`

``

`Max stream index = 102`

`Min stream index = 2`

`Stream: 2 Tributary 0: 6 Tributary 1: 4`

`Stream: 4 Tributary 0: 0 Tributary 1: 0`

`Stream: 6 Tributary 0: 14 Tributary 1: 8`

`Stream: 8 Tributary 0: 12 Tributary 1: 10`

`Stream: 10 Tributary 0: 0 Tributary 1: 0`

`Stream: 12 Tributary 0: 0 Tributary 1: 0`

`Stream: 14 Tributary 0: 22 Tributary 1: 16`

`Stream: 16 Tributary 0: 20 Tributary 1: 18`

`Stream: 18 Tributary 0: 0 Tributary 1: 0`

`Stream: 20 Tributary 0: 0 Tributary 1: 0`

`Stream: 22 Tributary 0: 102 Tributary 1: 24`

`Stream: 24 Tributary 0: 36 Tributary 1: 26`

`Stream: 26 Tributary 0: 30 Tributary 1: 28`

`Stream: 28 Tributary 0: 0 Tributary 1: 0`

`Stream: 30 Tributary 0: 34 Tributary 1: 32`

`Stream: 32 Tributary 0: 0 Tributary 1: 0`

`Stream: 34 Tributary 0: 0 Tributary 1: 0`

`Stream: 36 Tributary 0: 100 Tributary 1: 38`

`Stream: 38 Tributary 0: 42 Tributary 1: 40`

`Stream: 40 Tributary 0: 0 Tributary 1: 0`

`Stream: 42 Tributary 0: 98 Tributary 1: 44`

`Stream: 44 Tributary 0: 96 Tributary 1: 46`

`Stream: 46 Tributary 0: 50 Tributary 1: 48`

`Stream: 48 Tributary 0: 0 Tributary 1: 0`

`Stream: 50 Tributary 0: 58 Tributary 1: 52`

`Stream: 52 Tributary 0: 56 Tributary 1: 54`

`Stream: 54 Tributary 0: 0 Tributary 1: 0`

`Stream: 56 Tributary 0: 0 Tributary 1: 0`

`Stream: 58 Tributary 0: 90 Tributary 1: 60`

`Stream: 60 Tributary 0: 88 Tributary 1: 62`

`Stream: 62 Tributary 0: 66 Tributary 1: 64`

`Stream: 64 Tributary 0: 0 Tributary 1: 0`

`Stream: 66 Tributary 0: 74 Tributary 1: 68`

`Stream: 68 Tributary 0: 72 Tributary 1: 70`

`Stream: 70 Tributary 0: 0 Tributary 1: 0`

`Stream: 72 Tributary 0: 0 Tributary 1: 0`

`Stream: 74 Tributary 0: 82 Tributary 1: 76`

`Stream: 76 Tributary 0: 80 Tributary 1: 78`

`Stream: 78 Tributary 0: 0 Tributary 1: 0`

`Stream: 80 Tributary 0: 0 Tributary 1: 0`

`Stream: 82 Tributary 0: 86 Tributary 1: 84`

`Stream: 84 Tributary 0: 0 Tributary 1: 0`

`Stream: 86 Tributary 0: 0 Tributary 1: 0`

`Stream: 88 Tributary 0: 0 Tributary 1: 0`

`Stream: 90 Tributary 0: 94 Tributary 1: 92`

`Stream: 92 Tributary 0: 0 Tributary 1: 0`

`Stream: 94 Tributary 0: 0 Tributary 1: 0`

`Stream: 96 Tributary 0: 0 Tributary 1: 0`

`Stream: 98 Tributary 0: 0 Tributary 1: 0`

`Stream: 100 Tributary 0: 0 Tributary 1: 0`

`Stream: 102 Tributary 0: 0 Tributary 1: 0`

``

If the program is run in debugging mode, then the values of stream, accumulation and direction of the head of the stream and the surrounding 8 pixels are also printed. Set the value of `DB_FIND_TRIBS` (defined in the beginning of the routine `find_tribs)` to a value of `1` and recompile the program to activate debugging mode.

The second table that is output by r.tribs list the streams and their associated tributaries in their proper computational order. The program loops through the data listed above several times. First streams with no tributaries are listed (`LOOP 1`). Then streams with only first-order streams as tributaries are listed (`LOOP 2`). Then, streams with tributaries listed in previous loops are listed. The program continues until all streams have been listed. The variable "`Order`" is also listed, which can be interpreted as the computational order. This is the order in which programs, such as KINEROS, must consider streams in the network such that data for tributaries will be available when considering the listed stream. An example listing of the computational order of streams is given below.

``

``

`Computational Order of Stream Segments:`

``

`LOOP: 1`

`Order: 0 Stream: 4 Tributary 0: 0 Tributary 1: 0`

`Order: 1 Stream: 10 Tributary 0: 0 Tributary 1: 0`

`Order: 2 Stream: 12 Tributary 0: 0 Tributary 1: 0`

`Order: 3 Stream: 18 Tributary 0: 0 Tributary 1: 0`

`Order: 4 Stream: 20 Tributary 0: 0 Tributary 1: 0`

`Order: 5 Stream: 28 Tributary 0: 0 Tributary 1: 0`

`Order: 6 Stream: 32 Tributary 0: 0 Tributary 1: 0`

`Order: 7 Stream: 34 Tributary 0: 0 Tributary 1: 0`

`Order: 8 Stream: 40 Tributary 0: 0 Tributary 1: 0`

`Order: 9 Stream: 48 Tributary 0: 0 Tributary 1: 0`

`Order: 10 Stream: 54 Tributary 0: 0 Tributary 1: 0`

`Order: 11 Stream: 56 Tributary 0: 0 Tributary 1: 0`

`Order: 12 Stream: 64 Tributary 0: 0 Tributary 1: 0`

`Order: 13 Stream: 70 Tributary 0: 0 Tributary 1: 0`

`Order: 14 Stream: 72 Tributary 0: 0 Tributary 1: 0`

`Order: 15 Stream: 78 Tributary 0: 0 Tributary 1: 0`

`Order: 16 Stream: 80 Tributary 0: 0 Tributary 1: 0`

`Order: 17 Stream: 84 Tributary 0: 0 Tributary 1: 0`

`Order: 18 Stream: 86 Tributary 0: 0 Tributary 1: 0`

`Order: 19 Stream: 88 Tributary 0: 0 Tributary 1: 0`

`Order: 20 Stream: 92 Tributary 0: 0 Tributary 1: 0`

`Order: 21 Stream: 94 Tributary 0: 0 Tributary 1: 0`

`Order: 22 Stream: 96 Tributary 0: 0 Tributary 1: 0`

`Order: 23 Stream: 98 Tributary 0: 0 Tributary 1: 0`

`Order: 24 Stream: 100 Tributary 0: 0 Tributary 1: 0`

`Order: 25 Stream: 102 Tributary 0: 0 Tributary 1: 0`

`-----------------------------------------------------------------`

``

`LOOP: 2`

`Order: 26 Stream: 8 Tributary 0: 12 Tributary 1: 10`

`Order: 27 Stream: 16 Tributary 0: 20 Tributary 1: 18`

`Order: 28 Stream: 30 Tributary 0: 34 Tributary 1: 32`

`Order: 29 Stream: 52 Tributary 0: 56 Tributary 1: 54`

`Order: 30 Stream: 68 Tributary 0: 72 Tributary 1: 70`

`Order: 31 Stream: 76 Tributary 0: 80 Tributary 1: 78`

`Order: 32 Stream: 82 Tributary 0: 86 Tributary 1: 84`

`Order: 33 Stream: 90 Tributary 0: 94 Tributary 1: 92`

`-----------------------------------------------------------------`

``

`LOOP: 3`

`Order: 34 Stream: 26 Tributary 0: 30 Tributary 1: 28`

`Order: 35 Stream: 74 Tributary 0: 82 Tributary 1: 76`

`-----------------------------------------------------------------`

``

`LOOP: 4`

`Order: 36 Stream: 66 Tributary 0: 74 Tributary 1: 68`

`-----------------------------------------------------------------`

``

`LOOP: 5`

`Order: 37 Stream: 62 Tributary 0: 66 Tributary 1: 64`

`-----------------------------------------------------------------`

``

`LOOP: 6`

`Order: 38 Stream: 60 Tributary 0: 88 Tributary 1: 62`

`-----------------------------------------------------------------`

``

`LOOP: 7`

`Order: 39 Stream: 58 Tributary 0: 90 Tributary 1: 60`

`-----------------------------------------------------------------`

``

`LOOP: 8`

`Order: 40 Stream: 50 Tributary 0: 58 Tributary 1: 52`

`-----------------------------------------------------------------`

``

`LOOP: 9`

`Order: 41 Stream: 46 Tributary 0: 50 Tributary 1: 48`

`-----------------------------------------------------------------`

``

`LOOP: 10`

`Order: 42 Stream: 44 Tributary 0: 96 Tributary 1: 46`

`-----------------------------------------------------------------`

``

`LOOP: 11`

`Order: 43 Stream: 42 Tributary 0: 98 Tributary 1: 44`

`-----------------------------------------------------------------`

``

`LOOP: 12`

`Order: 44 Stream: 38 Tributary 0: 42 Tributary 1: 40`

`-----------------------------------------------------------------`

``

`LOOP: 13`

`Order: 45 Stream: 36 Tributary 0: 100 Tributary 1: 38`

`-----------------------------------------------------------------`

``

`LOOP: 14`

`Order: 46 Stream: 24 Tributary 0: 36 Tributary 1: 26`

`-----------------------------------------------------------------`

``

`LOOP: 15`

`Order: 47 Stream: 22 Tributary 0: 102 Tributary 1: 24`

`-----------------------------------------------------------------`

``

`LOOP: 16`

`Order: 48 Stream: 14 Tributary 0: 22 Tributary 1: 16`

`-----------------------------------------------------------------`

``

`LOOP: 17`

`Order: 49 Stream: 6 Tributary 0: 14 Tributary 1: 8`

`-----------------------------------------------------------------`

``

`LOOP: 18`

`Order: 50 Stream: 2 Tributary 0: 6 Tributary 1: 4`

``

``

``Obtaining r.tribs via FTP or email

The C programs that are required to generate r.tribs are available via anonymous ftp to moon.cecer.army.mil. I have printed out the main segment of the r.tribs program in Appendix I. This was done to illustrate how routines in the GRASS library are used to read in raster data. Routines from the GRASS library begin with `"G_`". Two utility programs are also used: imatrix and ivector. These are discussed by Press and others (1989), and are used to allocate space for arrays. Comments in the code discuss the details of each routine. Appendix II lists the makefile used to compile r.tribs. Note that the variable `GIS` must be edited so that the proper path to the grass directory is specified.

As previously mentioned, the program can be obtained by anonymous ftp to moon.cecer.army.mil. Change to the "incoming/r.tribs" directly to get the files. The program can also be obtained by contacting the author via email (jfs5@po.cwru.edu). The code has been commented to help the user to understand the program structure.

References

Press, W.H., Flannery, B.P., Teukolsky, S.A., and Vetterling, W.T., 1989, Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, 735 pp.

Shapiro, M., Westervelt, J., Gerdes, D., Larson, M., and Brownfield, K.R., 1992, GRASS 4.0 Programmer's Manual, U.S. Army Construction Engineering Research Laboratory, Champaign, Illinois, 292 pp.

Smith, R.E., Goodrich, D.C., Woolhiser, D.A., and Unkrich, C.L., in press, A KINematic Runoff and EROSion Model, in: V.P. Singh (Ed.), Computer Models of Watershed Hydrology, Water Resources Pub., Highlands Ranch, Colorado

United States Army Corps of Engineers, 1993, GRASS4.1 Users's Reference Manual, U.S. Army Construction Engineering Research Laboratories, Champaign, Illinois, 556 pp.

Woolhiser, D.A., Smith, R.E., Goodrich, D.C., 1990, KINEROS, A Kinematic Runoff and Erosion Model: Documentation and Users Manual, U.S. Department of Agriculture, Agricultural Research Service, ARS-77, 130 pp.

Appendix I: main.c code for r.tribs

`#include "gis.h"`

`#include <stdio.h>`

``

`int **imatrix();`

`int *ivector();`

``

`/*`

`* Program to determing the topology of a stream network. A`

`* table is generated that reports the tributaries that are`

`* at the head of each stream segement. The computational order`

`* of streams is also determinted.`

`*`

`* Written by:`

`* Dr. John F. Stamm`

`* Department of Geological Sciences`

`* Case Western Reserve University`

`* Cleveland, OH 44106-7216`

`* email: jfs5@po.cwru.edu`

`*/`

``

`main(argc,argv)`

`int argc;`

`char *argv[];`

`{`

`/*`

`* Matricies`

`*/`

`int **accum;`

`int **chann;`

`int **aspect;`

``

`CELL *cell;`

``

`char *chann_name;`

`char *accum_name;`

`char *aspect_name;`

`char *mapset;`

``

`int col;`

`int fd_accum;`

`int fd_chann;`

`int fd_aspect;`

`int ncols;`

`int nrows;`

`int row;`

``

`struct {`

`struct Option *accum ;`

`struct Option *chann ;`

`struct Option *aspect ;`

`} parm;`

``

`/*`

`* Allocate memory for the Option structure and return a`

`* pointer to this structure. Do this for the structured`

`* variables parm.accum and parm.chann.`

`*`

`* Set values for parm.accum`

`*/`

`parm.accum = G_define_option() ;`

`parm.accum->key = "accumulation";`

`parm.accum->type = TYPE_STRING;`

`parm.accum->required = YES;`

`parm.accum->gisprompt = "old,cell,raster" ;`

`parm.accum->description= "Name of the ACCUMULATION map" ;`

`/*`

`* Set values for parm.chann`

`*/`

`parm.chann = G_define_option() ;`

`parm.chann->key = "stream";`

`parm.chann->type = TYPE_STRING;`

`parm.chann->required = YES;`

`parm.chann->gisprompt = "old,cell,raster" ;`

`parm.chann->description= "Name of the STREAM map" ;`

`/*`

`* Set values for parm.aspect`

`*/`

`parm.aspect = G_define_option() ;`

`parm.aspect->key = "drainage";`

`parm.aspect->type = TYPE_STRING;`

`parm.aspect->required = YES;`

`parm.aspect->gisprompt = "old,cell,raster" ;`

`parm.aspect->description= "Name of the DRAINAGE DIRECTION`

`map";`

``

`/*`

`* Initailize GIS library for this program`

`*/`

``

`G_gisinit(argv[0]);`

``

`/*`

`* Parse values from the command line. If this is not`

`* successful, then display a usage statement and exit.`

`*/`

``

`if (G_parser(argc, argv))`

`exit (-1);`

`accum_name = parm.accum->answer;`

`chann_name = parm.chann->answer;`

`aspect_name = parm.aspect->answer;`

``

`/*`

`* Find the name of mapset that we are going to use.`

`*/`

``

`mapset = G_find_cell2 (accum_name, "");`

``

`if (mapset == NULL) {`

`char msg[100];`

`sprintf (msg, "%s: <%s> cellfile not found\n",`

`G_program_name(), accum_name);`

`G_fatal_error (msg);`

`exit(1);`

`}`

``

`/*`

`* Open the cell files in "mapset".`

`*/`

``

`fd_accum = G_open_cell_old (accum_name, mapset);`

`if (fd_accum < 0)`

`exit(1);`

``

`fd_chann = G_open_cell_old (chann_name, mapset);`

`if (fd_chann < 0)`

`exit(1);`

``

`fd_aspect = G_open_cell_old (aspect_name, mapset);`

`if (fd_aspect < 0)`

`exit(1);`

``

`/*`

`* Open up a vector that is just long enough to hold one`

`* row of data.`

`*/`

``

`cell = G_allocate_cell_buf();`

``

`/*`

`* Determine the number of rows and columns.`

`*/`

``

`nrows = G_window_rows();`

`ncols = G_window_cols();`

`printf ("\n", nrows);`

`printf ("nrows: %d\n", nrows);`

`printf ("ncols: %d\n", ncols);`

``

`/*`

`* Allocate memory for matricies.`

`*/`

``

`accum = imatrix(0,nrows,0,ncols);`

`chann = imatrix(0,nrows,0,ncols);`

`aspect = imatrix(0,nrows,0,ncols);`

``

`/*`

`* Process DEM and Channel files.`

`*/`

``

`for (row=(nrows-1); row>=0; row--) {`

`if(G_get_map_row (fd_accum, cell, row) < 0)`

`exit(1);`

`for (col = 0; col < ncols; col++) {`

`accum[row][col] = (int)cell[col];`

`}`

`if(G_get_map_row (fd_chann, cell, row) < 0)`

`exit(1);`

`for (col = 0; col < ncols; col++) {`

`chann[row][col] = (int)cell[col];`

`}`

`if(G_get_map_row (fd_aspect, cell, row) < 0)`

`exit(1);`

`for (col = 0; col < ncols; col++) {`

`aspect[row][col] = (int)cell[col];`

`}`

`}`

``

`/*`

`* Compute the topology of the network.`

`*/`

``

`(void)find_tribs(nrows, ncols, accum, chann, aspect);`

``

`exit(0);`

`}`

``

``

Last changed: $Date: 2002/01/25 05:45:34 $