Skip to content
Prev 7000 / 29559 Next

Overlay of SpatialGrid/Raster and polygons

Ned, here is how I have done that using the overlay function. Not sure 
if it would work with a rasterstack. 
Alan

xy.locs <- readOGR(poly.shapefile.name,drop_unsupported_fields=T)
img<-readGDAL(image.name)
y<-slot(img,"data")
n.bands <- ncol(y)
bnames <- paste(varname,"_b",1:n.bands,sep="") }
out.frame <- 
as.data.frame(matrix(NA,nrow=nrow(xy.locs),ncol=n.bands,dimnames=list(NULL,bnames)))   

for(m in 1:(n.poly <- length(slot(xy.locs,"polygons")))){
    gc()
    z<-overlay(img,xy.locs[m,])
    for(j in 1:n.bands){
        xx <- mean(y[z==1 & !is.na(z),j],na.rm=T)
        out.frame[m,j] <- xx
        } # end inner loop over bands #
    } # end loop over individual polygons #
Ned Horning wrote: