Search code examples
rraster

r raster statistics per grid cell zone

I would like to calculate grid cell statistics in zones that are within graticules (see plot).

enter image description here

Getting the cellstats for the shown grid cell is not difficult, but how could this be done for all grid cells? 

library(raster)

filename <- system.file("external/test.grd", package="raster")
r=raster(filename)

plot(r)
e=extent(180000,181000,330000,331000)
plot(e,add=T)
grid()

x=crop(r,e)
cellStats(x,mean)
Solution

  • Your reproducible example did not work for me as I am not sure why r[r>500]=1 and r[r<=500]=NA are used while your original plot is different?! Plus, I was not sure what rx is?! One way to do this as below: 

    library(raster)

filename <- system.file("external/test.grd", package="raster")
r=raster(filename)

# r[r>500]=1
# r[r<=500]=NA

plot(r)
xmin <- seq(178000, 181000, 1000)
xmax <- seq(179000, 182000, 1000)
ymin <- seq(329000, 333000, 1000)
ymax <- seq(330000, 334000, 1000)

zonal.mtx <- matrix(nrow=length(xmin), ncol=length(ymin))

library(spatialEco)
for (i in 1:length(xmin)){
  for (j in 1:length(ymin)){
    e=extent(xmin[i],xmax[i],ymin[j],ymax[j])
    plot(e,add=T)
    grid()
    p <- as(e, 'SpatialPolygons')  
    crs(p) <- crs(r)

    p$ID=j
    p <- SpatialPolygonsDataFrame(p,p@data) 

    zonal.mtx[i,j] <- zonal.stats(x=p, y=r, stat=sum, trace=TRUE, plot=TRUE)
    #print(zonal.mtx[i,j])
  }

}

    Note: this plots each portion of the data while calculating zonal stats. I am sure you already know that you can also define your own function for using in zonal.stats rather than sum, mean etc.