EDN: ESRI Developer Network
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.
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:
The following text discusses these five steps and provides Visual Basic examples demonstrating the step.
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.
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
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.
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
pEnv.SetAsNewDefaultEnvironment
' Return reference to the default environment setting
Set SetNewDefaultEnvironment = pEnv
Exit Function
ERH:
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).
pNewEnv.Reset
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.
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
ERH:
Set OpenRasterDataset = Nothing
MsgBox "Failed in Opening RasterDataset. " & Err.Description
End FunctionA 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.
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
ERH:
Set CreateRasterDescriptor = Nothing
MsgBox "Failed in Creating RasterDescriptor " & Err.Description
End FunctionYou can also create a RasterDescriptor or FeatureClassDescriptor using any desired field. See 'Fields' later in this document).
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).
' 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.
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.
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 pActView.Refresh pMxDoc.UpdateContents Exit Sub ERH: 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.
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.
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 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.
' 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)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.
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.
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
ERH:
Set UsingRasterDescriptor = Nothing
End SubThe 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.
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.
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 SubSome 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.
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
Else
GetSecondRaster = Nothing
End If
End SubTo 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'.
' 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.HeightIf 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.
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.
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
Case PT_COMPLEX
MsgBox "Operation not supported on complex data", _
VbOKOnly, "Neighborhood Notation"
Exit Sub
Case PT_DCOMPLEX
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)))
Else
vSafeArray(I, J) = CDbl(pOutputRasProps.NoDataValue)
End If
Else
vSafeArray(I, J) = CDbl(pOutputRasProps.NoDataValue)
End If
Next J
Next I
' Write out the result
pOutputRawPixel.Write pPnt, pOutputBlock
End SubThe 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.
' 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
ERH:
Set GetRasterFromMap = Nothing
MsgBox Err.Number & Err.Description, , "Error: GetRasterFromMap"
End FunctionAnother 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.
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
ERH:
Set OpenFeatureClassFromShapefile = Nothing
MsgBox "Failed in creating FeatureClass " & Err.Description
End FunctionAnother common non-raster input comes from ArcInfo coverages as demonstrated in the following code.
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
ERH:
Set OpenFeatureClassFromCoverage = Nothing
MsgBox "Failed in open featureclass from coverage " & Err.Description
End Function