R/detrendGrid.R

In SantanderMetGroup/transformeR: A climate4R package for general climate data manipulation and transformation

#     detrendGrid.R Linear detrending of a grid
#
#     Copyright (C) 2017 Santander Meteorology Group (http://www.meteo.unican.es)
#
#     This program is free software: you can redistribute it and/or modify
#     it under the terms of the GNU General Public License as published by
#     the Free Software Foundation, either version 3 of the License, or
#     (at your option) any later version.
# 
#     This program is distributed in the hope that it will be useful,
#     but WITHOUT ANY WARRANTY; without even the implied warranty of
#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#     GNU General Public License for more details.
# 
#     You should have received a copy of the GNU General Public License
#     along with this program.  If not, see <http://www.gnu.org/licenses/>.

#' @title Linear detrending
#' @description Perform a linear detrending along the time dimension of a grid
#' @param grid Input grid (possibly multimember)
#' @param grid2 Optional grid. If provided, the output is the detrended \code{grid2}
#'  using the regression coefficients calculated using \code{grid}. Default to \code{NULL}.
#' @template templateParallelParams
#' @return A detrended grid. 
#' @details  Performs a simple linear detrending by fitting a linear model and retaining the residuals.
#' An attribute indicating the linear detrending is added to the \code{Variable} component of the output grid.
#' 
#' In the presence of missing data in the time series, it operates by filtering them prior to linear model fitting. The
#' missing data positions are then restored back to the output detrended series.
#' 
#' The function uses the \code{\link{fastLm}} implementation from package \pkg{RcppEigen}, significantly speeding-up 
#' the linear model fitting.
#' 
#' @template templateParallel
#' @export
#' @importFrom parallel stopCluster parApply
#' @importFrom stats coef
#' @importFrom magrittr %>%
#' @importFrom RcppEigen fastLm
#' @export
#' @author J Bedia, J Fernandez, M.D. Frias
#' @examples  \donttest{
#' require(climate4R.datasets) 
#' data("NCEP_Iberia_ta850")
#' monthly <- aggregateGrid(NCEP_Iberia_ta850, aggr.m = list(FUN = "mean"))
#' plot(monthly$Data[,4,2], ty = 'l')
#' abline(reg = lm(monthly$Data[,4,2] ~ I(1:length(monthly$Data[,4,2]))), lty = 2)
#' det <- detrendGrid(monthly, parallel = FALSE)
#' # Detrended series in red
#' lines(det$Data[,4,2], col = "red")
#' abline(reg = lm(det$Data[,4,2] ~ I(1:length(det$Data[,4,2]))), col = "red", lty = 2)
#' legend("topright", c("Raw", "Detrended"), lty = 1, col = c(1,2))
#' ## example of forecast detrend
#' # Suppose 1991-2001 is the reference period. We use it to detrend data from year 2002
#' grid <- subsetGrid(monthly, years = 1991:2001)
#' grid2 <- subsetGrid(monthly, years = 2002)
#' det2 <- detrendGrid(grid, grid2, parallel = FALSE)
#' # det2 is the grid2 (year 2002) corrected according to the linear trend computed for 1991-2001
#' }


detrendGrid <- function(grid,
                        grid2 = NULL,
                        parallel = FALSE,
                        max.ncores = 16,
                        ncores = NULL) {
      grid <- redim(grid, var = TRUE, member = TRUE)
      parallel.pars <- parallelCheck(parallel, max.ncores, ncores)
      apply_fun <- selectPar.pplyFun(parallel.pars, .pplyFUN = "apply")
      if (parallel.pars$hasparallel) on.exit(parallel::stopCluster(parallel.pars$cl))
      arr <- grid$Data
      x <- as.numeric(as.Date(getRefDates(grid)))
      dimNames <- getDim(grid)
      mar <- grep("^time$", dimNames, invert = TRUE)
      message("[", Sys.time(), "] - Detrending...")
      if (is.null(grid2)) {
            arr <- unname(apply_fun(X = arr, MARGIN = mar, FUN = function(y) {
                  out <- rep(NA, length(x))
                  ind <- intersect(which(!is.na(y)), 1:length(x))
                  out[ind] <- tryCatch(expr = summary(RcppEigen::fastLm(y ~ I(x)))$resid + mean(y, na.rm = TRUE),
                                       error = function(err) {
                                             out
                                       })
                  return(out)
            }))
            dimNames <- getDim(grid) 
            perm <- match(dimNames, c("time", "var", "member", "lat", "lon"))
            arr <- aperm(arr, perm)
      } else {
            grid2 <- redim(grid2, var = TRUE, member = TRUE)
            dims <- getShape(grid2)
            x2 <- as.numeric(as.Date(getRefDates(grid2)))
            arr <- unname(apply_fun(X = arr, MARGIN = mar, FUN = function(y) {
                  tryCatch(expr = {
                        coefs <- coef(RcppEigen::fastLm(y ~ I(x)))
                        coefs[1] + x2 * coefs[2]
                  },
                  error = function(err) {
                        rep(NA, length(x2))
                  })
            }))
            # apply drops the time dimension if a singleton
            # (https://radfordneal.wordpress.com/2008/08/20/design-flaws-in-r-2-%E2%80%94-dropped-dimensions/)
            if (length(dims) != length(dim(arr))) arr <- unname(abind(arr, along = -1L))
            dimNames <- getDim(grid2) 
            perm <- match(dimNames, c("time", "var", "member", "lat", "lon"))
            arr <- aperm(arr, perm)
            resid <- grid2$Data - arr
            clim <- redim(grid, drop = TRUE) %>% climatology(by.member = FALSE) %>% redim(var = TRUE, member = TRUE) %>% getElement("Data")
            nmem <- getShape(grid2, "member")
            ntime <- getShape(grid2, "time")
            parallel.pars <- parallelCheck(parallel, max.ncores, ncores)
            lapply_fun <- selectPar.pplyFun(parallel.pars, .pplyFUN = "lapply")
            arr <- lapply_fun(1:nmem, function(x) {
                  aux <- asub(resid, idx = x, dims = grep("member", dimNames), drop = FALSE)
                  aux1 <- lapply(1:ntime, function(y) {
                        aux1 <- asub(aux, idx = y, dims = grep("^time", dimNames), drop = FALSE)
                        aux1 + clim
                  })
                  unname(do.call("abind", c(aux1, along = grep("^time", dimNames))))
            })
            arr <- unname(do.call("abind", c(arr, along = grep("member", dimNames))))
            grid <- grid2
      }
      message("[", Sys.time(), "] - Done.")
      grid$Data <- arr
      attr(grid$Data, "dimensions") <- dimNames
      grid <- redim(grid, drop = TRUE) 
      attr(grid$Variable, "detrended:method") <- "linear"
      return(grid)
}
# End