16 min read

In this article by Daniela Cristiana Docan, author of ArcGIS for Desktop Cookbook, we will learn that the ArcGIS Spatial Analyst extension offers a lot of great tools for geoprocessing raster data. Most of the Spatial Analyst tools generate a new raster output. Before starting a raster analysis session, it’s best practice to set the main analysis environment parameters settings (for example, scratch the workspace, extent, and cell size of the output raster). In this article, you will store all raster datasets in file geodatabase as file geodatabase raster datasets.

(For more resources related to this topic, see here.)

Analyzing surfaces

In this recipe, you will represent 3D surface data in a two-dimensional environment. To represent 3D surface data in the ArcMap 2D environment, you will use hillshades and contours. You can use the hillshade raster as a background for other raster or vector data in ArcMap.

Using the surface analysis tools, you can derive new surface data, such as slope and aspect or locations visibility.

Getting ready

In the surface analysis context:

  • The term slope refers to the steepness of raster cells
  • Aspect defines the orientation or compass direction of a cell
  • Visibility identifies which raster cells are visible from a surface location

In this recipe, you will prepare your data for analysis by creating an elevation surface named Elevation from vector data. The two feature classes involved are the PointElevation point feature class and the ContourLine polyline feature class. All other output raster datasets will derive from the Elevation raster.

How to do it…

