R statistics per zone using raster::zonal

Question

I would like to calculate zonal statistics in areas that align with the graticules.

After getting a list of extents, the idea is to assign block numbers, because the raster::zonal function requires a raster layer with codes representing the zones.

When I try to fill the extent with a zonal number, the filled area does not correspond with the extent (see plot). Why is that?

library(raster)
library(foreach)

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

xmin = seq(178000, 181000, 1000)
ymin = seq(329000, 333000, 1000)

e=foreach(j=ymin,.combine=c) %:% 
  foreach(i=xmin) %do% {
    e=extent(i,i+1000,j,j+1000)
  }


# this should be going into the foreach loop
n=length(r[e[[1]]])
r[e[[1]]]=rep(1,n)

plot(r)
plot(e[[1]],add=T)

Answer 1

Is this what you wanted? I've modified the extent of raster if you don't want to modify your fishnet.

library(raster)
library(foreach)

filename <- system.file("external/test.grd", package="raster")
r=raster(filename)
extent(r) <- c(178000,182000,329000,334000 )

xmin = seq(178000, 181000, 1000)
ymin = seq(329000, 333000, 1000)

e=foreach(j=ymin,.combine=c) %:% 
  foreach(i=xmin) %do% {
    e=extent(i,i+1000,j,j+1000)

    n=length(r[e])
    r[e]=rep((j-i),n)#zonal stat instead of j-i?
  }
plot(r);plot(raster(filename), add=T, alpha=0.5, legend=F)