R/downscale.R

In jhardenberg/rainfarmr: Stochastic Precipitation Downscaling with the RainFARM Method

#' Downscale a precipitation field
#'
#' @description Downscales an input precipitation matrix using a metagaussian
#' spectral field `f` previously generated with [initmetagauss()].
#' The target resolution is defined by the dimensions of `f`.
#' An optional weights array can be specified.
#' @author Jost von Hardenberg, \email{j.vonhardenberg@isac.cnr.it}
#' @param r matrix of precipitation data to downscale.
#' @param f matrix containing a complex spectrum generated by the
#' [initmetagauss()] function.
#' @param weights matrix with climatological weights generated with the
#' [rfweights()] function.
#' @param fglob logical to conserve global average over domain.
#' @param fsmooth logical to use smoothing for conservation.
#' If neither option is set precipitation is conserved over each coarse box.
#' @return The downscaled field, with the same dimensions as `f`.
#' @import stats
#' @export
#' @examples
#' # Make some sample synthetic rainfall data
#' r <- exp(rnorm(4 * 4))
#' dim(r) <- c(4, 4)
#' r
#' #           [,1]      [,2]      [,3]      [,4]
#' # [1,] 1.8459816 1.8536550 2.1600665 1.3102116
#' # [2,] 1.3851011 1.4647348 0.2708219 0.4571810
#' # [3,] 0.2492451 0.8227134 0.4790567 1.9320403
#' # [4,] 0.5985922 3.3065678 2.1282795 0.6849944
#' # Create help field `f` with logarithmic slope 1.7
#' # with `dim(f) = c(8 * 4 ,8 * 4)`
#' f <- initmetagauss(1.7, 8 * 4)
#' rd <- downscale(r, f, fsmooth=FALSE) 
#' # Verify that downscaled data maintained original box averages
#' agg(rd, 4) 
#' #           [,1]      [,2]      [,3]      [,4]
#' # [1,] 1.8459816 1.8536550 2.1600665 1.3102116
#' # [2,] 1.3851011 1.4647348 0.2708219 0.4571810
#' # [3,] 0.2492451 0.8227134 0.4790567 1.9320403
#' # [4,] 0.5985922 3.3065678 2.1282795 0.6849944
downscale <- function(r, f, weights = 1., fglob = FALSE, fsmooth = TRUE) {
  nas <- dim(r)[1]
  ns <- dim(f)[1]
  r[r < 0] <- 0.
  rg <- gaussianize(r)
  g <- metagauss(f)

  gm <- mergespec(rg, g, nas / 2)
  st <- sd(gm)
  if (st == 0) {
    st <- 1.
  }
  gm <- gm / st
  fm <- exp(gm)
  fm <- fm * weights
  fm[ !is.finite(r) ] <- NA

  if (fglob) {
    imask <- !is.na(fm)
    ri <- interpola(r, ns)
    fm <- fm * mean(ri[imask]) / mean(fm[imask])
  } else if (fsmooth) {
    imask <- is.na(fm)
    fmi <- interpola(agg(fm, nas), ns);
    fmi[imask] <- NA
    fma <- smoothconv(fmi, nas)
    ri <- interpola(r, ns)
    ri[imask] <- NA
    raa <- smoothconv(ri, nas);
    fm <- raa / fma * fm
  } else {
    raa <- agg(r, nas)
    fma <- agg(fm, nas)
    ca <- raa / fma
    cai <- interpola(ca, ns)
    fm <- cai * fm
  }
  return(fm)
}