Newer
Older
#'@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 wetday logical indicating whether to perform wet day correction
#'or not OR a numeric threshold below which all values are set to zero (by
#'default is set to 'FALSE').
#'@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 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))
#' # to use DynBiasCorrection
#' dynbias1 <- DynBiasCorrection(exp = expL$data, obs = obsL$data, proxy= "dim",
#' quanti = 0.6)
#' # to use CST_DynBiasCorrection
#' dynbias2 <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim",
CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', wetday=FALSE,
proxy = "dim", quanti,
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,
Verónica Torralba-Fernández
committed
wetday=wetday,
proxy = proxy, quanti = quanti, ncores = ncores)
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
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:
#'@param wetday logical indicating whether to perform wet day correction
#'or not OR a numeric threshold below which all values are set to zero (by
#'default is set to 'FALSE').
#'@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 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, method='QUANT',
Verónica Torralba-Fernández
committed
DynBiasCorrection<- function(exp, obs, method = 'QUANT',wetday=FALSE,
proxy = "dim", quanti, 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 (!(any(names(dim(exp)) %in% 'time'))){
if (any(names(dim(exp)) %in% 'sdate')) {
if (any(names(dim(exp)) %in% 'ftime')) {
exp <- MergeDims(exp, merge_dims = c('ftime', 'sdate'),
rename_dim = 'time')
}
}
}
if (!(any(names(dim(obs)) %in% 'time'))){
if (any(names(dim(obs)) %in% 'sdate')) {
if (any(names(dim(obs)) %in% 'ftime')) {
obs <- MergeDims(obs, merge_dims = c('ftime', 'sdate'),
rename_dim = 'time')
}
}
}
dim_exp <- dim(exp)
names_to_check <- names(dim_exp)[which(names(dim_exp) %in%
c('time', 'lat', 'lon', 'sdate') == FALSE)]
if (length(names_to_check) > 0) {
dim_obs <- dim(obs)
if (any(names(dim_obs) %in% names_to_check)) {
if (any(dim_obs[which(names(dim_obs) %in% names_to_check)] !=
dim_exp[which(names(dim_exp) %in% names_to_check)])) {
for (i in names_to_check) {
pos <- which(names(dim_obs) == i)
names(dim(obs))[pos] <- ifelse(dim_obs[pos] !=
dim_exp[which(names(dim_exp) == i)],
paste0('obs_', names(dim_obs[pos])),
names(dim(obs)[pos]))
}
warning("Common dimension names with different length are renamed.")
}
}
}
adjusted <- Apply(list(exp, obs), target_dims = 'time',
Verónica Torralba-Fernández
committed
fun = .dynbias, method, wetday,
predyn.exp = predyn.exp$pred.dim$pos.d,
predyn.obs = predyn.obs$pred.dim$pos.d,
ncores = ncores, output_dims = 'time')$output1
adjusted <- Apply(list(exp, obs), target_dims = 'time',
Verónica Torralba-Fernández
committed
fun = .dynbias, method, wetday,
predyn.exp = predyn.exp$pred.theta$pos.t,
predyn.obs = predyn.obs$pred.theta$pos.t,
ncores = ncores, output_dims = 'time')$output1
} else {
stop ("Parameter 'proxy' must be set as 'dim' or 'theta'.")
}
if(any(names(dim(adjusted)) %in% 'memberObs')){
if(dim(adjusted)['memberObs'] == 1){
adjusted <- Subset(adjusted,along='memberObs',indices=1,drop = 'selected')
}else{
print('Dimension member in obs changed to memberObs')
}
}
if(any(names(dim(adjusted)) %in% 'datasetObs')){
if(dim(adjusted)['datasetObs'] == 1){
adjusted <- Subset(adjusted,along='datasetObs',indices=1,drop = 'selected')
}else{
print('Dimension dataset in obs changed to datasetObs')
}
}
Verónica Torralba-Fernández
committed
.dynbias <- function(exp, obs, method, wetday,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]]]
Verónica Torralba-Fernández
committed
adjust <- .qbiascorrection(exp_sub, obs_sub, method,wetday)
result[predyn.exp[[x]]] <<- adjust
return(NULL)
})
return(result)
}
Verónica Torralba-Fernández
committed
.qbiascorrection <- function(expX, obsX, method,wetday) {
## functions fitQmap and doQmap
if (method == "PTF") {
qm.fit <- fitQmap(obsX, expX, method = "PTF", transfun = "expasympt",
Verónica Torralba-Fernández
committed
cost = "RSS", wet.day = wetday)
qmap <- doQmap(expX, qm.fit)
} else if (method == "QUANT") {
Verónica Torralba-Fernández
committed
qm.fit <- fitQmap(obsX, expX, method = "QUANT", qstep = 0.01, wet.day = wetday)
qmap <- doQmap(expX, qm.fit, type = "tricub")
} else if (method == "RQUANT") {
Verónica Torralba-Fernández
committed
qm.fit <- fitQmap(obsX, expX, method = "RQUANT", qstep = 0.01,wet.day = wetday)
qmap <- doQmap(expX, qm.fit, type = "linear")
} else if (method == "SSPLIN") {
Verónica Torralba-Fernández
committed
qm.fit <- fitQmap(obsX, expX, qstep = 0.01, method = "SSPLIN",wet.day = wetday)