Using Spatial Analyst to work with DEMs

Working with Digital Elevation Models in ArcGIS 8.3

This lab continues your introduction to Raster data, using the Spatial Analyst Extension and Toolbar. This is the primary tool that is available for working with Raster data (Grids) in ArcGIS. In the first lab you created contour lines and hillshading using the Spatial Analyst extension. These are two fundamental tools to help display elevation data. Many people are familiar with contour lines - which are a specific form of an isoline. An isoline is a line of equal value. In the case of a contour map, the lines are of equal elevation. You might have lines of equal rainfall, which may be an easy way to view rainfall data. Rainfall data is another continous variable that is best represented for modeling in raster format, although for visualization, isolines are easily understandable. We will work with DEMs in this lab.

Contouring

Let's look more closely at the contouring process. We will use the Digital Elevation Models (DEM) that you downloaded from the USGS website in the previous lab.

Open ArcMap, start a new project, and add the projected NED data (DEM) that you created.

When you made the contour lines, you saw a dialog like the one below. (To find the Contour window, click on Spatial Anaylst in its toolbar, then click on Surface Analysis then click on Contour.)

You may have changed the contour interval or you may have run the command (to create contours) as it was presented to you.

Create contours with the DEM that you downloaded. (not the land clasification data that you also downloaded). If you downloaded a very large area, you may see that the contours are quite dense and unreadable, as compared to the contours on the rainfall maps above. That is because the default value that is given is based on the elevation range in the data, not by what would be readable. You should try to find a more reasonable contour interval for the area you are looking at. Try taking the suggested interval and multiplying it by 5. That should give a more reasonable interval. If not, try other values until you get a set of contours that appear reasonable and readable.

A trick that you might want to try: sometimes you are interested in a single contour line, like zero elevation when your data includes a coastline. This way you can find the coastline, as defined in your elevation data. To do this, set the base contour as the level you want (0 in the case where you are trying to find the coastline) and then change your contour interval to a number larger than the highest elevation in your DEM. Check your map legend to see what your highest elevation is.

Remember that the contouring process assumes that the single elevation represented in each cell is the elevation of the center point of the cell. In this case, the DEM is treated as a lattice, with data only at the midpoint of cells, not data representing the entire cell.

Contours can be of any type of data, including rainfall data. The command's name in ArcGIS is narrowly focused but its use is not.

Setting the working window

Sometimes you want to work with only part of your DEM. You may want to change the cell size because you are working with another dataset that has a larger cell size. You may want to do this today, because the data you retrieved from the USGS Seamless data server is quite large. For purposes of learning, you can use a small dataset, which will speed up processing considerably, both in this lab and in the Hydrology lab in the near future.

Try this: First, zoom in to a small area. Make this area no larger than 15,000 meters on a side. You can use the measure tool to get distances (make sure that you have the units - meters in this case - set properly in your Data Frame Properties). You don't want to have too small an area so try to get it at least 10,000 meters on a side. The USGS 1:24,000 topographic maps are roughly 10,000 by 13,000 meters at 42 degrees north. Next, open up the Spatial Analyst Options window. Do this by clicking on Spatial Analyst on its toolbar. Then click on Options. Next, click on the Extent tab. You should see a window like this one:

You can change the extent by changing the analysis extent. In the case above, I changed it to "Same as Display". This means that any function you run in Spatial Analyst will yield results only for the area that you set in the Analysis Options. If you had a definite area of interest, you could enter the coordinates in the Options window directly, rather than relying on the display window. Another way to limit the area in which Spatial Analyst operates is to set an "Analysis Mask" in the General tab. Right-click on "Analysis Mask" to see a description of it's capabilities.

You can change the cell size by clicking on the Cell Size tab. You could change it by entering a new number or using another grid that you have - such as the Land Characterization data. But don't do this before completing the lab exercise.

Click on General to change the working folder to a folder on your H:\ drive.

Click OK to set the Spatial Analyst options. You may get an error indicating that "The Working Directory location can not have a space in the path name...". To fix this, go to the General tab in the options window, and select your H: drive Lab5 directory as the working directory. Or better yet, you can create a directory called "Work" or "Temp" under your 1.963 directory, and use that as the working directory (where Spatial Analyst keeps it's temporary files, and sometimes it's output files).

Now, make a new set of contour lines. (First, remove the contours you created earlier). You should see contours only for the area of display (which is what you set your extent to).

