Illustration of how to use the R sp, raster, and ncdf packages to output GCAM data in netCDF format

gcam2nc.R

library(sp)
library(rgdal)
library(raster)
library(ncdf)

##aez map
## extract data from shape file

GIS_DIR <- "./cropdata"
GIS_FILE <- "GCAM_region_AEZ"
setwd(GIS_DIR)

## read the shape files for GCAM AEZs, etc. This is a "SpatialPolygonsDataFrame" (SPDF)
gcam <- readOGR( dsn=".", layer=GIS_FILE )

## subset the geo data to china for this example
gc.chn <- gcam[gcam[["CNTRY_NAME"]]=="China","GRIDCODE"]
## "GRIDCODE" is the GCAM AEZ.  This is *not* the combination of AEZ and region
## that we sometimes call "AEZ", so it runs from 1..18, not 1..~230

## create a grid.  There won't be any data in it; it's just a template
## for later use
grid <- raster(extent(gc.chn), resolution=0.5) 

## read in our example GCAM data:  crop pct. area for china for 3 crops and 
## 5 SSPs
d <- read.csv("./prodbyaezarea.csv", stringsAsFactors=FALSE)

scenarios <- levels(as.factor(d$Scenario))
crops <- levels(as.factor(d$Crop))
years <- seq(2005,2050,5)

for (s in scenarios){

    for (crop in crops) {
        cat('proc crop= ', crop, ' for scenario= ', s, '\n')

        ## create an empty list to collect the raster layers we will
        ## be creating
        levels <- list()

        ## subset the GCAM data to the crop and scenario we're working
        data <- subset(d, Crop==crop & Scenario==s)
        
        for(year in years) {
            ## column names are X2005, X2010, etc.
            token <- paste('X',year, sep="")

            ## rasterize converts to a grid of pixels.  Use it like this:
            ## rasterize(SPDF, template.grid, data.vector)
            ## In this case the data vector will be constructed by finding the row with AEZ
            ## equal to the GRIDCODE and the column with the stringified year.
            levels <- c(levels,         # append next layer to list
                        rasterize(gc.chn, grid,
                                  100*data[match(gc.chn[["GRIDCODE"]], data$AEZ), token]))
        }
        ## Collect the individual layers into a 'brick' and set the z (i.e., time) values
        ## Note that setZ is misleadingly named.  It doesn't actually set the z values on 
        ## the brick it is given; it returns a new brick with those values set.
        outdata <- setZ(brick(levels), years, name='year')
        filename <- paste(s,crop,"2005-2050.nc",sep='-')
        vname <- paste(crop,"crop_pct_area",sep='_')

        ## writeRaster recognizes the '.nc' on the end of the filename
        ## and knows it's supposed to write a netCDF file.  We have to
        ## set the name and units of the data variable, and the name
        ## and units of the z-axis.
        writeRaster(outdata, filename, varname=vname, varunit='pct', zname='year',
                    zunit='', overwrite=T)
        cat('outfile is: ', filename, '\n')
    } 
}