Iterating over a function to combine many raster stacks into one

Question

Been stuck on this for a while now. Looked everywhere for an answer, but I can't seem to find anything on Stack. Any help you all can give that would be very appreciated.

My main issue is that I need to import many, many netcdf4 files, create raster bricks of each, then combine many bricks to make a "master brick" per variable. To give you a clearer example, I have 40 years (netcdf = 40) of many climate variables (n = 15) that are at a daily resolution. The goal is to aggregate to monthly, but first I have to get this function that reads all years of netcdf's for one variable in and into one large stack.

What I have now reads as follows:

# 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)

Two years of the data can be found here:

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

Any help would be most welcomed. Thank you all in advance, and if this is a duplicate post I apologize; I looked far and wide to no avail.

Answer 1

Your code is fine, you just need to return the loaded brick rbrck from your function, otherwise you'll get the extent.

As for loading and stacking, I'd suggest using lapply to apply the function to each datafile. This will give you a neat list with a year per item. There you could do some more processing and finally just call stack on the list to produce your "master brick".

Mind that I only did this with two files, so I'm not sure about the size of the whole thing when you do it with 40.

Here's your modified code:

# 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