GIS Fundamentalslesson10: Introduction to Raster Analysis

GIS FundamentalsLesson10: Introduction to Raster Analysis

Lesson 10: Raster Analyses

What You’ll Learn:Spatial analysis and modeling with raster data. You will estimate the access costs for all points on a landscape, based on slope and distance to roads. You’ll then apply a threshold to this access cost.

There is a mix of old and new functions used in this lab. We’ll explain the new, but you are expected to review old labs if needed.

You should read chapter 10 in the GIS Fundamentals textbook before performing this lab.

Dataare located in the \L10 subdirectory, including Valley3 and Valley9, three and 9 meter DEMS of a portion of southeastern Minnesota, mar_rd83.shp, a vector road layer, andmardem, a 30 meter resolution raster elevation grid, and Shasta, a 30m DEM for part of northern California. All data are in NAD83 UTM coordinates, Z units in meters, the Minnesota files in zone 15, and the California files in zone 11.

What You’ll Produce: Three maps. First you’ll produce mosaic of the Valley DEMs. Then you’ll fix the Shasta DEM, and produce the first map, a hillshade of Mt. Shasta, a California volcano, after removing spikes and pits. Then you’ll create the second map, a cost surface with an applied threshold, and a third map, some of the intermediate data used for the cost surface.

Background

Raster analysis is commonly applied when working with continuous data, e.g. elevation, slope, for areas of interest. We often have to evaluate and improve the DEM data before we perform our analysis, by fixing, combining, or updating data sets. In the first project we will apply a few simple tools for DEM conditioning, and produce one map.

In the second part of this exercise we will calculate an access cost surface based on raster and vector data layers. This is a highly simplified example, but introduces basic tools that are useful in a range of problems.

Project 1: DEM Filtering, and Updating

We often obtain DEM data with various errors, a cell size different that we wish to have, or at different resolutions for different parts of our study areas. We may use raster operations, or functions, to improve our DEM data.

In our first project, we will mosaic data of two different resolutions. We often have data from various sources, for example, DEMs at approximately 10 m resolution have been developed for nearly all of the lower 48 US states. In addition, higher resolution DEMs are currently under production, typically 1 to 3 m, for many portions of the US. We may have part of our study area in each of these zones.

Our first task here will be combining two DEMs, valley3, at 3 meter resolution (cell size), and valley9, with a 9 meter resolution. We want to use the higher resolution data where we have it, but use the lower resolution data elsewhere.

(Video: Resample)

Start ArcMap and add both data sets.

Your view should be similar to that at the right.

Calculate the hillshade for both data sets (ArcToolboxSpatial Analyst Tools – Surface – Hillshade).

Inspect these hillshades carefully, and

note the enhanced detail with the 3meter DEM, as shown below. The figure on the right is from the hillshade of valley3, on the left thevalley9. Note the greater definition of the small streams and streambanks.

Remove the two hillshades from the data view to reduce clutter.

We can use the raster calculator to join these two data sets. Note that we want to use the detailed valley3 data where we have it, and the courser valley9 data everywhere else.

To avoid confusing results we need to only perform raster operations on similar structured rasters. Therefore our first step is to convert the 9 meter data to a 3 meter cell size. This does not make the 9 meter data better. Basically, we are making a copy of the data at a finer resolution.

The process to convert 9 meter to 3 meter cell size is resample.

Use ArcToolBox – Data Management Tools – Raster – Raster Processing - Resample

In the resultant window, specify the valley9 as input, an output file, a 3 meter cell size, and a bilinear resampling (the textbook describes the differences among resampling methods). Name the output raster valley9to_3.

After the resampling is done, examine the valley9to_3 and verify that it has a 3 meter resolution ( Properties - Source).

We can now combine the two data sets.

Although there are many ways to do this, perhaps one of the simplest is with clever use of two raster functions – IsNull and Con. The textbook describes these in detail, but in short, the IsNull returns True whenever a cell comparison or value is Null. The con function takes three values, the first is a true/false test, and second is the value to assign to a grid if the test is true, and the second is the value to assign if the test is false.

Using the Raster calculator, we can nest these functions:

ArcToolBox  Spatial Analyst Tools  Map Algebra  Raster Calculator

Con(IsNull("valley3"),"valley9to_3","valley3")

Name the output combinedDEM

Note: Remember when using the Raster Calculator ALWAYS use your mouse to select the layers and operators. ONLY type in data and values when you have no other choice. The reason for this is that the spacing of the Raster Calculator is very sensitive. If you type in the correct information it is VERY possible it will not work due to small spacing errors.

This function first applies the test, is the cell of Valley3 equal to Null? This is true everywhere in our study area outside the data region of Valley3. When the value is Null, the con function assigns the value found in Valley9_to3 to the cell in the output data set. When the value is not Null, the con function assigns the value found in Valley3.

