This topic is relevant for the following:
Product(s): ArcGIS Desktop: All, ArcGIS Engine
Version(s): 9.0, 9.1, 9.2
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).
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 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.
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 Function
You 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 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.
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 Sub
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.
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 Sub
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'.
' 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
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.
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 Sub
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.
' 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 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.
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 Function
Another 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