Working with rasters  

Working with ArcGIS Spatial Analyst objects

This topic is relevant for the following:
Product(s): ArcGIS Desktop: All, ArcGIS Engine
Version(s): 9.0, 9.1, 9.2
Language(s): VB6
Experience level(s): Insert comma-separated list of experience levels here

The ArcObjects components for the ArcGIS Spatial Analyst extension provide an easy to customize environment with a rich set of developer tools allowing model development. The Spatial Analyst works closely with the Raster Data Objects (RDO). The RDO model is used to access, analyze, manage, convert, and visualize raster data. RDO provides some features to compliment the Spatial Analyst object model such as unified data access and interaction independent format, a high level nonproprietary language, and an environment that is easy to extend, upgrade, and customize.

The Spatial Analyst object model is composed of a collection of objects to perform raster analysis. There are two main types of objects: objects that support the overall Spatial Analyst environment (e.g., the analysis environment, selection, and toolbar), and objects that perform analysis operations (e.g., slope, aspect, and distance analysis). The objects that perform analysis operations are functionally grouped according to the type of analysis they perform. That is, CostDistance, CostAllocation, EuclideanDistance, and CostPath, as well as others are methods of the distance operator object, RasterDistanceOp.

This paper is divided into two sections. The first section provides a basic overview of the steps necessary for implementing a Spatial Analyst operation. The second section provides details on additional features and how the system works.

Section 1: Steps for using the Spatial Analyst objects

There are many objects in the Spatial Analyst object model and hundreds of methods to perform spatial modeling. Even though there are many objects and methods, there are five main steps when implementing any Spatial Analyst operation. They are:

  1. Obtaining a Spatial Analyst license
  2. Setting the analysis environment
  3. Accessing the input data
  4. Performing an operation
  5. Using the output

The following text discusses these five steps and provides Visual Basic examples demonstrating the step.

1. Obtaining a Spatial Analyst license

In order to access the Spatial Analyst objects you must obtain a Spatial Analyst license. You can call for a license prior to accessing each Spatial Analyst object, call it once in a single subroutine, or call once globally in the Main routine in a program. Since you only need to call the single subroutine once, or the license is automatically enabled if you call it globally, either of these two approaches are the preferred methods for obtaining a license. The following code shows you how to obtain a Spatial Analyst license through a single subroutine.

[Visual Basic 6.0]
Public Sub CheckOutSpatialAnalystLicense()
  ' This subroutine checks out the SpatialAnalyst license in a standalone VB application.

  Dim pAoInitialize As IAoInitialize
  Set pAoInitialize = New AoInitialize
  pAoInitialize.Initialize (esriLicenseProductCodeArcView)
  pAoInitialize.CheckOutExtension (esriLicenseExtensionCodeSpatialAnalyst)

End Sub

2. Setting the analysis environment

Certain operators, such as the RasterDensityOp and the RasterInterpolationOp, require that some properties of the analysis output be specified. Two properties of particular importance to set are the output cell size and the extent of the analysis.

The analysis environment (through the RasterAnalysis object) controls four primary properties when performing analysis, the cell size, extent, mask, and spatial reference. Even though these settings do nothing to the original data, the extent and mask do affect the locations where the operation will occur. The extent defines the area on which the operation will be performed. The mask defines the cell locations within the extent on which to perform the operation. And, the specified cell size controls the resolution of the output.

Each operator object has its own analysis environment settings. You can set the environment parameters individually for each operator object (see 'The analysis environment' in the second part of this document). When first created, the new operator object obtains its parameters from the settings in the RasterAnalysis object. The following function demonstrates how to change the default analysis environment settings that subsequent newly created operator objects will inherit.

[Visual Basic 6.0]
Public Function SetNewDefaultEnvironment(pExtent As IEnvelope, nCellSize As Double, _
                         pMask As IGeoDataset, pWS As IWorkspace, pPrj As ISpatialReference) _
                         As IRasterAnalysisEnvironment
  ' Creates new RasterAnalysis object and sets it as a new default settings

  On Error GoTo ERH

  ' Create a RasterAnalysis object
  Dim pEnv As IRasterAnalysisEnvironment
  Set pEnv = New RasterAnalysis

  ' Set the new default extent
  If Not pExtent Is Nothing Then pEnv.SetExtent esriRasterEnvValue, pExtent

  ' Set the new default cellsize
  If nCellSize > 0 Then pEnv.SetCellSize esriRasterEnvValue, nCellSize

  ' Set the new default mask for analysis
  If Not pMask Is Nothing Then Set pEnv.Mask = pMask

  ' Set the new default output workspace
  If Not pWS Is Nothing Then Set pEnv.OutWorkspace = pWS

  ' Set the new default output spatial reference
  If Not pPrj Is Nothing Then Set pEnv.OutSpatialReference = pPrj

  ' Set it as the default settings

  ' Return reference to the default environment setting
  Set SetNewDefaultEnvironment = pEnv
  Exit Function
  Set SetNewDefaultEnvironment = Nothing