Examine the combineDEM raster. Compute and inspect a hillshade, and verify that it has the higher detail contributed by the Valley3 data set.

Add a new empty Data Frame.

Remove the Data Frame that contains the valley3,..etc.

We’d now like to introduce filtering as a tool to fix “noisy” data. This is often used with interpolated surfaces, particularly LiDAR data, and similar tools are used near edges for mosaiced DEMs and other continuous surfaces. Filters are described in Chapter 10 of the GIS Fundamentals textbook.

Load the DEM named Shasta. Calculate the slope (ArcToolboxSpatial Analyst Tools – SurfaceSlope). Pick either percent or degrees with a z factor of 1.

You should get an output surface similar to that on the right.

Examine the slope surface. Notice the “speckled” appearance in the northwest and southeast, flatter portions of the DEM.

Now calculate a hillshade surface for the Shasta DEM. Leave the Azimuth as is, set the Altitude to 25, and check the model shadows box.

Again, inspect this, and notice the funny artifacts, in the locations of the anomalous slopes. These are both data “spikes,” white points with a long thin shadow trailing to the southeast, or “pits,” dark areas on the northwest with a white “edge” on the southeast.

How do we remove these?

As described in the text, we may use a low-pass filter to identify and get rid of this speckle.

We may apply a low-pass filter with the ArcToolBox – Spatial Analyst Tools –Neighborhood – Filter command. (see Video Lowpass)

Specify the Shasta DEM as the input, and a filter type of LOW, and specify an output name, like lowpass.

Run the filter, and inspect the output. It is probably easiest to see the effects by calculating the hillshade of the output, and looking at the areas that have spikesand pits. The figure at right, corresponding the area above, shows the reduction in the size of the spikes and pits, although they are still visible.

We could stop here, and just accept the filtered data layer. But if we look carefully at the filtered and unfiltered hillshade, we’ll see we pay a cost for filtering. We lose some of the fine detail, apparent in the image of the unfiltered hillshade to the left, below, compared to the filtered hillshade to the right.

We’d like to keep both, and we can. Video: Calculator

First we should subtract the filtered layer from the original Shasta layer.

ArcToolBoxSpatial Analyst ToolsMap AlgebraRaster Calculator

Select the following in the Raster Calculator "shasta" - "lowpass"

Name the output difference

We can then only replace the cells were the difference is large. Here, large is a relative term, but after a few trials, and looking at the difference histogram, a threshold of about 15 works fairly well.

We want to replace the cells in the Shasta image when they are more than 15 meters different from the filtered surface. Otherwise, we leave the Shasta surface alone, and hence, we don’t get any of the degradation in detail in otherwise good data.

We can open the Spatial Analyst Raster Calculator, and apply the following function, shown to the right. Name the output smoothed.

If the absolute value of the difference is greater than 15, we write the filtered value to the output. Otherwise, we write the original data value.

We can verify that this is helpful by viewing the hillshade of the output raster, and comparing it to the original, and the filtered surfaces. Note that we have removed most of the speckle, but maintained the detail.

We could apply the filter successively, and average the local points further, and apply the con function to a difference layer again. It would lead to an improved surface, and we would do this if we were using the data in the project. But since you won’t learn much new with repetition here, we’ll just produce a map and move on to project 2.

Calculate a hillshade of the smoothed raster, name it HillSmooth.

Use an altitude of 25 degrees and model shadows.

Add the shaded relief (that is the Hillsmooth), making it semi-transparent (50%) over the smoothed elevation data. Make sure you model shadows in the shaded relief surface (points will be deducted if you do not – we need it to evaluate how well you’ve applied the filtering).

Create a map as shown below, add a reasonable legend, name, north arrow, description, and scale bar, and export the map as a .pdf.

Project 2: A Cost Surface

We will create a cost surface for locating a building. Our cost surface will depend on slope and distance to existing roads. In our problem, we will assign a road construction cost of $25 per meter of road required. We have a vector data layer of roads, digitized from USGS maps, and we will use grid functions to convert this to a cost data layer.

Slope also affects access costs, because roads on steeper terrain are more expensive. The cost is nonlinear, increasing slowly at first for low slopes, then more rapidly at steeper slopes. We will derive slope from a DEM data layer. We will modify the tables associated with both the derived slope and distance layers to include a cost column. To reflect the nonlinearity in slope costs, we will apply the trigonometric sine function to model this increase in cost. We will then add these two cost layers. Finally, we wish to apply an upper threshold of $5,000 to consider only those areas that are within our budget.

Before we start the second project, we need to describe a difference between a permanent reclassification you’ll be doing today, and a display reclass you’ve done before and you’ll also do today.

