GCM data from the IPCC data repository

From SourceWiki
Revision as of 19:31, 6 March 2008 by Wbuytaert (talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

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. Here is an example of how to read the data, extract a region of interest and convert the result to a SpatialGridDataFrame:

region we are interested in (in degree_north and degree_east):

   region <- data.frame(N = 12.45, E = 290, S = -5.23, W = 278)
   library(RNetCDF)
   library(sp)
   nc <- open.nc("HADCM3_SRA1B_1_pr-change_2080-2099.nc")
   subset <- find.subregion(nc,region)
   lat  <- var.get.nc(nc,"latitude")
   long <- var.get.nc(nc,"longitude")

make a grid for the region of interest (the resolution is given manually as it is not constant in the NetCDF files)

   grid <- GridTopology(c(long[subset$W],lat[subset$S]), c(3.75,2.5), c((subset$E-subset$W+1),(subset$N-subset$S+1)))

extract the relevant data from the nc connection and take subset

   prec <- var.get.nc(nc, "precipitation_flux")
   HADCM3_prec_rel <- prec[subset$E:subset$W,subset$S:subset$N,]

averaging over all months:

   HADCM3_prec_rel_y <- as.data.frame(as.vector(apply(HADCM3_prec_rel,c(1,2),"mean")))

converting to SpatialGridDataFrame

   HADCM3_prec_rel_y <- SpatialGridDataFrame(grid, HADCM3_prec_rel_y, proj4string = CRS("+proj=longlat +datum=WGS84"))

Here is the function that extracts the index limits of the smallest subgrid that covers the region of interest:

   find.subregion <- function(nc,region) {
       subset <- data.frame(N = NA, E = NA, S = NA, W = NA)
       bounds <- var.get.nc(nc, "latitude_bounds")
       # search for N limit:
           dummy <- bounds[2,] - region$N
           dummy[dummy < 0] <- 1000
           subset$N <- which.min(dummy)
       # search for S limit:
           dummy <- bounds[2,] - region$S
           dummy[dummy < 0] <- 1000
           subset$S <- which.min(dummy)
       bounds <- var.get.nc(nc, "longitude_bounds")
       # search for E limit:
           dummy <- bounds[2,] - region$E
           dummy[dummy < 0] <- 1000
           subset$E <- which.min(dummy)
       # search for W limit:
           dummy <- bounds[2,] - region$W
           dummy[dummy < 0] <- 1000
           subset$W <- which.min(dummy)
       return(subset)
   }