Massachusetts Institute of Technology
Department of Urban Studies and Planning

11.520: A Workshop on Geographic Information Systems
11.188: Urban Planning and Social Science Laboratory

Lab Exercise 7: Raster Spatial Analysis

Distributed: October 27, 2008
Due: Monday, Nov. 3, 2008

Overview

The purpose of this lab exercise is to introduce spatial analysis methods using raster models of geospatial phenomena. Thus far, we have represented spatial phenomena as discrete features modeled in the GIS as points, lines, or polygons--i.e., so-called 'vector' models of geospatial features. Sometimes it is useful to think of spatial phenomena as 'fields' such as temperature, wind velocity, or elevation. The spatial variation of these 'fields' can be modeled in various ways including contour lines and raster grid cells. In this lab exercise, we shall focus on raster models and examine ArcGIS's 'Spatial Analyst' extension.

We shall use raster models to create a housing value 'surface' for Cambridge. A housing value 'surface' for Cambridge would show the high- and low-value neighborhoods much like an elevation map shows height. To create the 'surface' we will explore ArcGIS's tools for converting vector data sets into raster data sets--in particular, we will 'rasterize' the 1989 housing sales data for Cambridge and the 1990 Census data for Cambridge block groups.

The block group census data and the sales data contain relevant information about housing values, but the block group data may be too coarse and the sales data may be too sparse. One way to generate a smoother housing value surface is to interpolate the housing value at any particular location based on some combination of values observed for proximate housing sales or block groups. To experiment with such methods, we will use a so-called 'raster' data model and some of the ArcGIS Spatial Analyst's capabilities.

The computation needed to do such interpolations involve lots of proximity-dependent calculations that are much easier using a so-called 'raster' data model instead of the vector model that we have been using. Thus far, we have represented spatial features--such as Cambridge block group polygons--by the sequence of boundary points that need to be connected to enclose the border of each spatial object--for example, the contiguous collection of city blocks that make up each Census block group. A raster model would overlay a grid (of fixed cell size) over all of Cambridge and then assign a numeric value (such as the block group median housing value) to each grid cell depending upon, say, which block group contained the center of the grid cell. Depending upon the grid cell size that is chosen, such a raster model can be convenient but coarse-grained with jagged boundaries, or fine-grained but overwhelming in the number of cells that must be encoded.

In this exercise, we only have time for a few of the many types of spatial analyses that are possible using rasterize data sets. Remember that our immediate goal is to use the cmbbgrp and sales89 data to generate a housing-value 'surface' for the city of Cambridge. We'll do this by 'rasterizing' the block group and sales data and then taking advantage of the regular grid structure in the raster model so that we can easily do the computations that let us smooth out and interpolate the housing values.


I. Setting Up Your Work Environment

