GEOG 486
Cartography and Visualization

Part III: Interpolation, Contouring and Symbolization

PrintPrint

The Spatial Analyst extension gives you access to several interpolation methods used for calculating regularly-spaced height values for a surface from known values of irregularly-spaced sample data. The resulting interpolated dataset is in a raster format. If one of the desired products of your surface analysis is regularly-spaced isolines (vector data) that represent the surface, these regularly-spaced lattices of values can then be further interpolated linearly and contoured at regularly-spaced intervals.

Before you begin the interpolation process, review these concepts in the concept gallery.

Concept Gallery

Learn more about interpolation, IDW and Kriging in the Concept Gallery.

A. Bring in data from a text file

In this section you will convert spatial and attribute data, that is contained in a text file, into a point layer. (Many of you did this in GEOG 484, Lesson1, performing these steps with a .txt file that contained x-y coordinates and elevation values for benchmarks).

In this lesson you will convert data in a .csv file into a layer containing point features. The csv stands for Comma Separated Values. Each line in the text file is comparable to a row, or record, in a database table. The data in each line of the file is separated by commas. The commas denote the breaks between columns or fields in the database table. It is not too much different from a .txt file containing comma-delimited values. A .csv file was more than likely saved out of a spreadsheet application. If you have Microsoft Excel loaded on your computer the .csv file will most likely be associated with Excel by default. You can, however, also open the .csv file with a text editor. In the text editor (Notepad, for example) you will see the commas. In Excel the data will be separated into columns.

The point features in the layer you create represent meteorological data collection stations in the state of Oregon. In addition to longitude and latitude values the data file contains: the elevation of the station; precipitation amounts for the month of December 2003; and average precipitation for 1971-2000. Based on personal communication with the staff of the Oregon Climate Center, I can tell you that the datum associated with the lon-lat coordinates is WGS84.

Thanks to Mandy Matzke and the staff of the Oregon Climate Center for providing the precipitation data and the state of Oregon polygon coverage that we are using in this exercise.

Let's look at the contents of the .csv file.

  1. Open the oregon_ppt.csv file with Excel, if you have the Excel program.
  2. Now, open the oregon_ppt.csv file with Notepad.
    The first line, or record, in the file contains column headings.
    Note that there is a field named LongNeg. This is a field that I added. It contains the same information you see in the LONG field, except that there are negative signs in front of each value. Do you recall why this is necessary?
    When you convert the contents of the text file to a GIS layer you will indicate that the x-values will come from the LongNeg field (not from the LONG field), and that the y-values will come from the LAT field.
  3. In your oregon_lesson6.mdb geodatabase, create a feature class from the contents of oregon_ppt.csv. Name it oregon_ppt.
    I am not going to provide you with step-by-step instructions, but I will provide you with some hints:
    There are two ways to ingest text file-based x-y data. One method is performed in ArcMap (File > Add Data > Add XY Data...), the other in ArcCatalog (right-click on the csv file). If you have the .csv file open in Excel, close it.
    You will need to correctly define the coordinate system of the x-y data, one way or another, since it is not inherent in the data contained in the .csv file.  From the above information you know that the x-y coordinates are lon-lat and that the datum is WGS84.
    noteIllegal characters: The GIS is particular when it comes to which characters it will tolerate in the field heading names of a text file. If there are illegal characters in the field names, you will get an error message to that effect in ArcCatalog, but you will not receive an error message in ArcMap. In ArcMap, the process simply will not work, with no explanation.

B. Let's not do GIS in lon-lat coordinates

Your new point layer of precipitation gauging stations is in lon-lat coordinates. Again, in Lesson 7 we will discuss and see reasons for not performing GIS analysis when the coordinates of the spatial data are in degrees of longitude and latitude. For now, let's just take this to be true and get the points out of lon-lat coordinates.

  1. Create a version of the oregon_ppt layer that is in the same projection as the or_bound_lesson6 feature class. (HINT: you can use the Project tool and import the coordinate system information by selecting Import...  under the Add Coordinate System button.) Name the new feature class oregon_ppt_ea.
    If you do import the coordinate system parameters from or_bound_lesson6 you will notice WGS_1984_(ITRFOO)_To_NAD_1983 appears in the window below the Geographic Transformation section.  The software will automatically use this to transform your data from WGS84 to NAD83.

