Difference between revisions of "GCM data from the IPCC data repository"

From SourceWiki
Jump to navigation Jump to search
(No difference)

Revision as of 19:31, 6 March 2008

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)
   }