From 5127d274e9ad3a044ba1e163ad5bfae0d329dcef Mon Sep 17 00:00:00 2001 From: lpalma Date: Tue, 2 Mar 2021 10:59:02 +0100 Subject: [PATCH 1/5] 1st working version --- R/CST_BiasCorrection.R | 41 +++++++++++++++++++++++++++++------------ 1 file changed, 29 insertions(+), 12 deletions(-) diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index 9baf897b..bda531ba 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -90,20 +90,37 @@ BiasCorrection <- function (exp, obs , na.rm = FALSE) { 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) { + # defining forecast,hindcast and observation in cross-validation + fcst <- var_exp[ , t] + hcst <- var_exp[ , -t] + obs <- var_obs[-t] + + # 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) + + # bias corrected forecast + corrected[ , t] <- ((fcst - clim_exp) * (sd_obs / sd_exp)) + clim_obs + names(dim(corrected)) <- c('member', 'sdate') + } + } else { + # defning forecast,hindcast and observation in cross-validation + fcst <- var_cor + hcst <- var_exp + obs <- var_obs # parameters sd_obs <- sd(obs , na.rm = na.rm) @@ -112,8 +129,8 @@ BiasCorrection <- function (exp, obs , na.rm = FALSE) { clim_obs <- mean(obs , na.rm = na.rm) # bias corrected forecast - corrected[ , t] <- ((fcst - clim_exp) * (sd_obs / sd_exp)) + clim_obs + corrected <- ((fcst - clim_exp) * (sd_obs / sd_exp)) + clim_obs + names(dim(corrected)) <- c('member') } - names(dim(corrected)) <- c('member', 'sdate') return(corrected) } -- GitLab From e20b6e62b7c26a7dd6218838871ba83445f4aff7 Mon Sep 17 00:00:00 2001 From: lpalma Date: Tue, 30 Mar 2021 19:47:59 +0200 Subject: [PATCH 2/5] Removed redundant vars --- R/CST_BiasCorrection.R | 30 ++++++++++-------------------- 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index bda531ba..45db9d2c 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -101,35 +101,25 @@ BiasCorrection <- function (exp, obs , na.rm = FALSE) { if (is.null(var_cor)){ 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] - # 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[-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] <- ((fcst - clim_exp) * (sd_obs / sd_exp)) + clim_obs + corrected[ , t] <- ((var_exp[ , t] - clim_exp) * (sd_obs / sd_exp)) + clim_obs names(dim(corrected)) <- c('member', 'sdate') } } else { - # defning forecast,hindcast and observation in cross-validation - fcst <- var_cor - hcst <- var_exp - obs <- var_obs - # 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 <- ((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') } return(corrected) -- GitLab From 9a3029d13efe3bfa979773fad6a2df0adb072832 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 31 Mar 2021 17:17:52 +0200 Subject: [PATCH 3/5] FInclude exp_cor to be called by Apply and CST_ functions --- .Rbuildignore | 2 +- DESCRIPTION | 3 +- NAMESPACE | 1 + NEWS.md | 9 +++- R/CST_BiasCorrection.R | 87 +++++++++++++++++++++++++++++++-------- man/BiasCorrection.Rd | 41 ++++++++++++++++++ man/CST_BiasCorrection.Rd | 4 +- 7 files changed, 126 insertions(+), 21 deletions(-) create mode 100644 man/BiasCorrection.Rd diff --git a/.Rbuildignore b/.Rbuildignore index fa596e70..b2d8e5fc 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,4 +5,4 @@ ./.nc$ .*^(?!data)\.RData$ .*\.gitlab-ci.yml$ -^tests$ +#^tests$ diff --git a/DESCRIPTION b/DESCRIPTION index 7fe8d681..aecd89b3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: CSTools Title: Assessing Skill of Climate Forecasts on Seasonal-to-Decadal Timescales -Version: 4.0.0 +Version: 4.0.1 Authors@R: c( person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8568-3071")), person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "aut", comment = c(ORCID = "0000-0001-5221-0147")), @@ -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 7e313bc5..17befef4 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 23d19889..dd228936 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,12 @@ +### CSTools 4.0.1 +**Submission date to CRAN: 23-02-2021** + +- New features: + + CST_BiasCorrection and BiasCorrection allows to calibrate a forecast given the calibration in the hindcast by using parameter 'exp_cor'. + + ### 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 45db9d2c..b7bd3e1d 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,14 +126,24 @@ 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 , na.rm = na.rm)$output1 + } return(BiasCorrected) } -.sbc <- function(var_obs, var_exp ,var_cor=NULL, 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'))) { diff --git a/man/BiasCorrection.Rd b/man/BiasCorrection.Rd new file mode 100644 index 00000000..8b60a6de --- /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 55c325a2..adfc2798 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{ -- GitLab From 3c369d8440c3b4ad5399a8829df9440df5b122ad Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 14 May 2021 20:51:06 +0200 Subject: [PATCH 4/5] add output_dims to get correct dimnames --- R/CST_BiasCorrection.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index b7bd3e1d..0b837513 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -138,7 +138,8 @@ BiasCorrection <- function (exp, obs, exp_cor = NULL, na.rm = FALSE) { target_dims = list(target_dims_obs, c('member', 'sdate'), c('member', 'sdate')), - fun = .sbc , na.rm = na.rm)$output1 + fun = .sbc , output_dims = c('member','sdate'), + na.rm = na.rm)$output1 } return(BiasCorrected) } -- GitLab From 50a4fe8547554ba40c976ca0add3d5667f79c94c Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 25 May 2021 16:19:45 +0200 Subject: [PATCH 5/5] one automatic test to biascorrection --- tests/testthat/test-CST_BiasCorrection.R | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/testthat/test-CST_BiasCorrection.R b/tests/testthat/test-CST_BiasCorrection.R index c5525c57..fc6405ce 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))) }) -- GitLab