Try this with hillshading as well. To find the Hillshade window, click on Spatial Anaylst in its toolbar, then click on Surface Analysis then click on Hillshade. Zoom out to check that the results are only for the window that you selected.

Other Surface Analysis tools

You should use other tools that are included with Contouring and Hillshading.

Aspect tool.

What is aspect? From the help pages:
Aspect identifies the steepest downslope direction from each cell to its neighbors. It can be thought of as slope direction or the compass direction of a hill faces.

This is a very useful tool for the environmental engineer or planner. You can find north facing slopes where snow is likely to linger or south facing slopes that will allow greater solar heating of residences.

Use Spatial Analyst to create an Aspect map of your DEM. Look at the results. Notice that the results are numbers (in direction: 0 is north, 90 is east, 180 is south, and 270 is west). The results also groups these in the legend so that the numbers are easier to read. You will use the aspect results in your homework.

Slope tool.

What does the slope tool produce? From the help pages:
The Slope function calculates the maximum rate of change between each cell and its neighbors, for example, the steepest downhill descent for the cell (the maximum change in elevation over distance between the cell and its eight neighbors). Every cell in the output raster has a slope value. The lower the slope value, the flatter the terrain; the higher the slope value, the steeper the terrain. The output slope dataset can be calculated as percent slope or degree of slope.

Again, this is important information for the environmental engineer and planner. Create a Slope map for your DEM. Notice that the steep areas are easy to see, perhaps easier to see than with color and hillshading, and less cluttered than with contours. You will use the slope results in your homework.

Viewshed tool.

What does the viewshed tool produce? From the help pages:
Viewshed identifies the cells in an input raster that can be seen from one or more observation points or lines. Each cell in the output raster receives a value that indicates how many observer points can be seen from each location. If you have only one observer point, each cell that can see that observer point is given a value of 1. All cells that cannot see the observer point are given a value of 0. The observer points feature class can contain points or lines.

This is also a useful tool. The results can be seen in two ways - what can be seen from the observer location and what areas can see the location. The former might be useful for planning a scenic road and the latter might be useful for sighting a cell tower.

Since the Viewshed tool requires a shapefile as a starting point, you will have to edit shapefiles. Look at the addditional handout on this.. Please remember to add a projection to the new shapefile you are making! Make one new shapefile, containing a single point, somewhere in your area of interest (on your DEM). To demostrate the idea of a viewshed using a point shapefile, you should try to locate the point where you expect to see a sizable portion of the landscape. The side or top of a hill are good candidates. The bottom of a valley is not a good place.

Once you have created and saved your point shapefile, try the viewshed tool. Make sure that the "Input Surface" is set to your DEM, and not one of the other grids that you have created. Notice that you can submit one shapefile that contains observer points, which the help described as being points or lines. Try it with your new point shapefile.

Using the Raster Calculator

The raster calculator is both a calculator and query tool in ArcGIS. It can be used for finding all areas of a specific slope range or of certain aspects. Open the Raster Calculator by clicking on Spatial Analyst on it's toolbar, and selecting "Raster Calculator". The calculator looks like this:

You can use this to build expressions or to query datasets. We will use this to query the various datasets you have created. A word of warning: the syntax is best handled using the interface until you have it memorized. Rather than type in a grid name, click on it. Do the same with numbers, operators and bolean operators (And, Or, etc). The proper use of punctuation and spaces is very important, and easy to mess up if you try to type it yourself.

Find all areas of a certain elevation using this query (remember that you need to use your DEM name in place of "ashfield_nos", which you simply double click on):

[ashfield_nos] >= 200 & [ashfield_nos] < 300

The result of this query is a grid that has cells in your extent with one of two values, a 0 where the condition was not met and a 1 where the condition was met. You can use the raster calculator two calculate a new grid that has elevation values only where the above condition was met. You can do this by multiplying the original dem with the new grid. Note that the new grids are called Calculation or Calculation 2, etc. Your query might look like this:

[ashfield_nos] * [calculation]

Again, the result will be named calculation # where # is a number higher than the last calculated grid you made in the Raster Calculator.

