diff --git a/.Rbuildignore b/.Rbuildignore index fa596e707601da63df8c53cf4f087a70a953dbea..b2d8e5fcebca62bff5e5380a881580283874cd54 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,4 +5,4 @@ ./.nc$ .*^(?!data)\.RData$ .*\.gitlab-ci.yml$ -^tests$ +#^tests$ diff --git a/DESCRIPTION b/DESCRIPTION index 6bc5898fbaf0cccef6786eefbae6da7100b6096d..aecd89b3bf1497c6f1c057cb9cfb7d9392cc41f2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,6 +10,7 @@ Authors@R: c( person("Jost", "von Hardenberg", , email = c("j.vonhardenberg@isac.cnr.it", "jost.hardenberg@polito.it"), role = "aut", comment = c(ORCID = "0000-0002-5312-8070")), person("Llorenç", "LLedo", , "llledo@bsc.es", role = "aut"), person("Nicolau", "Manubens", , "nicolau.manubens@bsc.es", role = "aut"), + person("Lluís", "Palma", , "lluis.palma@bsc.es", role = "aut"), person("Eroteida", "Sanchez-Garcia", , "esanchezg@aemet.es", role = "aut"), person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "aut"), person("Veronica", "Torralba", , "veronica.torralba@bsc.es", role = "aut"), diff --git a/NAMESPACE b/NAMESPACE index 7e313bc5e0cdbdfbd3f2735c6958745477e20276..17befef4842581b9c5b473ea802b99fc1fe089e1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ export(AdamontAnalog) export(Analogs) export(BEI_PDFBest) export(BEI_Weights) +export(BiasCorrection) export(CST_Analogs) export(CST_AnalogsPredictors) export(CST_Anomaly) diff --git a/NEWS.md b/NEWS.md index 60d119cf3492af2c342d10bc03e7bef9909139a2..ae7955dd26f648927f46437124eddcebf62384f6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,13 +1,17 @@ ### CSTools 4.0.1 +<<<<<<< HEAD **Submission date to CRAN: XX-06-2021** +- New features: + + CST_BiasCorrection and BiasCorrection allows to calibrate a forecast given the calibration in the hindcast by using parameter 'exp_cor'. + - Fixes: + Calibration retains correlation absolute value + Calibration fixed when cal.methodi == rpc-based, apply_to == sign, eval.method == 'leave-one-out' and the correlation is not significant ### CSTools 4.0.0 -**Submission date to CRAN: XX-12-2020** +**Submission date to CRAN: 23-02-2021** - New features: + ADAMONT downscaling method: requires CST_AdamontAnalogs and CST_AdamontQQCor functions diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index 9baf897bb680eff59d5b3bd4b3169e4a880ee65b..0b8375130cb170345fd7e356214f577b3703662b 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -5,6 +5,7 @@ #' #'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data} #'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}. +#'@param exp_cor an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonl forecast experiment to be corrected. If it is NULL, the 'exp' forecast will be corrected. #'@param na.rm a logical value indicating whether missing values should be stripped before the computation proceeds, by default it is set to FALSE. #' #'@return an object of class \code{s2dv_cube} containing the bias corrected forecasts in the element called \code{$data} with the same dimensions of the experimental data. @@ -30,28 +31,70 @@ #'a <- CST_BiasCorrection(exp = exp, obs = obs) #'str(a) #'@export -CST_BiasCorrection <- function(exp, obs, na.rm = FALSE) { +CST_BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE) { if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") } - if (dim(obs$data)['member'] != 1) { stop("The length of the dimension 'member' in the component 'data' ", "of the parameter 'obs' must be equal to 1.") } - dimnames <- names(dim(exp$data)) - BiasCorrected <- BiasCorrection(exp = exp$data, obs = obs$data, na.rm = na.rm) - pos <- match(dimnames, names(dim(BiasCorrected))) - BiasCorrected <- aperm(BiasCorrected, pos) - names(dim(BiasCorrected)) <- dimnames - exp$data <- BiasCorrected - exp$Datasets <- c(exp$Datasets, obs$Datasets) - exp$source_files <- c(exp$source_files, obs$source_files) - return(exp) + if (is.null(exp_cor)) { + dimnames <- names(dim(exp$data)) + BiasCorrected <- BiasCorrection(exp = exp$data, obs = obs$data, na.rm = na.rm) + pos <- match(dimnames, names(dim(BiasCorrected))) + BiasCorrected <- aperm(BiasCorrected, pos) + names(dim(BiasCorrected)) <- dimnames + exp$data <- BiasCorrected + exp$Datasets <- c(exp$Datasets, obs$Datasets) + exp$source_files <- c(exp$source_files, obs$source_files) + return(exp) + } else { + if (!inherits(exp_cor, 's2dv_cube')) { + stop("Parameter 'var_cor' must be of the class 's2dv_cube'.") + } + dimnames <- names(dim(exp_cor$data)) + BiasCorrected <- BiasCorrection(exp = exp$data, obs = obs$data, + exp_cor = exp_cor$data, na.rm = na.rm) + pos <- match(dimnames, names(dim(BiasCorrected))) + BiasCorrected <- aperm(BiasCorrected, pos) + names(dim(BiasCorrected)) <- dimnames + exp_cor$data <- BiasCorrected + exp_cor$Datasets <- c(exp_cor$Datasets, exp$Datasets, obs$Datasets) + exp_cor$source_files <- c(exp_cor$source_files, + exp$source_files, obs$source_files) + return(exp_cor) + } } - -BiasCorrection <- function (exp, obs , na.rm = FALSE) { +#' Bias Correction based on the mean and standard deviation adjustment +#' +#'@author Verónica Torralba, \email{veronica.torralba@bsc.es} +#'@description This function applies the simple bias adjustment technique described in Torralba et al. (2017). The adjusted forecasts have an equivalent standard deviation and mean to that of the reference dataset. +#' +#'@param exp a multidimensional array with named dimensions containing the seasonal forecast experiment data with at least 'member' and 'sdate' dimensions. +#'@param obs a multidimensional array with named dimensions containing the observed data with at least 'sdate' dimension. +#'@param exp_cor a multidimensional array with named dimensions containing the seasonl forecast experiment to be corrected. If it is NULL, the 'exp' forecast will be corrected. +#'@param na.rm a logical value indicating whether missing values should be stripped before the computation proceeds, by default it is set to FALSE. +#' +#'@return an object of class \code{s2dv_cube} containing the bias corrected forecasts in the element called \code{$data} with the same dimensions of the experimental data. +#' +#'@references Torralba, V., F.J. Doblas-Reyes, D. MacLeod, I. Christel and M. Davis (2017). Seasonal climate prediction: a new source of information for the management of wind energy resources. Journal of Applied Meteorology and Climatology, 56, 1231-1247, doi:10.1175/JAMC-D-16-0204.1. (CLIM4ENERGY, EUPORIAS, NEWA, RESILIENCE, SPECS) +#' +#'@import multiApply +#'@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) +#'a <- BiasCorrection(exp = mod1, obs = obs1) +#'str(a) +#'@export +BiasCorrection <- function (exp, obs, exp_cor = NULL, na.rm = FALSE) { if (!all(c('member', 'sdate') %in% names(dim(exp)))) { stop("Parameter 'exp' must have the dimensions 'member' and 'sdate'.") @@ -83,37 +126,55 @@ BiasCorrection <- function (exp, obs , na.rm = FALSE) { if ('member' %in% names(dim(obs))) { target_dims_obs <- c('member', target_dims_obs) } - - BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp), - target_dims = list(target_dims_obs, c('member', 'sdate')), + if (is.null(exp_cor)) { + BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp), + target_dims = list(target_dims_obs, + c('member', 'sdate')), fun = .sbc , na.rm = na.rm)$output1 + } else { + BiasCorrected <- Apply(data = list(var_obs = obs, + var_exp = exp, + var_cor = exp_cor), + target_dims = list(target_dims_obs, + c('member', 'sdate'), + c('member', 'sdate')), + fun = .sbc , output_dims = c('member','sdate'), + na.rm = na.rm)$output1 + } return(BiasCorrected) } -.sbc <- function(var_obs, var_exp , na.rm = FALSE) { +.sbc <- function(var_obs, var_exp , var_cor = NULL, na.rm = FALSE) { nmembers <- dim(var_exp)['member'][] ntime <- dim(var_exp)['sdate'][] - if (all(names(dim(var_exp)) != c('member','sdate'))) { - var_exp <- t(var_exp) - } + #if (all(names(dim(var_exp)) != c('member','sdate'))) { + # var_exp <- t(var_exp) + #} corrected <- 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] - + + if (is.null(var_cor)){ + for (t in 1 : ntime) { + # parameters + sd_obs <- sd(var_obs[-t], na.rm = na.rm) + sd_exp <- sd(var_exp[ , -t], na.rm = na.rm) + clim_exp <- mean(var_exp[ , -t], na.rm = na.rm) + clim_obs <- mean(var_obs[-t], na.rm = na.rm) + + # bias corrected forecast + corrected[ , t] <- ((var_exp[ , t] - clim_exp) * (sd_obs / sd_exp)) + clim_obs + names(dim(corrected)) <- c('member', 'sdate') + } + } else { # parameters - sd_obs <- sd(obs , na.rm = na.rm) - sd_exp <- sd(hcst , na.rm = na.rm) - clim_exp <- mean(hcst , na.rm = na.rm) - clim_obs <- mean(obs , na.rm = na.rm) + sd_obs <- sd(var_obs, na.rm = na.rm) + sd_exp <- sd(var_exp, na.rm = na.rm) + clim_exp <- mean(var_exp, na.rm = na.rm) + clim_obs <- mean(var_obs, na.rm = na.rm) # bias corrected forecast - corrected[ , t] <- ((fcst - clim_exp) * (sd_obs / sd_exp)) + clim_obs + corrected <- ((var_cor - clim_exp) * (sd_obs / sd_exp)) + clim_obs + names(dim(corrected)) <- c('member') } - names(dim(corrected)) <- c('member', 'sdate') return(corrected) } diff --git a/man/BiasCorrection.Rd b/man/BiasCorrection.Rd new file mode 100644 index 0000000000000000000000000000000000000000..8b60a6de655fd4261a51bf72ad02f99ac2756498 --- /dev/null +++ b/man/BiasCorrection.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_BiasCorrection.R +\name{BiasCorrection} +\alias{BiasCorrection} +\title{Bias Correction based on the mean and standard deviation adjustment} +\usage{ +BiasCorrection(exp, obs, exp_cor = NULL, na.rm = FALSE) +} +\arguments{ +\item{exp}{a multidimensional array with named dimensions containing the seasonal forecast experiment data with at least 'member' and 'sdate' dimensions.} + +\item{obs}{a multidimensional array with named dimensions containing the observed data with at least 'sdate' dimension.} + +\item{exp_cor}{a multidimensional array with named dimensions containing the seasonl forecast experiment to be corrected. If it is NULL, the 'exp' forecast will be corrected.} + +\item{na.rm}{a logical value indicating whether missing values should be stripped before the computation proceeds, by default it is set to FALSE.} +} +\value{ +an object of class \code{s2dv_cube} containing the bias corrected forecasts in the element called \code{$data} with the same dimensions of the experimental data. +} +\description{ +This function applies the simple bias adjustment technique described in Torralba et al. (2017). The adjusted forecasts have an equivalent standard deviation and mean to that of the reference dataset. +} +\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) +a <- BiasCorrection(exp = mod1, obs = obs1) +str(a) +} +\references{ +Torralba, V., F.J. Doblas-Reyes, D. MacLeod, I. Christel and M. Davis (2017). Seasonal climate prediction: a new source of information for the management of wind energy resources. Journal of Applied Meteorology and Climatology, 56, 1231-1247, doi:10.1175/JAMC-D-16-0204.1. (CLIM4ENERGY, EUPORIAS, NEWA, RESILIENCE, SPECS) +} +\author{ +Verónica Torralba, \email{veronica.torralba@bsc.es} +} diff --git a/man/CST_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd index 55c325a2db6dc53f8b315c4c0afc718746f79b1a..adfc2798ac9c9ddfa3cb942a0c2faf69233bfaed 100644 --- a/man/CST_BiasCorrection.Rd +++ b/man/CST_BiasCorrection.Rd @@ -4,13 +4,15 @@ \alias{CST_BiasCorrection} \title{Bias Correction based on the mean and standard deviation adjustment} \usage{ -CST_BiasCorrection(exp, obs, na.rm = FALSE) +CST_BiasCorrection(exp, obs, exp_cor = NULL, na.rm = FALSE) } \arguments{ \item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}} \item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}.} +\item{exp_cor}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonl forecast experiment to be corrected. If it is NULL, the 'exp' forecast will be corrected.} + \item{na.rm}{a logical value indicating whether missing values should be stripped before the computation proceeds, by default it is set to FALSE.} } \value{ diff --git a/tests/testthat/test-CST_BiasCorrection.R b/tests/testthat/test-CST_BiasCorrection.R index c5525c57012ec0b92c46e7bc9e2284831a26a66e..fc6405ce7f9ee002c70aa7d27fbb8417bd854e35 100644 --- a/tests/testthat/test-CST_BiasCorrection.R +++ b/tests/testthat/test-CST_BiasCorrection.R @@ -42,4 +42,9 @@ test_that("Sanity checks", { expect_warning(CST_BiasCorrection(exp = exp2, obs = obs2), "Parameter 'obs' contains NA values", "Parameter 'exp' contains NA values.") + + hinc <- array(1:6, c(sdate = 3, member = 2)) + obs <- array(3:6, c(sdate = 3, member = 1)) + res <- round(BiasCorrection(exp = hinc, obs = obs, exp_cor = hinc), 2) + expect_equal(res, array(c(2.66, 4.27, 3.2, 4.8, 3.73, 5.34), c(member = 2, sdate = 3))) })