避免 raster::extract(rst,shp) 中的 for 循环

Question

我正在使用 R 提取某些建筑物 3 米缓冲区内的栅格的平均值和最大值。

为此，我创建了一个 for 循环，它遍历每个建筑物以提取这两个值。我当前的代码如下所示：

for (b in c(1:nrow(buildings_shp))){
    
    building <- buildings_shp[b,]
    
    buffered <- st_buffer(building, 3)
    
    raster_cropped <- crop(raster, extent(buffered))

    mean <- extract(depths_cropped, buffered, fun = mean, na.rm = TRUE)
    max <- extract(depths_cropped, buffered, fun = max, na.rm = TRUE)
    
    buildings_shp[b,"mean"] <- mean
    buildings_shp[b,"max"] <- max
    
  }

然而，这个循环需要相当长的时间（1500 座建筑物大约需要 17 分钟），而且似乎花费最多时间的步骤是两条提取线。我想知道是否有方法可以通过以下方式加快此过程：

a) 避免使用循环——这个循环的原因是我担心如果我在整个数据集上使用 st_buffer，那么当建筑物距离小于 3 米时我会生成重叠的几何图形，这可能会导致错误。

b) 并行化 for 循环（我尝试了栅格聚类功能，但它并没有加快进程，可能是因为它没有并行化循环本身，而是并行化提取函数）

c) 使用除 raster::extract 之外的其他函数。我看过一些帖子 recommending the velox package，但好像这个包已经从 CRAN 中删除了。

一些虚拟数据（从上面引用的问题中复制）

library(raster)
library(sf)

raster <- raster(ncol=1000, nrow=1000, xmn=2001476, xmx=11519096, ymn=9087279, ymx=17080719)
raster []=rtruncnorm(n=ncell(raster ),a=0, b=10, mean=5, sd=2)
crs(raster ) <- "+proj=utm +zone=51 ellps=WGS84"
    
x1 <- runif(100,2001476,11519096)
y1 <- runif(100, 9087279,17080719)

buildings_shp <- st_buffer(st_sfc(st_point(c(x1[1],y1[1]), dim="XY"),crs=32651),200000)

Answer 1

你不需要循环。来自 ?raster::extract

的示例数据

library(raster)
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))
buildings <- spPolygons(cds1, cds2)

获取缓冲区并提取。由于您要计算两个统计量，因此在这种情况下不使用汇总函数会更容易。

b <- buffer(buildings, width=3, dissolve=FALSE)
e <- extract(r, b)

现在计算统计数据

sapply(e, mean, na.rm=TRUE)
#[1] 379.4167 330.0741
sapply(e, max, na.rm=TRUE)
#[1] 507 498

terra

应该会更快

library(terra)
v <- vect(b)
x <- rast(r)
ee <- extract(x, v)

避免 raster::extract(rst,shp) 中的 for 循环

Avoid a for loop in raster::extract(rst,shp)

performance

r

raster

sf