End Function

To reset the global analysis environment back to the defaults you will run the following code snippet (here pNewEnv is the new default settings that is set by the above function).

[Visual Basic 6.0]

3. Accessing the input data

Once you have a license, and have optionally set the analysis environment, you then need to identify the input data on which to apply the Spatial Analyst operation. There are four acceptable raster object inputs into a Spatial Analyst operation that calls for raster input: RasterDataset, RasterBand, Raster, or RasterDescriptor. A RasterDataset represents an existing dataset stored on disk or in a database in a particular raster format (e.g., TIFF). RasterBand are the individual bands of a RasterDataset. Some RasterDatasets have multiple bands like a satellite image while others contain a single band such as an ESRI grid. You cannot alter the values or projection of a RasterDataset or the RasterBands. A Raster is a virtual representation of raster data derived from a data source on disk. A Raster may be a subset of the original RasterDataset, may be projected into a different projection, and can be composed of bands from different RasterDatasets. A RasterDescriptor is used to present a subset of the original Raster (e.g., through an attribute selection), and/or use an alternative value for each cell defined in a field in the datas attribute table. The raster data object model is discussed in detail in the DataSourcesRaster and DataSourcesRasterUI libraries.

In addition, two feature data objects are valid inputs, FeatureClass and FeatureClassDescriptor. The features will be converted to rasters using the existing environment settings. A FeatureClass is converted directly, while a subset of the original FeatureClass will be converted in a FeatureClassDescriptor.

As previously mentioned, a Spatial Analyst operation can be performed on a RasterDataset, a RasterBand, a Raster, a RasterDescriptor, a FeatureClass, or a FeatureClassDescriptor. When accessing a RasterDataset the RasterWorkspaceFactory must be used. The following function accesses a raster on disk.

[Visual Basic 6.0]
Public Function OpenRasterDataset(sPath As String, sFileName As String) As IRasterDataset
   ' Returns RasterDataset object given a file name and its directory
   ' sPath: path of the input raster dataset
   ' sFileName: name of the input raster dataset

  On Error GoTo ERH
  Dim pWSFact As IWorkspaceFactory
  Dim pRasterWS As IRasterWorkspace

  Set pWSFact = New RasterWorkspaceFactory
  If pWSFact.IsWorkspace(sPath) Then
    Set pRasterWS = pWSFact.OpenFromFile(sPath, 0)
    Set OpenRasterDataset = pRasterWS.OpenRasterDataset(sFileName)
  End If
    Exit Function
  Set OpenRasterDataset = Nothing
  MsgBox "Failed in Opening RasterDataset. " & Err.Description
End Function

A Raster object, once created or if already in existence, can be used directly in an operation.

To access a raster layer from ArcMap or a shapefile or coverage from disk, see the 'Using other ArcGIS objects with Spatial Analyst' section later in this document.

You may not want to perform a Spatial Analyst operation on the entire input dataset, but rather on a subset of the data. Or you may want an operation to use values defined in a certain field in the datas attribute table.

An attribute subset can be created from the input dataset identifying which cells or features on which to perform the operation. To create an attribute selection or subset of the cells of an input RasterDataset, you create a Raster object from the RasterDataset and then create a RasterDescriptor from the Raster. The RasterDescriptor identifies the qualifying selection. The following code uses a RasterDescriptor to create an attribute subset. Alternatively, the same thing can be done with feature data inputs by creating a FeatureClassDescriptor from a FeatureClass. Note that the following subroutine calls the OpenRasterDataset function, which was discussed in the previous section.

[Visual Basic 6.0]
Public Function CreateRasterDescriptor(sPath As String, sFileName As String, _
                          sWhereClause As String, sSubField As String, sInField As String) _
                          As IRasterDescriptor
  ' This function returns a RasterDescriptor object
  ' sPath: path of the raster dataset
  ' sFileName: name of the raster dataset
  ' sWhereClause: where clause for the QueryFilter
  ' sSubField: a list of field names. "" will use all fields
  ' sInField: a field name

  On Error GoTo ERH
  ' Create and set query filter
  Dim pFilt As IQueryFilter
  Set pFilt = New QueryFilter
  If Len(sWhereClause) > 0 Then pFilt.WhereClause = sWhereClause
  If Len(sSubField) > 0 Then pFilt.SubFields = sSubField

  ' Create raster from dataset
  Dim pRaster As IRaster
  Set pRaster = OpenRasterDataset(sPath, sFileName).CreateDefaultRaster

  ' Create RasterDescriptor
  Dim pDesc As IRasterDescriptor
  Set pDesc = New RasterDescriptor
  pDesc.Create pRaster, pFilt, sInField
  Set CreateRasterDescriptor = pDesc
  Exit Function

  Set CreateRasterDescriptor = Nothing
  MsgBox "Failed in Creating RasterDescriptor " & Err.Description
