问题描述
我有 9000 多个 netCDF 文件，我正在从中提取数据。瓶颈是将数据读入 raster::bricks（或堆栈）。
netCDF 只有 1 个时间步长（一天），但有 72 个高程层和 84 个变量。我正在为每一天、每个变量和第 72 层提取数据。
我的方法是使用 raster::brick，然后使用 extract()。我使用 lapply（或 purrr::map，这里没有太大区别）。我的问题是，经过几次迭代后，数据的读取速度要慢得多。最初，读入一个变量需要约 0.25 秒，但随后会减慢到约 6 秒。光栅砖本身并没有在每次迭代中被保存 - 没有其他应该增长的东西，所以我不知所措。我在 Mac 上，R402。 TIA！
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()
解决方法
您提供的代码并不完全清楚。大概 i 是一天（9000 天之一）？然后 lapply 将 US.extent 的 RasterBrick 放入列表 z 中，用于所有 varname。然后将这些层通过管道传输到 raster::brick 中。我假设在那之后您提取兴趣点的值，然后循环到第二天？
在这种情况下，您最多将保存 RasterBricks 两天（某些层为 168 层）的数据（直到一天的数据覆盖前一天之前；假设未使用的内存立即被释放）。
以下可能表现更好。在 brick 之后调用 lapply 比调用 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,...
#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,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,...
#time : 0
i <- grep("lev=72",names(r))
head(i)
#[1] 72 144 216 288 360 432
现在在文件名的循环中（如果您愿意，可以使用 lapply），创建一个 SpatRaster，带有 i 的子集，并使用 extract
r <- rast(f)[[i]]
extract(r,v)
其中 v 是 SpatVector。