как эффективно извлечь значения из (многих) больших растров для многих полигонов

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 здесь , похоже, не работает, так как пытается загрузить растр в память, которая выходит из строя. Я мог бы попытаться загрузить его по кусочкам, но тогда мне нужно было бы выяснить, какие точки на каком кусочке есть ...

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) перед растеризацией.