How to buffer raster pixels by their values?

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? • Did you try to convert raster to points, then buffer by field then convert back to raster? – user681 Feb 10 '14 at 17:46
• It helps to realize that this is a non-local operation, because that non-locality shows there are inherent limits on how the calculation can be carried out. For instance, your output would radically change almost everywhere if just a single isolated pixel in the input were to change to a large value. Thus, if you know of restrictions on the input values, then please share them, because that can lead to improved solutions. For instance, will all your input values always be in the set {2,3,4}? – whuber Feb 10 '14 at 17:46
• @Dan Patterson That is how I came up with the image to the right. However, I am trying to avoid vector operations altogether and avoid those steps. – Aaron Feb 10 '14 at 17:50
• @whuber This dataset represents trees with variable crown diameters. Given that, the tree crown radius measurements can realistically vary from 1-10. I should also add that the buffered output need only be 0's for crown absence and 1's for crown presence. – Aaron Feb 10 '14 at 17:54
• OK, thanks. It looks like you created the example output by unioning the 3-buffers of the points with value 3, the 4-buffers of the points with value 4, and the 5-buffers of the points with value 5. (You seem to have forgotten to process the points with value 2.) That process not only answers your question, but (I believe) it is the simplest solution using the tools available in Spatial Analyst. – whuber Feb 10 '14 at 18:07

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
mask = x**2 + y**2 <= radius**2
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

fig = plt.figure()
plt.imshow(C, interpolation='none')
plt.show()

Input: Output: • +1 for the dilation approach! It works with near points too. – Antonio Falciano Feb 13 '14 at 16:01
• This is a great example of why that old jet colourscheme was terrible. This looks much clearer with viridis. – naught101 May 1 at 1:11

Vector-based approach

This task can be done in three steps:

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.

• Thanks for the response afalciano. This method is how I created the image to the right and uses a vector approach--one I am trying to avoid. – Aaron Feb 10 '14 at 18:59
• Ok Aaron, here's a raster-based approach now. Hope this helps. – Antonio Falciano Feb 11 '14 at 9:29

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.

• This idea's the right one. The details can be improved: (1) the selection creates a binary (0,1) indicator of the trees of a given crown radius. (2) A focal sum of that selection--using a circular neighborhood of the given radius--is fast to compute using the FFT. (3) Adding the focal sums (pointwise) and comparing that to 0 gives the desired buffer. – whuber Feb 10 '14 at 22:57

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)

• +1 This is a creative approach. I will test to see if it is feasible to scale up using this approach. – Aaron Feb 10 '14 at 19:06
• The Euclidean distance approach won't work, because it computes only the distance to the nearest tree, which is not necessarily the distance to a tree whose crown covers the point. – whuber Feb 10 '14 at 22:52

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.

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.

• I know it is not a matter of the question but, do you want to know more about graphic primitives and scan line polygon fill? For cirles it is very easy. Convex combination is a buzzword and so on.... – huckfinn Feb 10 '14 at 21:50
• How would this be applied using basic raster operations? – whuber Feb 10 '14 at 22:54
• If you try to handle this in the raster space, determin the circle points, sorting them by the y or x and fill the space by a straight line (scan line) is on way to fill the circle. In the triangular approach, if you build the circle by an approximation of tringular sectors and try to fill the triangle you need a test if the point is inside or outside (convex combination) and is the other way. And in the "GIS" appoach, building polygons (clock wise oriented polygons) and make a union is the third (IMO the most computual expensive one). – huckfinn Feb 10 '14 at 23:03
• To be clear: And in the "GIS" appoach... make an algebraric operation like union, intersection, touch.... is the third IMO the most computual expensive one. – huckfinn Feb 10 '14 at 23:16