Difference between revisions of "GCM data from the IPCC data repository"
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) }