C. Interpolate using the Inverse Distance Weighted method

Inverse Distance Weighting is a much-used interpolation method. The height at a given location of any continuous surface is estimated as a distance-weighted sum of the height of known data points in some surrounding neighborhood. Determining the level of importance to assign to the distance (the power of the function) of neighboring points, and determining the nature and extent of that neighborhood are the challenges when it comes to implementing the IDW method of interpolation. Refer to the Concept Gallery item on Inverse Distance Weighting.

  1. In ArcMap open your gtopo30_equalarea map document with the rain gauge station location points on your Oregon map (i.e. with the oregon_ppt_ea feature class added to your ArcMap session).
  2. Before we interpolate, let's find out what the range in values is for the December 2003 precipitation data. Open the attribute table of the oregon_ppt_ea layer.
  3. Right-click on the PPT_DEC03 field label and,
  4. Choose one of the Sort choices. Take note of the range of the data values.
    Now, back to the interpolation.
  5. In ArcToolbox go to > Spatial Analyst Tools > Interpolation and start the IDW tool,
    • Be certain that oregon_ppt_ea is in the Input point features: window.
    • Set the Z value field to PPT_DEC03.
    • Name the resulting IDW raster and place it in your Lesson 6 directory using the folder icon next to the Output raster window.
    • Let the cell size default.
    • Set the remainder of the values to your liking, based on what you gleaned from the concept gallery, or simply experiment to see what happens.
  6. Click the OK button.
    By default the results of the interpolation will be displayed in 9 classes. Experiment with using Stretched or unclassed symbolization if you like.
    Zoom in, if necessary, in order to get a feel for the cell size of the resulting raster dataset.
    With the Oregon state boundary and the masked hillshade image visible you can see that because we do not have precipitation data that extend beyond the extent of the boundaries of the area in question, we fall short of the state boundary polygon. The cells making up the raster surface created by the interpolation routine extend only as far as the extent of the data points.

D. Contour the IDW surface

With the precipitation data interpolated to a regularly-spaced grid of values you are now able to interpolate contour lines. Go ahead and do so.

  1. In ArcToolbox go to > Spatial Analyst Tools > Surface and start the Contour tool,
    Supply values for the settings in the dialog window, being certain to contour the IDW results from the previous section, and saving the results to your Lesson6 folder. The output will be a feature class.
  2. Open the attribute table of the new contour layer. Note the CONTOUR field.

E. Interpolate based on the Kriging method

According to O'Sullivan and Unwin (2010), Kriging is a statistical interpolation method that is optimal in the sense that it makes best use of what can be inferred about the spatial structure in the surface to be interpolated from an analysis of the control point data. O'Sullivan, D. and D. Unwin. (2010). Geographic Information Analysis, 2nd Edition. Hoboken, NJ: John Wiley & Sons.

Spatial Analyst gives us access to basic Kriging functionality but, unfortunately, independent knowledge of which function is best to use in order to estimate the semivariogram is needed, as is the case with the Range, Sill and Nugget parameters.

Also from O'Sullivan and Unwin (2010): How best to estimate the experimental semivariogram and how to choose an appropriate mathematical function to model it is to some extent a "black art" that calls for careful analysis informed by good knowledge of the variable being analyzed. It is definitely not something that you should leave to be decided automatically by a simpleminded computer program, although this is precisely what some GIS programs attempt to do. Many experienced workers fit models by eye, others use standard but complex numerical approaches, and still others proceed by trial and error.

