GCM data from the IPCC data repository

Introduction

You can download summary GCM data from the IPCC data centre for free for research purposes. They come in NetCDF format, which you can read into R with the RNetCDF package. However, the GCMs run at different resolutions, which complicates the combination of them to calculate characteristics such as the predicted mean and range. The best method is to interpolate the maps to a common grid for the region of interest, and then calculate the statistics for each cell. The interpolation methods used here is the nearest neighbour.

Note: some IPCC maps do not come in a perfectly regular grid. Gridcells are often smaller at the boundaries, and for some reason the grid centers are not completely aligned. This creates problems for making a SpatialGridDataFrame. Interpolating to a self-defined regular grid avoids this issue.

Required packages

RNetCDF (CRAN) is used to import the NetCDF data. sp (CRAN) defines the spatial classes.

   library(RNetCDF)
   library(sp)

Procedure

Define the grid we are interested in. The example is for the northern Andes region in South America:

   grid <- GridTopology(c(-81.75,-5.25), c(0.5,0.5), c(24,36))

Open the NetCDF file and get the required variables. Here we are interested in the air temperature.

   nc <- open.nc("C:/GIAOM_SRA1B_1_tas-change_2080-2099.nc")
   lat <- var.get.nc(nc,"latitude")
   long <- var.get.nc(nc,"longitude")
   data <- var.get.nc(nc, "air_temperature_anomaly")

To make a SpatialPointsDataFrame we need a dataframe where each row consists of an entry of the form (x,y,a) where x and y are the coordinates and a is the value. So we need to reformat the lat, long and data objects to be in the right format, We also average the temperature to yearly values by taking the mean over the 12 months.

   latx  <- rep(lat,length(long))
   longx <- rep(long, length(lat))
   longx <- as.vector(t(matrix(longx,nrow=length(long))))
   data <- data.frame(var  = as.vector(t(apply(data,c(1,2),"mean"))))
   map <- SpatialPointsDataFrame(cbind(longx,latx), data, proj4string = CRS("+proj=latlong +datum=WGS84"))

Now interpolate using the closest grid cell of the original predictions. First we need to convert the grid to a SpatialGridDataFrame, which we initialise with NA, and then we look for the closest gridcell in map using spDistsN1().

   data <- data.frame(temperature = rep(NA,length(coordinates(grid)[,1])))
   grid <- SpatialGridDataFrame(grid, data, proj4string = CRS("+proj=longlat +datum=WGS84"))
   GridPoints<- coordinates(grid)
   MapPoints <- coordinates(map)
   for(i in 1:dim(GridPoints)[1]) {
       distance <- spDistsN1(MapPoints, GridPoints[i,], longlat = TRUE)
       grid$temperature[i] <- map$var[which.min(distance)]
   }

Finally we can plot the map:

   spplot(grid, scales=list(draw=T))