Exercise 2: Delineation of Watersheds and Streams

CE 397 GIS in Water Resources
University of Texas at Austin

Table of Contents

Goals of the Exercise

The purpose of this exercise is to acquaint you with the Arc/Info Grid commands that do watershed and stream network delineation using a digital elevation model of the land surface terrain. You will use a 3 arc-second DEM for the West portion of the Austin mapsheet which includes part of Lake Travis and the areas East and North of Lake Travis Once you've mastered this exercise, you can do similar work anywhere else in the US using 3 arc-second data downloaded via Internet, although we won't download data in this exercise.

Computer and Data Requirements

To carry out this exercise you need access to a workstation with Arc/Info version 6.1, or 7.0 installed, including the Grid module which is the focus of this exercise. Instructions are presented later for obtaining the required data.

Assignment

(1) Prepare a map showing watersheds and streams delineated from the DEM and a table with the properties of these watersheds.

(2) Answer other questions about the analysis that you perform as noted on the instructions in the following pages. Turn in these answers along with your map and give a brief comment on any difficulties that you had with this exercise.

Procedure

I am assuming that you have logged on to your account on the workstation. If you have not done that, do so now. Make a new directory for this exercise and change directories so that you are in it:

$ mkdir ex2
$cd ex2

1. Acquire the Digital Elevation Model

We have aquired the DEM file from the internet and placed it in an export format file called ausfil.e00 (e zero zero not e o o) in the directory /class/watershed/gisfiles. Yoy can also download the files via anonymous ftp from ftp.crwr.utexas.edu/pub/gisclass/ex2. Instructions on how to use anonymous ftp.

Alternatively you can aquire the Digital Elevation Model (DEM) over the internet from the U.S. Geological Survey. For this exercise you may work with any 3" DEM of the U.S.

You are now ready to run ARC/INFO. Type the following to start up ARC/INFO and import the grid you need.

$ arc

Import the grid ausfil as a directory of files

Arc: import grid ausfil ausfil

2. Display the DEM

Now you can take a look at this grid. To do so, start up grid and open an Xwindow

Arc: grid [this starts up the grid system]

Grid: display 9999 [this produces a graphics display window --- drag it to a convenient place on the screen]