Follow these steps to prepare your data for spatial analysis:

    1. Start ArcMap and open the existing map document, SurfaceAnalysis.mxd, from <drive>:PacktPublishingDataSpatialAnalyst. Go to Customize | Extensions and check the Spatial Analyst extension.
    2. Open ArcToolbox, right-click on the ArcToolbox toolbox, and select Environments. Set the geoprocessing environment as follows:
      • Workspace | Current Workspace: DataSpatialAnalystTOPO5000.gdb and Scratch Workspace: DataSpatialAnalystScratchTOPO5000.gdb.
      • Output Coordinates: Same as Input.
      • Raster Analysis | Cell Size: As Specified below: type 0.5 with unit as m.
      • Mask: SpatialAnalystTOPO5000.gdbTrapezoid5k.
      • Raster Storage | Pyramid: check Build pyramids and Pyramid levels: type 3. Click on OK.
    3. In ArcToolbox, expand Spatial Analyst Tools | Interpolation, and double-click on the Topo to Raster tool to open the dialog box. Click on Show Help to see the meaning of every parameter. Set the following parameters:
      • Input feature data:

        PointElevation Field: Elevation and Type: PointElevation

        ContourLine Field: Elevation and Type: Contour

        WatercourseA Type: Lake

      • Output surface raster: …ScratchTOPO5000.gdbElevation
      • Output extent (optional): ContourLine
      • Drainage enforcement (optional): NO_ENFORCE
    4. Accept the default values for all other parameters. Click on OK. The Elevation raster is a continuous thematic raster. The raster cells are arranged in 4,967 rows and 4,656 columns.
    5. Open Layer Properties | Source of the raster and explore the following properties: Data Type (File Geodatabase Raster Dataset), Cell Size (0.5 meters) or Spatial Reference (EPSG: 3844).
    6. In the Layer Properties window, click on the Symbology tab. Select the Stretched display method for the continuous raster cell values as follows: Show: Stretched and Color Ramp: Surface. Click on OK.
    7. Explore the cell values using the following two options:
      • Go to Layer Properties | Display and check Show MapTips
      • Add the Spatial Analyst toolbar, and from Customize | Commands, add the Pixel Inspector tool

      Let’s create a hillshade raster using the Elevation layer:

    8. Expand Spatial Analyst Tools | Interpolation and double-click on the Hillshade tool to open the dialog box. Set the following parameters:
      • Input raster: ScratchTOPO5000.gdbElevation
      • Output raster: ScratchTOPO5000.gdbHillshade
      • Azimuth (optional): 315 and Altitude (optional): 45
    9. Accept the default value for Z factor and leave the Model shadows option unchecked. Click on OK.

      From time to time, please ensure to save the map document as MySurfaceAnalysis.mxd at …DataSpatialAnalyst.

    10. The Hillshade raster is a discrete thematic raster that has an associated attribute table known as Value Attribute Table (VAT). Right-click on the Hillshade raster layer and select Open Attribute Table.
    11. The Value field stores the illumination values of the raster cells based on the position of the light source. The 0 value (black) means that 25406 cells are not illuminated by the sun, and 254 value (white) means that 992 cells are entirely illuminated. Close the table.
    12. In the Table Of Contents section, drag the Hillshade layer below the Elevation layer, and use the Effects | Transparency tool to add a transparency effect for the Elevation raster layer, as shown in the following screenshot:

      In the next step, you will derive a raster of slope and aspect from the Elevation layer.

    13. Expand Spatial Analyst Tools | Interpolation and double-click on the Slope tool to open the dialog box. Set the following parameters:
      • Input raster: Elevation
      • Output raster: ScratchTOPO5000.gdbSlopePercent
      • Output measurement (optional): PERCENT_RISE
    14. Click on OK. Symbolize the layer using the Classified method, as follows: Show: Classified. In the Classification section, click on Classify and select the Manual classification method. You will add seven classes. To add break values, right-click on the empty space of the Insert Break graph. To delete one, select the break value from the graph, and right-click to select Delete Break. Do not erase the last break value, which represents the maximum value. Secondly, in the Break Values section, edit the following six values: 5; 7; 15; 20; 60; 90, and leave unchanged the seventh value (496,6). Select Slope (green to red) for Color Ramp. Click on OK. The green areas represent flatter slopes, while the red areas represent steep slopes, as shown in the following screenshot:

    1. Expand Spatial Analyst Tools | Interpolation and double click on the Aspect tool to open the dialog box. Set the following parameters:
      • Input raster: Elevation
      • Output raster: ScratchTOPO5000.gdbAspect
    2. Click on OK. Symbolize the Aspect layer. For Classify, click on the Manual classification method. You will add five classes. To add or delete break values, right-click on the empty space of the graph, and select Insert / Delete Break. Secondly, edit the following four values: 0; 90; 180; 270, leaving unchanged the fifth value in the Break Values section. Click on OK.
    3. In the Symbology window, edit the labels of the five classes as shown in the following picture. Click on OK. In the Table Of Contents section, select the <VALUE> label, and type Slope Direction. The following screenshot is the result of this action:

      In the next step, you will create a raster of visibility between two geodetic points in order to plan some topographic measurements using an electronic theodolite. You will use the TriangulationPoint and Elevation layers:

    4. In the Table Of Contents section, turn on the TriangulationPoint layer, and open its attribute table to examine the fields.

      There are two geodetic points with the following supplementary fields: OffsetA and OffsetB. OffsetA is the proposed height of the instrument mounted on its tripod above stations 8 and 72. OffsetB is the proposed height of the reflector (or target) above the same points.

    5. Close the table. Expand Spatial Analyst Tools | Interpolation and double-click on the Visibility tool to open the dialog box. Click on Show Help to see the meaning of every parameter. Set the following parameters:
      • Input raster: Elevation
      • Input point or polyline observer features: TOPO5000.gdbGeodeticPointsTriangulationPoint
      • Output raster: ScratchTOPO5000.gdbVisibility
      • Analysis type (optional): OBSERVERS
      • Observer parameters | Surface offset (optional): OffsetB
      • Observer offset (optional): OffsetA
      • Outer radius (optional): For this, type 1600

      Notice that OffsetA and OffsetB were automatically assigned. The Outer radius parameter limits the search distance, and it is the rounded distance between the two geodetic points. All other cells beyond the 1,600-meter radius will be excluded from the visibility analysis.

    6. Click on OK. Open the attribute table of the Visibility layer to inspect the fields and values.

      The Value field stores the value of cells. Value 0 means that cells are not visible from the two points. Value 1 means that 6,608,948 cells are visible only from point 8 (first observer OBS1). Value 2 means that 1,813,578 cells are visible only from point 72 (second observer OBS2). Value 3 means that 4,351,861 cells are visible from both points. In conclusion, there is visibility between the two points if the height of the instrument and reflector is 1.5 meters. Close the table.

    7. Symbolize the Visibility layer, as follows: Show: Unique Values and Value Field: Value. Click on Add All Values and choose Color Scheme: Yellow-Green Bright. Select <Heading> and change the Label value to Height 1.5 meters. Double-click on the symbol for Value as 0, and select No Color. Click on OK. The Visibility layer is symbolized as shown in the following screenshot:

  1. Turn off all layers except the Visibility, TriangulationPoint, and Hillshade layers. Save your map as MySurfaceAnalysis.mxd and close ArcMap.

You can find the final results at <drive>:PacktPublishingDataSpatialAnalystSurfaceAnalysis.

How it works…

You have started the exercise by setting the geoprocessing environment. You will override those settings in the next recipes. At the application level, you chose to build pyramids. By creating pyramids, your raster will be displayed faster when you zoom out. The pyramid levels contain the copy of the original raster at a low resolution. The original raster will have a cell size of 0.5 meters. The pixel size will double at each level of the pyramid, so the first level will have a cell size of 1 meter; the second level will have a cell size of 2 meters; and the third level will have a cell size of 4 meters.

Even if the values of cells refer to heights measured above the local mean sea level (zero-level surface), you should consider the planimetric accuracy of the dataset. Please remember that TOPO5000.gdb refers to a product at the scale 1:5,000. This is the reason why you have chosen 0.5 meters for the raster cell size. At step 4, you used the PointElevation layer as supplementary data when you created the Elevation raster.

