如何在raster::extract()中传递多个函数

tvokkenx  于 2023-03-10  发布在  其他
关注(0)|答案(3)|浏览(150)

我正在使用一个栅格CHM,我必须从一个多边形shapefile中提取几个指标,现在我正在做这样的事情:

library(raster)
library(sp)

#from the help page of extract
r <- raster(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 <- spPolygons(cds1, cds2)

#metrics extraction
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...)

我想知道是否有一种方法可以将所有的汇总函数传递给extract()函数的fun=参数,以及是否有可能并行执行。
注:这是我在StackOverflow中的第一个问题,我为任何不当之处道歉

vs3odd8k

vs3odd8k1#

我一直在想,在一个提取调用中运行多个函数是否更快,而且确实如此。
下面是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
ht4b089n

ht4b089n2#

这个问题已经过去了很多年,我认为对“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方法比其他两种方法快得多

ni65a41a

ni65a41a3#

正如@dww在上面的评论中所建议的,下面是一个函数,它计算一些汇总统计数据,并将它们作为向量返回。它被传递给raster::extractfun参数。请注意,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

相关问题