End Function

You can also create a RasterDescriptor or FeatureClassDescriptor using any desired field. See 'Fields' later in this document).

4. Performing an operation

The majority of the objects in the Spatial Analyst object model are operator objects that contain the methods for spatial analysis. These operator objects contain the majority of the commands and functions that were available in ArcGRID and ArcView 3 Spatial Analyst. The methods are grouped functionally into objects based on the analysis they perform. There are two main steps for applying a Spatial Analyst operation once the data has been identified. You first query (which creates) the desired Op and second, you call the method. The following code queries the RasterSurfaceOp and invokes the Slope method on an input elevation raster (note that the OpenRasterDataset function described in the previous section is called).

[Visual Basic 6.0]
' Gets an elevation raster from disk
Dim pElevation As IRasterDataset
Set pElevation = OpenRasterDataset("D:\SpatialTest", "elevation")

' Queries the RasterSurfaceOp
Dim pSurfaceOp As ISurfaceOp
Set pSurfaceOp = New RasterSurfaceOp

' Creates a geodataset to store the results and runs the operation
Dim pSlope As IGeoDataset
Set pSlope = pSurfaceOp.Slope(pElevation, esriGeoAnalysisSlopeDegrees)

Even though there are many operator objects with hundreds of methods, to implement any of them the pattern of use is the same. The only things that change from this example is the operator queried, the method called, and the parameters specified.

5. Using the output

If raster data is produced as the output from an operator it will be a temporary. The data on disk will be deleted as soon as it is not referenced by any objects. This temporary output can be used as input into another operator, be saved as a permanent raster dataset, or displayed in ArcMap. To see how to make the output permanent refer the section 'Temporary outputs' later this document.

The code below demonstrates how to display the temporary output in ArcMap.

[Visual Basic 6.0]
Public Sub AddRasterLayer(pApp As IApplication, pRaster As IRaster)
  On Error GoTo ERH
  ' Adds a layer to an ArcMap session
  Dim pRasterLayer As IRasterLayer
  Set pRasterLayer = New RasterLayer
  pRasterLayer.CreateFromRaster pRaster
  Dim pMap As IBasicMap
  Dim pMxDoc As IMxDocument
  Dim pActView As IActiveView
  Set pMxDoc = pApp.Document
  Set pMap = pMxDoc.FocusMap
  Set pActView = pMxDoc.ActiveView
  pMap.AddLayer pRasterLayer
  Exit Sub
  MsgBox "Failed to Add Raster Layer." & Err.Description
End Sub

Even though there are many ways to qualify the input data, perform operations, and use the output, and many possible sequences to perform the desired analysis, the pattern of use for implementing any of the Spatial Analyst operator follows the same five steps identified above. The only things that change in any analysis are the input datasets and qualifiers, the operator, method, and its parameters, and how the output is used.

Section 2: Features of the ArcGIS Spatial Analyst object model

In this section, additional details are given about how the overall Spatial Analyst object model works and how the individual objects interact. Discussed are the issues to consider when using the objects in conjunction with one another to perform a desired analysis.

The functionality included in the ArcGIS Spatial Analyst objects provides approximately 90% of the complete functionality of the ArcGRID module and ArcView 3 Spatial Analyst. Some of the functions not present in the ArcGIS Spatial Analyst include the IF, WHILE, DOCELL, and the multivariate statistics functionality. IF, WHILE, and DOCELL can be simulated using the PixelBlock object (see below) within the Raster core objects. The multivariate statistics functions and a few others can only be accessed from the ArcGRID prompt at this time.

In addition to the standard operator objects in the ArcGIS Spatial Analyst, another way to apply raster operations is through the RasterMapAlgebraOp and RasterModel. These two objects allow you to directly submit GRID Map Algebra syntax command strings to be executed on your data. The use of GRID Map Algebra is discussed in detail in Appendix A of the book Using ArcGIS Spatial Analyst.

The properties of the output from an ArcGIS Spatial Analyst operator can be specified using the raster analysis environment. Not only can you control the analysis extent, cell size, and mask properties, with ArcGIS Spatial Analyst you have control of the spatial reference system the output dataset will be created in and the workspace in which the output will be placed. The default settings for most of these properties are the same as their corresponding settings in ArcGRID.

Referencing the Spatial Analyst objects and their libraries

The Spatial Analyst objects are divided among three different object libraries, each of which contains some, but not all, of the coclasses and interfaces used in the Spatial Analyst extension. The Spatial Analyst objects are divided in this manner because of different licensing schemes. Some of the objects are available with the core ArcGIS product, some are available with the 3D Analyst or Spatial Analyst, and others are available only with the Spatial Analyst.

The analysis environment

