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.

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?

share|improve this question
    
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
add comment

3 Answers 3

up vote 2 down vote accepted

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[0]) / gt[1]) #x pixel
        py = int((my - gt[3]) / gt[5]) #y pixel

        #Shift pixel coords
        cx=window[0]//2
        cy=window[1]//2
        spx=px-cx
        spy=py-cy

        #copy for modification
        shape=window[:]

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

        yield rb.ReadAsArray(spx,spx,shape[0],shape[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
share|improve this answer
    
Thanks! I think I can implement this into my code! –  ShaunLangley Mar 15 '13 at 17:34
add comment

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.

share|improve this answer
add comment

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
share|improve this answer
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.