You can also use this to run commands such as slope and aspect (but you won't get the nice legends and color scheme):

slope([ashfield_nos])

You can also use this to look at flow direction and flow accumulation. We will do so in a future lab.

To find all of the functionality that you can access from the raster calculator, find it in the on-line help. Open up help (ArcMap's Help menu >> ArcGIS Desktop Help) and go to the Contents tab. Go to Extensions >> Spatial Analyst >> Spatial Analyst Functional Reference.

Site Selection using Raster Data

For this part of the exercise, attach the 1.966 locker as the T drive (do this in the Command Prompt window, found in the Accessories, using the command: attach -DT 1.966). The data is in the rastersiteselection folder.

You will do a site selection exercise using raster processing. This involves constructing a query using the raster calculator. The base materials that you need are the slope, aspect from the DEM (dem2) and the Land Cover layer (landcover) that can be downloaded with the elevation data. You will need to rerun the commands to create the slope and aspect layers for this new DEM (which is in western Massachusetts).

You need to map out areas of southern exposure, areas of low slope, and of a certain land cover type. Imagine that you are working for an engineering firm that is assisting in finding a site for building a solar energy complex. You have certain parameters that need to be met:

slopes less than 5 degrees (for construction purposes)
aspect is generally southern, you are willing to accept any thing from 110 to 160 (SE to SW)
this is very fertile farm land, therefore you should avoid farm land. You should also avoid open water and commerical or residential development. Therefore, the best choices are areas that are forested (41, 42, and 43 in your NLCD layer - see handout or the USGS site for the complete code list).

You will need to construct a query that will yield a resulting grid which matches all of these constraints. A couple of clues:

your aspect grid will have to appear in the calculator twice, once for >= 110 and once for <= 160
your will need to connect each of the parts of the query with "And" (&)
the viewshed grid needs to be set so that the value(s) representing "Not Visible" are included.

An example might look like this, in part:

[Aspect of dem2] >= 110 & [Aspect of dem] <= 160

Note that there is not a "between" operator and if you want to pick a range, you need to have the upper and lower limit set as in the above example.

Try making a map that illustrates what you are trying to show - it should include contour lines with the land cover data. Change the text so that the land cover is a short text string rather than the obscure numbers.

Additional work

Try working with the DEM in the NeighborhoodExercise folder on the T drive (attach again, if you do this outside of the lab session) and the Neighborhood Statistics functions. You will want to find local low points. There are a couple of steps to this process:

Run the neighborhood function and search for the minimum area, using DEM as the Input data. For an initial neighborhood, try a rectangle that is 3 x 3 Cells. Your statistic type should be "Minimum". and the Neighborhood should be a rectangle.
Using the Raster Calculator, find the cells in the original DEM and in the output of the neighborhood function that are the same. This shows you all of the local low points. Be sure to use the interface to write the statement in the calculator. You will see tha the expression includes "= =" rather than a single "=".
Change the symbology of the results of the Raster Calculator so that the values of zero have no color. These are places where the values are not equal. The cells with a value of 1 is where the local low point is equal to the cell value is a local low point.
Use the orthoimages to check if these are wetlands (water appears as black in these images). Use the OrthoIndex Shapefile to find the Sheet_ID for the two ortho images in the NW (or upper left) corner of the DEM. The Ortho images can be downloaded from ftp://data.massgis.state.ma.us/pub/wave/coq2005hsid. Use IE, rather than Mozilla Firefox to download these, since you can drag and drop the files into a Windows Explorer window. Note that you need all files with the same sheet_id, not just the .sid file (sid is an compressed image format). Download the wetlands data (Massachusetts 1:12,000 data and possible and confirmed vernal pools databases from http://mass.gov/mgis/laylist.htm) to check if these have been recorded as wetlands. You will use all of the downloaded data on the field trip on 10/14. You may also want to try to download images of USGS maps for these same sheet ids. Find these on the MASSGIS data layers page.

Try this with a larger neighborhood (say 10 cells). Also try with a larger neighborhood and with a Neighborhood as a Circle, rather than the default rectangle. Are the results different?

Can you tell if the local low point is a depression or an error in the data? What neighborhood statistic would you use to help find the difference between the calculated local low point and the surrounding area? There are at least two possibilities. Try this.

Write a short (1-2 page) report on what you tried and how it worked. Detail differences, if any and why, between the results of different neighborhoods. Why would the results of a 3x3 rectangle be suspect? Were you successful in finding wetlands that were identified preiously?

Created by Daniel Sheehan 9/26/03
last updated 10/4/06 - Daniel Sheehan.