The analysis environment controls the properties of the analysis performed by almost every Spatial Analyst operation. It contains the cell size, extent, and mask properties ArcGRID users are familiar with, and provides a way to control the output spatial reference, the output workspace into which temporary results are placed, the verify environment used when over-writing previous results, and the prefix of the output dataset name.

The session default analysis environment controls the analysis environment that is assigned to each operator object when it is created. When you begin a new session of ArcGIS or your application using the Spatial Analyst objects, it contains these initial default settings:

The session default analysis environment can be overridden so each new operator object begins with a different set of analysis environment properties. This can be useful if you are performing an analysis that requires use of multiple operator objects and this analysis needs to use an analysis environment setting that is not the same as the default. To change the session default analysis environment, cocreate either the RasterAnalysis object or any Op objects, then query for the IRasterAnalysisEnvironment interface. Change the settings of the analysis environment to those you wish to make the default, then invoke the SetAsNewDefaultEnvironment method. The SetAsNewDefaultEnvironment method will not change the analysis environment associated with any existing operator objects, but each new operator object created will contain the settings that have been set as the default. See the step 'Setting the analysis environment' in the previous section for example code on how to set the default parameters.

Each operator object has its own analysis environment. This environment can be accessed for any operator object through the IRasterAnalysisEnvironment interface. The methods on this interface allow direct control of the analysis environment on that object. When an operator object is created, the initial analysis environment properties are taken from the session default analysis environment. The following code sets the desired analysis environment properties for a newly created operator object. Note that the example gives an explicit workspace path, cell size, extent, and uses the MakeConstant method from the RasterMakerOp. You can change these parameters and methods to your desired analysis.

[Visual Basic 6.0]
' Create a Op (RasterMaker operator)
Dim pRMakerOp As IRasterMakerOp
Set pRMakerOp = New RasterMakerOp

' Query Op for IRasterAnalysisEnvironment
Dim pEnv As IRasterAnalysisEnvironment
Set pEnv = pRMakerOp

' Set output workspace for the Op
Dim pWS As IWorkspace
Dim pWSF As IWorkspaceFactory
Set pWSF = New RasterWorkspaceFactory
Set pWS = pWSF.OpenFromFile("c:\temp", 0)
Set pEnv.OutWorkspace = pWS

' Set cell size for the Op
pEnv.SetCellSize esriRasterEnvValue, 30

' Set output extent for the Op
Dim pExt As IEnvelope
Set pExt = New Envelope
pExt.XMin = 0
pExt.YMin = 0
pExt.XMax = 3000
pExt.YMax = 3000
pEnv.SetExtent esriRasterEnvValue, pExt

' Perform spatial operation
Dim pOutRaster As IRaster
Set pOutRaster = pRMakerOp.MakeConstant(10, True)

Data integration

While virtually all of the methods in the Spatial Analyst perform their analysis on raster data, inputs can be specified in many vector formats and any raster format supported in ArcGIS. This data conversion is implicitly performed by the Spatial Analyst operator objects using properties specified in the analysis environment and is transparent to the user. If multiband raster data is input into a Spatial Analyst operations, most methods will use the first band of the dataset for analysis. There are a handful of exceptions to this rule (e.g., ILocalOp.Combine and ILocalOp.LocalStatistics). If more control over data conversion properties is desired, the RasterConversionOp can be used.

Only geographic data types are acceptable as inputs, and to provide compiler type checking, an interface supported by all appropriate input types must be used. The IGeoDataset interface is supported by the widest range of geographic data types. Objects that support this interface and are acceptable inputs to most Spatial Analyst operations include the Raster, RasterDataset, RasterBand, RasterDescriptor, FeatureClass, and FeatureClassDescriptor objects. However, many other data types that cannot be used as inputs to spatial operations also support this interface. These data types include any layer type including the RasterLayer and FeatureLayer, as well as Tin objects. If a dataset of a type not supported by Spatial Analyst is input to an operation, the error "Error During Conversion" or "Could not reference the GRID dataset" is returned. These errors may also occur if a vector dataset is input to a method with no cell size set, or when there are other problems with the input data.

The GRID engine that provides most of the functionality present in the Spatial Analyst can only accept raster input in the ESRI GRID format, while a few methods that use feature data as inputs can accept only shapefiles or coverages (e.g., IDW and Density). The Spatial Analyst objects allow as input all supported raster formats and many vector formats into these methods. This is accomplished by converting these datasets to the correct input format prior to performing the operation. Thus, it is not recommended to perform analysis on large, highly compressed rasters, such as CCITT Group 4, JPEG, or MrSID format. The conversion may be slow and consume significant system resources. If you wish to perform analysis on a dataset repeatedly, this data should be converted into ESRI grid or shapefile format to prevent this implicit conversion from occurring for each operation to minimize processing time.

Defaults for methods

