CE 547 - GIS in Water Resources Engineering

In this assignment, you will download DEMs from the RGIS site and mosaic the DEMs.  You’ll then choose a smaller area of the DEM for which to create watersheds.  All of the work is done in the ArcView licensing level of ArcGIS with the Spatial Analyst Extension.

Go to http://rgis.unm.edu (Note that there is a link to this webpage from the class webpage.  To get past the flash screen, press Enter.) Select along the left hand side of the screen.    Note that many of the folders are expandable if you click directly on the folder.  Click the Elevation folder to expand it and select 15' Quad DEM.

You may want to check out some of the other data RGIS has to offer.  If you have trouble with data that refers to quad names, click the Reference Map along the top bar.  Also note that on the first page for some data sources, you can do a text search.

Create a workspace and make the data usable.

Create a directory for this assignment. Remember that it will be used as an ArcInfo workspace.  You’ll be working with grids which will have information in an INFO folder. The directory structure is imperative.

 The files you downloaded from RGIS are ArcInfo Export files that have been zipped.  You can use a WinZip program or the operating system to unzip the file.  To use the operating system, highlight the file to be unzipped and go to File > Compressed (zipped) Folders.  You need to unzip all three files and place them in your working directory (or workspace).

 
Once all three files are un-zipped, we need to import the DEMs, using ArcToolbox

In the ArcToolbox, navigate to Coverage Tools > Conversion > To Coverage > Import From Interchange File.  In ArcGIS 9, the capability of importing to coverages has been limited to an ArcInfo license.  (However, the ArcView 8x toolbar can be added to ArcCatalog to allow one with an ArcView license to import interchange files.  Simply, click the View menu, point to Toolbars, and click ArcView 8x Tools.)  Either way, ArcGIS will give you an error if there are any spaces in the pathname.  Make sure that none of the directories in the path have spaces.

We’re using DEMs (which are grids, not coverages) so the To Coverage is not necessarily the tool you’d expect to use but grids and coverages have the same type of file structure.  Once you click Import From Interchange File, you will get the following dialogue box requesting the name of the input file and the output dataset.  It should be fine to leave the Feature Type to AUTO.  If you click the Show Help button you will see information for each item to which you are responding. Use the Browse buttons to navigate to the correct Input file. 

Then use the Browse button to name the Output dataset.  Click OK and you should see a box similar to that shown below. 

Once you save, you can look in a Windows explorer window to confirm that a grid (and an info) folder have been created.

Repeat for all three files.  Your workspace will contain 4 folders:  e177, e178, e179, and info.

Begin viewing the data

Start using ArcMap with a New Empty Map.

 
Press the Add Data button to bring in the DEMs.

  

 

 

 

Use the shift key to select all three DEMs.  By default, they will appear in black and white.  You should recognize the shape of the Sandia Mountains in e179.  Each DEM has a different legend, so the entire DEM does not appear to be continuous.

We must use the Spatial Analyst Extension to mosaic the DEMs.  Make the Extension active by going to Tools > Extensions.

Make sure the Spatial Analyst is checked.  Next we need to get the Spatial Analyst Toolbar by going to Tools > Customize.

  Alternatively, you can access toolbars by right-clicking in the empty space at the top of the program.

Once you’ve checked the Spatial Analyst toolbar, it will appear somewhere in ArcMap.  You can move toolbars wherever they are most convenient for you.

 Typically, your first step in the Spatial Analyst will be to set the analysis extent.  Press the Spatial Analyst pulldown menu on the Spatial Analyst toolbar.  Select Options.

There are three tabs, all of which you need to access.  Type in your working directory in the General tab.

Set your analysis extent to be the Same as Display.

Set Analysis cell size to be Maximum of Inputs.

Click OK.  Now we are ready to use the Spatial Analyst Raster Calculator which is accessed from the Spatial Analyst pulldown menu.

