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?
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.Sign up to join this community
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]]
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