The graphics display window sometimes takes a little time to come up. (If you get a response to the effect "Xserver not defined" it means that you haven't specified where the display is to be produced. Quit from grid and arc and use the "export DISPLAY=mymachine:0" to define where the graphics are to be displayed (on mymachine), restart arc and grid and continue.)

First display the grid that you will be using:

Grid: mapextent ausfil

Grid: gridpaint ausfil

You'll see something with rather lurid colors appear. This is the default option for gridpaint. You can get a gray shaded image using

Grid: clear [clear the display]

Grid: gridpaint ausfil value linear nowrap gray

The white areas are higher elevations and the gray areas are lower elevations.

Any time you get a response "map extent is undefined" you have to redefine it again as the geographic extent of one of your active grids, eg mape ausfil.

You can query the display to find out what the elevation values (m) are at different points on the grid.

Grid: cellvalue ausfil *

When the cursor is on the graphics window a pair of crosslines appears. Click the mouse anywhere in the display window and the cell value (elevation) corresponding to the intersection of the crosslines will be reported.

(Pressing "9" on the keyboard while the mouse is located in the display window will exit this mode). If you get in a jam, type "Control-C" to get back to grid.

You can get some descriptive statistics about the grid by typing

Grid: describe ausfil

To be turned in:

(1) Make a note of the maximum, minimum and average elevation of the cells in the grid

(2) How many rows and columns does the grid have? What is the cellsize? What range of latitude and longitude does this grid cover?

Note: You can open a Notepad window, use the cursor to highlight the area of the Decterm text screen whose contents you want to save, use the center key of the three keys on the alpha mouse to paste text highlighted in the Decterm text display window into the Notepad window for later reference and printing. Save your Notepad file each time you do this so that you have a record of what you have done.

(3) Selecting a region for Analysis

If you have lots of time, you can analyze the whole of the ausfil grid, but to limit the amount of time each operation takes, it is better to select an area within the ausfil grid for analysis. If you have the ausfil grid displayed on the screen, proceed to the setwindow step; otherwise, repeat the commands you typed earlier:

Grid: display 9999

Grid: mape ausfil

Grid: gridpaint ausfil value linear nowrap gray

You are now going to set an analysis window within the ausfil grid

Grid: setwindow * [this allows a user defined box to be selected within the grid]

Move the cursor onto the display, click the button in one corner of your desired area, then move the cursor to the opposite corner and click again. You should see a box appear on the screen as you do this. If you don't see a box, do the setwindow * again and repeat the exercise. Choose and area that is about 1/4 as high and 1/4 as wide as the original grid.

Grid: ausdem = ausfil [creates a new grid called ausdem containing only the portion of ausfil that is within the analysis window]

Grid: clear [clears display]

Grid: gridpaint ausdem value linear nowrap gray [shows where ausdem is in relation to ausfil]

Grid: mape ausdem [resizes display to extent of the smaller grid]

Grid: gridpaint ausdem value linear nowrap gray [gives a closer look at ausdem]

Grid: describe ausdem

Check that you have a grid with between 200 and 300 rows and columns. Larger than this slows down the processing; smaller than this means there isn't much to see.

***************************************************************************

If your grid is the wrong size, kill it and start over by defining a new setwindow

Grid: setwindow ausfil

Grid: kill ausdem all

Grid: setwindow * [ etc, etc as above]

***************************************************************************

(4) Projecting the Grid onto a Flat Map Surface

We are going to use a standard US albers projection to convert the grid points defined on a lat-long mesh to points defined on a cartesian (x,y) mesh. This particular set of parameters for map projection is a standard set used by the US Geological Survey and other agencies to present national data sets for the United States.

When typing grid commands, always leave a space around any operator like (, =, #, 100, project, etc). If the commands don't work, try again leaving spaces liberally around each item in the command line. This helps the system interpret what you have typed.

Grid: ausalb = project ( ausdem, #, #, 100 ) [initiates a projection dialog that follows]

The characteristics of the original grid (ausdem) were defined when it was created from ausfil, but you must specify the characteristics of the new grid (ausalb). You will prompted for this dialog:

Project: output

Project: projection albers

Project: datum NAD83

Project: units meters

Project: parameters

1st standard parallel [0 0 0.000]: 29 30 00

2nd standard parallel [0 0 0.000]: 45 30 00

central meridian [0 0 0.000]: -96 00 00

latitude of projection origin [0 0 0.000]: 23 00 00

false easting (meters) [0.00000]: 0

false northing (meters) [0.00000]: 0

Project: end

Grid: [seeing the grid prompt appear tells you that the projection is done].

Grid: listgrids [lists the existing grids --- ausfil, ausdem, ausalb]

This will create a grid in albers projection with a cellsize of 100 m. The # allows for default options to be chosen. You have now described how the output should look. The dimensions of your grid are now in meters referenced to (0,0) at the intersection of the central meridian and the latitude of the projection origin.

If you don't know the syntax of an arc command, just type the command and you'll see it echoed back to you. If you need help, type help and an online help facility will come up.

Grid: describe ausalb

To be turned in: the cellsize of the new grid, its number of rows and columns and the minimum, maximum and average of the values shown in the grid. What area in km2 does this grid cover?

Use the text editor on the workstation and the keys for cutting pasting data to save the statistics of ausalb for printing as part of the material you turn in for this exercise

5. Filling in the Pits

In order that water can "flow" across the landscape, any spurious pits have to be filled in. First you have to change the setwindow analysis area to the new set of projected coordinates:

Grid: setwindow ausalb

Grid: fill ausalb albfil SINK 100 #

As you watch the screen, you'll see some things about number of sinks here. If you have a big grid, filling pits will take a while. This command fills all sinks that are less than 100 m deep. I've had some problems with this function. If necessary, you can copy a complete filled and projected grid in the file /res2/maidment/gisdata/ausfill.e00 and import it as you did earlier for the ausfil grid. To check whether you have problems:

Grid: describe albfil

If your grid has 0 for min, max and mean you have a problem! In that event, copy the filled grid ausfill and select a window from it for analysis and proceed on.

6. Flowdirection Grid

The cells flow to their nearest neighbour along 1 of 8 compass directions labelled as East = 1, SE = 2, S = 4, SW=8, W=16, NW=32, N=64, NE=128

Grid: ausfdr = flowdirection ( albfil )

The flowdirection function in grid assigns to each cell a number corresponding to which of the 8 neighbouring cells lies on the path of steepest descent.

You can get a sense of flow direction in the watershed by typing:

Grid: clear

Grid: mape ausfdr

Grid: gridpaint ausfdr

(numbers 32, 64, 128 all show up as black in this plot).

You can again check the individual values of flow direction by typing

Grid: cellvalue ausfdr *

and querying individual locations. You will get "value" which is the number in the cell and "count" which is the number of cells possessing the same value. Again, type 9 while the cursor is over the graphic display to get out of this querying mode and get the grid prompt back again. Make sure your cursor is within your text display window before continuing to type grid commands.

The different flowdirection areas in grid are called "zones" and grid carries an attribute table for zones like the arc and polygon attribute tables. It is called the "value attribute table" (vat). To see the value attribute table for this grid type

Grid: list ausfdr.vat

Under "value" are listed the 8 direction numbers, and under "count" the number of cells which have that value.

To be turned in: Note down the number of cells corresponding to each flow direction. What are the principal flow directions on this watershed?

(you can copy the vat to the Text Editor for later reference)

7. Flow Accumulation Grid

Now, lets get the flowaccumulation grid. This grid counts all the cells upstream of a given cell:

Grid: ausfac = flowaccumulation ( ausfdr ) [spaces before and after the =]

Lets display the flowaccumulation grid as a grayshaded image:

Grid: clear

Grid: gridpaint ausfac value linear nowrap gray

Just from this picture you can see the approximate locations of the streams. Pretty neat!! Using the cellvalue command to query the flow accumulation. You'll see big numbers on the streams (sometimes its hard to click on them just right) and low values on the land surface. Cells with flowaccumulation of zero are on ridge lines.

Grid: cellvalue ausfac *

and 9 to quit with the cursor on the graphics screen.

To get a closer view of some portion of the screen, click the cursor on the "Pan/Zoom" button on the top left corner of the graphics screen, drag the pulldown menu to "Create" and release the mouse button. A new smaller graphics window will appear. Repeat the gridpaint command given above to display the flow accumulation grid in both windows. Use the zoomed in window to click on cells and check their flow accumulation.

Describe the statistics of the grid using

Grid: describe ausfac

To be turned in: The maximum and the mean of the flow accumulation grid.

8. Stream Network Delineation.

Now we will create a grid of the streams. You will get a different result depending on the threshold value of flow accumulation used. Each cell is 100x100 m = 10,000 m2 = 0.01 km2 in area, so 1000 cells = 10 km2 in drainage area, and 10,000 cells = 100 km2 in area. Lets delineate the streams for a 10,000 cell flow accumulation threshold (this means the minimum drainage area upstream necessary to define a stream).

Grid: str10 = con ( ausfac > 10000, 1 ) [The 1 means a stream is labeled with 1]

[Remember to leave spaces before and after the = and >]

Running... 100%

The con function is a boolean query which assigns the value 1 where the test ausfac > 10000 is true and NODATA where the test is false. You can verify this by displaying the streams and using cellvalue:

Grid: clear

Grid: gridpaint str10

Grid: cellvalue str10 *

and 9 to quit as usual. Use the Pan/Zoom to get a closer look if necessary.

Lets repeat the above exercise for flow accumulations of 1000 cells (10 km2) to see how the stream network grows as we lower the threshold flow accumulation.

Grid: str1 = con ( ausfac > 1000, 1 )

Running... 100%

Grid: clear

Grid: gridpaint str10

Grid: gridpaint str1

You can see how the stream network extends.

Use the cellvalue command and the Grid cursor to query individual locations on the grid:

Lets convert the grids to coverages so we can use them later:

Grid: covstr10 = gridline ( str10 )

Grid: covstr1 = gridline ( str1 )

To view the streams thus created against a backdrop of the original dem:

Grid: mape ausalb

Grid: gridpaint ausalb value linear nowrap gray

Grid: linecolor 5

Grid: arcs covstr1

Grid: linecolor 4

Grid: arcs covstr10

You can see the main drainage network in dark blue and the smaller streams in light blue, both lying the DEM valleys. Pretty neat!!

You can check the number of cells in the stream network by using

Grid: list str10.vat

Grid: list str1.vat

To be turned in: What is the number of cells in the 10,000 cell threshold stream network and the 1000 cell threshold network?

Now, we're going to separate each reach of the stream network into a different zone of cells:

Grid: lnk10 = streamlink ( str10, ausfdr )

Grid: gridpaint lnk10

Grid: cellvalue lnk10 *

Notice how each reach has a different color and a different number. All cells in this reach or "zone" have the same number.

9. Watershed delineation.

The first thing that we need to find on the watershed delineation is the location of the outlet cell at the bottom end of each watershed. We'll delineate watersheds from each intersection of stream reaches. These outlet cells are located as follows:

Grid: acc10 = zonalmax ( lnk10, ausfac ) [puts the max flow accumulation in each zone as the value of all cells in that zone]

Grid: out10 = con ( acc10 == ausfac, lnk10 ) [picks the outlet cells as those for which the flow accumulation is maximum in each reach i.e the downstream end]

Note the == operator here. This is "really, really" equals!

Now, we'll find the watershed draining to this point using the watershed function which uses the flowdirection grid and the outlet grid as inputs:

Grid: shd10 = watershed ( ausfdr, out10 )

Delineating watershed...

Now, lets take a look at the watershed and its stream network

Grid: clear

Grid: gridpaint shd10

Grid: linecolor 10 [changes the linecolor to #10 (black) for arc display]

Grid: arcs covstr10

You are actually viewing an arc coverage on top of a grid here. It is necessary to do this because if you try to draw a grid of stream on top of a grid of the watershed, the grid of the watershed is obscured.

Lets determine the number of cells in each watershed:

Grid: list shd10.vat

To be turned in: A listing of the number of cells in each subwatershed . What is the percentage of the total watershed area occupied by each subwatershed?

Grid: list lnk10.vat

Grid: list out10.vat

Notice how there are the same number of subwatersheds as there are stream links and each watershed has an outlet at the downstream end of its link.

We'll convert the watersheds to a polygon coverage so they can be displayed easily in Arcview:

Grid: covshd10 = gridpoly ( shd10 )

Grid: list covshd10.pat

The polygon attribute table (pat) has an attribute grid-code for the zone number of grid which produced each polygon. Check that the polygons are consistent with the grids that produced them:

Grid: arcs covshd10

10. Cleaning Up

In this exercise you have created 12 grids, 3 coverages and a file. Within grid you can get a listing of all the active grids by using the listgrids command. Delete all of these grids from your workspace using the kill command.

Grid: kill ausfil all

Grid: kill ausdem all

Grid: kill ausalb all

Grid: kill albfil all

Grid: kill ausfac all

Grid: kill ausfdr all

Grid: kill str1 all

Grid: kill str10 all

Grid: kill lnk10 all

Grid: kill acc10 all

Grid: kill out10 all

Grid: kill shd10 all

Type quit to exit grid. Use the listcoverages command to list the active coverages:

Arc: listcoverages [covstr1, covstr10 and covshed10 --- we'll use these in Arcview so don't kill them yet!]

and type quit again to exit arc. At the $ prompt remove the file you copied:

$ rm ausfil.e00

11. Making a Watershed and Stream Map in Arcview

at the $ prompt, fire up arcview

$arcview & [the & returns the $ prompt to you in the Decterm window so you can do other things while Arcview is on display if you want to]

Open a new view, add covstr10, covstr1, and covshd10 to it as new themes. Color them in appropriately using the legend editor. To label the map with numbers on each watershed, highlight the watershed theme in the View legend by clicking on it, go to the Theme button on the top and scroll down to Properties, choose Labels (2nd one down), and with the button at the top right hand side, scroll to covshd-10id as the item to label the watersheds by, then go to the Theme button and scroll down to Autolabel. If you want text of a different size, highlight by clicking on the View screen and under Window get the Show Symbol Palette tool and choose the button with letters on it to get access to a tool for resizing the lettering.

With View highlighted, hit the open tables button on the 2nd row of the scroll bar to get the "attributes of covshd10" table. With View highlighted, Scroll down from the View button on the top scroll bar to Layout, create a layout with the view in it and use the tool on the right hand end of the bottom row of icons to create a box for a table in the layout. However your table is displayed in the table view is the way it will show up in the Layout, so manipulate it so that it looks good.

Print the Layout and save the view so you can retrieve it and correct it if your printout is not ok.

To Be Turned In: a copy of your watershed and stream network.

Quit your workstation. Whew, its done!!


These materials may be used for study, research, and education, but please credit the authors and the Center for Research in Water Resources, The University of Texas at Austin. All commercial rights reserved. Copyright 1997 Center for Research in Water Resources.