We will use Grid functions in the Raster Calculator.  Type the expression seen in the figure below. It may be useful to open the  Help menu when typing in raster functions so that you better understand what they mean.

 

We’ll use the Raster Calculator with some grid functions to fill in the NO DATA cells.  Type the expression seen in the figure below.  When entering expressions into the raster calculator, syntax will be very important. Alternatively, you could use a series of calculation rather than this nested statement. 

Now our DEM looks really good!  You may want to change the way you look at it using the Symbology tab of the Layer Properties.  (more details in Chapter 5 of your text)

We’d like this continuous DEM with NO DATA cells to be a grid instead of a “calcuation” so right-click on the Calculation2 in the Table of Contents and select Data >Make Permanent.  Add your new permanent data set that you saved as an ESRI GRID.  I named the permanent grid Albuq.

We really don’t need any of the other grids or calculations so they can be removed by right-clicking on their names, then clicking Remove. 

 Click the Add Data button and add your continuous DEM (Albuq).  We’ll then use a number of functions in the Raster Calculator again.  If you'd like to use Fill from the ArcToolbox, you should do it at this time.  By filling the DEM, your stream network will be more continuous.

 Under options make extent and cell size “the same as Albuq”

Next you will create several calculations using the Raster Calculator.  I’ve shown the expressions below rather than the entire Raster Calculator window each time.  If you have removed all of your previous calculations, the program will start renaming the calculations as shown below.

 If you create all of these calculations as quickly as possible, you will probably have no idea what you are doing.  I recommend that you look up these grid functions in the ArcGIS help menu as you go.

 Let’s look at a few of the calculation layers you created.  Bring your flow direction calculation to the top.  Edit the symbology such that just the eight major pour point directions are shown.

Now when you view the layer you should be able to note prominent flow directions in different parts of the City.  Now make the calculation with the flow accumulation visible by bringing it to the top.  You should be able to see a general stream network.  The symbology was changed to have 5 classifications with natural breaks for the view below.  There are places where the stream network is not continuous, this is because of pits in the DEM.  In class, we will discuss filling DEMs and altering DEMs to obtain a better drainage network.

After the flow direction and flow accumulation calculations are created, we will decrease our area of analysis.  In the Data Frame, zoom in on an area.  I chose the Rio Puerco because the stream network was more continuous in this area.  With the Data Frame showing a smaller area we can change the analysis extent to just what is visible.  From the Spatial Analyst pulldown menu, select Options to change the Analysis Extent. 

Your next calculation creates a stream network. Each cell is approximately 60x60 m = 3,600 m² = 0.0036 km² in area, so 278 cells = 1 km² in drainage area.

After creating the stream network, it may be useful to convert this raster into a line coverage.  Use the Spatial Analyst pulldown menu and select Convert > Raster to Features.

Make sure you change the Output geometry type to Polyline.

Proceed with the other calculations as shown in the above table for your analysis area that is displayed.  Remember to work with the symbology for your layers to show what you want.  The layout shown below was created following the steps of this assignment.  Create a final layout showing your streams, watersheds, and area of analysis.

Delineation of Watersheds and Streams using GRID

1. Import the DEM's

Open Arc and create a workspace for this assignment.

Arc: workspace J:

WARNING: New location is not a workspace.

Arc: workspace

Current location: j:\

Arc: createworkspace hw04

Copy the 'Arc export' files (e177.e00, e178.e00, e179.e00) into your workspace. They are available in the cegis account on aix.unm.edu. From Arc, you can confirm the files are in your workspace using &system which allows you to use operating system commands.

Arc: workspace hw04

Arc: &system dir

Volume in drive J is GIS workspace

Volume Serial Number is A420-F71D

Directory of j:\hw04

02/10/99 10:41a <DIR> .

02/10/99 10:41a <DIR> ..

08/05/96 11:35p 1,212,081 E177.e00

08/05/96 11:35p 1,208,262 E178.e00

08/05/96 11:35p 1,204,661 E179.e00