In ArcGRID, most function arguments have default values. No input is needed if the default value is to be used for a specific argument. In the Spatial Analyst, there are no default values, so these parameters have been made optional. These optional arguments are not strongly typed, thus, optional arguments are specified as Variant. A Variant argument can accept any data type and therefore the compiler cannot check to make sure the type of the input argument is correct. However, an operator is expecting a specific data type for the argument, and errors may occur if the wrong type of argument is input. Thus, to obtain the desired results, care must be given to the arguments of the methods.


Any numeric field in your input data can be specified as the analysis field for most operations and many operations also support analysis on fields containing character strings. Furthermore, almost all Spatial Analyst operations allow analysis to be performed on a subset of the input data. This subset can be specified using a logical query on any numeric or string field in the input data. Examples of this include a pH field in a soil raster that was created from soil type, a string field containing landuse type in a landuse raster, or a string field of road type in a roads feature dataset. The following VBA code demonstrates how to create a Raster Descriptor. The subroutine allows an attribute subset to be created from the original input raster based on a specified field and where clause. For instance, pInRaster is a Raster object created from a landuse raster dataset, sFieldName is the Field Name for Land use type, and sValue is one of the land use types such as Forest.

[Visual Basic 6.0]
Public Function UsingRasterDescriptor (pInRaster as IRaster, sFieldName as String, sValue as_
                          String) As IRasterDescriptor
  ' Creates a RasterDescriptor from an input raster and applies a where clause
  On Error GoTo ERH

  ' Create the RasterDecriptor
  Dim pRDescr As IRasterDescriptor
  Set pRDescr = New RasterDescriptor

  ' Creates the query
  Dim pQFL as IQueryFilter
  Set pQFL = New QueryFilter
  pQFL.WhereClause = sFieldName + " = " +  sValue

  ' Performs the selection
  pRDescr.Create pInRaster, pQFL, sFieldName
  set UsingRasterDescriptor = pRDescr
  Exit Function
  Set UsingRasterDescriptor = Nothing
End Sub

The resulting raster descriptor can be input into any method. Any method that is applied to the descriptor will only be applied to cells containing the Forest land use type.

Temporary outputs

All raster outputs from Spatial Analyst operations are temporary. A temporary raster is created in the output workspace specified by the IRasterAnalysisEnvironment::OutWorkspace method. If an output workspace is not specified then it goes to a folder specified by the TMP environment variable. In Windows XP, the TMP variable usually points to \Documents and Settings\username\Local Settings\Temp folder unless set to a different location. If you wish to keep an output dataset after objects referencing it go out of scope, you must explicitly make it permanent. To make the output raster permanent, you must first get the RasterDataset associated with the Raster object returned from a Spatial Analyst operation. From the RasterDataset object the ITemporaryDataset interface can be obtained and used to make the dataset permanent. Vector outputs are always permanent.

The temporary dataset created by the Spatial Analyst has a name that is the prefix specified by the analysis environment followed by an integer. This integer is appended to the prefix so each raster will be uniquely named. Each time an output is created, this integer sequentially increases.

The following example shows how to make a temporary output raster permanent.

[Visual Basic 6.0]
Public Sub MakePermanent(pResultOfSpatialOp As IRaster)
  ' Query the output (a Raster object) for IRasterBandCollection
  Dim pRasBandC As IRasterBandCollection
  Set pRasBandC = pResultOfSpatialOp

  ' Get the dataset from the first band
  Dim pRasterDS As IRasterDataset
  Set pRasterDS = pRasBandC.Item(0).RasterDataset

  ' Query the dataset for ITemporaryDataset
  Dim pTemp As ITemporaryDataset
  Set pTemp = pRasterDS

  Dim pWSF As IWorkspaceFactory
  Dim pWS As IWorkspace
  Set pWSF = New RasterWorkspaceFactory
  Set pWS = pWSF.OpenFromFile("c:\temp", 0)

  ' Format can be "GRID", "TIFF", or "IMAGINE Image"
  ' Note: format string is Case sensitive
   pTemp.MakePermanentAs "Result2.img", pWS, "IMAGINE Image"
End Sub

MultiBand output

Some methods, such as Curvature, EucDistanceFull, CostDistanceFull, and FlowDirection, allow additional optional outputs to be returned from the operation. In these cases, the return is a Raster object that is a collection of bands. The collection of bands consists of bands from the different output datasets. To separate them, obtain the IRasterBandCollection interface on the output Raster and individually access the bands. Once you have the individual bands, you are able to get to the datasets to make them permanent, or whatever. The first band will be the default output of the method, and the bands will be placed in the same order as they appear in the method signature. The following code demonstrates how to access the second band, the profile curve, from the results of the Curvature method.

[Visual Basic 6.0]
Sub GetSecondRaster(pResultOfCurvature As IRaster) as IRaster
  ' Gets the second band of a multiband raster
  Dim pRasBandC As IRasterBandCollection
  Dim pRasBand as IRasterBand
  Dim pOutRasterDS As IRasterDataSet
  Set pRasBandC = pResultOfCurvature
  If pRasBandC.Count > 1 Then
    Set pRasBand = pRasBandC.Item(1)
    Set pOutRasterDS = pRasBand.RasterDataset
    Set GetSecondRaster = pOutRasterDS.CreateDefaultRaster
    GetSecondRaster = Nothing
  End If
