##' Bias correction of climate model output using KDDM.
##'
##' This is the description of the function. It's the first full
##' paragraph, and should briefly describe what the function does.
##'
##' Third and subsequent paragraphs are details: a long section shown
##' after the argument that goes into detail about how the function
##' works.
##'
##' @param bcdata The data to be bias-corrected. A list with elements
##' named "obs", "cur", and "fut".
##'
##' @param norm The type of normalization to use. See
##' \code{\link{normalize}} for more details.
##'
##' @param minobs The minimum number of data points to attempt a bias
##' correction. If there are fewer that this many values in the
##' observations or in the current period model output (e.g., for
##' precipitation in a very dry location), the bias correction falls
##' back on a simple scaling by the ratio of obs mean to cur mean.
##' (If the ratio is undefined, future values are passed through
##' unchanged.)
##'
##' @param dmap Logical; if TRUE, returns the \code{\link{distmap}}
##' object generated during bias correction as an element named "dmap"
##' in the returned list. N.B.: This option defaults to FALSE because
##' the distmap object is significantly larger than the inputs and
##' outputs of this function combined.
##'
##' @param ... Arguments to pass to \code{\link{distmap}}.
##'
##' @return The return value of the function.
##'
##' @examples
##' x <- "example code"
##' y <- "more example code"
##' paste(x, y)
##'
##' @importFrom stats predict
##'
##' @export
biascorrect <- function(bcdata, norm="zscore", minobs=250, dmap=FALSE, ...){
## If there's too little data to get a good estimate of the PDF for
## either obs or cur, it's not possible to bias-correct using KDDM;
## fall back to simple ratio scaling. Also set an uncorrectable
## attribute, in case knowing where this has happened is of
## interest.
if(sum(is.finite(unlist(bcdata$obs))) < minobs |
sum(is.finite(unlist(bcdata$cur))) < minobs ){
muobs <- mean(unlist(bcdata$obs), na.rm=TRUE)
mucur <- mean(unlist(bcdata$cur), na.rm=TRUE)
ratio <- muobs / mucur
if(!is.finite(ratio)){
muobs <- 1
mucur <- 1
ratio <- 1
}
unfixable <- function(x){
x@uncorrectable <- is.finite(x)
x * ratio
}
result <- rapply(bcdata, unfixable, how="replace")
result$obs <- bcdata$obs
if(dmap){ result$distmap <- distmap(c(0,mucur),c(0,muobs)) }
return(result)
}
## normalize the three data components
if(norm=="boxcox"){
## You can't really apply an adjustment to the denormalization
## for Box-Cox; you need to use the same parameter values you
## used for the obs or things go wonky. This means that you
## also want to use the same params over all times in both cur
## and fut, so that you're not discarding any climate change
## signal when you normalize. Given that, it makes sense to use
## a single value for all three, so we pool data to fit gamma to
## get the best compromise value.
gamma <- megamma(unlist(bcdata, use.names=FALSE))
nbcd <- rapply(bcdata, normalize, how="replace", norm=norm, gamma=gamma)
} else {
nbcd <- rapply(bcdata, normalize, how="replace", norm=norm)
}
## construct distribution mapping
mapping <- distmap(unlist(nbcd$cur), unlist(nbcd$obs), ...)
## apply KDDM transfer function
## note: could skip predict(obs), but need its norm atts for next step
fixed <- rapply(nbcd, function(x){predict(mapping, x)}, how="replace")
## extract & average normalization attributes
adj <- rapply(fixed, attributes, how="replace")
adj <- lapply(adj, renest)
adj <- lapply(adj, function(x){lapply(x, function(y){unlist(y)})})
adj <- lapply(adj, function(x){x$norm<-NULL; x})
adj <- rapply(adj, mean, how="replace")
## Probably fix for entire month missing; need to test.
# adj <- rapply(adj, mean, how="replace", na.rm=TRUE)
adj <- renest(adj)
## calculate adjustments for denormalization
shift <- NA
scale <- NA
pscale <- NA
if(norm == "zscore"){
shift <- adj$mu$obs - adj$mu$cur
scale <- adj$sigma$obs / adj$sigma$cur
}
if(norm == "boxcox"){
shift <- 0
pscale <- 1
}
if(norm == "scale"){
scale <- adj$mu2$obs / adj$mu2$cur
}
## no adjustments for range or identity norms
## denormalize bias-corrected data
result <- list()
result$obs <- bcdata$obs
result$cur <- rapply(fixed$cur, denormalize, how="replace",
shift=shift, scale=scale, pscale=pscale)
result$fut <- rapply(fixed$fut, denormalize, how="replace",
shift=shift, scale=scale, pscale=pscale)
if(dmap){ result$distmap <- mapping }
return(result)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.