02/10/99 10:40a <DIR> info

6 File(s) 3,625,004 bytes

557,891,584 bytes free

Files that end with extension e00 (i.e. *.e00) are Arc export files. These files cannot readily be used in Arc; however the files take a lot less space and there is no directory structure to maintain. Use the import command at the arc prompt to make the export file a grid. (You will need to repeat the import for all 3 export files.) Alternatively, you can use the Import71 program that comes with ArcView.

Arc: import grid e177 e177

Importing e177 from interchange file e177.e00...

Arc: import grid e178 e178

Importing e178 from interchange file e178.e00...

Arc: import grid e179 e179

Importing e179 from interchange file e179.e00...

You should now have 3 grids in your workspace. Each grid is a 2-arc second DEM from RGIS Clearinghouse Resource Data. The USGS makes 3-arc second DEM's available for the entire country and you can download them from the internet. (Keep this option in mind for your class project.)

Access the GRID module by typing grid at the Arc prompt. You can then use the following successive commands to view one of the DEM's.

Arc: grid

Copyright (C) 1982-1998 Environmental Systems Research Institute, Inc.

All rights reserved.

GRID Version 7.2.1 (Thu Apr 2 15:59:38 PST 1998)

Grid: display 9999

Grid: mapextent e177

Grid: gridpaint e177

(The coloring is random.)



Grid: clear

Grid: gridpaint e177 value linear nowrap gray

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



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

Grid: cellvalue e177 *

<9 to END>

The cell containing point (-106.855,35.219) has value 1851.000

The cell containing point (-106.895,35.144) has value 1713.000

The cell containing point (-106.937,35.060) has value 1613.000

 

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 may be interested in a region that extends across two different DEM's. You can combine the grids into one grid using Mosaic. View the combined grid as before.

Grid: albdem = mosaic(e177, e178, e179)

Mosaicing grids ... 100%

Grid: mape albdem

Grid: gridpaint albdem value linear nowrap gray



You are looking at the Albuquerque area with the Rio Grande in the darkest gray and the Sandia mountains in white.

Note that there may be a dashed line where two of your grids do not match properly. This is a problem with the data source. The black dashes have a value of NO DATA.

You can get some descriptive statistics about the grid by typing

Grid: describe albdem

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?

2. Select a region for Analysis

Because GRID operations can take a lot of time (and a lot of space) you should choose a smaller region for analysis. With the albdem grid displayed on the screen, we can set a window that will define our region of analysis.

Grid: setwindow *

Define the box

 

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. (Do not choose a window which includes a portion where the grids did not match properly.)

Grid: arroydem = albdem

Running... 100%

Grid: clear

Grid: mape arroydem

Grid: gridpaint arroydem value linear nowrap gray



 

Since you have zoomed in on a potion of the original grid, your new grid will appear more "rastery". Describe the new grid.

Grid: describe arroydem

3. Project 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. We will, however, used a different central meridian.

When typing grid functions, always leave a space around operators.

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

