parallel_rast_extract.md

service <- file.path(
  "https://s3-us-west-2.amazonaws.com", 
  "mrlc"
)
layers <- paste0(
  service, 
  "/Annual_NLCD_LndCov_", 
  years, 
  "_CU_C1V0.tif"
)

layers <- sprintf("/vsicurl/%s", layers)
layers

## write to an actual VRT file
vrt(layers, myvrt <- "layers.vrt", options = "-separate", overwrite = F)

r <- rast(myvrt)

## make polygons from points in the crs of the raster
pt <- structure(c(-1955363.12995246, -1727153.93819334, -1942684.84152139, 
-1714475.64976228, -1955363.12995246, -1562336.18858954, -1955363.12995246, 
-1473588.16957211, -1777867.09191759, -1403857.58320127, -1372161.86212361, 
-928421.76703645, -383255.364500792, 123876.17274168, 688060.00792393, 
973321.497622821, 459850.816164818, -85315.5863708397, -459325.095087163, 
-858691.18066561, 3026651.27575277, 2969598.97781299, 2722372.35340729, 
2620946.0459588, 2329345.41204437, 2253275.681458, 1942657.61489699, 
1853909.59587956, 1638378.6925515, 1613022.11568938, 1302404.04912837, 
1346778.05863708, 1346778.05863708, 966429.40570523, 1048838.28050713, 
1479900.08716323, 1898283.60538827, 2196223.38351822, 2449789.15213946, 
2614606.90174326), dim = c(20L, 2L), dimnames = list(NULL, c("x", 
"y")))

p <- buffer(vect(pt, crs = crs(r)), 4000)

library(terra)
## we get a matrix of values for each polygon
#fun <- function(vrt, p) {
#  extract(crop(rast(vrt), p), p)
#}
## parallelize this loop
#l <- vector("list", nrow(p))
#for (i in seq_along(l)) l[[i]] <- fun(myvrt, p[i, ])

#options(parallelly.fork.enable = TRUE, future.rng.onMisuse = "ignore")
library(furrr); plan(multicore)

## put the VRT on each polygon, so each row is self-contained 
p$vrt <- myvrt
listof <- split(p, 1:nrow(p))
fun1 <- function(p) {
  extract(crop(rast(p$vrt), p), p, raw = TRUE)
}
a <- future_map(listof, fun1)
plan(sequential)