迭代函數許多光柵堆棧組合成一個

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