R/cslice.R

In sethmcg/climod: Bias correction and other tools for climate model output

##' Climatological slicing of dates
##'
##' Given a sequence of dates, creates an object that indexes the sets
##' of dates falling into each of a series of climatological moving
##' windows.
##'
##' \code{cslice} generates indices for a set of inner and outer
##' windows.  The inner windows are contiguous and non-overlapping;
##' the outer windows overlap and are centered on the inner windows.
##' Applying an operation like the mean or the variance to the outer
##' window gives a moving-window statistic for the set of inner
##' windows.  For example, if the inner windows are 1 day long and the
##' outer windows are 30 days long, averaging over the outer windows
##' generates a daily climatology with a 30-day filter, which is
##' smoother and more robust than simply averaging over the Nth day of
##' each year.
##'
##' The windows generated by \code{cslice} are all of uniform size,
##' and there are the same number of inner and outer windows.  The
##' \code{slice} and \code{unslice} functions (q.v.) use a
##' \code{cslice} object to subset data for climatological operations.
##'
##' Typically, when people generate climatologies they either don't
##' use an inner window (e.g., when generating a monthly climatology
##' from daily data) or they use an inner window with the same
##' frequency as the data.  However, allowing the length of the inner
##' window to vary is useful for smoothing out irregularities and
##' differences in the calendar.  Using this approach to generate
##' daily climatologies with 360 inner windows for multiple different
##' datasets seamlessly accommodates model data with a 360-day
##' calendar, model data with a noleap (365-day) calendar, and real
##' world data with leap years.  Likewise, it will also handle
##' irregularly-spaced and missing timesteps.
##'
##' Caveats: the function operates on dates encoded as time elapsed
##' since an epoch (start time).  When generating climatologies from
##' multiple sources, be sure that the epochs match, or the
##' climatologies may not align with one another.  Using an inner
##' window that's shorter than the sampling frequency of the data
##' could result in some output values being based on many fewer
##' inputs than others; think carefully about window sizes and the
##' temporal structure of the data.
##'
##' If the indexes are split into contiguous segments by year using
##' the \code{split} argument, any outer window segment that doesn't
##' contain any inner window indexes will be discarded.  This happens
##' when the first (last) outer window overlaps the start (end) of the
##' dataset; discarding these segments ensures a one-to-one
##' correspondence between inner and outer window segments and
##' prevents some of the intermediate statistics created in KDDM
##' bias-correction from being distorted by segments containing only a
##' few points.
##' 
##' The size and number of the inner and outer windows can be
##' specified using any two of the four parameters \code{num},
##' \code{ratio}, \code{inner}, and \code{outer}.  By default,
##' \code{inner} and \code{outer} are assumed to have units of days,
##' but other units should work as long as they are consistent with
##' the value of \code{year}.
##'
##' See the \link[=climod/doc/kddm-bc.html]{KDDM bias-correction
##' vignette} for examples.
##'
##' @param time A vector of times, encoded as time elapsed since some
##' start date.
##'
##' @param num The number of inner windows per year.
##'
##' @param ratio The length of the outer window as a multiple of the
##' length of the inner window.
##'
##' @param inner The length of the inner window.
##'
##' @param outer The length of the outer window.
##'
##' @param year The length of the year.
##'
##' @param names (optional) A vector of names for the windows (e.g.,
##' names of months).
##' 
##' @param split Logical: whether to split the indexes for each window
##' into contiguous segments (i.e., grouped by year).  Defaults to
##' TRUE.
##'
##' @return An object of class 'cslice' containing lists of indices
##' associated with moving windows across multiple years and a list of
##' the parameters specifying the object.
##'
##' @seealso \code{\link{slice}}
##' 
##' @export

# add plot, print, str? methods


cslice <- function(time,
                   num = ceiling(year / inner),
                   ratio = outer / inner,
                   inner = year / num,
                   outer = inner * ratio,
                   year  = yearlength(time),
                   names = NULL,
                   split = TRUE
                   ){

    for(a in c(num, ratio, inner, outer, year)){
      if(a < 0 | !is.finite(a)){
        stop("slicing parameters must be finite and positive")
      }
    }
    if(outer < inner){ stop("outer must be >= inner")}
    if(num < 1){ stop("num must be > 1")}
    if(num != as.integer(num)){ warning("non-integer number of inner windows")}
    
    N <- length(time)
    
    span <- diff(range(time))
    cycles <- span / year
  
    cs <- list(
        time      = time,
        inner     = vector(mode="list", length=num),
        outer     = vector(mode="list", length=num),
        params    = namelist(span, cycles, year, ratio, inner, outer, split)
        )

    class(cs) <- "cslice"
    
    names(cs$inner)    <- names
    names(cs$outer)    <- names
    
    for(i in 1:num){
      innerind <- which( (time - inner*(i-1)) %% year < inner)
      outerind <- which( (time - inner*(i-1/2) + outer/2 ) %% year < outer)
      
      if(!split) {
        cs$inner[[i]] <- innerind
        cs$outer[[i]] <- outerind
      } else {
        ## find segment boundaries
        ibdy <- which(diff(innerind) > 1)
        obdy <- which(diff(outerind) > 1)
        
        ## split vectors into segments
        isegs <- mapply(function(a,b){innerind[a:b]}, SIMPLIFY=FALSE,
                        c(0,ibdy)+1, c(ibdy,length(innerind)))
        osegs <- mapply(function(a,b){outerind[a:b]}, SIMPLIFY=FALSE,
                        c(0,obdy)+1, c(obdy,length(outerind)))
          
        ## discard outer segments not containing an inner segment
        osegs <- osegs[sapply(osegs, function(x){any(innerind %in% x)})]
        
        cs$inner[[i]] <- isegs
        cs$outer[[i]] <- osegs
      }
    }    
  return(cs)
}