If one of your ArcToolbox tools fails to execute or you have obtained an empty raster output, you have some options here:

    • Open the Results dialog from the Geoprocessing menu to explore the error report. This will help you to identify the parameter errors. Right-click on the previous execution of the tool and choose Open (step 1). Change the parameters and click on OK to run the tool. Choose Re Run if you want to run the tool with the parameters unchanged (step 2) as shown in the following screenshot:

  • Run the ArcToolbox tool from the ArcCatalog application. Before running the tool, check the geoprocessing environment in ArcCatalog by navigating to Geoprocessing | Environments.

There’s more…

What if you have a model with all previous steps?

Open ArcCatalog, and go to …DataSpatialAnalystModelBuilder. In the ModelBuilder folder, you have a toolbox named MyToolbox, which contains the Surface Analysis model. Right-click on the model and select Properties. Take your time to study the information from the General, Parameters, and Environments tabs. The output (derived data) will be saved in Scratch Workspace: ModelBuilder ScratchTopo5000.gdb. Click on OK to close the Surface Analysis Properties window. Running the entire model will take you around 25 minutes. You have two options:

  • Tool dialog option: Right-click on the Surface Analysis model and select Open. Notice the model parameters that you can modify and read the Help information. Click on OK to run the model.
  • Edit mode: Right-click on the Surface Analysis model and select Edit. The colored model elements are in the second state—they are ready to run the Surface Analysis model by using one of those two options:
    • To run the entire model at the same time, select Run Entire Model from the Model menu.
    • To run the tools (yellow rounded rectangle) one by one, select the Topo to Raster tool with the Select tool, and click on the Run tool from the Standard toolbar. Please remember that a shadow behind a tool means that the model element has already been run.

You used the Visibility tool to check the visibility between two points with 1.5 meters for the Observer offset and Surface offset parameters. Try yourself to see what happens if the offset value is less than 1.5 meters. To again run the Visibility tool in the Edit mode, right-click on the tool, and select Open. For Surface offset and Observer offset, type 0.5 meters and click on OK to run the tool. Repeat these steps for a 1 meter offset.

Interpolating data

Spatial interpolation is the process of estimating an unknown value between two known values taking into account Tobler’s First Law:

“Everything is related to everything else, but near things are more related than distant things.”

This recipe does not undertake to teach you the advanced concept of interpolation because it is too complex for this book. Instead, this recipe will guide you to create a terrain surface using the following:

  • A feature class with sample elevation points
  • Two interpolation methods: Inverse Distance Weighted (IDW) and Spline

For further research, please refer to: Geographic Information Analysis, David O’Sullivan and David Unwin, John Wiley & Sons, Inc., 2003, specifically the 8.3 Spatial interpolation recipe of Chapter 8, Describing and Analyzing Fields, pp.220-234.

Getting ready

In this recipe, you will create a terrain surface stored as a raster using the PointElevation sample points. Your sample data has the following characteristics:

  • The average distance between points is 150 meters
  • The density of sample points is not the same on the entire area of interest
  • There are not enough points to define the cliffs and the depressions
  • There are not extreme differences in elevation values

How to do it…

Follow these steps to create a terrain surface using the IDW tool:

  1. Start ArcMap and open an existing map document Interpolation.mxd from <drive>:PacktPublishingDataSpatialAnalyst.
  2. Set the geoprocessing environment, as follows:
    • Workspace | Current Workspace: DataSpatialAnalystTOPO5000.gdb and Scratch Workspace: DataSpatialAnalystScratchTOPO5000.gdb
    • Output Coordinates: Same as the PointElevation layer
    • Raster Analysis | Cell Size: As Specified Below: 1
    • Mask: DataSpatialAnalystTOPO5000.gdbTrapezoid5k

