[R-sig-Geo] neighbourhood analysis with R spatial objects

[Alexandre Villers]

Dear list,

I am struggling a bit with this issue and was hoping that someone could 
indicate me a solution.

I have a raster image describing the landscape (with several habitat 
classes) at 10 meters resolution. I would like to compute the different 
proportions of each habitat type, but on a coarser grid (let's say 250 
meters), each cell getting the proportions of habitat classes from 
different buffer radii, taking the centre of each grid cell as 
representative of the grid (might not be very accurate for small 
buffers, but with the radii increasing, it should quickly be OK).

I used to do similar tasks with the r.neighbors function in GRASS but 
the problem is that this time, I need to run it on the 10 meters 
resolution grid (because some habitats are very scarce but a priori 
relevant to the questions my colleague and I want to answer), meaning a 
long computation time (already several days at 25 meters resolution), 
most of which being useless (as r.neighbors in GRASS compute the 
neighbourhood on all the pixels).

For the moment, I came up with this the code below, but this is not 
efficient at all (57 seconds for a final grid of 400 cells as in the 
example, mine containing ~20000 pixels...). Moreover, as written 
earlier, I have 11 buffer radii and 16 images on which to batch.
Basically, I create a regular grid of points at 250 meters over the 
SpatialGridDataFrame at 10 meters resolution, and then, loop over each 
points (create a buffer, overlay with the raster to get the number of 
pixels of each habitat class, fill the final raster with the proportion 
after computations and save the object).

I'm not sure if this can be done more efficiently with GRASS (computing 
the neighbourhood just for a set of locations, not for all pixels) and I 
was hoping to be able to do that with R somehow. I couldn't find 
anything in the discussion  lists (with keywords "neighbourhood 
analysis", "buffer raster", etc.) but maybe I missed something (with te 
different possible approaches in sp, raster, spatstat, it is highly 
probable).

Any link towards a potential solution that could speed up the 
computations would be highly appreciated


Cheers

Alex

P.S.: the code below should give you an idea of what I would like to 
achieve (but much faster if possible !)


library(sp)


#Create the environmental covariate
extent<-5000 #extent of the raster
sp.res<-10 #resolution of the raster

bb<-matrix(c(0,0,extent,extent), ncol=2, dimnames=list(NULL, 
c("min","max")))
cs<-c(sp.res,sp.res) #the resolution of our grid of habitat
cc<-bb[,1] + (cs/2)
cd<-ceiling(diff(t(bb))/cs)
landGT.temp<-GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
land.spg<-SpatialGrid(landGT.temp)
Habitat<-rpois(length(coordinates(land.spg)[,1]), 5) #create dummy 
variable of integers for landscape category
land.spgdf<-SpatialGridDataFrame(land.spg, data.frame(Habitat=Habitat))


#create the final grid on which we want the data
sp.res.p<-250 #resolution of the raster

bb.p<-matrix(c(0,0,extent,extent), ncol=2, dimnames=list(NULL, 
c("min","max")))
cs.p<-c(sp.res.p,sp.res.p) #the resolution of our grid of habitat
cc.p<-bb.p[,1] + (cs.p/2)
cd.p<-ceiling(diff(t(bb.p))/cs.p)
fgridGT.temp<-GridTopology(cellcentre.offset=cc.p, cellsize=cs.p, 
cells.dim=cd.p)
fgrid.spg<-SpatialGrid(fgridGT.temp)
fgrid.spp<-as(fgrid.spg, "SpatialPoints")




#the idea is to get the proportion of the different habitats described 
by land.spgdf in buffer of increasing radii
# such as in the example

image(land.spgdf)
plot(fgrid.spp[210,], col="black", add=T)
plot(gBuffer(fgrid.spp[210,], width=1000), add=T)
plot(gBuffer(fgrid.spp[210,], width=250),  add=T)



n.points<-length(coordinates(fgrid.spp)[,1])
buffers<-c(250,1000)
n.hab.class<-length(unique(Habitat))

system.time(for (j in 1:length(buffers)){ #
          temp.df<-as.data.frame(matrix(data=NA, nrow=n.points, 
ncol=n.hab.class))
         spgdf.res<-SpatialGridDataFrame(fgrid.spg, temp.df)
         res.temp<-NULL

         for (k in 1:n.points){ # n.points
             buff.temp<-gBuffer(fgrid.spp[k,], width=buffers[j]) #create 
the buffer for the kth point with radius buffers[j]
             assol.buff<-land.spgdf$Habitat[!is.na(over(land.spgdf, 
buff.temp))] #overlay and remove the useless NAs
             n.pix<-tapply(assol.buff, assol.buff, 
length)[as.character(sort(unique(Habitat)))][] #number of pixel per 
habitat class (sorted)
             res.temp<-rbind(res.temp,n.pix) #store the data in a 
temporary array
                     }
         spgdf.res at data[1:n.points,]<-res.temp; 
spgdf.res at data[is.na(spgdf.res at data[,1:21])]<- 0 #replace the NAs 
(absence of the habitat class) with zeros
         spgdf.prop<-spgdf.res; spgdf.prop at data[,1:21]<-NA; 
spgdf.prop at data[,1:21]<-t(apply(spgdf.res at data, 1, function(x) 
x/sum(x))) #the last to get the proportion of habitat
         name.temp<-paste("prop.",buffers[j],"m",sep="")
         assign(name.temp, spgdf.prop) #create the final R object (that 
could saved as well as GTiff or any other raster format)
         rm(spgdf.res, name.temp, spgdf.prop, temp.df)
           }
         )

#An example
spplot(prop.250m, "V5")