End Sub

Deferred evaluation

To increase system performance, raster outputs from Spatial Analyst methods are not computed until they are used. This deferred evaluation of output datasets can result in outputs never being calculated if they are never used. Any use of the output, including creating a layer, examining the height and width of the output, or making the dataset permanent, will trigger evaluation of the expression to create the dataset.

Global functions, such as CostDistance, are not deferred and are evaluated when called whether the output is referred to or not. For all other operations the evaluation of the output is deferred until the output dataset is used. The following code forces the evaluation of the Slope method by calling the GetProperty function that queries the height of the output dataset. The following code calls the OpenRasterDataset function which is provided in the first section of this document in the step 'Accessing the input data'.

[Visual Basic 6.0]
' Gets an elevation raster from disk
Dim pElevation As IRasterDataset
Set pElevation = OpenRasterDataset("D:\SpatialTest", "elevation")

' Queries the RasterSurfaceOp
Dim pSurfaceOp As ISurfaceOp
Set pSurfaceOp = New RasterSurfaceOp

' Creates a geodataset to store the results of the operation
Dim pOutput As IGeoDataset
Set pOutput = pSurfaceOp.Slope(pElevation, esriGeoAnalysisSlopeDegrees)

' Forces the evaluation by querying the Height property
Dim pProp as IRasterProps
Set pProp = pOutput
Msgbox pProp.Height

Spatial reference handling

If all of your data is in the same coordinate system, it is best to perform analysis in that coordinate system. However, the output spatial reference system of any Spatial Analyst object can be specified using the analysis environment by setting the OutSpatialReference property. If the output spatial reference system is not specified, a Spatial Analyst operation will use the spatial reference of the first raster input that has a spatial reference other than Unknown. If no raster inputs have spatial reference information, the spatial reference will be taken from the first feature input with a defined spatial reference. If no input dataset has a defined spatial reference, the output will have an unknown coordinate system. Any other input with a spatial reference different from the output will be projected on the fly into the output coordinate system before the operation begins.

For vector data, this projection on the fly is very accurate, however, raster data is projected using an algorithm that produces faster results than a true projection at the cost of high accuracy. For raster datasets at low latitudes (<50 north or south) and smaller than a few degrees on each side this projection is fairly accurate, but this projection is inappropriate for larger or polar raster datasets. In these cases use ArcGRID to project your data using a cell-based projection for highest accuracy.

Using the PixelBlock for custom analysis

The DOCELL, IF, and WHILE commands or ArcGRID are not available within ArcGIS Spatial Analyst. Each of these constructs loops through all cells in a raster dataset and performs some computation on each cell. The core raster data objects provide the ability to perform direct pixel reading and writing through the PixelBlock object, allowing you to replicate the these tools and create new tools of your own

The PixelBlock allows you to read and write pixel data from a band of a raster dataset on disk or from a Raster. The following code shows how to access the individual cell values within an input raster dataset, perform a specified function, and then write the output to a new raster using the PixelBlock. The example not only demonstrates how to access the individual values for each cell but also shows through neighborhood notation (as discussed in the ArcGIS Desktop Help topic 'An overview of Neighborhood tools' within the Spatial Analyst functional reference) how to access the neighboring cell values around the processing cell (which is accomplished by specifying the neighboring cell address such as (-1, -1) is the bottom left diagonal cell from the processing cell). To access the value of the processing cell (not its neighbors) in the following example use pInputBlock.GetVal(0, I, J). The example is long, but demonstrates a powerful capability in ArcGIS. By slightly changing this example, any number of new functions can be created or any filters specified.

