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.
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.
(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.
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
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.
Import the grid ausfil as a directory of files
Arc: import grid ausfil ausfil
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.
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]
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: projection albers
Project: datum NAD83
Project: units meters
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
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
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.
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: 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)
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: 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.
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 >]
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: 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 )
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.
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 )
Now, lets take a look at the watershed and its stream network
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
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
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.