# query raster cell values by point

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.

Ideas?

• the specific language used to solve the problem isn't that important. – ShaunLangley Mar 14 '13 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 Mar 14 '13 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. – ShaunLangley Mar 14 '13 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? – ShaunLangley Mar 14 '13 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 Mar 15 '13 at 0:01

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()
ncols=ds.RasterXSize
nrows=ds.RasterYSize

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) / gt) #x pixel
py = int((my - gt) / gt) #y pixel

#Shift pixel coords
cx=window//2
cy=window//2
spx=px-cx
spy=py-cy

#copy for modification
shape=window[:]

#Handle edge/corner cases by shrinking the kernel.
if spx < 0:
shape+=spx
spx=0
elif px > ncols-cx-1:
shape=ncols-spx-1
if spy < 0:
shape+=spy
spy=0
elif py > nrows-cy-1:
shape=nrows-spy-1

ds = gdal.Open( someraster )
band=1
window=[3,3]
coords=([x1,y1],[x2,y2],...,[xn,yn])
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