[Visual Basic 6.0]
Public Sub NeighborhoodNotation(sInputPath As String, sInputFileName As String, _
                  sOutputPath As String, sOutputFileName As String)
  ' sInputPath: path of the input raster dataset
  ' sInputFileName: name of input raster dataset
  ' sOutputPath: path of the output raster dataset
  ' sOutputFileName: name of the output raster dataset

  ' This example reproduces the result you would get by
  ' by running the following at the ArcGRID prompt:
  ' Grid: outgrid = ingrid + (ingrid(-1,-1) - ingrid(1,1))
  ' It demonstrates how one would use the PixelBlock
  ' object to perform neighborhood notation.

  ' Open input raster dataset
  Dim pRWS As IRasterWorkspace2
  Dim pWSF As IWorkspaceFactory
  Set pWSF = New RasterWorkspaceFactory
  If Not pWSF.IsWorkspace(sInputPath) Then Exit Sub
  Set pRWS = pWSF.OpenFromFile(sInputPath, 0)
  Dim pInputRasterDS As IRasterDataset
  Set pInputRasterDS = pRWS.OpenRasterDataset(sInputFileName)

  ' Get RasterBand from the input raster and QI raster
  ' properties interface from input raster dataset
  Dim pInputBand As IRasterBand
  Dim pInputBandCol As IRasterBandCollection
  Set pInputBandCol = pInputRasterDS
  Set pInputBand = pInputBandCol.Item(0)
  Dim pInputRasProps As IRasterProps
  Set pInputRasProps = pInputBand

  ' Create a default raster from input raster dataset
  Dim pInputRaster As IRaster
  Set pInputRaster = pInputRasterDS.CreateDefaultRaster

  ' Create output raster dataset
  If Not pWSF.IsWorkspace(sOutputPath) Then Exit Sub
  Set pRWS = pWSF.OpenFromFile(sOutputPath, 0)
  Dim pOrigin As IPoint
  Set pOrigin = New Point
  pOrigin.X = pInputRasProps.Extent.XMin
  pOrigin.Y = pInputRasProps.Extent.YMin
  Dim pOutputRasterDS As IRasterDataset
  Dim outputPixelType As rstPixelType
  Select Case pInputRasProps.PixelType
    Case PT_CHAR
      outputPixelType = PT_LONG
      MsgBox "Operation not supported on complex data", _
          VbOKOnly, "Neighborhood Notation"
      Exit Sub
       MsgBox "Operation not supported on complex data", _
           vbOKOnly, "Neighborhood Notation"
       Exit Sub
    Case PT_DOUBLE
      outputPixelType = PT_FLOAT
    Case PT_FLOAT
      outputPixelType = PT_FLOAT
    Case PT_LONG
      outputPixelType = PT_LONG
    Case PT_SHORT
      outputPixelType = PT_LONG
    Case PT_U1
      outputPixelType = PT_LONG
    Case PT_U2
      outputPixelType = PT_LONG
    Case PT_U4
      outputPixelType = PT_LONG
    Case PT_UCHAR
      outputPixelType = PT_LONG
    Case PT_ULONG
      outputPixelType = PT_LONG
    Case PT_USHORT
      outputPixelType = PT_LONG
  End Select
  Set pOutputRasterDS = pRWS.CreateRasterDataset(sOutputFileName, "GRID", pOrigin, _
				  pInputRasProps.Width, pInputRasProps.Height, _
        				  pInputRasProps.MeanCellSize.X, _
                                         pInputRasProps.MeanCellSize.Y, 1, outputPixelType, _
                                         pInputRasProps.SpatialReference, True)

  ' Create a default raster from output raster dataset
  Dim pOutputRaster As IRaster
  Set pOutputRaster = pOutputRasterDS.CreateDefaultRaster

  ' Get RasterBand from the output raster
  Dim pOutputBand As IRasterBand
  Dim pOutputBandCol As IRasterBandCollection
  Set pOutputBandCol = pOutputRasterDS
  Set pOutputBand = pOutputBandCol.Item(0)

  ' QI RawPixel interface for input
  Dim pInputRawPixel As IRawPixels
  Set pInputRawPixel = pInputBand

  ' QI RawPixel interface for output
  Dim pOutputRawPixel As IRawPixels
  Set pOutputRawPixel = pOutputBand

  ' Create a DblPnt to hold the PixelBlock size
  Dim pPnt As IPnt
  Set pPnt = New DblPnt
  pPnt.SetCoords pInputRasProps.Width, pInputRasProps.Height

  ' Creates PixelBlock of input with defined size
  Dim pInputBlock As IPixelBlock
  Set pInputBlock = pInputRawPixel.CreatePixelBlock(pPnt)

  ' Creates PixelBlock of output with defined size
  Dim pOutputBlock As IPixelBlock
  Set pOutputBlock = pOutputRawPixel.CreatePixelBlock(pPnt)

  ' Read input PixelBlock
  pPnt.X = 0
  pPnt.Y = 0
  pInputRawPixel.Read pPnt, pInputBlock

  ' Get the SafeArray associated with the first band
  ' of output
  Dim vSafeArray As Variant
  vSafeArray = pOutputBlock.SafeArray(0)

  ' QI RasterProps for output for NoData handling
  Dim pOutputRasProps As IRasterProps
  Set pOutputRasProps = pOutputBand

  ' Loop through the SafeArray and calculate each pixel value
  ' using equation given above
  Dim I, J As Long
  For I = 0 To pInputRasProps.Width - 1
    For J = 0 To pInputRasProps.Height - 1
      If I - 1 >= 0 And I + 1 <= pInputRasProps.Width - 1 And J - 1 >= 0 _
         And J + 1 <= pInputRasProps.Height - 1 Then
        If Not CDbl(pInputBlock.GetVal(0, I, J)) = CDbl(pInputRasProps.NoDataValue) _
           And Not CDbl(pInputBlock.GetVal(0, I - 1, J - 1)) = _
           CDbl(pInputRasProps.NoDataValue) And Not _
           CDbl(pInputBlock.GetVal(0, I + 1, J + 1)) = _
           CDbl(pInputRasProps.NoDataValue) Then
          vSafeArray(I, J) = CDbl(pInputBlock.GetVal(0, I, J)) + _
         			            (CDbl(pInputBlock.GetVal(0, I - 1, J - 1)) - _
         			            CDbl(pInputBlock.GetVal(0, I + 1, J + 1)))
          vSafeArray(I, J) = CDbl(pOutputRasProps.NoDataValue)
       End If
        vSafeArray(I, J) = CDbl(pOutputRasProps.NoDataValue)
     End If
   Next J
 Next I

  ' Write out the result
  pOutputRawPixel.Write pPnt, pOutputBlock

