Newer
Older
#' Calibration of a CSTools object based on the ensemble inflation of a forecast following Doblas-Reyes et al. (2005)
#'
#'@author Verónica Torralba, \email{veronica.torralba@bsc.es}
#'@description This function applies a variance inflation technique described in Doblas-Reyes et al. (2005) in leave-one-out cross-validation. This bias adjustment method produces calibrated forecasts with equivalent mean and variance to that of the reference dataset, but at the same time preserve reliability.
#'
#'@param data a CSTools object (an s2dverification object as output by the \code{Load} function from the s2dverification package).
#'
#'@return a CSTools object (s2dverification object) with the calibrated forecasts (provided in $mod) with the same dimensions as data$mod.
#'
#'@import s2dverification
#'@examples
#'
#'# Example
#'# Creation of sample s2dverification objects. These are not complete
#'# s2dverification objects though. The Load function returns complete objects.
#'mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
#'dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
#'obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
#'dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
#'lon <- seq(0, 30, 5)
#'lat <- seq(0, 25, 5)
#'data1 <- list(mod = mod1, obs = obs1, lat = lat, lon = lon)
#'str(a)
#'@export
CST_Calibration <- function(data) {
is_s2dv_object <- function(x) {
if (all(c('mod', 'obs', 'lon', 'lat') %in% names(x))) { #&& length(x) == 11) {
TRUE
} else {
FALSE
}
}
wrong_input <- FALSE
if (!is.list(data)) {
wrong_input <- TRUE
}
if (!is_s2dv_object(data)) {
wrong_input <- TRUE
}
if (wrong_input) {
stop("Parameter 'data' must be an object of s2dverification as returned by ",
"the s2dverification::Load function.")
}
if (dim(data$obs)[which(names(dim(data$obs))=='member')]!=1) {
stop("The length of the dimension 'member' in 'data$obs must be equal to 1")
}
Calibrated <- Calibration(data$mod, data$obs)
Calibrated <- aperm(Calibrated, c(3, 1, 2, 4, 5, 6))
data$mod <- Calibrated
data$obs <- NULL
return(data)
}
Calibration <- function(exp, obs) {
if (!all(c('member', 'sdate') %in% names(dim(exp)))) {
stop("Parameter 'exp' must have the dimensions 'member' and 'sdate'.")
}
if (!all(c('sdate') %in% names(dim(obs)))) {
stop("Parameter 'obs' must have the dimension 'sdate'.")
}
if (any(is.na(exp))) {
warning('There are NA in exp.')
}
if (any(is.na(obs))) {
warning('There are NA in obs.')
}
target_dims_obs <- 'sdate'
if ('member' %in% names(dim(obs))) {
target_dims_obs <- c('member', target_dims_obs)
}
Calibrated <- Apply(data = list(var_obs = obs, var_exp = exp),
target_dims = list(target_dims_obs, c('member', 'sdate')),
fun = '.cal')$output1
return(Calibrated)
}
.cal <- function(var_obs, var_exp) {
ntime <- dim(var_exp)[which(names(dim(var_exp))=='sdate')][]
nmembers<-dim(var_exp)[which(names(dim(var_exp))=='member')][]
if (all(names(dim(var_exp))!=c('member','sdate'))) {
var_exp <- t(var_exp)
}
climObs <- NA * var_obs
climPred <- NA * var_obs
for (t in 1:length(var_obs))
{
climObs[t] <- mean(var_obs[,-t])
climPred[t] <- mean(var_exp[,-t])
}
var_obs <- var_obs - climObs
calibrated <- NA * var_exp
for (t in 1:ntime) {
# defining forecast,hindcast and observation in cross-validation
fcst<-NA*var_exp[,t]
hcst<-NA*var_exp[,-t]
for (i in 1:nmembers) {
fcst[i] <- var_exp[i,t] - climPred[t]
hcst[i,]<- var_exp[i,-t]- climPred[t]
}
obs <- var_obs[-t]
#coefficients
em_fcst <- mean(fcst)
em_hcst <- apply(hcst, c(2), mean)
corr <- cor(em_hcst, obs)
sd_obs <- sd(obs)
sd_em_hcst <- sd(em_hcst)
fcst_diff <- fcst - em_fcst
hcst_diff <- NA * hcst
for (n in 1:nmembers) {
hcst_diff[n,] <- hcst[n,] - em_hcst
}
sd_hcst_diff <- sd(hcst_diff)
a <- corr * (sd_obs / sd_em_hcst)
b <- (sd_obs / sd_hcst_diff) * sqrt(1 - (corr ^ 2))
calibrated[, t] <- (a * em_fcst) + (b * fcst_diff) + climObs[t]
}
names(dim(calibrated)) <- c('member', 'sdate')
return(calibrated)
}