diff --git a/NAMESPACE b/NAMESPACE index 7bdea4c7b791bb7d6214f22653b95ad94121e895..297790529456371eb825624945ffc27b322c86de 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +export(CST_BiasCorrection) export(CST_Calibration) export(CST_MultiMetric) export(CST_MultivarRMSE) diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R new file mode 100644 index 0000000000000000000000000000000000000000..3dfb18ac8deaae7ea41ee82e754d948ea52da087 --- /dev/null +++ b/R/CST_BiasCorrection.R @@ -0,0 +1,121 @@ +#' 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 data CSTools object (an s2dverification object as output by the \code{Load} function from the s2dverification package). +#' +#'@return a CSTools object (s2dverification object) with the bias corrected forecasts (provided in $mod) with the same dimensions as data$mod. +#' +#'@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 s2dverification +#'@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) +#'lon <- seq(0, 30, 5) +#'lat <- seq(0, 25, 5) +#'data1 <- list(mod = mod1, obs = obs1, lat = lat, lon = lon) +#'a <- CST_BiasCorrection(data1) +#'str(a) +#'@export +CST_BiasCorrection <- 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 as output of Load function + from s2dverification package") + } + + if (dim(data$obs)['member']!=1) { + stop("The length of the dimension 'member' in label 'obs' of parameter 'data' must be equal to 1.") + } + + BiasCorrected <- BiasCorrection(data$mod, data$obs) + dimnames <- names(dim(BiasCorrected)) + BiasCorrected <- aperm(BiasCorrected, c(3, 1, 2, 4, 5, 6)) + names(dim(BiasCorrected)) <- dimnames[c(3, 1, 2, 4, 5, 6)] + data$biascorrection <- BiasCorrected + data$obs <- NULL + data$mod <- NULL + data <- data[c("biascorrection", names(data)[-which(names(data) == "biascorrection")])] + return(data) +} + +BiasCorrection <- 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 (any(is.na(exp))) { + warning("Parameter 'exp' contains NA values.") + } + + if (any(is.na(obs))) { + warning("Parameter 'obs' contains NA values.") + } + + target_dims_obs <- 'sdate' + 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')), + fun = .sbc)$output1 + return(BiasCorrected) +} + +.sbc <- function(var_obs, var_exp) { + nmembers <- dim(var_exp)['member'][] + ntime <- dim(var_exp)['sdate'][] + 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] + + # parameters + sd_obs <- sd(obs) + sd_exp <- sd(hcst) + clim_exp <- mean(hcst) + clim_obs <- mean(obs) + + # bias corrected forecast + corrected[ , t] <- ((fcst - clim_exp) * (sd_obs / sd_exp)) + clim_obs + } + names(dim(corrected)) <- c('member', 'sdate') + return(corrected) +} diff --git a/man/CST_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd new file mode 100644 index 0000000000000000000000000000000000000000..d187ba1243517e25c90898ea49707d24ec47918b --- /dev/null +++ b/man/CST_BiasCorrection.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_BiasCorrection.R +\name{CST_BiasCorrection} +\alias{CST_BiasCorrection} +\title{Bias Correction based on the mean and standard deviation adjustment} +\usage{ +CST_BiasCorrection(data) +} +\arguments{ +\item{data}{CSTools object (an s2dverification object as output by the \code{Load} function from the s2dverification package).} +} +\value{ +a CSTools object (s2dverification object) with the bias corrected forecasts (provided in $mod) with the same dimensions as data$mod. +} +\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) +lon <- seq(0, 30, 5) +lat <- seq(0, 25, 5) +data1 <- list(mod = mod1, obs = obs1, lat = lat, lon = lon) +a <- CST_BiasCorrection(data1) +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/tests/testthat/test-CST_BiasCorrection.R b/tests/testthat/test-CST_BiasCorrection.R new file mode 100644 index 0000000000000000000000000000000000000000..acc7425d51c29fa2958145400f223e6520ff4c18 --- /dev/null +++ b/tests/testthat/test-CST_BiasCorrection.R @@ -0,0 +1,44 @@ +context("Generic tests") +test_that("Sanity checks", { + expect_error( + CST_BiasCorrection(data = 1), + "Parameter 'data' must be a list as output of Load function + from s2dverification package") + + mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) + obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7) + dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 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) + data <- list(mod = mod1, obs = obs1, lat = lat, lon = lon) + + expect_equal(length(CST_BiasCorrection(data = data)), 3) + expect_equal(dim(CST_BiasCorrection(data = data)$biascorrection), + c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)) + expect_equal(CST_BiasCorrection(data = data)$lat, lat) + expect_equal(CST_BiasCorrection(data = data)$lon, lon) + + data <- list(mod = mod1, obs = mod1, lat = lat, lon = lon) + expect_error(CST_BiasCorrection(data = data), + "The length of the dimension 'member' in label 'obs' of parameter 'data' must be equal to 1.") + + mod2 <- mod1 + mod2[1, 2, 1, 1, 1, 1] <- NA + data <- list(mod = mod2, obs = obs1, lat = lat, lon = lon) + expect_warning( + CST_BiasCorrection(data = data), + "Parameter 'exp' contains NA values.") + + obs2 <- obs1 + obs2[1, 1, 2, 1, 1, 1] <- NA + data <- list(mod = mod1, obs = obs2, lat = lat, lon = lon) + expect_warning( + CST_BiasCorrection(data = data), + "Parameter 'obs' contains NA values.") + + data <- list( mod = mod2, obs = obs2, lat = lat, lon = lon) + expect_warning( + CST_BiasCorrection(data = data), + "Parameter 'obs' contains NA values", "Parameter 'exp' contains NA values.") +}) \ No newline at end of file