1. Launch the ArcGIS and add five data layers listed below.

  • M:\data\cambbgrp_point.shp
    		• Census 1990 block group' centroids for Cambridge
  • M:\data\cambbgrp.shp
    		• Census 1990 block group polygons for Cambridge
  • M:\data\cambtigr (coverage - arcs)
    		•  U.S. Census 1990 TIGER file for Cambridge
  • M:\data\sales89 (coverage - points)
    		•  Cambridge Housing Sales Data
  • M:\data\camborder (coverage - polygons)
    		•  Cambridge polygon

    2. Set Display unit = meter. In this exercise you will use "Meter" instead of using "Mile"

    II. Spatial Analyst Setup

    ArcGIS's raster manipulation tools are bundled with its Spatial Analyst extension. It's a big bundle so lets open ArcGIS's help system first to find out more about the tools. Open the ArcGIS help page by clicking Help > ArcGIS Desktop help from the menu bar. Click the index tab and type "Spatial analyst'. During the exercise, you'll find these online help pages helpful in clarifying the choices and reasoning behind a number of the steps that we will explore. Be sure, at some point, to take a look at the Overview section. (You can also open ArcGIS Help directly from the Windows operating system via: Start / Programs / ArcGIS / ArcGIS-Desktop-Help.)

    The Spatial Analyst module is an extension, so it must be loaded into ArcGIS separately.
    To load the Spatial Analyst extension:

    III. Interpolating Housing Values Using SALES89

    This part of the lab will demonstrate some techniques for filling in missing values in your data using interpolation methods. In this case, we will explore different ways to estimate housing values for Cambridge. Keep in mind that there is no perfect way to determine the value of a property.

    A city assessor's database of all properties in the city would generally be considered a good estimate of housing values because the data set is complete and maintained by an agency which has strong motivation to keep it accurate. This database does have drawbacks, though. It is updated at most every three years, people lobby for the lowest assessment possible for their property, and its values often lag behind market values by several years.

    Recent sales are another way to get at the question. On the one hand, their numbers are believable because the price should reflect an informed negotiation between a buyer and a seller that results in the 'market value' of the property being revealed (if you are a believer in the economic market-clearing model). However, the accuracy of such data sets are susceptible to short-lived boom or bust trends, not all sales are 'arms length' sales that reflect market value and, since individual houses (and lots) might be bigger or smaller than those typical of their neighborhood, individual sale prices may or may not be representative of housing prices in their neighborhood.

    Finally, the census presents us with yet another estimate of housing value--the median housing values aggregated to the block group level. This data set is vulnerable to criticism from many angles. The numbers are self-reported and only a sample of the population is asked to report. The benefit of census data is that they are cheaply available and they cover the entire country.

    We will use sales89 and cambbgrp to explore some of these ideas. Let's begin with sales89.

     Fig. 3. Interpolation without mask


    The interpolated surface is shown thematically by shading each cell dark or light depending upon whether that cell is estimated to have a lower housing value (lighter shades) or higher housing value (darker shades). Based on the parameters we set, the cell value is an inverse-distance weighted average of the 12 closest sales. Since the power factor was set to the default (2), the weights are proportional to the square of the distance. This interpolation heuristic seems reasonable, but the surface extends far beyond the Cambridge borders (all the way to the rectangular bounding box that covers Cambridge). We can prevent the interpolation from computing values outside of the Cambridge boundary by 'masking' off those cells that fall outside of Cambridge. Do this by adding a mask to the Analysis Properties as follows:
     


    With this analysis mask set, interpolate the Realprice values in sales89 once again and save it as sales89_pw2-2. The sales89_pw2-2 should look like this:
    (When using Quantile, 9 classes)
    Fig. 4. Interpolation with mask

    All the values inside Cambridge are the same as before, but the cells outside Cambridge are 'masked off'. [Note: When you 'interpolate to rater' using inverse distance weighted (IDW) in ArcMap, the thematic map display does a further smoothing of the computed values for each cell so the map in Fig. 4 looks more like a contour plot then the display of the discrete values for each grid cell - zoom in to see what I mean. But this is for display purposes only and the grid cell values saved in the attribute table are a single value for each grid cell as calculated via the IDW averaging.)
    To get some idea of how the interpolation method will affect the result, redo the interpolation (using the same mask) with the power set to 1 instead of 2. Label this surface 'Sales89_pw1'. Use the identify tool to explore the differences between the values for the point data set, sales89, and the two raster grids that you interpolated. You will notice that the grid cell has slightly different values than the realprice in sales89, even if there is only one sale falling within a grid cell. This is because the interpolation process looks at the city as a continuous value surface with the sale points being sample data that gives an insight into the local housing value.  The estimate assigned to any particular grid cell is a weighted average (with distant sales counting less) of the 12 closest sales (including any within the grid cell). In principle, this might be a better estimate of typical values for that cell than an estimate based only on the few sales that might have occurred within the cell. On your lab assignment sheet, write down the original and interpolated values for the grid cell in the upper left (NorthWest part of Cambridge) that contains the most expensive Realprice value in the original sale89 data set. Do you understand why the interpolated value using the power=1 model is considerably lower than the interpolated value using the power=2 model? There was only one sale in this cell and it is the most expensive 1989 sale in Cambridge. Averaging it with its 11 closest neighbors (all costing less) will yield a smaller number. Weighting cases by the square of the inverse-distance-from-cell (power=2) gives less weight to the neighbors and more to the expensive local sale compared with the case where the inverse distance weights are not squared (power=1).

    IV. Interpolating Housing Values Using CAMBBGRP

    Finally, let's smooth this new grid layer using the Spatial Analyst > Neighborhood Statistics option.  Let's recalculate each cell value to be the average of all the neighboring cells - in this case we'll use the 9 cells (a 3x3 matrix) in and around each cell. To do this, choose the following settings: (they are the defaults)
    Many other variations on these interpolations are possible. For example, we know that med_hvalue is zero for several block groups--presumably those around Harvard Square and MIT where campus and commercial/industrial activities results in no households residing in the block group. Perhaps we should exclude these cells from our interpolations--not only to keep the 'zero' value cells from being displayed, but also to keep them from being included in the neighborhood statistics averages. Copy and paste the cambbgrp.shp layer and use the query tools in the Layer Properties > Definition Query tab to exclude all block groups with med_hvalue = 0 (which means include all block groups with med_hvalue > 0). Now, recompute the polygon-based interpolation (including the neighborhood averaging) and call this grid layer 'hvalue_non0'.  Select the same color scheme as before. In the data window, turn off all layers except the original camborder layer (displayed in a non-grayscale color like blue) and the new hvalue_non0 layer that you just computed. The resulting view window should look something like this:
    (When using Quantile, 9 classes)
    Fig. 8. hvalue_non0

    Notice the no-data cambordergd cells sticking out from under the new surface and notice that the interpolated values don't fall off close to the no-data cells as rapidly as they did before (e.g., near Harvard Square). You'll also notice that the low-value categories begin above $100,000 rather than at 0 the way they did before. This surface is about as good an interpolation as we are going to get using the block group data.

    Comment briefly on some of the characteristics of this interpolated surface of med_hvalue compared with the ones derived from the sales89 data.  Are the hot-spots more concentrated or diffuse? Does one or another approach lead to a broader range of spatial variability?

    V. Combining Grid Layers Using the Map Calculator

    Finally, let us consider combining the interpolated housing value surfaces computed using the sales89 and med_hvalue methods. ArcGIS provides a 'Raster calculator' option that allows you to create a new grid layer based on a user-specified combination of the values of two or more grid cell layers. Let's compute the simple arithmetic average of the sale89_pw2-2 grid layer and the med_hvalue_non0 layer. Select Spatial Analyst > Raster Calculator and enter this formula:
    ([hvalue_non0] + [sale89_pw2-2]) / 2
    and click Evaluate.
    The result is a new grid which is the average of the two estimates and looks something like this:
    (When using Quantile, 9 classes)
    Fig. 9. Raster Calculation
    The map calculator is a powerful and flexible tool. For example, if you felt the sales data was more important than the census data, you could assign it a higher weight with a formula such as:
    ([hvalue_non0] * 0.7 + ( [Sales_Price] * 1.3) / 2

    The possibilities are endless--and many of them won't be too meaningful! Think about the reasons why one or another interpolation method might be misleading, inaccurate, or particularly appropriate. For example, you might want to compare the mean and standard deviation of the interpolated cell values for each method and make some normalization adjustments before combining the two estimates using a simple average. For the lab assignment, however, all you need do at this point is determine the final interpolated value (using the first map-calculator formula) for the cell containing the highest price sales89 house (in the Northwest corner of Cambridge).  Write this value on the assignment sheet.

    VI. Combining Grid Layers Using Weighted Sum in ModelBuilder

    Next, we can 'automate' some of these raster analyses by creating a simple ModelBuilder model of one or another of the above steps. We can do an identical analysis, or we can refine things somewhat using additional capabilities of the "Weighted Overlay Tool"

    Setup

    To create a new ModelBuilder model, first make sure that you know where ArcGIS is going to save your model. The default location can be found by going under the main menubar to Tools->Options, and in the resulting dialog box, looking at the "Geoprocessing" tab. The third item there allows you to specify the location of your 'My Toolboxes' folder. The default location on MIT's network is H:\WinData\Application Data\ESRI\ArcToolbox." This should allow you to access the same model from different workstations, but if you want to share (or back up!) a model, it is good to know where it is located. You can't directly "open" or "save" ModelBuilder models. However you can move them arround in ArcCatalog. In general, if you are creating models specific to a particular project, consider putting them in a directory near the project data.

    Geoprocessing Options Dialog Box

     

    To Create a New ModelBuilder Model, first create a "Toolbox" to contain it. Open the main ArcGIS toolbox (red tools icon). Then right mouse on ArcToolbox to get the contextual menu shown, and select New Toolbox.

    New Toolbox

    Once you have a container located in a writeable folder, variously called Toolbox (the default) or a name you choose (such as 'RasterAnalysis'), then you can create a new Model within it.

     

    New Model

    This will bring up and empty Model diagram window:

    Empty Model

    You should then be able to drag and drop your two grids (called camb_sale_idw and camb_val_idw in the figure below) from the main map with "Table of Contents" into the Model diagram window. You can also select layers to add from disk using the standard ArcGIS "yellow Plus" icon.

    One you have your input data, you need to select the geoprocessing operator you want to use. In this case, it is the "Weighted Sum " (under Spatial Analyst Tools, Overlay in the toolbox). Remember that if you cannot find a tool, there is an index and a search available. Drag and drop the weighted overlay operator onto the ModelBuilder window.

    Weighted Sum

    After dragging and dropping the 'weighted sum' operator, the model window will look something like this:

    Connecting Data to Operators

     

    Next, use the "link tool" (second from the right) to link from source data to the operator. When your model is correctly connected, the operator will turn yellow to indicate that it has all required inputs. Once you make the connection (see below), double click on the Weighted Overlay operator. This brings up a dialog box in which you can specify the weights to use for each input grid cell value. Click the triangular 'run' button to run the entire model.

     

    Connected

     

    Try to create a weighted overlay in which the med_hvalue is weighted at 70% influence where non-0, and the sales89_pw2-2 layer interpolated from the sales_89 point has a 30% influence. When the census-derived data is 0, you want the weighted overlay tool to fill in with sales_89 estimates. Getting the weights to be handled correctly for grid cells that are missing in one of the layers can be tricky and may require adding additional steps to the model. Also, using map algebra to combine raster layers can get tricky when you need to rescale grid values or convert them to integer values before doing the map algebra that you want. (ArcMap provides additional Math functions to help in this regard. See, for example, the Int() function to convert floating point grid cell values to integers.) For this exercise, you do not need to turn in any results from your use of Model Builder. However, we will make further use of Model Builder in Homework Set #3.

    -----

    We have only scratched the surface of all the raster-based interpolation and analysis tools that are available. If you have extra time, review the help files regarding the Spatial Analyst extension and work on those parts of the homework assignment that ask you to compute a population density surface for youths. Another suggestion is to use model builder to cut out non-residential areas and reallocate population to the residential portions of the block groups - then compute population density.

    Lab Assignment


    Created by Raj Singh and Joseph Ferreira.
    Modified for 1999-2008 by Joseph Ferreira, Thomas H. Grayson, Jeeseong Chung, Jinhua Zhao, Xiongjiu Liao, Diao Mi, Michael Flaxman, and Yang Chen.
    Last Modified: October 27, 2008 by Joe Ferreira.

    Back to the 11.520 Home Page.
    Back to the CRON Home Page.