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. The calibrated forecasts have an equivalent 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)
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#'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 a list of s2dverification objects as returned by ",
"the s2dverification::Load function.")
}
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 (length(which(is.na(exp))) > 0) {
warning('There are NA in exp.')
}
if (length(which(is.na(obs))) > 0) {
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 <- length(var_obs)
if (dim(var_exp)[1] != dim(var_exp)[2]) {
nmembers <- dim(var_exp)[-which(dim(var_exp) == length(var_obs))]
} else{
nmembers <- dim(var_exp)[1]
}
if (!all(dim(var_exp) == c(nmembers, ntime))) {
var_exp <- t(var_exp)
}
if (!all(dim(var_exp) == c(nmembers, ntime))) {
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)
climPred[t] <- mean(var_exp[,-t])
}
var_obs <- var_obs - climObs
for (i in 1:nmembers) {
var_exp[i,] <- var_exp[i,] - climPred
}
calibrated <- NA * var_exp
for (t in 1:ntime) {
# defining forecast,hindcast and observation in cross-validation
fcst <- var_exp[, t]
hcst <- var_exp[,-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)
}