I have a binary raster (0s and 1s) that I need to query given a set of point coordinates. This is easy enough to do with gdal (gdallocationinfo), but where it get's complicated is that I also need to return the sum of the 8-cell neighborhood. Is there an easy way to do this without processing the entire map?

In GRASS the only way I can think to do this is to use r.mapcalc to return a sum of the cell I'm interested in and then change the value of all other cells to 0. This allows me to compute an aggregate sum for the map with r.sum. But this requires a great deal of processing and if I need to do this for hundreds of points, the task becomes enormous.


  • the specific language used to solve the problem isn't that important. Commented Mar 14, 2013 at 19:15
  • I'm not sure I completely understand what you want, but wouldn't precomputing a focal sum do the trick? Then all you have to do is query the focal sum grid.
    – whuber
    Commented Mar 14, 2013 at 21:01
  • Yes, it would. However, I need to be able to repeat this process thousands of times, with map data changing frequently. If I'm only interested in one point, it seems silly to compute focal sums for the whole map just to get at one 8 cell neighborhood. Commented Mar 14, 2013 at 23:08
  • It was suggested to me by someone else that I look at v.random.cover. Have you had experience with this @whuber? Commented Mar 14, 2013 at 23:09
  • Could you explain further how you have map data that are "changing frequently"? Incidentally, if you really just need to poke once at each grid, why not extract the neighborhood and compute the sum? That's almost no processing at all--it's essentially I/O-bound.
    – whuber
    Commented Mar 15, 2013 at 0:01

3 Answers 3


Here's some gdal python code to get NxN numpy arrays from point coordinates:

def getwindow(ds,band,coords,window=[3,3]):
    rb = ds.GetRasterBand(band)
    gt = ds.GetGeoTransform()

    for mx,my in (coords):

        #Convert from map to pixel coordinates.
        #Only works for georansforms with no rotation.
        #Does no checking that coordinates will fall within raster extent so can return negative values or values > ncols or nrows
        px = int((mx - gt[0]) / gt[1]) #x pixel
        py = int((my - gt[3]) / gt[5]) #y pixel

        #Shift pixel coords

        #copy for modification

        #Handle edge/corner cases by shrinking the kernel. 
        #No fancy reflection/padding/wrapping.
        if spx < 0:
        elif px > ncols-cx-1:
        if spy < 0:
        elif py > nrows-cy-1:

        yield rb.ReadAsArray(spx,spx,shape[0],shape[1])

ds = gdal.Open( someraster )
for arr in getwindow(ds,band,coords,window):
    print arr.sum()
    print arr.shape #Note: can be smaller than NxN if it is a corner/edge case

In GRASS GIS, use r.neighbors which calculates new values as a function of the category values assigned to the cells around it, and stores new cell values in an output raster map layer. It has a sum operator (focal sum). Then query the resulting map.


If you are familiar with R, use:

column <- extract(RASTER, Samplingsites.spdf, buffer=...)
SamplingSites.df <- cbind(SamplingSites.df, column)

Where RASTER is your binary raster, and the SpatialPointsDataFrame are your set of point coordinates.

You can even do this from within grass with spgrass6 to speed up the process.

If you want to stay within grass I would assume you use a workflow with:

  1. r.neighbors with your set of points(it calculates your sum with a moving window, size:3x3 cells ... so your point and 8 additional points)
  2. v.what.rast to do a spatial query

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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