使用 terra 包 R 并行处理提取多个多边形的大时间序列

Question

我正在研究北美大片区域的 NETCDF 数据。时间序列文件超过1000个，空间分辨率相同，时间分辨率不同（年、月、日、次日）。

在我的实际代码中，我实现了一个例程，从每个 NETCDF 文件中读取时间分辨率，并使用它来确定哪些每日文件具有 360、365 和 366 天的日历。

接下来，我使用

terra 

包从所有 NETCDF 文件中提取数据，一次一个文件，一次一个时间步，针对 30 个多边形。

这是一个简单的可重现示例（我无法重现数据文件的确切结构）：

`library(sf)
library(terra)
library(dplyr)
library(raster)
library(sp)


r <- raster(ncol = 100, nrow = 364)
r[] <- runif(ncell(r), -1, 1)
r <- stack(r)
for(i in 1:365){
  x=r[[1]]; x[] <- runif(ncell(x), -1, 1)
  r <- addLayer(r, x)  
}
names(r) <- paste0("evi_", 1:466)     
polys <- spPolygons(rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)), 
                    rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)))
plot(r[[1]])
plot(polys, add=TRUE)

r.med <- terra::extract(r, polys, mean)     ` 

任何想法都值得赞赏。

问题：提取单个NETCDF文件的数据需要6个多小时。如何实现并行处理以加快数据提取速度？

我使用的是超级计算服务器，内存为 140GB，每个节点有 20 个核心。

Answer 1

您没有使用“terra”中的

extract

方法；您正在使用“栅格”中的方法。

extract

是一个通用函数。使用的实现取决于第一个参数的类。在您的情况下，这是一个 RasterStack，因此您正在使用

extract

的“光栅”实现。指定

terra::

不会改变这一点。

这里是使用 terra 方法的示例（通过创建 SpatRaster）

library(terra)
r <- rast(ncol = 100, nrow = 364, nlyr=365)
values(r) <- runif(size(r), -1, 1)
names(r) <- paste0("evi_", 1:nlyr(r))     

polys <- vect(rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)), 
             rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)), type="polygons")
             
plot(r, 1); lines(polys)
r.med <- terra::extract(r, polys, mean)         
r.med[, 1:5]
#  ID         evi_1       evi_2       evi_3      evi_4
#1  1 -0.0008853134 -0.01698036 0.006471589 -0.0138343

你说你

从所有 NETCDF 文件中提取数据，一次一个文件，一次一个时间步，针对 30 个多边形。

这是一个禁忌。一步提取所有时间步的值会更有效。

您当然可以通过循环遍历文件或多边形或两者来并行化它。如何做到这一点取决于您的 HPC 的设置