如何有效地从（许多）大型栅格中提取多个多边形的值

Question

我有一组点，我想从几个大栅格中提取值，以获得围绕这些点的缓冲区。栅格太大而无法存放在内存中（> 1e10个单元格）。我将在下面说明我目前的方法，但如果有更快的方法，我会感兴趣。

library(maps)
library(sf)
library(raster)
library(dplyr)
library(parallel)

# sf object with polygones for which we want values
crs <- "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
map <- sf::st_as_sf( maps::map(regions = "Sweden", plot = FALSE, fill = TRUE))
map <- st_transform(map, crs = crs)
sf_points <- st_sfc(st_sample(map, 100))

sf_points <- 
  data.frame(A = 1:length(sf_points)) %>% 
  st_set_geometry(sf_points) 

# raster too large to fit in memory
# the raster(s) I am working on has 10m resolution

r <- raster(extent(map), nrow = 15000, ncol = 7000,
       crs = crs)

values(r) <- rep(sample(1:10, 77, replace = T),  length.out = ncell(r))

#use the parallel package for parallel processing
cluster <- makeCluster(4)
clusterExport(cluster, c("r","sf_points", "as_Spatial"))

List_points <- 
  sf_points %>% 
  mutate(split = rep(1:ceiling(n()/4), each=4, length.out=n())) %>% # 4 cores
  split(f = .$split) %>% 
  parLapply(cl = cluster, X = ., function(x) raster::extract(r, y = as_Spatial(x), buffer = 5000)) %>% 
  unlist(recursive = F)

我重复每个栅格的提取。随着值的排序，我可以汇总列表中的像素值。我不能（轻松地）创建栅格堆栈，因为栅格具有不同的范围。

使用velox包有sugested here似乎没有工作，因为它试图将光栅加载到失败的内存中。我可以尝试加载它，然后我需要找出哪个点在哪个块...

Answer 1

也许你可以使用多边形来加快速度，而不需要聚合（溶解）缓冲区

library(raster)
swe <- getData("GADM", country="SWE", level=0)
set.seed(0)
pts <- spsample(swe, 100, "regular")
r <- raster(swe, nrow = 15000, ncol=7000)
values(r) <- rep(sample(1:10, 77, replace = T),  length.out = ncell(r))

b1 <- buffer(pts, 5000, dissolve=FALSE)
b2 <- buffer(pts, 5000, dissolve=TRUE)

system.time(e1 <- extract(r, pts, buffer=5000))
#   user  system elapsed 
#   1.39    0.02    1.40 
system.time(e2 <- extract(r, b1))
#   user  system elapsed 
#   0.88    0.00    0.88 
system.time(e3 <- extract(r, b2))
#   user  system elapsed 
#  26.34   25.02   51.52 

显然b1比b2表现更好;但并不比第一种方法快得多。

你说你不能制作RasterStack，因为栅格具有不同的范围。但是，如果（并且只有！）它们具有相同的原点和分辨率，您可以先将所有区域转换为xy坐标，然后使用它们。

像这样的东西：

z <- rasterize(b, r)
pts <- rasterToPoints(z, xy=TRUE)

以上需要时间，但之后

system.time(a <- extract(r, zz[,1:2]))
   user  system elapsed 
   0.04    0.00    0.04 

对每个点执行此操作可能更快，并且在栅格化之前使用crop(raster(r), polygon)。