0

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?

9
  • I can help you do this in numpy if that would work for you
    – Paul H
    May 7 '19 at 5:30
  • Don't really have a lot of python experience, if you think you can walk me through it though it would be greatly appreciated! May 7 '19 at 5:32
  • 1
    Start by making a true/false raster for agricultural and forest where true is any positive value and false is nodata (SetNull might help here) then you can use a Con Agricultural > 0 and Focal Statistics (rectangle 3 by 3, statistic type of maximum) > 0 - this will give you any cell that is Agricultural type adjacent to a forest type assuming the cell sizes are the same. May 7 '19 at 5:53
  • @MichaelStimson Not really sure what you are describing here, can you clarify? The landuse data is all in one raster, so the cell size is not a concern. You said to start by making a true/false, but said true is any positive value and false is no data, why? Did you mean to say where True is Forested and False is Agricultural? Not really following on raster calc expression you gave either. May 7 '19 at 6:15
  • No, you'll need two rasters, both of which are true/false. Do you have only one value for agricultural and one value for forest? Most land use data I have used splits forest into native/plantation/orchard etc. and agricultural as cropping/grazing etc. You can create a true false raster for each by using Con from your landuse raster then using both determine if the agricultural is positive and one of the surrounding cells in the forest raster is also positive. May 7 '19 at 6:20
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]]
0
1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.