## Copyright 2009-2011 Pacific Climate Impacts Consortium 
## This program is distributed under the terms of the GNU General Public License 
## See http://www.gnu.org/licenses/

## update the following source path to where you have saved the R script with the functions
## you must have installed the libraries: rgdal, ncdf, raster
source("/home/haileye/code/darcs/r_code/functions.for.narccap.group.r")

shp <- readOGR('/home/', 'columbia_boundary', stringsAsFactors=F)
nc.fname <- ('/home/narccap_hrm3-sresa2-tasmin-hadcm3.nc')
nc <- open.ncdf(nc.fname, readunlim=FALSE)

# define which RCM it is & get the projection
rcm <- "hrm3"
proj.string <- get.rcm.projection(rcm)

# define WGS84 coordinate system
wgs84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

# Read in the grid and transform it to a SpatialGridDataFrame
nc.grid <- get.netcdf.grid(nc) 

## Special case for HRM3
## Round the xc & yc fields to 2 decimal degrees so they are a perfect grid.
## Define projection as wgs84 (even though it's _actually_ rotated pole)
if (rcm == "hrm3") {
  # define the coordinates & projection & then it will be a SpatialPointsDataFrame object from the rgdal library
  coordinates(nc.grid) <- c('xc', 'yc')
  proj4string(nc.grid) <- wgs84
} else {
  # define the coordinates & projection & then it will be a SpatialPointsDataFrame object from the rgdal library
  coordinates(nc.grid) <- c('xc', 'yc')
  proj4string(nc.grid) <- CRS(proj.string)
}

## here is the magic step, get the polygon into hrm3 rotated pole projection:
if (rcm == "hrm3") {
  shp <- get.hrm3.vector.coordinates(shp, proj.string)
} else {
  shp <- spTransform(shp, CRS(proj.string))
}

i <-  which(overlay(nc.grid, shp) != 'NA')

clipped <- nc.grid[i,]
coord.list <- clipped@data
## this function returns the clipped indicies in a format that is compatible with get.var.ncdf
contig <- find.continuous.columns(coord.list)
data <- NULL
for (region in contig) {
  x0 <- coord.list[region,'xi'][1]
  y0 <- coord.list[region,'yi'][1]
  x.count <- length(region)
  ## this next step, make sure that the third arguement in the start & count parameters are the time start & count you want.
  new.data <- get.var.ncdf(nc, start=c(x0, y0, 1), count=c(x.count, 1, 1))
  data <- rbind(data, new.data)
  rm(new.data)
}

# voila, the variable data now stores your clipped netcdf data with space & time dimensions.