End Sub

Using other ArcGIS objects with Spatial Analyst

The Spatial Analyst objects rely on obtaining input data from the Raster Data Object model as well as other object models available in ArcGIS. In Spatial Analyst methods, a RasterDataset, RasterBand, Raster, and RasterDescriptor are valid inputs as well as FeatureClass and FeatureClassDescriptor. Any of these inputs should provide acceptable output behavior. In the 'Accessing the Input Data' section in the steps for implementing a Spatial Analyst operation discussed above, it was demonstrated how to obtain a raster from disk to create a RasterDataset using the RasterWorkspace object.

Besides creating a RasterDataset from data on disk, the second way of getting input for Spatial Analyst operators is to get a Raster object from a RasterLayer. The following VBA code shows how to access a RasterLayer from a data frame in ArcMap and obtain a Raster from it that can be used in a Spatial Analyst operator.

[Visual Basic 6.0]
' This function will select the top raster layer in a map document
' If no raster layers are found in the map, function will return empty
Function GetRasterFromMap() As IRaster
  On Error GoTo ERH
  Dim mpMxDoc As IMxDocument
  Dim lCntr As Long
  Dim pLayer As ILayer
  Dim pRasLyr as IRasterLayer
  Set mpMxDoc = ThisDocument
  For lCntr = 0 To mpMxDoc.FocusMap.LayerCount - 1
    Set pLayer = mpMxDoc.FocusMap.Layer(lCntr)
    If TypeOf pLayer Is IRasterLayer Then
      Set pRasLyr = pLayer
      Set GetRasterFromMap = pRasLyr.Raster
      Exit Function
    End If
  Next lCntr
  Set GetRasterFromMap = Nothing
  Exit Function
  Set GetRasterFromMap = Nothing
  MsgBox Err.Number & Err.Description, , "Error: GetRasterFromMap"
End Function

Another input to most operators are feature class data. In some operators the feature class data is used directly as in the RasterIntepolationOp and in most others the feature class data will be converted automatically to a raster using the analysis environment settings. The following code accesses a shapefile dataset using the Geodatabase object model in ArcGIS.

[Visual Basic 6.0]
Public Function OpenFeatureClassFromShapefile(sPath As String, _
        sShapeName As String) As IFeatureClass
  ' Returns a FeatureClass object given a shapefile's name and it's path
  On Error GoTo ERH
  Dim pWSFact As IWorkspaceFactory
  Dim pFWS As IFeatureWorkspace

  Set pWSFact = New ShapefileWorkspaceFactory
  If pWSFact.IsWorkspace(sPath) Then
    Set pFWS = pWSFact.OpenFromFile(sPath, 0)
    Set OpenFeatureClassFromShapefile = pFWS.OpenFeatureClass(sShapeName)
  End If
  Exit Function
  Set OpenFeatureClassFromShapefile = Nothing
  MsgBox "Failed in creating FeatureClass " & Err.Description
End Function

Another common non-raster input comes from ArcInfo coverages as demonstrated in the following code.

[Visual Basic 6.0]
Public Function OpenFeatureClassFromCoverage(sPath As String, sCovName As String, _
     ClassID As Integer) As IFeatureClass
  ' Open a feature class in the specified coverage by giving a class ID
  On Error GoTo ERH
  ' Open a feature dataset
  Dim pFeatDs As IFeatureDataset
  Dim pFeatWs As IFeatureWorkspace
  Dim pWSF As IWorkspaceFactory
  Set pWSF = New ArcInfoWorkspaceFactory
  If pWSF.IsWorkspace(sPath) Then
    Set pFeatWs = pWSF.OpenFromFile(sPath, 0)
    Set pFeatDs = pFeatWs.OpenFeatureDataset(sCovName)
    ' Get a feature class from dataset
    Dim pFC As IFeatureClassContainer
    Set pFC = pFeatDs
    Set OpenFeatureClassFromCoverage = pFC.Class(ClassID)
  End If
  Exit Function
  Set OpenFeatureClassFromCoverage = Nothing
  MsgBox "Failed in open featureclass from coverage " & Err.Description
End Function