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

First, 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, we interpolate on a user defined grid using the nearest neighbour approach. Define a grid, convert it to a SpatialGridDataFrame, which we initialise with NA, and then we look for the closest gridcell in map using spDistsN1().

mymap <- GridTopology(c(-81.75,-5.25), c(0.5,0.5), c(24,36))
    
mymap <- cbind(coordinates(mymap), temperature = NA)
coordinates(mymap) <- ~s1 + s2
gridded(mymap) <- T
proj4string(mymap) <- CRS("+proj=latlong +datum=WGS84")
    
GridPoints <- coordinates(mymap)
MapPoints  <- coordinates(map)
    
for(i in 1:dim(GridPoints)[1]) {
  distance <- spDistsN1(MapPoints, GridPoints[i,], longlat = TRUE)
  mymap$temperature[i] <- map$var[which.min(distance)]
}

Finally we can plot the map:

spplot(grid, scales=list(draw=T))
 
guides/tutorials/hydrological_data_analysis/gcm_data.txt · Last modified: 2009/12/08
 
Recent changes RSS feed R Wiki powered by Driven by DokuWiki and optimized for Firefox Creative Commons License