Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
#'@rdname CST_DynBiasCorrection
#'@title Performing a Bias Correction conditioned by the dynamical
#'properties of the data.
#'
#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it}
#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it}
#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr}
#'
#'@description This function perform a bias correction conditioned by the
#'dynamical properties of the dataset. This function internally uses the functions
#''Predictability' to divide in terciles the two dynamical proxies
#'computed with 'CST_ProxiesAttractor'. A bias correction
#'between the model and the observations is performed using the division into
#'terciles of the local dimension 'dim' and inverse of the persistence 'theta'.
#'For instance, model values with lower 'dim' will be corrected with observed
#'values with lower 'dim', and the same for theta. The function gives two options
#'of bias correction: one for 'dim' and/or one for 'theta'
#'
#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D.,
#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large
#'scale atmospheric predictability.Nature Communications, 10(1), 1316.
#'DOI = https://doi.org/10.1038/s41467-019-09305-8 "
#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017).
#' Dynamical proxies of North Atlantic predictability and extremes.
#' Scientific Reports, 7-41278, 2017.
#'
#'@param exp an s2v_cube object with the experiment data
#'@param obs an s2dv_cube object with the reference data
#'@param method a character string indicating the method to apply bias correction among these ones:
#'"PTF","RQUANT","QUANT","SSPLIN"
#'@param proxy a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method.
#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta
#'@param time_dim a character string indicating the name of the temporal dimension
#'@param ncores The number of cores to use in parallel computation
#'
#'@return dynbias an s2dvcube object with a bias correction performed
#'conditioned by local dimension 'dim' or inverse of persistence 'theta'
#'
#'@examples
#'# example 1: simple data s2dvcube style
#'set.seed(1)
#'expL <- rnorm(1:2000)
#'dim (expL) <- c(time =100,lat = 4, lon = 5)
#'obsL <- c(rnorm(1:1980),expL[1,,]*1.2)
#'dim (obsL) <- c(time = 100,lat = 4, lon = 5)
#'time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-")
#'time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-")
#'lon <- seq(-1,5,1.5)
#'lat <- seq(30,35,1.5)
# qm=0.98 # too high for this short dataset, it is possible that doesn't
#'# get the requirement, in that case it would be necessary select a lower qm
#'# for instance qm=0.60
#'expL <- s2dv_cube(data = expL, lat = lat, lon = lon,
#' Dates = list(start = time_expL, end = time_expL))
#'obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon,
#' Dates = list(start = time_obsL, end = time_obsL))
#'dynbias <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim",
#' quanti = 0.6, time_dim = 'time')
#'
#'@export
CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT',
proxy = "dim", quanti, time_dim = 'ftime',
ncores = NULL) {
if (!inherits(obs, 's2dv_cube')) {
stop("Parameter 'obs' must be of the class 's2dv_cube', ",
"as output by CSTools::CST_Load.")
}
if (!inherits(exp, 's2dv_cube')) {
stop("Parameter 'exp' must be of the class 's2dv_cube', ",
"as output by CSTools::CST_Load.")
}
exp$data <- DynBiasCorrection(exp = exp$data, obs = obs$data, method = method,
proxy = proxy, quanti = quanti,
time_dim = time_dim, ncores = ncores)
return(exp)
}
#'@rdname DynBiasCorrection
#'@title Performing a Bias Correction conditioned by the dynamical
#'properties of the data.
#'
#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it}
#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it}
#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr}
#'
#'@description This function perform a bias correction conditioned by the
#'dynamical properties of the dataset. This function used the functions
#''CST_Predictability' to divide in terciles the two dynamical proxies
#'computed with 'CST_ProxiesAttractor'. A bias correction
#'between the model and the observations is performed using the division into
#'terciles of the local dimension 'dim' and inverse of the persistence 'theta'.
#'For instance, model values with lower 'dim' will be corrected with observed
#'values with lower 'dim', and the same for theta. The function gives two options
#'of bias correction: one for 'dim' and/or one for 'theta'
#'
#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D.,
#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large
#'scale atmospheric predictability.Nature Communications, 10(1), 1316.
#'DOI = https://doi.org/10.1038/s41467-019-09305-8 "
#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017).
#' Dynamical proxies of North Atlantic predictability and extremes.
#' Scientific Reports, 7-41278, 2017.
#'
#'@param exp a multidimensional array with named dimensions with the experiment data
#'@param obs a multidimensional array with named dimensions with the observation data
#'@param method a character string indicating the method to apply bias correction among these ones:
#'"PTF","RQUANT","QUANT","SSPLIN"
#'@param proxy a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method.
#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta
#'@param time_dim a character string indicating the name of the temporal dimension
#'@param ncores The number of cores to use in parallel computation
#'
#'@return a multidimensional array with named dimensions with a bias correction performed conditioned by local dimension 'dim' or inverse of persistence 'theta'
#'
#'@import multiApply
#'@importFrom s2dverification Subset
#'@import qmap
#'@examples
#'expL <- rnorm(1:2000)
#'dim (expL) <- c(time =100,lat = 4, lon = 5)
#'obsL <- c(rnorm(1:1980),expL[1,,]*1.2)
#'dim (obsL) <- c(time = 100,lat = 4, lon = 5)
#'dynbias <- DynBiasCorrection(exp = expL, obs = obsL,
#' proxy= "dim", quanti = 0.6)
#'@export
DynBiasCorrection<- function(exp, obs, method = 'QUANT',
proxy = "dim", quanti,
time_dim = 'time', ncores = NULL){
if (is.null(obs)) {
stop("Parameter 'obs' cannot be NULL.")
}
if (is.null(exp)) {
stop("Parameter 'exp' cannot be NULL.")
}
if (is.null(method)) {
stop("Parameter 'method' cannot be NULL.")
}
if (is.null(quanti)) {
stop("Parameter 'quanti' cannot be NULL.")
}
if (is.null(proxy)) {
stop("Parameter 'proxy' cannot be NULL.")
}
dims <- dim(exp)
attractor.obs <- ProxiesAttractor(data = obs, quanti = quanti)
predyn.obs <- Predictability(dim = attractor.obs$dim,
theta = attractor.obs$theta)
attractor.exp <- ProxiesAttractor(data = exp, quanti = quanti)
predyn.exp <- Predictability(dim = attractor.exp$dim,
theta = attractor.exp$theta)
if (proxy == "dim") {
adjusted <- Apply(list(exp, obs), target_dims = time_dim,
fun = .dynbias, method,
predyn.exp = predyn.exp$pred.dim$pos.d,
predyn.obs = predyn.obs$pred.dim$pos.d,
ncores = ncores, output_dims = time_dim)$output1
} else if (proxy == "theta") {
adjusted <- Apply(list(exp, obs), target_dims = time_dim,
fun = .dynbias, method,
predyn.exp = predyn.exp$pred.theta$pos.t,
predyn.obs = predyn.obs$pred.theta$pos.t,
ncores = ncores, output_dims = time_dim)$output1
} else {
stop ("Parameter 'proxy' must be set as 'dim' or 'theta'.")
}
return(adjusted)
}
.dynbias <- function(exp, obs, method, predyn.exp, predyn.obs) {
result <- array(rep(NA, length(exp)))
res <- lapply(1:3, function(x) {
exp_sub <- exp[predyn.exp[[x]]]
obs_sub <- obs[predyn.obs[[x]]]
adjust <- .qbiascorrection(exp_sub, obs_sub, method)
result[predyn.exp[[x]]] <<- adjust
return(NULL)
})
return(result)
}
.qbiascorrection <- function(expX, obsX, method) {
## functions fitQmap and doQmap
if (method == "PTF") {
qm.fit <- fitQmap(obsX, expX, method = "PTF", transfun = "expasympt",
cost = "RSS", wett.day = TRUE)
qmap <- doQmap(expX, qm.fit)
} else if (method == "QUANT") {
qm.fit <- fitQmap(obsX, expX, method = "QUANT", qstep = 0.01)
qmap <- doQmap(expX, qm.fit, type = "tricub")
} else if (method == "RQUANT") {
qm.fit <- fitQmap(obsX, expX, method = "RQUANT", qstep = 0.01)
qmap <- doQmap(expX, qm.fit, type = "linear")
} else if (method == "SSPLIN") {
qm.fit <- fitQmap(obsX, expX, qstep = 0.01, method = "SSPLIN")
qmap <- doQmap(expX, qm.fit)
} else {
stop ("Parameter 'method' doesn't match any of the available methods.")
}
return(qmap)
}