## 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/

library(rgdal)
library(raster)
library(ncdf)

# rcm proj4 string
get.rcm.projection <- function(rcm) {
  switch(rcm,
         crcm = "+proj=stere +lat_0=90.0 +lat_ts=60.0 +lon_0=263.0 +y_0=7475000.0 +x_0=3475000.0 +a=6371000 +b=6371000 +units=m",
         ecp2 = "+proj=stere +lat_0=90.0 +y_0=8400000 +x_0=4700000 +lat_ts=60.0 +lon_0=263.0 +a=6371205 +b=6371205 +units=m",
         hrm3 = "+proj=ob_tran +o_proj=latlon +o_lon_p=83 +o_lat_p=42.5 +lon_0=180", 
         mm5i = "+proj=lcc +lat_0=47.5 +lat_1=30.0 +lat_2=60.0 +lon_0=-97.0 +y_0=3200000.0 +x_0=3825000.0 +a=6370000 +b=6370000 +units=m",
         rcm3 = "+proj=somerc +lat_0=47.5 +y_0=3150000.0 +x_0=3900000.0 +lon_0=-97.0 +a=6371229.0 +b=6371229.0 +units=m", 
         wrfg = "+proj=lcc +lat_0=47.5 +lat_1=30.0 +y_0=2700000.0 +lat_2=60.0 +x_0=3325000.0 +lon_0=-97 +a=6370000 +b=6370000 +units=m",
         stop(paste("Projection for RCM", rcm, "not supported"))
         )
}


# a function to create a grid from the input netcdf files
get.netcdf.grid <- function(nc, degrees.east=T) {
  lat <- get.var.ncdf(nc, "lat")
  lon <- get.var.ncdf(nc, "lon")
  if (degrees.east) {
	#convert grid from eliptical longitude (0-360) to -180 to 180 longitude
    lon <- ((lon + 180) %% 360) - 180
  }
  ## For NARCCAP like grids
  if ("xc" %in% names(nc$dim)) {
    i <- expand.grid(seq(nc$dim$xc$len), seq(nc$dim$yc$len))
    names(i) <- c("xi", "yi")
    xc <- as.vector(nc$dim$xc$vals[i$xi])
    yc <- as.vector(nc$dim$yc$vals[i$yi])
    ll <- as.vector(apply(i, 1, function(x) {lat[x[1], x[2]]}))
    ln <- as.vector(apply(i, 1, function(x) {lon[x[1], x[2]]}))
    return(data.frame(lon=ln, lat=ll, xc=xc, yc=yc, i))
  }
  ## For GCM grids
  else if ("lon" %in% names(nc$dim)) {
    i <- expand.grid(seq(nc$dim$lon$len), seq(nc$dim$lat$len))
    names(i) <- c("xi", "yi")
    lon <- as.vector(nc$dim$lon$vals[i$xi])
    i <- expand.grid(seq(nc$dim$lon$len), seq(nc$dim$lat$len))
    names(i) <- c("xi", "yi")
    lon <- as.vector(nc$dim$lon$vals[i$xi])
    if (degrees.east) {
      lon <- ((lon + 180) %% 360) - 180
    }
    lat <- as.vector(nc$dim$lat$vals[i$yi])
    return(data.frame(lon=lon, lat=lat, i))
  }
  ## For other grids with projected coordinates (ex. NARR)
  else if ("x" %in% names(nc$dim)) {
    i <- expand.grid(seq(nc$dim$x$len), seq(nc$dim$y$len))
    names(i) <- c("xi", "yi")
    xc <- as.vector(nc$dim$x$vals[i$xi])
    yc <- as.vector(nc$dim$y$vals[i$yi])
    ll <- as.vector(apply(i, 1, function(x) {lat[x[1], x[2]]}))
    ln <- as.vector(apply(i, 1, function(x) {lon[x[1], x[2]]}))
    return(data.frame(lon=ln, lat=ll, xc=xc, yc=yc, i))
  }
}

## hrm3 projection is not very compliant. Must bend it to my will...
## There is no known proj4 string for hrm3... to convert it we have to:
## 1) set the CRS to be wgs84
## 2) reproject the vector to wgs84
## 3) convert the vector coordinates to radians
## 4) define vector projection to be oblique transverse (rotated pole) which is actually hrm3's projection
## 5) reproject vector to wgs84
## 6) Now the vector is in rotated pole coordinate space (WIN!)
get.hrm3.vector.coordinates <- function(vect, proj.string){
  cls <- class(vect)
  wgs84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
  rad <- 360 / (2 * pi)
  vect <- spTransform(vect, wgs84)
  if (cls == "SpatialPolygonsDataFrame") {
    
    for (i in seq(1:length(vect@polygons))){
      for (j in seq(1:length(vect@polygons[[i]]@Polygons))){
        vect@polygons[[i]]@Polygons[[j]]@coords <- vect@polygons[[i]]@Polygons[[j]]@coords/rad
        vect@polygons[[i]]@Polygons[[j]]@labpt <- vect@polygons[[i]]@Polygons[[j]]@labpt/rad
      }
    }
    #vect@bbox <- vect@bbox/rad
  } else if (cls == "SpatialLinesDataFrame"){
    for (i in seq(1:length(vect@lines))){
      for (j in seq(length(vect@lines[[i]]@Lines))){
        vect@lines[[i]]@Lines[[j]]@coords <- vect@lines[[i]]@Lines[[j]]@coords/rad
      }
    }
    #vect@bbox <- vect@bbox/rad
  } else if (grepl("SpatialPoints", cls)){
    vect@coords <- vect@coords/rad
    #vect@bbox <- vect@bbox/rad
  } 
  if (class(proj.string) == "CRS"){
    proj4string(vect) <- proj.string
  } else {
    proj4string(vect) <- CRS(proj.string)
  }
  vect <- spTransform(vect, wgs84)
  return(vect)
}

# Returns a list of rows which for which all x locations are contiguous
# df should be a data.frame with columns xi and yi
# The purpose of this is to limit the number of calls to get.var.ncdf.  Knowing
# contiguous regions allows the user to only make one call per contiguous region
#
# An improvement on this would be to return the maximum sized rectangles, but this
# is simple and will do for now
find.continuous.columns <- function(df) {
  breaks <- which(! diff(df$xi) == 1)
  lengths <- diff(c(0, breaks, nrow(df)))
  mapply(seq, from=c(0, breaks)+1, length.out=lengths)
}