If you’ve read through the Introducing the Raster Cell Iterator blog, you’ll be familiar with the basic concepts of the Raster Cell Iterator (RCI), a new ArcPy class added in ArcGIS Pro 2.5 in the Spatial Analyst module that allows you to iterate over raster cells, query and modify individual cell values.
In this blog, let’s explore how RCI can be used to customize your analysis on problems that cannot be easily solved using the existing out-of-the-box geoprocessing tool. In the first example, we will use RCI to create a new raster with a specific spatial pattern – checkerboard. And in two other examples, we will use RCI to customize two focal operations: 1) count the number of neighboring cells that have a different value from the center cell, and 2) calculate the minimum slope from a DEM.
Create a new raster with the checkerboard pattern
In atmospheric modeling, a heterogenous surface raster represented by a checkerboard pattern can be used as an initial condition for idealized model simulations. However, there is not an existing geoprocessing tool with which you can create such a raster directly. In earlier ArcGIS releases, you can accomplish this task by using NumPy, but it would require a multi-step workflow to do so. In ArcGIS Pro 2.5 with RCI, this task has become much easier! You simply need to assign a binary value to each cell based on the current row and column while iterating through a raster. Following is a simple code sample which does exactly that.
import json
from arcpy.sa import RasterInfo, Raster
myRasterInfo = RasterInfo()
myRasterInfoData = {
'bandCount': 1,
'extent': {
'xmin': -107.0,
'ymin': 38.0,
'xmax': -104.0,
'ymax': 40.0,
'spatialReference': {'wkid': 4326},
},
'pixelSizeX': 0.01,
'pixelSizeY': 0.01,
'pixelType': 'U8',
}
myRasterInfo.fromJSONString(json.dumps(myRasterInfoData))
myRaster = Raster(myRasterInfo)
for (r, c) in myRaster:
if r % 20 < 10 and c % 20 < 10 or r % 20 >= 10 and c % 20 >= 10:
myRaster[r, c] = 1
else:
myRaster[r, c] = 0
myRaster.save(r 'C:\output\checkerboard.tif')
From the code snippet above, lines 8-23 show the steps to construct a RasterInfo object from scratch. Line 26 creates an empty raster based on the raster info specified in the previous step. Lines 28-33 show the main logic of creating a checkerboard pattern raster by assigning binary values to the cell based on the current row and column. An implicit way of using Raster Cell Iterator is demonstrated by calling a simple For loop to go through all the cells. Since there is only one raster involved in this analysis, you don’t need to create a RasterCellIterators object explicitly. It is recommended to follow this pattern when you are working with only one raster. Finally, at line 35, we are saving the output raster. Figure 1 shows the output raster with the checkerboard pattern.
Figure 1. Output raster with the checkerboard pattern
Custom Focal Operations
Focal operations produce an output raster dataset in which the output value at each cell is a function of the input cell value at that location and its neighboring cell values. With RCI, you can customize these kinds of focal operations with more flexibility.
Count the number of neighboring cells that have a different value from the center cell
Stathakis & Tsilimigkas (2015) proposed several metrics to calculate the compactness of cities using land use adjacency information. To derive the compactness metric, one of the pre-processing steps is to iterate through the land use data and count the number of neighboring cells that have a value different to the center cell. This focal operation process is illustrated in Figure 2.
Figure 2. Examples of input and output with the operation counting the number of neighboring cells that have a value different to the center cell.
RCI makes it much easier to iterate over all cells and compare each cell value with its neighbors using a simple if logic. Here’s how you can do it yourself:
from math import isnan
from arcpy.sa import Raster, RasterCellIterator
landu = Raster(r"C:\sample_data\landuse.tif")
output = Raster(landu.getRasterInfo())
with RasterCellIterator({'rasters':[landu, output]}) as rci:
for i,j in rci:
count = 0
if isnan(landu[i,j]):
output[i,j] = math.nan
continue
for x in [-1,0,1]:
for y in [-1,0,1]:
if not isnan(landu[i+x,j+y]) and landu[i+x,j+y] != landu[i,j]:
count+=1
output[i,j] = count
output.save(r"C:\output\landuse_diff_count.tif")
In line 6, we are creating a new empty raster using the land use raster as an input template. Line 8 creates a RasterCellIterator object explicitly using input and output. In lines 9-21, we are looping through all cells in the input raster, and for each cell comparing its value to each of its immediate neighbors based on its relative position.
While performing the analysis, it is required to check for the NoData cells in the input, so that an appropriate value can be assigned to the output raster. In lines 12 and 19, we are using isnan() function from the math module to check if the cell value in the input raster is NoData. Line 13 uses nan from the same math module to represent NoData. In this example, NoData values are not considered while iterating through the cells and they are excluded from the neighbor cells when counting. The input land use data and the final output raster used in the sample code are shown in Figure 3.
Figure 3. The figure on the left shows the input land use raster dataset. The figure on the right shows the output raster with values ranging from 0 to 8, which counts the number of neighboring cells that have a different value from the center cell.
Calculate minimum slope from a DEM
For the Slope geoprocessing tool, the slope value of each cell is determined by the rates of change of the elevation in the horizontal and vertical directions from the center cell (How slope works). In some cases, you might be interested in calculating the slope using a different equation. For an example, when building a fish trawlability model, an analyst might want to use a slope of the seabed that is based on the minimum change of elevation across each cell. By creating a custom focal operation with RCI, the following example shows how to calculate the minimum slope for your bathymetric DEM:
from arcpy.sa import Raster, RasterCellIterator
dem = Raster(r'C:\sample_data\demtif.tif')
cell_x = dem.meanCellHeight
cell_y = dem.meanCellWidth
raster_info = dem.getRasterInfo()
raster_info.setPixelType('F32')
min_slope = Raster(raster_info)
with RasterCellIterator({'rasters':[dem, min_slope]}) as rci:
for r,c in rci:
slopes = []
for x,y in [(-1,-1),(-1,0),(1,0),(0,-1),(0,1),(1,1),(-1,1),(1,-1)]:
slope = abs(dem[r,c]-dem[r+x,c+y])/math.sqrt(cell_x**2+cell_y**2)
slopes.append(slope)
min_slope[r,c] = min(slopes)
min_slope.save(r'C:\output\min_slope.tif')
In the code sample, lines 13-21 implement the main logic for calculating the minimum slope. For each cell from DEM, we iterate through its eight neighbors and calculate the rate of change. Then we find the minimum slope and assign that value to the output cell. Note that in line 9, the pixel type is set as ‘F32’, which means that the output is expected to be a 32-bit floating raster. In this example, the raster info of the slope output is inherited from the input DEM data. If the pixel type of your DEM data is integer, you need this step to make sure that the slope value of each cell is calculated as a floating-point value instead of an integer.
Summary
In this blog, we went through three different examples of using Raster Cell Iterator in ArcPy, each of them showing how easy it is to perform custom raster analysis. We hope you start using this simple but powerful capability to extend your own raster analysis capabilities. Please let us know if you have any questions or want to see other examples and use cases. You can also pass along any cool applications you come up with, and we can share your solutions with the community at large!
Reference:
Stathakis, D., & Tsilimigkas, G. (2015). Measuring the compactness of European medium-sized cities by spatial metrics based on fused data sets. International Journal of Image and Data Fusion, 6(1), 42-64.
The original blog post Unleash the power of RasterCellIterator to perform custom raster analysis has been written by Hao Hu and Yu Wang.