Remember, a reclassification is a conversion from one set of numbers to another. We do this in a raster GIS through a reclass table. This table has a column for input values (Old values in the figure at right) and a column for output values (New values in the figure at right). Each cell value is examined, and input value matched to an entry in the table, and the corresponding output value reassigned according to the table. For example, the table at right specifies that all Old values between 3.385998 and 7.336329 are assigned a new value of 2.

In a permanent reclassification, each output value is saved to a new raster. In a display reclassification, the value is used only to assign symbols for display. No data are changed in the source file, nor are new files saved. In previous lessons we have only performed reclassifications for display. Today we will perform a permanent reclassification. It is easy to get confused, because the classify menus for applying these classifications are similar.

Project 1: Cost Surface

Start ArcGIS - ArcMAP

• Create a new map project, add the raster mardem to the view, and inspect it.

Use the cursor and the layer Properties –Source tab. What are the units of the elevation? What are the highest and lowest elevation values? Does it make sense?

• Derive the slope for mardem. Select ArcToolboxSpatial Analyst ToolsSurface AnalysisSlope. Specifydegrees units for slope. Name the output filemar_slope. Video: Slope & Reclass

• To keep the view uncluttered, remove the mardem grid from the map.

• Examine the slope layer. They should be values from 0 to about 33 degrees.

• Select ArcToolboxSpatial Analyst ToolsReclassReclassify. You’ll get a popup menu, with a reclassification table, similar to the table in the figure on the previous page.

• At this popup menu, click on the Classify button. This will open a classification window that is exactly equal to the window you see when changing the symbology for a layer. However, in this case you are using to change the assignment or classification table.

• Here,

• Select aDefined Interval classification with an interval size of 1.
• Left click on OK

Note this returns you back to the Reclassify menu, but that the reclassification table has been changed. Now the Old values to New values list should reflect the reclassification you specified, as illustrated in the figure below.

Also note that you can save this reclassification table, for example, if you wish to use it again in the future, and you can load saved tables.

Also note you can specify how missing data are assigned, and you must specify the name and the directory of your output file

Name the output grid something likeSlpcls.

During the subsequent steps we walk you through you may wonder, why work on all these derived layers, why don’t we just work on the slope layer? You could easily skip several of these steps but the first time it is best that you see the details of each step.

• Now remove the original slope layer, it just clutters your display.

Next we will apply a formula that determines the cost of building on slopes.

• Left click on ArcToolboxSpatial Analyst ToolsMap AlgebraRaster Calculator (Video: Calculator)

• Use the mouse to select the following operators and values for the center window: Sin(“Slpcls”/57.2958) * 200.

Verify that the cost layer makes sense, and that they are highest where slopes are steepest.

Next, we need to display and generate our distance costs from the roads

layer.

Add the Mar_rd83.shp file.

Make sure Slope_Costis the target Layer for the Spatial Analyst.

Then leftclickArcToolboxSpatial Analyst ToolsDistance Euclidean Distance.

• Name the output rasterDistance.

• Examine the result layer, and make sure it is reasonable.

Now, multiply the distance layer you just produced by the cost per unit distance to estimate distance cost.

• Left click ArcToolboxSpatial Analyst ToolsMap AlgebraRaster Calculator and enter the equation as shown below. Name the output raster Dist_Cost

Our next step is to combine the two sets of costs. Open the Raster map calculator again (ArcToolboxSpatial Analyst ToolsMap AlgebraRaster Calculator), and add the two cost layers, as below:

Name the output Total_Cost.

Examine the Total_Cost layer and make sure it makes sense.

Think a minute about what you’ve just done. You first calculated a slope, and then a cost associated with building a road per unit distance across the slope. Then you calculated a distance, and then a cost associated with building a road to that distance from an existing road. Both of these were calculated for every grid cell in your study area. You then added these two together for an estimated total cost to build a road to any portion of the mapped area. A real problem would include many other factors, like soils, surface vegetation, slope constraints over minimum segments, etc., but this would only lengthen the analysis, and not change the basic way you are applying the tools.

Our job now is to select those areas below the $5,000 threshold. We will do this by creating a mask grid. This grid will have 1 at all locations where the costs are below $5000 and 0 where the costs are above $5,000. We will then multiply this with our total cost grid, to zero out those areas we don’t wish to consider.

Reclassify Total_Cost by ArcToolboxSpatial Analyst ToolsReclassifyClassify.

In the reclassification window (not shown), set the

Method to Equal Interval,

Classes to 2,

Type in 5000 in place of the firstBreak Value,

Leave the 2ndBreak Value at the current (maximum) level.

Left click OK

This should result in a reclassification table as shown below.

Make sure you have NoData for the New value of the 5000 to 24753.49024 category.