1

一直停留在这个问题上。到处寻找答案，但我似乎在 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

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

4

1 回答 1

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

于 2017-07-21T09:08:26.827 回答