我一直在想，在一个提取调用中运行多个函数是否更快，而且确实如此。
下面是Elia的示例和qdread使用terra包和microbenchmark评估的答案。第二个选项几乎快三倍：

library(terra)
library(microbenchmark)
# library(sp)

#from the help page of extract
r <- rast(ncol=36, nrow=18, vals=1:(18*36))
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- terra::vect(rbind(cds1, cds2), type="polygons")

#metrics extraction
microbenchmark(A={
mean <- extract(r, polys,mean,df=T)
min<-extract(r, polys,min,df=T)
max<-extract(r, polys,max,df=T)
#and so on for other summary functions (like sd, mode, median, sum etc...)
},
B={
my_summary <- function(x, na.rm = T) c(mean = mean(x, na.rm=na.rm), min = min(x, na.rm=na.rm), max = max(x, na.rm=na.rm))
r_summary <- extract(r, polys, fun = my_summary)
names(r_summary) <- c('ID', 'mean', 'min', 'max')
})

评价：

Unit: milliseconds
 expr    min      lq     mean  median      uq     max neval cld
    A 5.8892 5.98180 6.429374 6.12040 6.25440 19.2841   100   b
    B 2.1281 2.16045 2.357640 2.18765 2.30115 14.4645   100  a

这个问题已经过去了很多年，我认为对“Gerald”的答案进行一些扩展会有所帮助。在写这篇文章的时候，到目前为止，最好的解决方案是exactextractr包中的exact_extract函数。该函数可以接受函数的字符向量作为fun参数。exact_extract考虑了每个单元的面积被多边形覆盖的部分，并且其被优化以与sf对象和大SpatRaster对象以及Raster*对象一起工作。
下面是一个示例，显示了栅格比原始示例更大时的可伸缩性

sf_polys <- sf::st_as_sf(polys)
r_summary <- exact_extract(r,sf_polys,fun=c("mean","min","max"))

具有更大光栅的基准：

r <- rast(ncol=3600, nrow=1800, vals=1:(1800*3600))
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))

microbenchmark(A = {
  mean <- extract(r, polys, mean, df = T)
  min <- extract(r, polys, min, df = T)
  max <- extract(r, polys, max, df = T)
  #and so on for other summary functions (like sd, mode, median, sum etc...)
},
B = {
  my_summary <-
    function(x, na.rm = T)
      c(
        mean = mean(x, na.rm = na.rm),
        min = min(x, na.rm = na.rm),
        max = max(x, na.rm = na.rm)
      )
  r_summary <- extract(r, polys, fun = my_summary)
  names(r_summary) <- c('ID', 'mean', 'min', 'max')
},
C = {
  r_summary <- exact_extract(r, sf_polys, fun = 'mean')
})

结果如下：

Unit: milliseconds
 expr       min        lq      mean    median       uq       max neval cld
    A 6309.1327 6362.1017 6489.0986 6459.2223 6632.241 6710.0625    100   c
    B 2100.8945 2116.3198 2156.2030 2137.7992 2162.850 2283.2969    100  b 
    C  130.3108  146.3663  184.8943  151.0666  157.701  505.5154    100 a

可以看出，C方法比其他两种方法快得多

正如@dww在上面的评论中所建议的，下面是一个函数，它计算一些汇总统计数据，并将它们作为向量返回。它被传递给raster::extract的fun参数。请注意，raster::extract的文档说明该函数必须接受na.rm参数。我无法更改extract命名数据框输出列的默认行为，因此我在以后手动设置名称。

代码

my_summary <- function(x, na.rm) c(mean = mean(x, na.rm=na.rm), min = min(x, na.rm=na.rm), max = max(x, na.rm=na.rm))
r_summary <- extract(r, polys, fun = my_summary, df = TRUE)
names(r_summary) <- c('ID', 'mean', 'min', 'max')

输出

ID     mean min max
1  1 387.8158 326 507
2  2 321.0800 172 498