Grid: arroyalb = project(arroydem, #, #, 60)

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

* The INPUT projection has been defined. *

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

Use OUTPUT to define the output projection and END

to finish.

Project: output

Project: projection albers

Project: units meter

Project: parameters

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

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

central meridian [ 0 0 0.000 ]: -106 0 0

latitude of projection's origin [ 0 0 0.000 ]: 23 0 0

false easting (meters) [ 0.00000 ]: 0

false northing (meters) [ 0.00000 ]: 0

Project: end

Inversibility is 100 per cent

Project...

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

Grid: listgrids

This will create a grid in an Albers projection with a cellsize of 60 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.

Grid: describe arroyalb

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?

4. Fill 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 arroyalb

Grid: fill arroyalb arroyfil 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. Describe the grid and see if it has the statistics you expect.

Grid: describe arroyfil

5. Create Flowdirection and Flowaccumulation Grids

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

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

Grid: arroyfdr = flowdirection(arroyfil)

Computing flow direction...

You can view predominant flow directions by looking at the grid.

Grid: clear

Grid: mape arroyfdr

Grid: gridpaint arroyfdr

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

You can also check the individual values of the flow direction grid.

Grid: cellvalue arroyfdr *

<9 to END>

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 arroyfdr.vat

Record VALUE COUNT

1 1 2041

2 2 1103

3 4 1848

4 8 3908

5 16 9649

6 32 2925

7 64 1532

8 128 630

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

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

Next, determine the flowaccumulation grid and view it. This grid counts all the cells upstream of a given cell.

Grid: arroyfac = flowaccumulation(arroyfdr)

Computing flow accumulation...

Grid: clear

Grid: gridpaint arroyfac value linear nowrap gray

From the picture you should see the approximate locations of the streams. Use the cellvalue command to query the flow accumulation. You'll see big numbers on the streams (sometimes it's hard to click on them just right) and lower values on the land surface. Cells with flowaccumulation of zero are on ridge lines.

Grid: cellvalue arroyfac *

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.

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

 

6. Delineate the Stream Network

Next create a grid of the streams. You will get a different result depending on the threshold value of flow accumulation used. Each cell is 60x60 m = 3,600 m² = 0.0036 km² in area, so 278 cells = 1 km² in drainage area, and 2,780 cells = 10 km² in area. Lets delineate the streams for a 1 km² and a 10km² flow accumulation threshold (this means the minimum drainage area upstream necessary to define a stream).

Grid: str1km = con(arroyfac > 278, 1)

Running... 100%

Grid: str10km = con(arroyfac > 2780, 1)

Running... 100%

The 1 means a stream is labeled with 1.

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

You can see how the stream network extends. Use the cellvalue command and the Grid cursor to query individual locations on the grid.

Grid: gridpaint str10km

Grid: gridpaint str1km

Convert the grids to coverages so we can use them later.

Grid: covstr10 = gridline(str10km)

Grid: covstr1 = gridline(str1km)

Now view everything to see what you've done.

Grid: mape arroyalb

Grid: gridpaint arroyalb 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 in the DEM valleys. You can check the number of cells in the stream network .

Grid: list str10km.vat

Grid: list str1km.vat

 

To be turned in: What is the number of cells in the 2780 cell threshold stream network and the 278 cell threshold network?

Separate each reach of the stream network into a different zone of cells and view the results.

Grid: lnk10 = streamlink(str10km, arroyfdr)

Labeling stream links...

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.

 

7. Delineate the Watershed

The first thing that we need to find on the watershed delineation is the location of the outlet cell at the lowest point of each watershed. We'll delineate watersheds from each intersection of stream reaches.

Put the maximum flow accumulation in each zone as the value of all cells in that zone.

Grid: acc10 = zonalmax(lnk10, arroyfac)

Percentage of Cells Processed: 100%

Pick the outlet cells as those for which the flow accumulation is maximum in each reach i.e the downstream end.

Grid: out10 = con(acc10 == arroyfac, lnk10)

Running... 100%

Note the == operator here. This operator is used within the argument of the grid function.

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(arroyfdr, out10)

Delineating watershed...

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

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?

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.

Look at the watershed and its stream network.

Grid: clear

Grid: gridpaint shd10

Grid: linecolor 10

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.

8. Make a Watershed and Stream Map in Arcview

Be creative.

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

 

9. Clean Up

In this exercise you have created numerous 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 e177 all

Grid: kill e178 all

Grid: kill e178 all

Grid: kill arroyfil all

Grid: kill albdem all

Grid: kill arroyalb all

Grid: kill arroyfac all

Grid: kill arroyfdr all

Grid: kill str1km all

Grid: kill str10km 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. Use "kill cover all" to get rid of these coverages. Type quit again to exit arc.

 Adapted by Coonrod from an assignment by Maidment.