How to buffer raster pixels by their values?

Question

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?

Answer 1

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:

Output:

Answer 2

Vector-based approach

This task can be done in three steps:

• Raster To Point;
• Buffer (using the VALUE field as buffer field);
• Feature To Raster.

Note: using the buffer field avoids the calculation of a buffer for each crown radius value.

---

Raster-based approach

Avoiding the vector-based solution, this problem suggests to use a kind of Cellular Automata based on the nearest neighbours. Assuming that all the black pixels are zeros, the pixels are squared and their size is equal to 1 (or, alternatively, are opportunely scaled), the rules to adopt are very simple:

1. If the pixel value (VALUE) is greater than 1, its value becomes VALUE-1 and then consider its surrounding pixels. If their values are less than VALUE-1, these pixels born or grow and their value becomes VALUE-1. Otherwise, these pixels survive and are left unchanged.
2. If VALUE<=1, do nothing (the pixel is dead!).

These rules have to be applied until all the pixels are dead, i.e. their values is equal to 0 or 1. So N-1 times, where N is the maximum value you have in the input raster. This approach can be quite easily implemented with a bit of Python and numpy.

Answer 3

Another option would be to create separate rasters for each pixel value, in this case 4 rasters, with a condition. Then expand the rasters by a pixel count corresponding to the raster's value (by possibly iterating over a value list). Lastly, join the rasters (either algebraic or spatially), to create one binary raster for the tree crowns.

Answer 4

It is a challenging question to do this in raster because you don't have the opportunity to use the value of the pixel for defining the size of the buffer. Therefore you would need to do the focal filter for each value, as you already said.

Here is a possible answer to do it with only 3 filters (I couldn't find less), but not perfectly as mentioned by Whuber : your buffers will be truncated when trees are close to each others.

1) EDIT : Euclidian allocation (this does not completely solve the issue, as it cuts the buffers in the vicinity of smaller trees, but it is better than the artefacts of my first solution).

2) euclidian distance around each pixel

3) raster calculator (map algebra) with a conditional statement

Con("allocation_result" > ("distance_result" / pixel_size) , 1 , 0)

Note that you can adjust the condtion depending on your needs in terms of radius (with or withot the central pixel)

Answer 5

Wonder why don't you use ArcGIS's expand tool?

import arcpy
from arcpy.sa import *

raster_in  = r'c:\test.tif'
raster_out = r'c:\test_out.tif'

outExpand1 = Expand(raster_in, 2, 2)
outExpand2 = Expand(outExpand1, 3, 3)
outExpand3 = Expand(outExpand3, 4, 4)
outExpand4 = Expand(outExpand4, 5, 5)

outExpand4.save(raster_out)

In case of overlaping: the latest expand command will cover the previous ones.

Answer 6

If you have the pixel position, the radius and the Midpoint circle algorithm (a variant of the Bresenham Alg.) gives you a clue. IMO it is easy to create a polygon from this approach and I think it easy to implement this in Python. A union of this set of polygons gives you the covering area.