It is not within the scope of this course to teach the statistics of Kriging. What we can do is provide you with a sense of the nature of the method, via the Concepts in Lesson 6, and make you aware of the interface provided in ArcGIS.

  1. In ArcToolbox go to > Spatial Analyst Tools > Interpolation and start the Kriging tool
    • Be certain that oregon_ppt_ea is in the Input point features: window.
    • Set the Z value field to PPT_DEC03.
    • Name the resulting Kriging raster and place it in your Lesson 6 directory using the folder icon next to the Output surface raster window.
    • Let the cell size default.
    • Based on.... whim, provide values for the other settings. Of course, if you are already schooled in the use of Kriging and can come up with an educated estimate of the semivariogram, please do so, and let the rest of the class know what you come up with. Refer to the Concept Gallery item on Kriging.

F. Contour the Kriging results

If you do create a Kriging surface, please go ahead and contour it. The steps are the same as those outlined above in Section D.

G. Render the results - control line symbolization, add hypsometric tinting, control the effect of the shade relief image

When it comes to displaying all of the information we have been creating there are many symbology variables to take into consideration. We will discuss and illustrate a few of them.

  1. Arrange your layers, top to bottom, as follows:
    • or_bound_lesson6
    • oregon_mask
    • one of the precipitation contour layers
    • a hillshaded layer (you may have more than one)
    • the interpolation raster that corresponds with the precipitation contours
    Spend some time setting and experimenting with the following for examples of the sort of thing you can end up with.
    The display of the created information, showing a range of blue colors that have been assigned to each of the classes of grid data.
    Figure 6.3.1 The interpolated precipitation surface grid has been classified into ranges that correspond to certain of the contour lines. A range of blue colors has been assigned to each of the classes of grid data. The hillshaded image has a transparency setting of 35%.

    A view of the image without the contour lines, so that the blue shading is easily seen.
    Figure 6.3.2 At times you may want to let the tinting stand on its own, without the contour lines.
    These are not stellar examples. You should strive to make yours better.
  2. Thickness of the contour lines: Weighting, or varying the thickness of, isolines can help the visual communication of the data. With this example we will increase the weight of certain index contours.
    • Open the Properties of a contour line layer.
    • Click the Symbology tab.
    • In the Show: pane choose Categories, then choose Unique.
    • Set the Value Field to CONTOUR.
    • Click the Add Values... button and Ctrl-click certain values that you would like to appear as index contours, depending on the intervals you contoured. For example, if you created a contour for every inch of precipitation you could select every 5 inches to represent a slightly heavier weight index contour.

      From here you can double-click the individual line symbols and change their appearance. You can also shift-click a group of line symbols, then right-click and choose Properties For Selected Symbol(s) in order to give all symbols selected the same characteristics. Then you can give <all other values>, which would be all the contours that are not the index contours you added earlier, a thinner line, for example.
  3. Tinting: You can emphasize the contoured information by putting the cells of the interpolated height grid into categories that correspond to the contour interval.
    • Open the Properties of an interpolated layer (the results of using the IDW method, for example).
    • Click the Symbology tab.
    • In the Show: pane choose Classified.
    • Change the number of Classes: to however many contour intervals (or ranges of intervals) you want to colorize.
    • Hit the Classify... button.
    • In the next window that appears, in the Break Values pane, type in the values that correspond to your contour intervals. This will make the classes of the IDW correspond to the contour intervals.
    • Choose a Color Ramp, or right-click on individual color symbols and assign colors to the ranges of cell values.
    • Click on the Label column header and select Format Labels... to edit general appearance of the labels for the classes.
  4. Controlling the dominance of the shaded relief image:
    • Open the Properties of your hillshaded layer.
    • Click the Display tab.
    • Modify the value in the Transparency slot.

That is it for Part III, ....and for Lesson 6

Before you leave this lesson I am going to make a plug for a course. It is called Geography 586: Geographic Information Analysis. It was originally developed by Dr. David O'Sullivan, coauthor of the book: Geographic Information Analysis, 2nd Edition. Hoboken, NJ: John Wiley & Sons. (2010) by O'Sullivan, D. and D. Unwin. The course will go into greater depth as far as things like Kriging are concerned, and it will take advantage of the Geostatistical Analyst extension.

Contact the instructor if you have difficulty viewing this image
Figure 6.3.3 A glimpse of the kriging interface that is available in the Geostatistical Analyst extension.