R/extend.global.data.R

In jjvanderwal/climates: Methods for working with weather & climate

#' Extend Extents of Global Data to Interpolate Spherical Data
#' 
#'  a function to extend the extents of a GCM in preparation for bilinear downscaling
#' 
#'  inputs are from retrieve.nc, append.nc or extract.monthly.averages
#' 
#' @param nc  ~~Describe \code{nc} here~~
#' @return ...
#' @note  ~~further notes~~
#' @author  ~~who you are~~
#' @seealso  ~~objects to See Also as \code{\link{help}}, ~~~
#' @references  ~put references to the literature/web site here ~
#' @keywords ~kwd1 ~kwd2
#' @export
#' @examples
#' \dontrun{
#' ##---- Should be DIRECTLY executable !! ----
#' ##-- ==>  Define data, use random,
#' ##--	or do  help(data=index)  for the standard data sets.
#' 
#' ## The function is currently defined as
#' function(nc) {
#' 	#ensure nc1 & nc2 are from retrieve.nc
#' 	if (all(c('dat','tim','lat', 'lon') %in% names(nc))==FALSE) stop('nc must have 
#'   objects named dat, lat, lon and tim as from retrieve.nc of clim.pact package, 
#'   append.nc or extract.monthly.averages')
#' 	#do the work
#' 	#assuming this is a global model, we need to buffer the edges to interpolate
#' 	#adjust the latitude
#' 	nc$lat = c(nc$lat[1]-mean(diff(nc$lat),na.rm=T),
#'   nc$lat,nc$lat[length(nc$lat)]+mean(diff(nc$lat),na.rm=T))
#' 	if(min(nc$lat,na.rm=T)>-90) nc$lat[nc$lat==min(nc$lat,na.rm=T)] = -90
#' 	if(max(nc$lat,na.rm=T)<90) nc$lat[nc$lat==max(nc$lat,na.rm=T)] = 90
#' 	#adjust the longitude
#' 	nc$lon = c(nc$lon[1]-mean(diff(nc$lon),na.rm=T),
#'   nc$lon,nc$lon[length(nc$lon)]+mean(diff(nc$lon),na.rm=T))
#' 	if(min(nc$lon,na.rm=T)>-180) nc$lon[nc$lon==min(nc$lon,na.rm=T)] = -180
#' 	if(max(nc$lon,na.rm=T)<180) nc$lon[nc$lon==max(nc$lon,na.rm=T)] = 180
#' 	#populate a new set of outa with new dim of lat & lon
#' 	out = array(0,dim=c(length(nc$tim),length(nc$lat),length(nc$lon)))
#' 	out[1:dim(nc$dat)[1],1:dim(nc$dat)[2]+1,1:dim(nc$dat)[3]+1] = nc$dat[,,]
#' 	#adjust the latitude data
#' 	tmid = ceiling(dim(nc$dat)[3] / 2)
#' 	out[,1,] = cbind(NA,nc$dat[,1,tmid:dim(nc$dat)[3]],nc$dat[,1,1:(tmid-1)],NA)
#' 	out[,dim(out)[2],] = cbind(NA,nc$dat[,dim(nc$dat)[2],tmid:dim(nc$dat)[3]],
#'   nc$dat[,dim(nc$dat)[2],1:(tmid-1)],NA)
#' 	#adjust the longitude information... appending ends from other side of globe
#' 	out[,,1] = out[,,dim(out)[3]-1]
#' 	out[,,dim(out)[3]] = out[,,2]
#' 	#return the extended outa
#' 	return(list(dat=out,lon=nc$lon,lat=nc$lat,tim=nc$tim))
#'   }
#' }
extend.global.data = function(nc) {
	#ensure nc1 & nc2 are from retrieve.nc
	if (all(c('dat','tim','lat', 'lon') %in% names(nc))==FALSE) stop('nc must have objects named dat, lat, lon and tim as from retrieve.nc of clim.pact package, append.nc or extract.monthly.averages')
	#do the work
	#assuming this is a global model, we need to buffer the edges to interpolate
	#adjust the latitude
	nc$lat = c(nc$lat[1]-mean(diff(nc$lat),na.rm=T),nc$lat,nc$lat[length(nc$lat)]+mean(diff(nc$lat),na.rm=T))
	if(min(nc$lat,na.rm=T)>-90) nc$lat[nc$lat==min(nc$lat,na.rm=T)] = -90
	if(max(nc$lat,na.rm=T)<90) nc$lat[nc$lat==max(nc$lat,na.rm=T)] = 90
	#adjust the longitude
	nc$lon = c(nc$lon[1]-mean(diff(nc$lon),na.rm=T),nc$lon,nc$lon[length(nc$lon)]+mean(diff(nc$lon),na.rm=T))
	if(min(nc$lon,na.rm=T)>-180) nc$lon[nc$lon==min(nc$lon,na.rm=T)] = -180
	if(max(nc$lon,na.rm=T)<180) nc$lon[nc$lon==max(nc$lon,na.rm=T)] = 180
	#populate a new set of outa with new dim of lat & lon
	out = array(0,dim=c(length(nc$tim),length(nc$lat),length(nc$lon)))
	out[1:dim(nc$dat)[1],1:dim(nc$dat)[2]+1,1:dim(nc$dat)[3]+1] = nc$dat[,,]
	#adjust the latitude data
	tmid = ceiling(dim(nc$dat)[3] / 2)
	out[,1,] = cbind(NA,nc$dat[,1,tmid:dim(nc$dat)[3]],nc$dat[,1,1:(tmid-1)],NA)
	out[,dim(out)[2],] = cbind(NA,nc$dat[,dim(nc$dat)[2],tmid:dim(nc$dat)[3]],nc$dat[,dim(nc$dat)[2],1:(tmid-1)],NA)
	#adjust the longitude information... appending ends from other side of globe
	out[,,1] = out[,,dim(out)[3]-1]
	out[,,dim(out)[3]] = out[,,2]
	#return the extended outa
	out = list(dat=out,lon=nc$lon,lat=nc$lat,tim=nc$tim)
	class(out) = unique(c(class(out),'nc'))
	return(out)
}