迭代函数以将许多栅格堆栈合并为一个

Question

在这个问题上卡住了一段时间。到处寻找答案，但我似乎在 Stack 上找不到任何东西。大家可以提供的任何帮助将不胜感激。

我的主要问题是我需要导入很多很多 netcdf4 文件，为每个文件创建光栅积木，然后将许多积木组合成每个变量的“主积木”。为了给你一个更清楚的例子，我有 40 年 (netcdf = 40) 的许多气候变量 (n = 15) 每天都有分辨率。目标是每月汇总一次，但首先我必须获得这个函数，它可以读取所有年份的 netcdf 的一个变量，并将其放入一个大堆栈中。

我现在的内容如下：

# Libraries --------------------------------------------------------------
library(raster)
library(ncdf4)

# Directories -------------------------------------------------------------

tmp_dl <- list.files("/Users/NateM", pattern = "nc",
                 full.names = TRUE)
# Functions ---------------------------------------------------------------
rstlist = stack()

netcdf_import <- function(file) {
   nc <- nc_open(file)
   nc_att <- attributes(nc$var)$names
   ncvar <- ncvar_get(nc, nc_att)
   rm(nc)
   proj <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
   rbrck <- brick(ncvar, crs= proj)
   rm(ncvar)
   extent(rbrck) <- c(-124.772163391113, -67.06383005778, 25.0626894632975, 
                       49.3960227966309) 
   }      

t <- for(i in 1:length(tmp_dl)) {
         x <- netcdf_import(tmp_dl[i])
         rstlist <- stack(rstlist, x)
      }

allyears <- stack(t)

两年的数据可以在这里找到：

https://www.northwestknowledge.net/metdata/data/pdsi_2016.nc https://www.northwestknowledge.net/metdata/data/pdsi_2015.nc

非常欢迎任何帮助。提前谢谢大家，如果这是重复的帖子，我深表歉意；我四处张望无济于事。

Answer 1

你的代码没问题，你只需要从你的函数中return加载的砖块rbrck，否则你会得到范围。

至于加载和堆叠，我建议使用lapply 将函数应用于每个数据文件。这将为您提供一份整洁的清单，其中每个项目一年。在那里您可以进行更多处理，最后只需调用列表中的stack 即可生成您的“主积木”。

请注意，我只对两个文件执行此操作，所以当您使用 40 个文件时，我不确定整个文件的大小。

这是您修改后的代码：

# Libraries --------------------------------------------------------------
library(raster)
library(ncdf4)

# Directories -------------------------------------------------------------

tmp_dl <- list.files("/Users/NateM", pattern = "nc",
                     full.names = TRUE)
# Functions ---------------------------------------------------------------

netcdf_import <- function(file) {
  nc <- nc_open(file)
  nc_att <- attributes(nc$var)$names
  ncvar <- ncvar_get(nc, nc_att)
  rm(nc)
  proj <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
  rbrck <- brick(ncvar, crs= proj)
  rm(ncvar)
  extent(rbrck) <- c(-124.772163391113, -67.06383005778, 25.0626894632975, 
                     49.3960227966309) 
  return(rbrck)
}      

# apply function
allyrs <- lapply(tmp_dl,netcdf_import)

# stack to master brick
allyrs <- do.call(stack,allyrs)

HTH