Creating a list of raster bricks from a multidimensional netCDF file

Question

I work with spatial data RCP (representative concentration path). This is a good grid dataset in netCDF format. How can I get a list of bricks where each element represents a single variable from a multidimensional netCDF file (by variable I do not mean lat, lon, time, depth ... etc.). This is what Ive tried to do. I cannot post sample data, but I configured the script below to be reproducible if you want to look at it. Obviously, questions are welcome ... Perhaps I could not clearly express the language associated with the code. Greetings.

A: Package Requirements

library(sp) library(maptools) library(raster) library(ncdf) library(rgdal) library(rasterVis) library(latticeExtra) 

B: collect data and see the structure of the netCDF file

  td <- tempdir() tf <- tempfile(pattern = "fileZ") download.file("http://tntcat.iiasa.ac.at:8787/RcpDb/download/R85_NOX.zip", tf , mode = 'wb' ) nc <- unzip( tf , exdir = td ) list.files(td) ## Take a look at the netCDF file structure, beyond this I don't use the ncdf package directly ncFile <- open.ncdf(nc) print(ncFile) vars <- names(ncFile$var)[1:12] # I'll try to use these variable names later to make a list of bricks 

C: Create a raster brick for one variable. Levels correspond to years

  r85NOXene <- brick(nc, lvar = 3, varname = "emiss_ene") NAvalue(r85NOXene) <- 0 dim(r85NOXene) # [1] 360 720 12 

D: names for faces

  data(wrld_simpl) # in maptools worldPolys <- SpatialPolygons( wrld_simpl@polygons ) cTheme <- rasterTheme(region = rev(heat.colors(20))) levelplot(r85NOXene,layers = 4,zscaleLog = 10,main = "2020 NOx Emissions From Power Plants", margin = FALSE, par.settings = cTheme) + layer(sp.polygons(worldPolys)) 

E: sum all grid cells for each year with one variable "emis_ene", I want to do this for every variable of the netCDF file that I work with.

  gVals <- getValues(r85NOXene) dim(gVals) r85NOXeneA <- sapply(1:12,function(x){ mat <- matrix(gVals[,x],nrow=360) matfun <- sum(mat, na.rm = TRUE) # Other conversions are needed, but not for the question return(matfun) }) 

F: Meet and meet. See what E looks like

  library(ggplot2) # loaded here because of masking issues with latticeExtra years <- c(2000,2005,seq(2010,2100,by=10)) usNOxDat <- data.frame(years=years,NOx=r85NOXeneA) ggplot(data=usNOxDat,aes(x=years,y=(NOx))) + geom_line() # names to faces again detach(package:ggplot2, unload=TRUE) 

G: try creating a list of bricks. List of objects created in Part C

  brickLst <- lapply(1:12,function(x){ tmpBrk <- brick(nc, lvar = 3, varname = vars[x]) NAvalue(tmpBrk) <- 0 return(tmpBrk) # I thought a list of bricks would be a good structure to do (E) for each netCDF variable. # This doesn't break but, returns all variables in each element of the list. # I want one variable in each element of the list. # with brick() you can ask for one variable from a netCDF file as I did in (C) # Why can't I loop through the variable names and return on variable for each list element. }) 

H: Get rid of the trash you might have loaded ... Sorry

  file.remove(dir(td, pattern = "^fileZ",full.names = TRUE)) file.remove(dir(td, pattern = "^R85",full.names = TRUE)) close(ncFile) 

Answer 1

Step (E) can be simplified using cellStats .

 foo <- function(x){ b <- brick(nc, lvar = 3, varname = x) NAvalue(b) <- 0 cellStats(b, 'sum') } sumLayers <- sapply(vars, foo) 

sumLayers is the result you are looking for if I understood your question correctly.

Alternatively, you can use the zoo package because you are dealing with time series.

 library(zoo) tt <- getZ(r85NOXene) z <- zoo(sumLayers, tt) xyplot(z)