几次迭代后，在 R 中读取多个栅格/netCDF 数据会变慢

Question

我有 9000 多个 netCDF 文件，我正在从中提取数据。瓶颈是将数据读入 raster::bricks（或堆栈）。

netCDF 只有 1 个时间步（单天），但有 72 个高程层和 84 个变量。我每天都在为每个变量和第 72 层提取数据。

我的方法是使用 raster::brick 后跟 extract()。我使用 lapply （或 purrr::map，这里没有太大区别）。我的问题是，经过几次迭代后，数据的读取速度要慢得多。最初，读取一个变量需要约 0.25 秒，但随后会减慢至约 6 秒。光栅砖本身并没有在每次迭代中保存 - 没有其他东西应该增长，所以我很茫然。我在 Mac 上，R402。蒂亚！

https://portal.nccs.nasa.gov/datashare/merra2_gmi/

示例代码：我只是将 print(x) 用于跟踪，但在问题解决后最终会删除。

z <- lapply(1:length(varnames), 
    function(x) { 
         print(x)
         raster::brick(MERRA2.GMI.filenames.read[i],varname = varnames[x])[[72]] 
         %>% crop(US.extent)
     }
   ) 
%>% raster::brick()

Answer 1

您提供的代码并不完全清楚。大概i 是那天（9000 之一）？然后lapply 将US.extent 的RasterBrick 放入列表z 中，用于所有变量名。然后将这些层通过管道传输到 raster::brick。我假设之后您提取兴趣点的值，然后循环到第二天？

在这种情况下，您最多将 RasterBricks 保存两天（168 层）数据（直到一天的数据覆盖前一天；假设未使用的内存被立即释放）。

以下可能表现更好。在lapply 之后调用brick 比调用stack 更昂贵，因为后者实际上还没有读取任何值，它只指向文件。它应该（但你永远不知道）只使用一次crop 也会更快；并且这样做不会在内存中创建包含栅格数据的列表（仅引用文件，直到调用 crop）。

z <- lapply(1:length(varnames), 
    function(x) { 
         raster::brick(filenames[i], varname = varnames[x])[[72]] 
     }
   ) 
 s <- stack(z) %>% crop(US.extent)

我下载了一个文件，我建议使用以下工作流程，使用 terra 通常更快，特别是如果您使用多边形进行提取。但在这种情况下，terra 还可以更轻松地选择感兴趣的变量，完全避免使用 lapply。

对于一个文件，我创建了一个 SpatRasterDataSet 来了解子数据集结构

library(terra)
f <- "MERRA2_GMI.inst0_3d_ovp_Nv.19900202_1200z.nc4"
s <- sds(f)
s
#class       : SpatRasterDataset 
#subdatasets : 74 
#dimensions  : 361, 576 (nrow, ncol)
#nlyr        : 72, 72, 72, 72, 72, 72, 72, 72, 72, ... 
#resolution  : 0.625, 0.5  (x, y)
#extent      : -180.3125, 179.6875, -90.25, 90.25  (xmin, xmax, ymin, ymax)
#coord. ref. :  
#source(s)   : MERRA2_GMI.inst0_3d_ovp_Nv.19900202_1200z.nc4 
#names       : OVP10_AGE, OVP10_CH2O, OVP10_CH4, OVP10_CO, OVP10_EM_LGTNO, OVP10_EPV, OVP10_FCLD, OVP10_GOCART_SO2_VMR, OVP10_GOCART_SO2v_VMR, OVP10_NO, OVP10_NO2, OVP10_O3, OVP10_PL, OVP10_PPBL, OVP10_PS, OVP10_QLTOT, OVP10_QV_VMR, OVP10_T, OVP10_TAUTT, OVP10_TOTEXTTAU, OVP10_TROPP, OVP10_U10, OVP10_V10, OVP10_stO3, OVP14_AGE, OVP14_Br, OVP14_BrCl, OVP14_BrO, OVP14_BrONO2, OVP14_CH2O, OVP14_CH3Cl, OVP14_CH4, OVP14_CO, OVP14_Cl, OVP14_Cl2, OVP14_Cl2O2, OVP14_ClO, OVP14_ClONO2, OVP14_EM_LGTNO, OVP14_EPV, OVP14_FCLD, OVP14_GOCART_NH3_VMR, OVP14_GOCART_SO2_VMR, OVP14_GOCART_SO2v_VMR, OVP14_HBr, OVP14_HCOOH, OVP14_HCl, OVP14_HNO3, OVP14_HNO3COND, OVP14_HO2, OVP14_HOBr, OVP14_HOCl, OVP14_ISOP, OVP14_MOH, OVP14_N2O, OVP14_NO, OVP14_NO2, OVP14_O3, OVP14_OClO, OVP14_OH, OVP14_PAN, OVP14_PL, OVP14_PPBL, OVP14_PS, OVP14_QLTOT, OVP14_QV_VMR, OVP14_R4N2, OVP14_T, OVP14_TAUTT, OVP14_TOTEXTTAU, OVP14_TROPP, OVP14_U10, OVP14_V10, OVP14_stO3 

选择感兴趣的图层。有些变量只有一层，所以我将第 72 层用于那些有这么多的变量。

n <- nlyr(s)
i <- cumsum(n)
i <- i[n==72]

或者您可以执行以下操作，注意rast(f) 将所有子数据集组合成一个多层数据集

r <- rast(f)
#class       : SpatRaster 
#dimensions  : 361, 576, 4334  (nrow, ncol, nlyr)
#resolution  : 0.625, 0.5  (x, y)
#extent      : -180.3125, 179.6875, -90.25, 90.25  (xmin, xmax, ymin, ymax)
#coord. ref. :  
#sources     : MERRA2_GMI.inst0_3d_ovp_Nv.19900202_1200z.nc4:OVP10_AGE  (72 layers) 
#              MERRA2_GMI.inst0_3d_ovp_Nv.19900202_1200z.nc4:OVP10_CH2O  (72 layers) 
#              MERRA2_GMI.inst0_3d_ovp_Nv.19900202_1200z.nc4:OVP10_CH4  (72 layers) 
#              ... and 71 more source(s)
#varnames    : OVP10_AGE (Age of air (uniform source) tracer_10am_local) 
#              OVP10_CH2O (Formaldehyde_10am_local) 
#              OVP10_CH4 (Methane_10am_local) 
#              ...
#names       : OVP10~lev=1, OVP10~lev=2, OVP10~lev=3, OVP10~lev=4, OVP10~lev=5, OVP10~lev=6, ... 
#unit        :        days,        days,        days,        days,        days,        days, ... 
#time        : 0 
 
i <- grep("lev=72", names(r))
head(i)
#[1]  72 144 216 288 360 432

现在循环（如果您愿意，可以使用）而不是文件名，创建一个 SpatRaster，使用 i 的子集，并使用 extract

r <- rast(f)[[i]]
extract(r, v)

v 是一个 SpatVector。