In the next two steps, you will use the IDW tool. Running IDW with barrier polyline features will take you around 15 minutes:

  1. In ArcToolbox, go to Spatial Analyst Tools | Interpolation, and double-click on the IDW tool. Click on Show Help to see the meaning of every parameter. Set the following parameters:
    • Input point features: PointElevation
    • Z value field: Elevation
    • Output raster: ScratchTOPO5000.gdbIDW_1
    • Power (optional): 0.5
    • Search radius (optional): Variable
    • Search Radius Settings | Number of points: 6
    • Maximum distance: 500
    • Input barrier polyline features (optional): TOPO5000.gdbHydrographyWatercourseL
  2. Accept the default value of Output cell size (optional). Click on OK.
  3. Repeat step 3 by setting the following parameters:
    • Input point features: PointElevation
    • Z value field: Elevation
    • Output raster: ScratchTOPO5000.gdbIDW_2
    • Power (optional): 2
  4. The rest of the parameters are the same as in step 3. Click on OK. Symbolize the IDW_1 and IDW_2 layers as follows: Show: Classified; Classification: Equal Interval: 10 classes; Color Scheme: Surface. Click on OK. You should obtain the following results:

    In the following steps, you will use the Spline tool to generate the terrain surface:

  5. In ArcToolbox, go to Spatial Analyst Tools | Interpolation, and double-click on the Spline tool. Set the following parameters:
    • Input point features: PointElevation
    • Z value field: Elevation
    • Output raster: ScratchTOPO5000.gdbSpline_Regular
    • Spline type (optional): REGULARIZED
    • Weight (optional): 0.1 and Number of points (optional): 6
  6. Accept the default value of Output cell size (optional). Click on OK.
  7. Run again the Spline tool with the following parameters:
    • Input point features: PointElevation
    • Z value field: Elevation
    • Output raster: ScratchTOPO5000.gdbSpline_Tension
    • Spline type (optional): TENSION
    • Weight (optional): 0.1 and Number of points (optional): 6
  8. Accept the default value of Output cell size (optional). Click on OK. Symbolize the Spline_Regular and Spline_Tension raster layers using the Equal Interval method classification with 10 classes and the Surface color ramp:

    In the next steps, you will use the Spline with Barriers tool to generate a terrain surface using an increased number of sample points. You will transform the ContourLine layer in a point feature class. You will combine those new points with features from the PointElevation layer:

  9. In ArcToolbox, go to Data Management Tools | Features, and double-click on the Feature vertices to Points tool. Set the following parameters:
    • Input features: ContourLine
    • Output Feature Class: TOPO5000.gdbRelief ContourLine_FeatureVertices
    • Point type (optional): ALL
  10. Click on OK. Inspect the attribute table of the newly created layer. In the Catalog window, go to …TOPO5000.gdbRelief and create a copy of the PointElevation feature class. Rename the new feature class as ContourAndPoint. Right-click on ContourAndPoint and select Load | Load Data. Set the following parameters from the second and fourth panels:
    • Input data: ContourLine_FeatureVertices
    • Target Field: Elevation
    • Matching Source Field: Elevation
  11. Accept the default values for the rest of the parameters and click on Finish.
  12. In ArcToolbox, go to Spatial Analyst Tools | Interpolation, and double-click on the Spline with Barriers tool. Set the following parameters:
    • Input point features: ContourAndPoint
    • Z value field: Elevation
    • Input barrier features (optional): TOPO5000.gdbHydrographyWatercourseA
    • Output raster: ScratchTOPO5000.gdbSpline_WaterA
    • Smoothing Factor (optional): 0
  13. Accept the default value of Output cell size (optional). Click on OK. You should obtain a similar terrain surface to what’s shown here:

    Explore the results by comparing the similarities or differences of the terrain surface between interpolated raster layers and the ContourLine vector layer. The IDW method works well with a proper density of sample points. Try to create a new surface using the IDW tool and the ContourAndPoint layer as sample points.

  14. Save your map as MyInterpolation.mxd and close ArcMap.

You can find the final results at <drive>:PacktPublishingDataSpatialAnalyst Interpolation.

How it works…

The IDW method generated an average surface that will not cross through the known point elevation values and will not estimate the values below the minimum or above the maximum given point values. The IDW tool allows you to define polyline barriers or limits in searching sample points for interpolation. Even if the WatercourseL polyline feature classes do not have elevation values, river features can be used to interrupt the continuity of interpolated surfaces.

To obtain fewer averaged estimated values (reduce the IDW smoother effect) you have to:

  • Reduce the sample size to 6 points
  • Choose a variable search radius
  • Increase the power to 2

The Power option defines the influence of sample point values. This value increases with the distance. There is a disadvantage because around a few sample points, there are small areas raised above the surrounding surface or small hollows below the surrounding surface.

The Spline method has generated a surface that crosses through all the known point elevation values and estimates the values below the minimum or above the maximum sample point values. Because the density of points is quite low, we reduced the sample size to 6 points and defined a variable search radius of 500 meters in order to reduce the smoothening effect.

The Regularized option estimates the hills or depressions that are not cached by the sample point values. The Tension option will force the interpolated values to stay closer to the sample point values.

Starting from step 12, we increased the number of sample points in order to better estimate the surface. At step 14, notice that the Spline with Barriers tool allows you to use the polygon feature class as breaks or barriers in searching sample points for interpolation.

Summary

In this article, we learned about the ArcGIS Spatial Analyst extension and its tools.

Resources for Article:

 

Further resources on this subject:

Packt

LEAVE A REPLY

Please enter your comment!
Please enter your name here