Using ArcGIS Pro raster calculator to give cell value based on cells in its neighborhood?

Question

I have a land cover dataset, I want to find forested cells that are adjacent (touch) agricultural cells.

Is there a way I can achieve this with the ArcGIS Pro raster calculator?

Answer 1

First, let's start with a function that will take your land cover dataset (an N x M) raster, and make it an N x M x 9 raster where the last axis contains the neighbors of each original cell:

import numpy

def stack_neighbors(padded, radius=1):
    # new rows and cols count
    M, N = padded.shape

    # width is typically 3 -- length of each of
    # block that defines the neighbors
    width = radius * 2 + 1
    row_length = N - width + 1
    col_length = M - width + 1

    # Linear indices for the starting width-by-width block
    idx1 = numpy.arange(width)[:, None] * N + numpy.arange(width)

    # Offset (from the starting block indices) linear indices for all the blocks
    idx2 = numpy.arange(col_length)[:, None] * N + numpy.arange(row_length)

    # Finally, get the linear indices for all blocks
    all_idx = idx1.ravel()[None, None, :] + idx2[:, :, None]

    # Index into padded for the final output
    blocked = padded.ravel()[all_idx]

    return blocked

OK, now let's make a fake land cover dataset, where 8 is the code for forest, 7 is agg, then pad it in "edge" mode such that the edges and corners all still have nine neighbors. in this example, I expect all of the 7's to make the cut except for those in the upper right corner and in the first column, the third row from the top as displayed here.

land_cover = numpy.array([
    [5, 0, 3, 3,  7],
    [9, 7, 5, 2,  4],
    [7, 7, 8, 8,  7],
    [1, 7, 8, 8,  8],
])

padded = numpy.pad(land_cover, pad_width=1, mode='edge')
print(padded)

And that gives us:

[[5 5 0 3 3 7 7]
 [5 5 0 3 3 7 7]
 [9 9 7 5 2 4 4]
 [7 7 7 8 8 7 7]
 [1 1 7 8 8 8 8]
 [1 1 7 8 8 8 8]

So now let's stack the neighbors with our function:

blocked = _stack_neighbors(padded)
print(blocked.shape)
# (4, 5, 9)

Now we can use numpy's efficient array operations again to check the status of each cell and its neighbors:

is_agg = land_cover == 7
touches_forest = (stacked == 8).any(axis=-1)
agg_next_to_forest = is_agg & touches_forest
print(agg_next_to_forest.astype(int)) # only convert to int b/c easier to "see"

And I get:

[[0 0 0 0 0]
 [0 1 0 0 0]
 [0 1 0 0 1]
 [0 1 0 0 0]]

Answer 2

Not the most efficient, but the easiest way is to

• first, run the expand tool on your agricultural cells

• second, run the raster calculator to find pixels of the original raster that overlap agricultural pixels in the expanded raster.

In commands, that would give something like this (assuming that your integer values for agri = a1 or a2 and forest = f1 or f2 )

OutRas = Expand(InRas, 1, [a1,a2])

finalRas = Con( ( ((OutRas==a1)or (OutRas==a2)) and ((InRas==f1)or (InRas==f2)) ), 1,0)

You could put everything in one line if you do not store OutRas