[GIS] How to buffer raster pixels by their values

arcgis-desktoparcpybufferpythonraster

The pixels to the left represent tree locations and their associated crown radii (i.e. pixel values ranging from 2 – 5). I would like to buffer these raster pixels by their crown radius value. The image to the right is what I am hoping to accomplish using only raster processing methods.

I would initially think to use a circular focal sum in ArcGIS, although the neighborhood setting is a fixed value, which would not take into account the variable sized crown radius.

What is a good method to "buffer" pixels by their values?

enter image description here

Best Answer

Here is a pure raster solution in Python 2.7 using numpy and scipy:

import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt

#create tree location matrix with values indicating crown radius
A = np.zeros((120,320))
A[60,40] = 1
A[60,80] = 2
A[60,120] = 3
A[60,160] = 4
A[60,200] = 5
A[60,240] = 6
A[60,280] = 7

#plot tree locations
fig = plt.figure()
plt.imshow(A, interpolation='none')
plt.colorbar()

#find unique values
unique_vals = np.unique(A)
unique_vals = unique_vals[unique_vals > 0]

# create circular kernel
def createKernel(radius):
    kernel = np.zeros((2*radius+1, 2*radius+1))
    y,x = np.ogrid[-radius:radius+1, -radius:radius+1]
    mask = x**2 + y**2 <= radius**2
    kernel[mask] = 1
    return kernel

#apply binary dilation sequentially to each unique crown radius value 
C = np.zeros(A.shape).astype(bool)   
for k, radius in enumerate(unique_vals):  
    B = ndimage.morphology.binary_dilation(A == unique_vals[k], structure=createKernel(radius))
    C = C | B #combine masks

#plot resulting mask   
fig = plt.figure()
plt.imshow(C, interpolation='none')
plt.show()

Input: enter image description here

Output: enter image description here

Related Question