Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
1  
Did you try to convert raster to points, then buffer by field then convert back to raster? –  Dan Patterson Feb 10 at 17:46
2  
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 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 at 17:50
2  
@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 at 17:54
1  
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 at 18:07
add comment

5 Answers

up vote 3 down vote accepted

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

share|improve this answer
    
Wow...great solution! Thanks for the help. –  Aaron Feb 12 at 17:17
1  
+1 for the dilation approach! It works with near points too. –  afalciano Feb 13 at 16:01
add comment

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)

share|improve this answer
    
+1 This is a creative approach. I will test to see if it is feasible to scale up using this approach. –  Aaron Feb 10 at 19:06
1  
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 at 22:52
add comment

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.

share|improve this answer
1  
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 at 22:57
add comment

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.

share|improve this answer
1  
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 at 18:59
1  
Ok Aaron, here's a raster-based approach now. Hope this helps. –  afalciano Feb 11 at 9:29
2  
+1 The CA idea is very nice. –  whuber Feb 11 at 19:59
add comment

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. Bye Huck

share|improve this answer
    
+1 Another creative approach. –  Aaron Feb 10 at 20:56
    
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 at 21:50
    
How would this be applied using basic raster operations? –  whuber Feb 10 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 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 at 23:16
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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