Zonal Statistics R (raster/polygon)

spatialEco::zonal.stats uses exactextractr (I have not checked the code, but it told me to install it to be able to use zonal.stats) which should be more exact if you are considering polygons (the raster package turns them into rasters first, see zonal below). However, the below example (it is only one case) suggest that spatialEco is less precise.

Example (avoid random numbers, but if you do use them, use set.seed). I start with very large grid cells.

library(raster)
library(spatialEco)

ras <- raster(nrows=4, ncols=4, xmn=0, xmx=1000, ymn=0, ymx=800)
values(ras) <- 1:ncell(ras)
set.seed(1)
xy <- cbind(runif(3,0,1000), runif(3,0,800))
xy <- rbind(xy, xy[1,])
sp <- spPolygons(xy, attr=data.frame(x=1))


### zonal statistics using "spatialECO"
zonal.stats(sp, ras, stats="mean")
#  mean.layer
#1          7
### zonal statistics using "raster"
extract(ras, sp, fun=mean)
#     [,1]
#[1,]    6
### same as 
# x <- rasterize(sp, ras)
# zonal(ras, x, "mean")

With raster you can also get a more precise estimate like this

e <- extract(ras, sp, weights=T)[[1]]
weighted.mean(e[,1], e[,2])
#[1] 5.269565

To see how many cells are used

zonal.stats(sp, ras, stats="counter")
#  counter.layer
#1             6
extract(ras, sp, fun=function(x,...)length(x))
#     [,1]
#[1,]    3

One way to look at this is to create higher resolution raster data.

10x higher resolution, same values

ras <- disaggregate(ras, 10)
zonal.stats(sp, ras, stats="mean")
#  mean.layer
#1        5.5
extract(ras, sp, fun=mean)
#         [,1]
#[1,] 5.245614

zonal.stats(sp, ras, stats="counter")
#  counter.layer
#1           218
extract(ras, sp, fun=function(x,...)length(x))
#     [,1]
#[1,]  171

100x higher resolution, same values

ras <- disaggregate(ras, 10)
zonal.stats(sp, ras, stats="mean")
#mean.layer
#1   5.299915
extract(ras, sp, small=TRUE, fun=mean)
#         [,1]
#[1,] 5.271039


zonal.stats(sp, ras, stats="counter")
# counter.layer
#1         17695

extract(ras, sp, fun=function(x,...)length(x))
#      [,1]
#[1,] 17289

At the highest resolution the mean values are similar (and the relative difference in the number of cells is small); but raster was closer to the correct value (whatever that is, exactly) at a lower resolution (and also with the weighted mean). That is unexpected.

For better speed, there is now also the terra package

library(terra)    
r <- rast(ras)
v <- vect(sp)
extract(r, v, "mean")    
#     ID    layer
#[1,]  1 5.271039

Each algorithm is using a different number of pixels to calculate the zonal stats; thus, the difference is probably caused by that (z1 is higher than z2). According to this tiny example I would infer that zonal.stats is less restrictive than extract. Thus, probably zonal.stats takes into account every value of the raster that falls inside the polygon; however, extract only takes into account pixels whose center is found inside the polygon (check the function's documentation).

# Create a function to count the number of pixels used to calculate the zonal stats
counter <- function(x, na.rm = T) { 
  length(x)
} 

#z1 zonal statistics using "spatialECO"
z1 <- zonal.stats(spdf, ras, stats="counter")

#z2 zonal statistics using "raster"
z2 <- extract(ras, spdf, fun=counter)

Thanks @Robert Hijmans for your insights. Based on your example I've done some further comparisions of 1) precision and 2) computing time by calculating zonal means for different functions and resolutions:

library(raster)
library(spatialEco)
library(terra)
library(exactextractr)
library(sf)

ras <- raster(nrows=4, ncols=4, xmn=0, xmx=1000, ymn=0, ymx=800)
values(ras) <- 1:ncell(ras)
set.seed(1)
xy <- cbind(runif(3,0,1000), runif(3,0,800))
xy <- rbind(xy, xy[1,])
sp <- spPolygons(xy, attr=data.frame(x=1))

mn <- data.frame(matrix(ncol=6, nrow=4))
colnames(mn) <- c("disagr", "raster", "raster_weight", "spatialEco", "exactextractr", "terra")
mn[,1] <- c(2,10,50,250)
tm <- mn

for (i in 1:5){
  d <- mn[i,1]
  rasd <- disaggregate(ras, d)
  
  on <- Sys.time()
  mn[i,2] <- raster::extract(rasd, sp, fun=mean)
  off <- Sys.time()
  tm[i,2] <- off - on
  
  on <- Sys.time()
  mn[i,3] <- raster::extract(rasd, sp, fun=mean, weights=T)[[1]]
  off <- Sys.time()
  tm[i,3] <- off - on
  
  on <- Sys.time()
  mn[i,4] <- spatialEco::zonal.stats(sp, rasd, stats="mean")
  off <- Sys.time()
  tm[i,4] <- off - on
  
  on <- Sys.time()
  mn[i,5] <- exactextractr::exact_extract(rasd, st_as_sf(sp), fun="mean")
  off <- Sys.time()
  tm[i,5] <- off - on

  on <- Sys.time()
  val <- terra::extract(rast(rasd), vect(sp))
  mn[i,6] <- mean(val[,2])
  off <- Sys.time()
  tm[i,6] <- off - on
  print(i)
}
mn  # arithmetic mean
  disagr   raster raster_weight spatialEco exactextractr    terra
1      2 5.333333      5.269565   6.647059      5.271303 5.333333
2     10 5.245614      5.271039   5.500000      5.271303 5.245614
3     50 5.272370      5.271328   5.325525      5.271303 5.272370
4    250 5.271314      5.271303   5.282827      5.271303 5.271314
tm   # computing time in seconds
  disagr     raster raster_weight  spatialEco exactextractr       terra
1      2 0.03998685    0.03598809 0.003000021   0.002999067 0.008996964
2     10 0.09783196    0.04598618 0.003997803   0.003000021 0.008984089
3     50 0.37189507    0.40886688 0.004998922   0.003998041 0.021993160
4    250 4.10671687    8.29134583 0.035988092   0.019991875 0.336881876

Based on this example and when using polygons as zones exact_extract is the preferable choice.