From c67314a17ed1223e2a0219bd8d3e03f84b4750e0 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 19 Jul 2022 18:13:17 +0200 Subject: [PATCH 01/14] cross-validation --- R/CST_QuantileMapping.R | 475 ++++++++-------------------------------- 1 file changed, 93 insertions(+), 382 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index b27cbd4c..0a1c5d3b 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -1,248 +1,81 @@ -#'Quantiles Mapping for seasonal or decadal forecast data +#'Quaintiles Mapping for seasonal or decadal forecast data #' -#'@description This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. The quantile mapping adjustment between an experiment, tipically a hindcast, and observations is applied to the experiment itself or to a provided forecast. +#'@description This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. #' #'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} #'@param exp an object of class \code{s2dv_cube} #'@param obs an object of class \code{s2dv_cube} -#'@param exp_cor an object of class \code{s2dv_cube} in which the quantile mapping correction will be applied. If it is not specified, the correction is applied in object \code{exp}. -#'@param sample_dims a character vector indicating the dimensions that can be used as sample for the same distribution -#'@param sample_length a numeric value indicating the length of the timeseries window to be used as sample for the sample distribution and correction. By default, NULL, the total length of the timeseries will be used. +#'@parma exp_cor an object of class \code{s2dv_cube} in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'. +#'@param sampledims a character vector indicating the dimensions that can be used as sample for the probabilit #'@param method a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used. -#'@param ncores an integer indicating the number of parallel processes to spawn for the use for parallel computation in multiple cores. -#'@param ... additional arguments passed to the method specified by \code{method}. +#'@param wetday logical indicating whether to perform wet day correction or not. OR a numeric +#'@param qstep a numeric value between 0 and 1. The quantile mapping is fitted only for the quantiles defined by quantile(0,1,probs=seq(0,1,by=qstep). +#'@param type type of interpolation between the fitted transformed values. See details. #' -#'@details The different methods are: -#'\itemize{ -#'\item{'PTF'} {fits a parametric transformations to the quantile-quantile relation of observed and modelled values. See \code{?qmap::fitQmapPTF}.} -#' \item{'DIST'} {fits a theoretical distribution to observed and to modelled time series. See \code{?qmap::fitQmapDIST}.} -#'\item{'RQUANT'} {estimates the values of the quantile-quantile relation of observed and modelled time series for regularly spaced quantiles using local linear least square regression. See \code{?qmap::fitQmapRQUANT}.} -#'\item{'QUANT'} {estimates values of the empirical cumulative distribution function of observed and modelled time series for regularly spaced quantiles. See \code{?qmap::fitQmapQUANT}.} -#'\item{'SSPLIN'} {fits a smoothing spline to the quantile-quantile plot of observed and modelled time series. See \code{?qmap::fitQmapSSPLIN}.}} -#'All methods accepts some common arguments: -#'\itemize{ -#'\item{wet.day} {logical indicating whether to perform wet day correction or not.(Not available in 'DIS' method)} -#'\item{qstep} {NULL or a numeric value between 0 and 1.}} -#' When providing a forecast to be corrected through the pararmeter \code{exp_cor}, some inputs might need to be modified. The quantile correction is compute by comparing objects passed through 'exp' and 'obs' parameters, this correction will be later applied to the forecast provided in 'exp_cor'. Imaging the case of 'exp' and 'obs' having several start dates, stored using a dimension e.g. 'sdate', 'sample_dims' include this dimension 'sdate' and 'exp_cor' has forecasts for several sdates but different from the ones in 'exp'. In this case, the correction computed with 'exp' and 'obs' would be applied for each 'sdate' of 'exp_cor' separately. This example corresponds to a case of split a dataset in training set and validation set. -#' #'@return an oject of class \code{s2dv_cube} containing the experimental data after applyingthe quantile mapping correction. #') <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , #'@import qmap #'@import multiApply -#'@import abind #' -#'@seealso \code{qmap::fitQmap} and \code{qmap::doQmap} +#'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} #'@examples -#'library(qmap) -#'exp <- 1 : (1 * 5 * 10 * 6 * 2 * 3) -#'dim(exp) <- c(dataset = 1, member = 10, sdate = 5, ftime = 6 , -#' lat = 2, lon = 3) -#'exp <- list(data = exp) +#'exp$data <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) +#'dim(exp$data) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , +#' lat = 6, lon = 7) #'class(exp) <- 's2dv_cube' -#'obs <- 101 : (100 + 1 * 1 * 5 * 6 * 2 * 3) -#'dim(obs) <- c(dataset = 1, member = 1, sdate = 5, ftime = 6 , -#' lat = 2, lon = 3) -#'obs <- list(data = obs) +#'obs$data <- 101 : c(100 + 1 * 1 * 20 * 60 * 6 * 7) +#'dim(obs$data) <- c(dataset = 1, member = 1, sdate = 20, ftime = 60 , +#' lat = 6, lon = 7) #'class(obs) <- 's2dv_cube' -#'res <- CST_QuantileMapping(exp, obs, method = 'RQUANT') -#'\donttest{ -#'exp <- lonlat_data$exp -#'obs <- lonlat_data$obs #'res <- CST_QuantileMapping(exp, obs) -#' -#'exp_cor <- exp -#'exp_cor$data <- exp_cor$data[,,1,,,] -#'dim(exp_cor$data) <- c(dataset = 1, member = 15, sdate = 1, ftime = 3, -#' lat = 22, lon = 53) -#'res <- CST_QuantileMapping(exp, obs, exp_cor, -#' sample_dims = c('sdate', 'ftime', 'member')) -#'res <- CST_QuantileMapping(exp, obs, exp_cor, -#' sample_dims = c('ftime', 'member')) -#'data(obsprecip) -#'data(modprecip) -#'exp <- modprecip$MOSS[1:10000] -#'dim(exp) <- c(time = length(exp)) -#'exp <- list(data = exp) -#'class(exp) <- 's2dv_cube' -#'obs <- obsprecip$MOSS[1:10000] -#'dim(obs) <- c(time = length(obs)) -#'obs <- list(data = obs) -#'class(obs) <- 's2dv_cube' -#'res <- CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time', -#' method = 'DIST') -#'# Example using different lenght of members and sdates: #'exp <- lonlat_data$exp -#'exp$data <- exp$data[,,1:4,,,] -#'dim(exp$data) <- c(dataset = 1, member = 15, sdate = 4, ftime = 3, -#' lat = 22, lon = 53) #'obs <- lonlat_data$obs -#'obs$data <- obs$data[,,1:4, ,,] -#'dim(obs$data) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3, -#' lat = 22, lon = 53) -#'exp_cor <- lonlat_data$exp -#'exp_cor$data <- exp_cor$data[,1:5,5:6,,,] -#'dim(exp_cor$data) <- c(dataset = 1, member = 5, sdate = 2, ftime = 3, -#' lat = 22, lon = 53) -#'res <- CST_QuantileMapping(exp, obs, exp_cor, -#' sample_dims = c('sdate', 'ftime', 'member')) -#'exp_cor <- lonlat_data$exp -#'exp_cor$data <- exp_cor$data[,,5:6,,,] -#'dim(exp_cor$data) <- c(dataset = 1, member = 15, sdate = 2, ftime = 3, -#' lat = 22, lon = 53) -#'res <- CST_QuantileMapping(exp, obs, exp_cor, -#' sample_dims = c('sdate', 'ftime', 'member')) -#'} +#'res <- CST_QuantileMapping(exp, obs) #'@export -CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, - sample_dims = c('sdate', 'ftime', 'member'), - sample_length = NULL, method = 'QUANT', ncores = NULL, ...) { - if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { +CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', memb_dim = 'member', + window_dim = NULL, + method = 'QUANT', na.rm = FALSE, ncores = NULL, ...) { + if (!inherits(exp, 's2dv_cube') || !inherits(exp, 's2dv_cube')) { stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") } - if (!is.null(exp_cor)) { - if (!inherits(exp_cor, 's2dv_cube')) { - stop("Parameter 'exp_cor' must be of the class 's2dv_cube', ", - "as output by CSTools::CST_Load.") - } - } - if (!(method %in% c('PTF','DIST','RQUANT','QUANT','SSPLIN'))) { - stop("Parameter 'method' must be one of the following methods: ", - "'PTF','DIST','RQUANT','QUANT','SSPLIN'.") - } + dimnames <- names(dim(exp$data)) QMapped <- QuantileMapping(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, - sample_dims = sample_dims, sample_length = sample_length, - method = method, ncores = ncores, ...) - if (is.null(exp_cor)) { - exp$data <- QMapped - exp$source_files <- c(exp$source_files, obs$source_files) - } else { - exp_cor$data <- QMapped - exp_cor$source_files <- c(exp$source_files, obs$source_files, exp_cor$source_files) - exp <- exp_cor - } + sdate_dim = sdate_dim, memb_dim = memb_dim, window_dim = window_dim, + method = method, na.rm = na.rm, ncores = ncores, ...) + exp$data <- QMapped + exp$Datasets <- c(exp$Datasets, obs$Datasets) + exp$source_files <- c(exp$source_files, obs$source_files) return(exp) } -#'Quantiles Mapping for seasonal or decadal forecast data -#' -#'@description This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. The quantile mapping adjustment between an experiment, tipically a hindcast, and observations is applied to the experiment itself or to a provided forecast. -#' -#'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} -#'@param exp a multi-dimensional array with named dimensions containing the hindcast. -#'@param obs a multi-dimensional array with named dimensions (the same as the provided in 'exp') containing the reference dataset. -#'@param exp_cor a multi-dimensional array with named dimensions in which the quantile mapping correction will be applied. If it is not specified, the correction is applied in object \code{exp}. -#'@param sample_dims a character vector indicating the dimensions that can be used as sample for the same distribution -#'@param sample_length a numeric value indicating the length of the timeseries window to be used as sample for the sample distribution and correction. By default, NULL, the total length of the timeseries will be used. -#'@param method a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used. -#'@param ncores an integer indicating the number of parallel processes to spawn for the use for parallel computation in multiple cores. -#'@param ... additional arguments passed to the method specified by \code{method}. -#' -#'@details The different methods are: -#'\itemize{ -#'\item{'PTF'} {fits a parametric transformations to the quantile-quantile relation of observed and modelled values. See \code{?qmap::fitQmapPTF}.} -#' \item{'DIST'} {fits a theoretical distribution to observed and to modelled time series. See \code{?qmap::fitQmapDIST}.} -#'\item{'RQUANT'} {estimates the values of the quantile-quantile relation of observed and modelled time series for regularly spaced quantiles using local linear least square regression. See \code{?qmap::fitQmapRQUANT}.} -#'\item{'QUANT'} {estimates values of the empirical cumulative distribution function of observed and modelled time series for regularly spaced quantiles. See \code{?qmap::fitQmapQUANT}.} -#'\item{'SSPLIN'} {fits a smoothing spline to the quantile-quantile plot of observed and modelled time series. See \code{?qmap::fitQmapSSPLIN}.}} -#'All methods accepts some common arguments: -#'\itemize{ -#'\item{wet.day} {logical indicating whether to perform wet day correction or not.(Not available in 'DIS' method)} -#'\item{qstep} {NULL or a numeric value between 0 and 1.}} -#'@return an oject of class \code{s2dv_cube} containing the experimental data after applyingthe quantile mapping correction. -#') <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , -#'@import qmap -#'@import multiApply -#'@import abind -#' -#'@seealso \code{qmap::fitQmap} and \code{qmap::doQmap} -#'@export -QuantileMapping <- function(exp, obs, exp_cor = NULL, sample_dims = 'ftime', - sample_length = NULL, method = 'QUANT', ncores = NULL, ...) { +QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', memb_dim = 'member', + window_dim = NULL, method = 'QUANT', na.rm = FALSE, ncores = NULL) { obsdims <- names(dim(obs)) expdims <- names(dim(exp)) if (is.null(expdims)) { - stop("Parameter 'exp' must have dimension names.") + stop("Parameter 'exp' musth have dimension names.") } if (is.null(obsdims)) { - stop("Parameter 'obs' must have dimension names.") - } - if (any(is.na(exp))) { - warning("Parameter 'exp' contains NA values.") - } - if (any(is.na(obs))) { - warning("Parameter 'obs' contains NA values.") + stop("Parameter 'obs' musth have dimension names.") } if (!is.null(exp_cor)) { exp_cordims <- names(dim(exp_cor)) if (is.null(exp_cordims)) { - stop("Parameter 'exp_cor' must have dimension names.") + stop("Parameter 'exp_cor' musth have dimension names.") } } - if (!all(sample_dims %in% expdims)) { - stop("Parameter 'sample_dims' must be a vector of string character ", - "indicating names of exiting dimension in parameter 'exp'.") - } - ## The sample_dims could be of any length (even different between exp and obs) - ## but the repeated dims that aren't in sample_dims should be drop: - commondims <- obsdims[obsdims %in% expdims] - commondims <- names(which(unlist(lapply(commondims, function(x) { - dim(obs)[obsdims == x] != dim(exp)[expdims == x]})))) - if (any(!(commondims %in% sample_dims))) { - todrop <- commondims[!(commondims %in% sample_dims)] - todrop <- match(todrop, obsdims) - if (all(dim(obs)[todrop] != 1)) { - stop("Review parameter 'sample_dims' or the data dimensions", - "since multiple dimensions with different length have ", - "being found in the data inputs that don't match with ", - "'sample_dims' parameter.") - } else { - obs <- adrop(obs, drop = todrop) - } - } - if (!all(sample_dims %in% obsdims)) { - newobsdims <- sample_dims[!sample_dims %in% obsdims] - dim(obs) <- c(dim(obs), 1 : length(newobsdims)) - names(dim(obs))[-c(1:length(obsdims))] <- newobsdims - } - - if (!is.null(exp_cor)) { - commondims <- exp_cordims[exp_cordims %in% expdims] - commondims <- names(which(unlist(lapply(commondims, function(x) { - dim(exp_cor)[exp_cordims == x] != dim(exp)[expdims == x]})))) - if (any(commondims %in% sample_dims)) { - todrop <- commondims[(commondims %in% sample_dims)] - todroppos <- match(todrop, sample_dims) - if (all(dim(exp_cor)[todrop] != 1)) { - warning(paste("The sample_dims", paste(todrop, collapse = " "), - "are not used when applying the", - "correction to 'exp_cor'")) - sample_dims <- list(sample_dims, sample_dims, sample_dims[-todroppos]) - } else { - exp_cor <- adrop(exp_cor, drop = todroppos) - } - } else { - todrop <- commondims[!(commondims %in% sample_dims)] - todrop <- match(todrop, obsdims) - if (all(dim(exp_cor)[todrop] != 1)) { - stop("Review parameter 'sample_dims' or the data dimensions ", - "since multiple dimensions with different length have ", - "being found in the data inputs that don't match with ", - "'sample_dims' parameter.") - } else { - exp_cor <- adrop(exp_cor, drop = todrop) - } + if (!is.null(window_dim)) { + if (!(window_dim %in% obsdims)) { + stop("Dimension 'window_dim' not found in 'obs'.") } - } - - if (!is.null(sample_length) & !is.numeric(sample_length)) { - warning("Parameter 'sample_length' has not been correctly defined and ", - "the whole length of the timeseries will be used.") - sample_length <- NULL + obs <- CSTools::MergeDims(obs, c(memb_dim, window_dim)) } - if (length(sample_length) > 1) { - warning("Parameter 'sample_length' has length > 1 and only the first ", - "element will be used.") - sample_length <- sample_length[1] + sample_dims <- c(memb_dim, sdate_dim) + if (!all(memb_dim %in% obsdims)) { + obs <- InsertDim(obs, posdim = 1, lendim = 1, name = memb_dim[!(memb_dim %in% obsdims)]) + } + if (!all(sample_dims %in% expdims)) { + stop("Parameter 'exp' requires 'sdate_dim' and 'member_dim' dimensions.") } if (!is.character(method)) { warning("Parameter 'method' must be a character string indicating ", @@ -255,183 +88,61 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sample_dims = 'ftime', " is used.") method <- method[1] } - -# qmaped <- Apply(list(exp, obs), target_dims = sample_dims, fun = qmapcor, ..., -# exp_cor = exp_cor, sample_length = sample_length, - if (is.null(exp_cor)) { - qmaped <- Apply(list(exp, obs), target_dims = sample_dims, - fun = qmapcor, ..., sample_length = sample_length, - method = method, ncores = ncores)$output1 - } else { +print(dim(exp)) +print(dim(obs)) +print(sample_dims) + if (!is.null(exp_cor)) { qmaped <- Apply(list(exp, obs, exp_cor), target_dims = sample_dims, - fun = qmapcor, ..., sample_length = sample_length, - method = method, ncores = ncores)$output1 - } - pos <- match(expdims, names(dim(qmaped))) - out_names <- names(dim(exp)) - if (length(pos) < length(dim(qmaped))) { - toadd <- length(dim(qmaped)) - length(pos) - toadd <- seq(max(pos) + 1, max(pos) + toadd, 1) - pos <- c(pos, toadd) - new <- names(dim(qmaped))[names(dim(qmaped)) %in% out_names == FALSE] - out_names <- c(out_names, new) - } - qmaped <- aperm(qmaped, pos) - names(dim(qmaped)) <- out_names + fun = qmapcor, method = method, na.rm = na.rm, + ncores = ncores)$output1 + } else { + qmaped <- Apply(list(exp, obs), target_dims = sample_dims, + fun = qmapcor, exp_cor = NULL, method = method, na.rm = na.rm, + ncores = ncores)$output1 + } return(qmaped) } -qmapcor <- function(exp, obs, exp_cor = NULL, sample_length = NULL, method = 'QUANT', - ...) { - dimnames_exp <- names(dim(exp)) - dimnames_obs <- names(dim(obs)) - if (length(dimnames_exp) != length(dimnames_obs)) { - stop("Parameters 'exp' and 'obs' must have the same number of dimensions.") - } - if (!all(dimnames_exp %in% dimnames_obs)) { - stop("Parameters 'exp' and 'obs' must have the same dimension names.") - } - if (!(any(names(dim(exp)) %in% 'ftime') | any(names(dim(exp)) %in% 'time') | - any(names(dim(exp)) %in% 'sdate'))) { - stop("Parameters 'exp' and 'obs' must have a temporal dimension named ", - "'time', 'ftime' or 'sdate'.") - } - - dimensions <- dim(exp) - if (!is.null(exp_cor)) { - dimensions <- dim(exp_cor) - } - if (any(names(dim(exp)) %in% 'ftime') | any(names(dim(exp)) %in% 'time')) { - time_dim <- which(names(dim(exp)) == 'ftime' | names(dim(exp)) %in% 'time') - } else { - time_dim <- which(names(dim(exp)) == 'sdate') - } - if (any(names(dim(obs)) %in% 'ftime') | any(names(dim(obs)) %in% 'time')) { - time_dim_obs <- which(names(dim(obs)) == 'ftime' | names(dim(obs)) %in% 'time') - } else { - time_dim_obs <- which(names(dim(obs)) == 'sdate') - } - if (is.null(sample_length)) { - sample_length <- dim(exp)[time_dim] - } - nsamples <- dim(exp)[time_dim]/sample_length - if (nsamples %% 1 != 0) { - # add NA to complete the last sample - nsamples <- ceiling(nsamples) - fillsample1D <- rep(NA, nsamples * sample_length - dim(exp)[time_dim]) - if (length(dim(exp)) > 1) { - fillsample <- rep(fillsample1D, prod(dim(exp)[-time_dim])) - dims <- dim(exp) - exp <- c(exp, fillsample) - dim(exp) <- c(dims[-time_dim], sample = sample_length, - ceiling(dims[time_dim]/sample_length)) - fillsample <- rep(fillsample1D, prod(dim(obs)[-time_dim_obs])) - dims <- dim(obs) - obs <- c(obs, fillsample) - dim(obs) <- c(dims[-time_dim_obs], sample = sample_length, - ceiling(dims[time_dim_obs]/sample_length)) - } else { - exp <- abind(exp, fillsample1D, along = time_dim) - names(dim(exp)) <- dimnames_exp - obs <- abind(obs, fillsample1D, along = time_dim_obs) - names(dim(obs)) <- dimnames_obs - dim(exp) <- c(dim(exp)[-time_dim], sample = sample_length, - dim(exp)[time_dim]/sample_length) - dim(obs) <- c(dim(obs)[-time_dim_obs], sample = sample_length, - dim(obs)[time_dim_obs]/sample_length) - } - } else { - dim(exp) <- c(dim(exp)[-time_dim], sample = sample_length, - dim(exp)[time_dim]/sample_length) - dim(obs) <- c(dim(obs)[-time_dim_obs], sample = sample_length, - dim(obs)[time_dim_obs]/sample_length) - } - if (any(names(dim(exp)) %in% 'ftime') | any(names(dim(exp)) %in% 'time')) { - new_time_dim_exp <- which(names(dim(exp)) == 'ftime' | names(dim(exp)) %in% 'time') - } else { - new_time_dim_exp <- which(names(dim(exp)) == 'sdate') - } - if (any(names(dim(obs)) %in% 'ftime') | any(names(dim(obs)) %in% 'time')) { - new_time_dim_obs <- which(names(dim(obs)) == 'ftime' | names(dim(obs)) %in% 'time') - } else { - new_time_dim_obs <- which(names(dim(obs)) == 'sdate') - } - - if (!is.null(exp_cor)) { - if (any(names(dim(exp_cor)) %in% 'ftime') | any(names(dim(exp_cor)) %in% 'time')) { - time_dim_cor <- which(names(dim(exp_cor)) == 'ftime' | names(dim(exp_cor)) %in% 'time') - } else { - time_dim_cor <- which(names(dim(exp_cor)) == 'sdate') - } - - nsamples <- dimensions[time_dim_cor]/sample_length - if (nsamples %% 1 != 0) { - nsamples <- ceiling(nsamples) - fillsample1D <- rep(NA, nsamples * sample_length - dimensions[time_dim_cor]) - if (length(dimensions) > 1) { - fillsample <- rep(fillsample1D, prod(dimensions[-time_dim_cor])) - exp_cor <- c(exp_cor, fillsample) - dim(exp_cor) <- c(dim(exp_cor)[-time_dim_cor], sample = sample_length, - ceiling(dim(exp_cor)[time_dim_cor]/sample_length)) - } else { - exp_cor <- abind(exp_cor, fillsample1D, along = time_dim_cor) - names(dim(exp_cor)) <- names(dimensions) - } - } - dim(exp_cor) <- c(dim(exp_cor)[-time_dim_cor], sample = sample_length, - dim(exp_cor)[time_dim_cor]/sample_length) - if (any(names(dim(exp_cor)) %in% 'ftime') | any(names(dim(exp_cor)) %in% 'time')) { - new_time_dim_cor <- which(names(dim(exp_cor)) == 'ftime' | - names(dim(exp_cor)) %in% 'time') - } else { - new_time_dim_cor <- which(names(dim(exp_cor)) == 'sdate') - } - - } else { - time_dim_cor <- time_dim - exp_cor <- exp - new_time_dim_cor <- new_time_dim_exp - } - applied <- NULL - for (i in 1 : nsamples) { - if (i <= dim(obs)[new_time_dim_obs]) { - sample_obs <- as.vector(asub(obs, idx = i, dims = new_time_dim_obs)) - sample_exp <- as.vector(asub(exp, idx = i, dims = new_time_dim_exp)) - } else { - sample_obs <- as.vector(asub(obs, idx = dim(obs)[new_time_dim_obs], - dims = new_time_dim_obs)) - sample_exp <- as.vector(asub(exp, idx = dim(exp)[new_time_dim_exp], - dims = new_time_dim_exp)) - } - if (i >= dim(obs)[new_time_dim_obs]) { - sample_obs <- sample_obs[!is.na(sample_obs)] - sample_exp <- sample_exp[!is.na(sample_exp)] - } - if (sum(sample_obs) == 0) { - warning("The total sum of observed data sample in the sample number ", - i, ", is zero and the function may crash.") - } - if (sum(sample_exp) == 0) { - warning("The total sum of experimental data sample in the sample number ", - i, ", is zero and the function may crash.") - } - if (length(sample_exp) < sample_length) { - warning("The length of the sample used, ", length(sample_exp), - ", in the sample number ", i, - ", is smaller than the defined in parameter 'sample_length'.") - } - adjust <- fitQmap(sample_obs, sample_exp, method = method, - ...) - sample_cor <- as.vector(asub(exp_cor, idx = i, dims = new_time_dim_cor)) - if (i == nsamples) { - sample_cor <- sample_cor[!is.na(sample_cor)] - } - applied <- c(applied, doQmap(x = sample_cor, fobj = adjust, ...)) - } - if (any(is.na(exp_cor))) { - pos <- which(!is.na(exp_cor)) - exp_cor[pos] <- applied - applied <- exp_cor +qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUAN', na.rm = FALSE, ...) { + if (is.null(exp_cor)) { + applied <- exp * NA + for (sd in 1:dim(exp)['sdate']) { + if (na.rm) { + # select start date for cross-val + nas_pos <- which(!is.na(exp[,sd])) + obs2 <- as.vector(obs[,-sd]) + exp2 <- as.vector(exp[,-sd]) + exp_cor2 <- as.vector(exp[,sd]) + # remove NAs + obs2 <- obs2[!is.na(obs2)] + exp2 <- exp2[!is.na(exp2)] + exp_cor2 <- exp_cor2[!is.na(exp_cor2)] + tryCatch({ + adjust <- fitQmap(obs2, exp2, method = method, ...) + applied[nas_pos, sd] <- doQmap(exp_cor2, adjust, ...) + }, + error = function(error_message) { + return(applied[,sd]) + }) + } else { + adjust <- fitQmap(as.vector(obs[,-sd]), as.vector(exp[,-sd]), method = method, ...) + applied[,sd] <- doQmap(as.vector(exp[,sd]), adjust, ...) + } } - dim(applied) <- dimensions - return(applied) + } else { + applied <- exp_cor * NA + if (na.rm) { + tryCatch({ + adjust <- fitQmap(obs[!is.na(obs)], exp[!is.na(exp)], method = method, ...) + applied[!is.na(exp_cor)] <- doQmap(exp_cor[!is.na(exp_cor)], adjust, ...) + }, + error = function(error_message) { + return(applied) + }) + } else { + adjust <- fitQmap(as.vector(obs), as.vector(exp), method = method, ...) + applied <- doQmap(as.vector(exp_cor), adjust, ...) + } + } + return(applied) } + -- GitLab From f17f2c394cbae410a1f60a420d5be11ef1c73d83 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 20 Jul 2022 14:36:09 +0200 Subject: [PATCH 02/14] formatted code --- R/CST_QuantileMapping.R | 56 +++++++++++++++++++++++++---------------- 1 file changed, 34 insertions(+), 22 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 0a1c5d3b..fa0acc4a 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -6,11 +6,10 @@ #'@param exp an object of class \code{s2dv_cube} #'@param obs an object of class \code{s2dv_cube} #'@parma exp_cor an object of class \code{s2dv_cube} in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'. -#'@param sampledims a character vector indicating the dimensions that can be used as sample for the probabilit +#'@param sdate_dim +#'@param memb_dim +#'@param window_dim #'@param method a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used. -#'@param wetday logical indicating whether to perform wet day correction or not. OR a numeric -#'@param qstep a numeric value between 0 and 1. The quantile mapping is fitted only for the quantiles defined by quantile(0,1,probs=seq(0,1,by=qstep). -#'@param type type of interpolation between the fitted transformed values. See details. #' #'@return an oject of class \code{s2dv_cube} containing the experimental data after applyingthe quantile mapping correction. #') <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , @@ -32,24 +31,29 @@ #'obs <- lonlat_data$obs #'res <- CST_QuantileMapping(exp, obs) #'@export -CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', memb_dim = 'member', - window_dim = NULL, - method = 'QUANT', na.rm = FALSE, ncores = NULL, ...) { +CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', + memb_dim = 'member', window_dim = NULL, + method = 'QUANT', na.rm = FALSE, + ncores = NULL, ...) { if (!inherits(exp, 's2dv_cube') || !inherits(exp, 's2dv_cube')) { stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") } dimnames <- names(dim(exp$data)) - QMapped <- QuantileMapping(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, - sdate_dim = sdate_dim, memb_dim = memb_dim, window_dim = window_dim, - method = method, na.rm = na.rm, ncores = ncores, ...) + QMapped <- QuantileMapping(exp = exp$data, obs = obs$data, + exp_cor = exp_cor$data, + sdate_dim = sdate_dim, memb_dim = memb_dim, + window_dim = window_dim, method = method, + na.rm = na.rm, ncores = ncores, ...) exp$data <- QMapped exp$Datasets <- c(exp$Datasets, obs$Datasets) exp$source_files <- c(exp$source_files, obs$source_files) return(exp) } -QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', memb_dim = 'member', - window_dim = NULL, method = 'QUANT', na.rm = FALSE, ncores = NULL) { +QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', + memb_dim = 'member', window_dim = NULL, + method = 'QUANT', + na.rm = FALSE, ncores = NULL, ...) { obsdims <- names(dim(obs)) expdims <- names(dim(exp)) if (is.null(expdims)) { @@ -69,6 +73,10 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', memb_ stop("Dimension 'window_dim' not found in 'obs'.") } obs <- CSTools::MergeDims(obs, c(memb_dim, window_dim)) + if (window_dim %in% expdims) { + exp <- CSTools::MergeDims(exp, c(memb_dim, window_dim)) + warning("window_dim found in exp and it is merged to memb_dim.") + } } sample_dims <- c(memb_dim, sdate_dim) if (!all(memb_dim %in% obsdims)) { @@ -88,21 +96,22 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', memb_ " is used.") method <- method[1] } -print(dim(exp)) -print(dim(obs)) -print(sample_dims) if (!is.null(exp_cor)) { qmaped <- Apply(list(exp, obs, exp_cor), target_dims = sample_dims, - fun = qmapcor, method = method, na.rm = na.rm, + fun = qmapcor, method = method, na.rm = na.rm, ..., ncores = ncores)$output1 } else { - qmaped <- Apply(list(exp, obs), target_dims = sample_dims, - fun = qmapcor, exp_cor = NULL, method = method, na.rm = na.rm, + qmaped <- Apply(list(exp, obs), target_dims = sample_dims, + fun = qmapcor, exp_cor = NULL, method = method, + na.rm = na.rm, ..., ncores = ncores)$output1 } return(qmaped) } -qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUAN', na.rm = FALSE, ...) { +qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUAN', na.rm = FALSE, + ...) { + # exp[memb, sdate] + # obs[window, sdate] if (is.null(exp_cor)) { applied <- exp * NA for (sd in 1:dim(exp)['sdate']) { @@ -124,7 +133,8 @@ qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUAN', na.rm = FALSE, .. return(applied[,sd]) }) } else { - adjust <- fitQmap(as.vector(obs[,-sd]), as.vector(exp[,-sd]), method = method, ...) + adjust <- fitQmap(as.vector(obs[,-sd]), as.vector(exp[,-sd]), + method = method, ...) applied[,sd] <- doQmap(as.vector(exp[,sd]), adjust, ...) } } @@ -132,8 +142,10 @@ qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUAN', na.rm = FALSE, .. applied <- exp_cor * NA if (na.rm) { tryCatch({ - adjust <- fitQmap(obs[!is.na(obs)], exp[!is.na(exp)], method = method, ...) - applied[!is.na(exp_cor)] <- doQmap(exp_cor[!is.na(exp_cor)], adjust, ...) + adjust <- fitQmap(obs[!is.na(obs)], exp[!is.na(exp)], + method = method, ...) + applied[!is.na(exp_cor)] <- doQmap(exp_cor[!is.na(exp_cor)], + adjust, ...) }, error = function(error_message) { return(applied) -- GitLab From 5d20ae487a6a39b80eebad47fac0ff01a2c30222 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 20 Jul 2022 17:26:50 +0200 Subject: [PATCH 03/14] qmapcor default method fixed --- R/CST_QuantileMapping.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index fa0acc4a..641aa9ce 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -108,7 +108,7 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', } return(qmaped) } -qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUAN', na.rm = FALSE, +qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUANT', na.rm = FALSE, ...) { # exp[memb, sdate] # obs[window, sdate] -- GitLab From 8e11709b20bc97b4bbcb45ce6db959c6275f6c4b Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 22 Jul 2022 13:24:26 +0200 Subject: [PATCH 04/14] order lines for obs dims merging --- R/CST_QuantileMapping.R | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 641aa9ce..099ac7b1 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -15,6 +15,7 @@ #') <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , #'@import qmap #'@import multiApply +#'@import s2dv #' #'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} #'@examples @@ -68,6 +69,10 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', stop("Parameter 'exp_cor' musth have dimension names.") } } + if (!all(memb_dim %in% obsdims)) { + obs <- InsertDim(obs, posdim = 1, lendim = 1, + name = memb_dim[!(memb_dim %in% obsdims)]) + } if (!is.null(window_dim)) { if (!(window_dim %in% obsdims)) { stop("Dimension 'window_dim' not found in 'obs'.") @@ -79,9 +84,6 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', } } sample_dims <- c(memb_dim, sdate_dim) - if (!all(memb_dim %in% obsdims)) { - obs <- InsertDim(obs, posdim = 1, lendim = 1, name = memb_dim[!(memb_dim %in% obsdims)]) - } if (!all(sample_dims %in% expdims)) { stop("Parameter 'exp' requires 'sdate_dim' and 'member_dim' dimensions.") } -- GitLab From 94b9f0efb8458beb8629d15aa90663f79bc353e6 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Tue, 30 Aug 2022 10:59:36 +0200 Subject: [PATCH 05/14] corrected three typos --- R/CST_QuantileMapping.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 099ac7b1..b0d63cbd 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -58,15 +58,15 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', obsdims <- names(dim(obs)) expdims <- names(dim(exp)) if (is.null(expdims)) { - stop("Parameter 'exp' musth have dimension names.") + stop("Parameter 'exp' must have dimension names.") } if (is.null(obsdims)) { - stop("Parameter 'obs' musth have dimension names.") + stop("Parameter 'obs' must have dimension names.") } if (!is.null(exp_cor)) { exp_cordims <- names(dim(exp_cor)) if (is.null(exp_cordims)) { - stop("Parameter 'exp_cor' musth have dimension names.") + stop("Parameter 'exp_cor' must have dimension names.") } } if (!all(memb_dim %in% obsdims)) { -- GitLab From 522ef580a76002c98efcc731a47047ba6c85acae Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 30 Aug 2022 15:49:38 +0200 Subject: [PATCH 06/14] Fix typos and na.rm = F behaviour --- NAMESPACE | 1 + R/CST_QuantileMapping.R | 87 +++++++++++++++--- man/CST_QuantileMapping.Rd | 102 +++++----------------- man/QuantileMapping.Rd | 56 ++++++------ tests/testthat/test-CST_QuantileMapping.R | 67 +++++++------- 5 files changed, 166 insertions(+), 147 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 34851f66..157042d1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -61,6 +61,7 @@ import(multiApply) import(ncdf4) import(qmap) import(rainfarmr) +import(s2dv) import(stats) importFrom(ClimProjDiags,SelBox) importFrom(RColorBrewer,brewer.pal) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 099ac7b1..f4c4e55b 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -5,11 +5,14 @@ #'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} #'@param exp an object of class \code{s2dv_cube} #'@param obs an object of class \code{s2dv_cube} -#'@parma exp_cor an object of class \code{s2dv_cube} in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'. -#'@param sdate_dim -#'@param memb_dim -#'@param window_dim +#'@param exp_cor an object of class \code{s2dv_cube} in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'. +#'@param sdate_dim a character string indicating the dimension name in which cross-validation would be applied when exp_cor is not provided. 'sdate' by default. +#'@param memb_dim a character string indicating the dimension name where ensemble members are stored in the experimental arrays. 'member' by default. +#'@param window_dim a character string indicating the dimension name where samples has been stored. It can be NULL (default) in case all samples are used. #'@param method a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used. +#'@param na.rm logical indicating ig missing values should be removed (FALSE by default). +#'@param ncores an integer that indicates the number of cores for parallel computations using multiApply function. The default value is NULL (1). +#'@param ... additional parameters to be used by the method choosen. See qmap package. #' #'@return an oject of class \code{s2dv_cube} containing the experimental data after applyingthe quantile mapping correction. #') <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , @@ -19,6 +22,7 @@ #' #'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} #'@examples +#'exp <- NULL #'exp$data <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) #'dim(exp$data) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , #' lat = 6, lon = 7) @@ -36,10 +40,16 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', memb_dim = 'member', window_dim = NULL, method = 'QUANT', na.rm = FALSE, ncores = NULL, ...) { - if (!inherits(exp, 's2dv_cube') || !inherits(exp, 's2dv_cube')) { + 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 (!is.null(exp_cor)) { + if (!inherits(exp_cor, 's2dv_cube')) { + stop("Parameter 'exp_cor' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + } dimnames <- names(dim(exp$data)) QMapped <- QuantileMapping(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, @@ -51,6 +61,43 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', exp$source_files <- c(exp$source_files, obs$source_files) return(exp) } +#'Quaintiles Mapping for seasonal or decadal forecast data +#' +#'@description This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. +#' +#'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +#'@param exp a multidimensional array with named dimensions. +#'@param obs a multidimensional array with named dimensions. +#'@param exp_cor a multidimensional array with named dimensions in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'. +#'@param sdate_dim a character string indicating the dimension name in which cross-validation would be applied when exp_cor is not provided. 'sdate' by default. +#'@param memb_dim a character string indicating the dimension name where ensemble members are stored in the experimental arrays. 'member' by default. +#'@param window_dim a character string indicating the dimension name where samples has been stored. It can be NULL (default) in case all samples are used. +#'@param method a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used. +#'@param na.rm logical indicating ig missing values should be removed (FALSE by default). +#'@param ncores an integer that indicates the number of cores for parallel computations using multiApply function. The default value is NULL (1). +#'@param ... additional parameters to be used by the method choosen. See qmap package. +#' +#'@return an array containing the experimental data after applyingthe quantile mapping correction. +#') <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , +#'@import qmap +#'@import multiApply +#'@import s2dv +#' +#'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} +#'@examples +#'exp <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) +#'dim(exp) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , +#' lat = 6, lon = 7) +#'class(exp) <- 's2dv_cube' +#'obs <- 101 : c(100 + 1 * 1 * 20 * 60 * 6 * 7) +#'dim(obs) <- c(dataset = 1, member = 1, sdate = 20, ftime = 60 , +#' lat = 6, lon = 7) +#'res <- QuantileMapping(exp, obs) +#'exp <- lonlat_data$exp$data +#'obs <- lonlat_data$obs$data +#'res <- QuantileMapping(exp, obs) +#'@export + QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', memb_dim = 'member', window_dim = NULL, method = 'QUANT', @@ -58,17 +105,21 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', obsdims <- names(dim(obs)) expdims <- names(dim(exp)) if (is.null(expdims)) { - stop("Parameter 'exp' musth have dimension names.") + stop("Parameter 'exp' must have dimension names.") } if (is.null(obsdims)) { - stop("Parameter 'obs' musth have dimension names.") + stop("Parameter 'obs' must have dimension names.") } if (!is.null(exp_cor)) { exp_cordims <- names(dim(exp_cor)) if (is.null(exp_cordims)) { - stop("Parameter 'exp_cor' musth have dimension names.") + stop("Parameter 'exp_cor' must have dimension names.") } } + if (!(method %in% c('PTF','DIST','RQUANT','QUANT','SSPLIN'))){ + stop("Parameter 'method' must be one of the following methods: ", + "'PTF','DIST','RQUANT','QUANT','SSPLIN'.") + } if (!all(memb_dim %in% obsdims)) { obs <- InsertDim(obs, posdim = 1, lendim = 1, name = memb_dim[!(memb_dim %in% obsdims)]) @@ -135,9 +186,23 @@ qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUANT', na.rm = FALSE, return(applied[,sd]) }) } else { - adjust <- fitQmap(as.vector(obs[,-sd]), as.vector(exp[,-sd]), - method = method, ...) - applied[,sd] <- doQmap(as.vector(exp[,sd]), adjust, ...) + # na.rm = FALSE shouldn't fail, just return NA + if(any(is.na(obs[,-sd]) || is.na(exp[,-sd]))) { + applied[,sd] <- NA + } else { + adjust <- fitQmap(as.vector(obs[,-sd]), as.vector(exp[,-sd]), + method = method, ...) + exp2 <- exp[,sd] + if (sum(is.na(exp2)) >= 1) { + app <- rep(NA, length(exp2)) + nas_pos <- which(is.na(exp2)) + exp2 <- exp2[!is.na(exp2)] + app[-nas_pos]<- doQmap(as.vector(exp2), adjust, ...) + } else { + app <- doQmap(as.vector(exp2), adjust, ...) + } + applied[,sd] <- app + } } } } else { diff --git a/man/CST_QuantileMapping.Rd b/man/CST_QuantileMapping.Rd index ec5fc8a3..736ea7d9 100644 --- a/man/CST_QuantileMapping.Rd +++ b/man/CST_QuantileMapping.Rd @@ -2,15 +2,17 @@ % Please edit documentation in R/CST_QuantileMapping.R \name{CST_QuantileMapping} \alias{CST_QuantileMapping} -\title{Quantiles Mapping for seasonal or decadal forecast data} +\title{Quaintiles Mapping for seasonal or decadal forecast data} \usage{ CST_QuantileMapping( exp, obs, exp_cor = NULL, - sample_dims = c("sdate", "ftime", "member"), - sample_length = NULL, + sdate_dim = "sdate", + memb_dim = "member", + window_dim = NULL, method = "QUANT", + na.rm = FALSE, ncores = NULL, ... ) @@ -20,102 +22,46 @@ CST_QuantileMapping( \item{obs}{an object of class \code{s2dv_cube}} -\item{exp_cor}{an object of class \code{s2dv_cube} in which the quantile mapping correction will be applied. If it is not specified, the correction is applied in object \code{exp}.} +\item{exp_cor}{an object of class \code{s2dv_cube} in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'.} -\item{sample_dims}{a character vector indicating the dimensions that can be used as sample for the same distribution} +\item{sdate_dim}{a character string indicating the dimension name in which cross-validation would be applied when exp_cor is not provided. 'sdate' by default.} -\item{sample_length}{a numeric value indicating the length of the timeseries window to be used as sample for the sample distribution and correction. By default, NULL, the total length of the timeseries will be used.} +\item{memb_dim}{a character string indicating the dimension name where ensemble members are stored in the experimental arrays. 'member' by default.} + +\item{window_dim}{a character string indicating the dimension name where samples has been stored. It can be NULL (default) in case all samples are used.} \item{method}{a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used.} -\item{ncores}{an integer indicating the number of parallel processes to spawn for the use for parallel computation in multiple cores.} +\item{na.rm}{logical indicating ig missing values should be removed (FALSE by default).} + +\item{ncores}{an integer that indicates the number of cores for parallel computations using multiApply function. The default value is NULL (1).} -\item{...}{additional arguments passed to the method specified by \code{method}.} +\item{...}{additional parameters to be used by the method choosen. See qmap package.} } \value{ an oject of class \code{s2dv_cube} containing the experimental data after applyingthe quantile mapping correction. ) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , } \description{ -This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. The quantile mapping adjustment between an experiment, tipically a hindcast, and observations is applied to the experiment itself or to a provided forecast. -} -\details{ -The different methods are: -\itemize{ -\item{'PTF'} {fits a parametric transformations to the quantile-quantile relation of observed and modelled values. See \code{?qmap::fitQmapPTF}.} -\item{'DIST'} {fits a theoretical distribution to observed and to modelled time series. See \code{?qmap::fitQmapDIST}.} -\item{'RQUANT'} {estimates the values of the quantile-quantile relation of observed and modelled time series for regularly spaced quantiles using local linear least square regression. See \code{?qmap::fitQmapRQUANT}.} -\item{'QUANT'} {estimates values of the empirical cumulative distribution function of observed and modelled time series for regularly spaced quantiles. See \code{?qmap::fitQmapQUANT}.} -\item{'SSPLIN'} {fits a smoothing spline to the quantile-quantile plot of observed and modelled time series. See \code{?qmap::fitQmapSSPLIN}.}} -All methods accepts some common arguments: -\itemize{ -\item{wet.day} {logical indicating whether to perform wet day correction or not.(Not available in 'DIS' method)} -\item{qstep} {NULL or a numeric value between 0 and 1.}} -When providing a forecast to be corrected through the pararmeter \code{exp_cor}, some inputs might need to be modified. The quantile correction is compute by comparing objects passed through 'exp' and 'obs' parameters, this correction will be later applied to the forecast provided in 'exp_cor'. Imaging the case of 'exp' and 'obs' having several start dates, stored using a dimension e.g. 'sdate', 'sample_dims' include this dimension 'sdate' and 'exp_cor' has forecasts for several sdates but different from the ones in 'exp'. In this case, the correction computed with 'exp' and 'obs' would be applied for each 'sdate' of 'exp_cor' separately. This example corresponds to a case of split a dataset in training set and validation set. +This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. } \examples{ -library(qmap) -exp <- 1 : (1 * 5 * 10 * 6 * 2 * 3) -dim(exp) <- c(dataset = 1, member = 10, sdate = 5, ftime = 6 , - lat = 2, lon = 3) -exp <- list(data = exp) +exp <- NULL +exp$data <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) +dim(exp$data) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , + lat = 6, lon = 7) class(exp) <- 's2dv_cube' -obs <- 101 : (100 + 1 * 1 * 5 * 6 * 2 * 3) -dim(obs) <- c(dataset = 1, member = 1, sdate = 5, ftime = 6 , - lat = 2, lon = 3) -obs <- list(data = obs) +obs$data <- 101 : c(100 + 1 * 1 * 20 * 60 * 6 * 7) +dim(obs$data) <- c(dataset = 1, member = 1, sdate = 20, ftime = 60 , + lat = 6, lon = 7) class(obs) <- 's2dv_cube' -res <- CST_QuantileMapping(exp, obs, method = 'RQUANT') -\donttest{ -exp <- lonlat_data$exp -obs <- lonlat_data$obs res <- CST_QuantileMapping(exp, obs) - -exp_cor <- exp -exp_cor$data <- exp_cor$data[,,1,,,] -dim(exp_cor$data) <- c(dataset = 1, member = 15, sdate = 1, ftime = 3, - lat = 22, lon = 53) -res <- CST_QuantileMapping(exp, obs, exp_cor, - sample_dims = c('sdate', 'ftime', 'member')) -res <- CST_QuantileMapping(exp, obs, exp_cor, - sample_dims = c('ftime', 'member')) -data(obsprecip) -data(modprecip) -exp <- modprecip$MOSS[1:10000] -dim(exp) <- c(time = length(exp)) -exp <- list(data = exp) -class(exp) <- 's2dv_cube' -obs <- obsprecip$MOSS[1:10000] -dim(obs) <- c(time = length(obs)) -obs <- list(data = obs) -class(obs) <- 's2dv_cube' -res <- CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time', - method = 'DIST') -# Example using different lenght of members and sdates: exp <- lonlat_data$exp -exp$data <- exp$data[,,1:4,,,] -dim(exp$data) <- c(dataset = 1, member = 15, sdate = 4, ftime = 3, - lat = 22, lon = 53) obs <- lonlat_data$obs -obs$data <- obs$data[,,1:4, ,,] -dim(obs$data) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3, - lat = 22, lon = 53) -exp_cor <- lonlat_data$exp -exp_cor$data <- exp_cor$data[,1:5,5:6,,,] -dim(exp_cor$data) <- c(dataset = 1, member = 5, sdate = 2, ftime = 3, - lat = 22, lon = 53) -res <- CST_QuantileMapping(exp, obs, exp_cor, - sample_dims = c('sdate', 'ftime', 'member')) -exp_cor <- lonlat_data$exp -exp_cor$data <- exp_cor$data[,,5:6,,,] -dim(exp_cor$data) <- c(dataset = 1, member = 15, sdate = 2, ftime = 3, - lat = 22, lon = 53) -res <- CST_QuantileMapping(exp, obs, exp_cor, - sample_dims = c('sdate', 'ftime', 'member')) -} +res <- CST_QuantileMapping(exp, obs) } \seealso{ -\code{qmap::fitQmap} and \code{qmap::doQmap} +\code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} } \author{ Nuria Perez-Zanon, \email{nuria.perez@bsc.es} diff --git a/man/QuantileMapping.Rd b/man/QuantileMapping.Rd index 8771e549..1d797e5a 100644 --- a/man/QuantileMapping.Rd +++ b/man/QuantileMapping.Rd @@ -2,58 +2,64 @@ % Please edit documentation in R/CST_QuantileMapping.R \name{QuantileMapping} \alias{QuantileMapping} -\title{Quantiles Mapping for seasonal or decadal forecast data} +\title{Quaintiles Mapping for seasonal or decadal forecast data} \usage{ QuantileMapping( exp, obs, exp_cor = NULL, - sample_dims = "ftime", - sample_length = NULL, + sdate_dim = "sdate", + memb_dim = "member", + window_dim = NULL, method = "QUANT", + na.rm = FALSE, ncores = NULL, ... ) } \arguments{ -\item{exp}{a multi-dimensional array with named dimensions containing the hindcast.} +\item{exp}{a multidimensional array with named dimensions.} -\item{obs}{a multi-dimensional array with named dimensions (the same as the provided in 'exp') containing the reference dataset.} +\item{obs}{a multidimensional array with named dimensions.} -\item{exp_cor}{a multi-dimensional array with named dimensions in which the quantile mapping correction will be applied. If it is not specified, the correction is applied in object \code{exp}.} +\item{exp_cor}{a multidimensional array with named dimensions in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'.} -\item{sample_dims}{a character vector indicating the dimensions that can be used as sample for the same distribution} +\item{sdate_dim}{a character string indicating the dimension name in which cross-validation would be applied when exp_cor is not provided. 'sdate' by default.} -\item{sample_length}{a numeric value indicating the length of the timeseries window to be used as sample for the sample distribution and correction. By default, NULL, the total length of the timeseries will be used.} +\item{memb_dim}{a character string indicating the dimension name where ensemble members are stored in the experimental arrays. 'member' by default.} + +\item{window_dim}{a character string indicating the dimension name where samples has been stored. It can be NULL (default) in case all samples are used.} \item{method}{a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used.} -\item{ncores}{an integer indicating the number of parallel processes to spawn for the use for parallel computation in multiple cores.} +\item{na.rm}{logical indicating ig missing values should be removed (FALSE by default).} + +\item{ncores}{an integer that indicates the number of cores for parallel computations using multiApply function. The default value is NULL (1).} -\item{...}{additional arguments passed to the method specified by \code{method}.} +\item{...}{additional parameters to be used by the method choosen. See qmap package.} } \value{ -an oject of class \code{s2dv_cube} containing the experimental data after applyingthe quantile mapping correction. +an array containing the experimental data after applyingthe quantile mapping correction. ) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , } \description{ -This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. The quantile mapping adjustment between an experiment, tipically a hindcast, and observations is applied to the experiment itself or to a provided forecast. +This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. } -\details{ -The different methods are: -\itemize{ -\item{'PTF'} {fits a parametric transformations to the quantile-quantile relation of observed and modelled values. See \code{?qmap::fitQmapPTF}.} -\item{'DIST'} {fits a theoretical distribution to observed and to modelled time series. See \code{?qmap::fitQmapDIST}.} -\item{'RQUANT'} {estimates the values of the quantile-quantile relation of observed and modelled time series for regularly spaced quantiles using local linear least square regression. See \code{?qmap::fitQmapRQUANT}.} -\item{'QUANT'} {estimates values of the empirical cumulative distribution function of observed and modelled time series for regularly spaced quantiles. See \code{?qmap::fitQmapQUANT}.} -\item{'SSPLIN'} {fits a smoothing spline to the quantile-quantile plot of observed and modelled time series. See \code{?qmap::fitQmapSSPLIN}.}} -All methods accepts some common arguments: -\itemize{ -\item{wet.day} {logical indicating whether to perform wet day correction or not.(Not available in 'DIS' method)} -\item{qstep} {NULL or a numeric value between 0 and 1.}} +\examples{ +exp <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) +dim(exp) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , + lat = 6, lon = 7) +class(exp) <- 's2dv_cube' +obs <- 101 : c(100 + 1 * 1 * 20 * 60 * 6 * 7) +dim(obs) <- c(dataset = 1, member = 1, sdate = 20, ftime = 60 , + lat = 6, lon = 7) +res <- QuantileMapping(exp, obs) +exp <- lonlat_data$exp$data +obs <- lonlat_data$obs$data +res <- QuantileMapping(exp, obs) } \seealso{ -\code{qmap::fitQmap} and \code{qmap::doQmap} +\code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} } \author{ Nuria Perez-Zanon, \email{nuria.perez@bsc.es} diff --git a/tests/testthat/test-CST_QuantileMapping.R b/tests/testthat/test-CST_QuantileMapping.R index f8482967..967f13d7 100644 --- a/tests/testthat/test-CST_QuantileMapping.R +++ b/tests/testthat/test-CST_QuantileMapping.R @@ -59,55 +59,56 @@ test_that("Sanity checks", { "Parameter 'exp_cor' must have dimension names.") expect_error( CST_QuantileMapping(exp = exp, obs = obs), - paste0("Parameter 'sample_dims' must be a vector of string character ", - "indicating names of exiting dimension in parameter 'exp'.")) + paste0("Parameter 'exp' requires 'sdate_dim' and 'member_dim' dimensions.")) + + dim(exp$data) <- c(sdate = 5, member = 4) + dim(obs$data) <- c(sdate = 5, time = 4) expect_error( - CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time', method = 'x'), + CST_QuantileMapping(exp = exp, obs = obs, method = 'x'), paste0("Parameter 'method' must be one of the following methods: ", "'PTF','DIST','RQUANT','QUANT','SSPLIN'.")) - expect_warning( - CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time', sample_length = "month"), - paste0("Parameter 'sample_length' has not been correctly defined and ", - "the whole length of the timeseries will be used.")) - + expect_error( + CST_QuantileMapping(exp = exp, obs = obs, window_dim = 'ftime'), + paste0("Dimension 'window_dim' not found in 'obs'.")) + res <- NULL + res$data <- t(array(c(0, 2:20), c(sdate = 5, member = 4))) + names(dim(res$data)) <- c('member', 'sdate') + class(res) <- 's2dv_cube' + expect_equal(CST_QuantileMapping(exp = exp, obs = obs, window_dim = 'time'), + res) + exp$data[1] <- NA - expect_warning( - CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time'), - "Parameter 'exp' contains NA values.") + obs$data[1] <- NA + res$data[1,c(1,2)] <- c(NA, 0) + expect_equal( + CST_QuantileMapping(exp = exp, obs = obs, window_dim = 'time', na.rm = T), + res) + + res$data <- array(c(NA, 6, 11, 16, rep(NA, 16)), c(member = 4, sdate = 5)) + expect_equal( + CST_QuantileMapping(exp = exp, obs = obs, window_dim = 'time'), + res) exp$data[1] <- 1 obs$data[1] <- NA - expect_warning( - CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time'), - "Parameter 'obs' contains NA values.") + res$data[1] <- 0 + expect_equal( + CST_QuantileMapping(exp = exp, obs = obs, window_dim = 'time'), + res) obs$data[1] <- 1 - expect_equal(CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time'), exp) - expect_equal(CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time', - method = 'PTF'), exp) - expect_equal(CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time', - method = 'RQUANT'), exp) - expect_equal(CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time', - method = 'SSPLIN'), exp) - library(CSTools) - expect_error(CST_QuantileMapping(exp = lonlat_data$exp, obs = lonlat_data$obs, - exp_cor = lonlat_data$exp), - paste0("Review parameter 'sample_dims' or the data dimensions ", - "since multiple dimensions with different length have being ", - "found in the data inputs that don't match with 'sample_dims' parameter.")) - exp <- lonlat_data$exp + exp <- CSTools::lonlat_data$exp exp$data <- exp$data[,,1:4,,,] dim(exp$data) <- c(dataset = 1, member = 15, sdate = 4, ftime = 3, lat = 22, lon = 53) - obs <- lonlat_data$obs + obs <- CSTools::lonlat_data$obs obs$data <- obs$data[,,1:4, ,,] dim(obs$data) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3, lat = 22, lon = 53) - exp_cor <- lonlat_data$exp + exp_cor <- CSTools::lonlat_data$exp exp_cor$data <- exp_cor$data[,,5:6,,,] dim(exp_cor$data) <- c(dataset = 1, member = 15, sdate = 2, ftime = 3, lat = 22, lon = 53) - expect_warning(CST_QuantileMapping(exp, obs, exp_cor, - sample_dims = c('sdate', 'ftime', 'member')), - "The sample_dims sdate are not used when applying the correction to 'exp_cor'") + expect_equal(length(CST_QuantileMapping(exp, obs, exp_cor)), + 9) }) -- GitLab From d615e957d9b48d200a45d05e83861ebfe7df4c59 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Fri, 2 Sep 2022 15:22:55 +0200 Subject: [PATCH 07/14] fixed condition for na.rm=FALSE --- R/CST_QuantileMapping.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index f4c4e55b..0c4bb4d4 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -187,7 +187,7 @@ qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUANT', na.rm = FALSE, }) } else { # na.rm = FALSE shouldn't fail, just return NA - if(any(is.na(obs[,-sd]) || is.na(exp[,-sd]))) { + if (anyNA(obs[,-sd]) | anyNA(exp[,-sd])) { applied[,sd] <- NA } else { adjust <- fitQmap(as.vector(obs[,-sd]), as.vector(exp[,-sd]), -- GitLab From dcf5374b7cb07a2debdecfcc53a3e4fabaeb3f10 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 16 Sep 2022 18:01:04 +0200 Subject: [PATCH 08/14] output dimension names --- R/CST_QuantileMapping.R | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 0c4bb4d4..73768a60 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -151,17 +151,17 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', } if (!is.null(exp_cor)) { qmaped <- Apply(list(exp, obs, exp_cor), target_dims = sample_dims, - fun = qmapcor, method = method, na.rm = na.rm, ..., + fun = .qmapcor, method = method, na.rm = na.rm, ..., ncores = ncores)$output1 } else { qmaped <- Apply(list(exp, obs), target_dims = sample_dims, - fun = qmapcor, exp_cor = NULL, method = method, + fun = .qmapcor, exp_cor = NULL, method = method, na.rm = na.rm, ..., ncores = ncores)$output1 } return(qmaped) } -qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUANT', na.rm = FALSE, +.qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUANT', na.rm = FALSE, ...) { # exp[memb, sdate] # obs[window, sdate] @@ -205,6 +205,7 @@ qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUANT', na.rm = FALSE, } } } + names(dim(applied)) <- names(dim(exp)) } else { applied <- exp_cor * NA if (na.rm) { @@ -221,7 +222,10 @@ qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUANT', na.rm = FALSE, adjust <- fitQmap(as.vector(obs), as.vector(exp), method = method, ...) applied <- doQmap(as.vector(exp_cor), adjust, ...) } + dim(applied) <- length(applied) + names(dim(applied)) <- names(dim(exp_cor))[1] } return(applied) } + -- GitLab From 9a1978e99c3984d4aeadfe2d62804c92950a970e Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Fri, 7 Oct 2022 09:53:12 +0200 Subject: [PATCH 09/14] Correct format, documentation and output dimensions. --- DESCRIPTION | 2 +- NAMESPACE | 2 +- R/CST_QuantileMapping.R | 244 +++++++++++------- man/CST_QuantileMapping.Rd | 42 ++-- man/QuantileMapping.Rd | 42 ++-- tests/testthat/test-CST_QuantileMapping.R | 293 +++++++++++++++------- 6 files changed, 412 insertions(+), 213 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 52999f2a..87c57324 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -85,4 +85,4 @@ VignetteBuilder: knitr License: Apache License 2.0 Encoding: UTF-8 LazyData: true -RoxygenNote: 7.0.2 +RoxygenNote: 7.2.0 diff --git a/NAMESPACE b/NAMESPACE index 157042d1..638f39b2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -121,4 +121,4 @@ importFrom(utils,head) importFrom(utils,read.table) importFrom(utils,tail) importFrom(verification,verify) -useDynLib(CSTools, .registration = TRUE) +useDynLib(CSTools) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 73768a60..9e042131 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -1,26 +1,36 @@ -#'Quaintiles Mapping for seasonal or decadal forecast data +#'Quantile Mapping for seasonal or decadal forecast data #' -#'@description This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. +#'@description This function is a wrapper from fitQmap and doQmap from package +#''qmap'to be applied in CSTools objects of class 's2dv_cube'. #' #'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} -#'@param exp an object of class \code{s2dv_cube} -#'@param obs an object of class \code{s2dv_cube} -#'@param exp_cor an object of class \code{s2dv_cube} in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'. -#'@param sdate_dim a character string indicating the dimension name in which cross-validation would be applied when exp_cor is not provided. 'sdate' by default. -#'@param memb_dim a character string indicating the dimension name where ensemble members are stored in the experimental arrays. 'member' by default. -#'@param window_dim a character string indicating the dimension name where samples has been stored. It can be NULL (default) in case all samples are used. -#'@param method a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used. -#'@param na.rm logical indicating ig missing values should be removed (FALSE by default). -#'@param ncores an integer that indicates the number of cores for parallel computations using multiApply function. The default value is NULL (1). -#'@param ... additional parameters to be used by the method choosen. See qmap package. +#'@param exp An object of class \code{s2dv_cube}. +#'@param obs An object of class \code{s2dv_cube}. +#'@param exp_cor An object of class \code{s2dv_cube} in which the quantile +#' mapping correction should be applied. If it is not specified, the correction +#' is applied in object 'exp'. +#'@param sdate_dim A character string indicating the dimension name in which +#' cross-validation would be applied when exp_cor is not provided. 'sdate' by +#' default. +#'@param memb_dim A character string indicating the dimension name where +#' ensemble members are stored in the experimental arrays. 'member' by default. +#'@param window_dim A character string indicating the dimension name where +#' samples have been stored. It can be NULL (default) in case all samples are +#' used. +#'@param method A character string indicating the method to be used:'PTF', +#' 'DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping +#' 'QUANT' is used. +#'@param na.rm A logical value indicating ig missing values should be removed +#' (FALSE by default). +#'@param ncores An integer indicating the number of cores for parallel +#' computation using multiApply function. The default value is NULL (1). +#'@param ... Additional parameters to be used by the method choosen. See qmap +#' package. #' -#'@return an oject of class \code{s2dv_cube} containing the experimental data after applyingthe quantile mapping correction. -#') <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , -#'@import qmap -#'@import multiApply -#'@import s2dv -#' -#'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} +#'@return An oject of class \code{s2dv_cube} containing the experimental data +#'after applyingthe quantile mapping correction. ) <- c(dataset = 1, member = +#'10, sdate = 20, ftime = 60 , +#'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} #'@examples #'exp <- NULL #'exp$data <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) @@ -35,6 +45,9 @@ #'exp <- lonlat_data$exp #'obs <- lonlat_data$obs #'res <- CST_QuantileMapping(exp, obs) +#'@import qmap +#'@import multiApply +#'@import s2dv #'@export CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', memb_dim = 'member', window_dim = NULL, @@ -50,40 +63,58 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', "as output by CSTools::CST_Load.") } } - dimnames <- names(dim(exp$data)) + QMapped <- QuantileMapping(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, sdate_dim = sdate_dim, memb_dim = memb_dim, window_dim = window_dim, method = method, na.rm = na.rm, ncores = ncores, ...) - exp$data <- QMapped - exp$Datasets <- c(exp$Datasets, obs$Datasets) - exp$source_files <- c(exp$source_files, obs$source_files) - return(exp) + if (is.null(exp_cor)) { + exp$data <- QMapped + return(exp) + + } else { + exp_cor$data <- QMapped + return(exp_cor) + } + + } -#'Quaintiles Mapping for seasonal or decadal forecast data + +#'Quantile Mapping for seasonal or decadal forecast data #' -#'@description This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. +#'@description This function is a wrapper from fitQmap and doQmap from package +#''qmap'to be applied in CSTools objects of class 's2dv_cube'. #' #'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} -#'@param exp a multidimensional array with named dimensions. -#'@param obs a multidimensional array with named dimensions. -#'@param exp_cor a multidimensional array with named dimensions in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'. -#'@param sdate_dim a character string indicating the dimension name in which cross-validation would be applied when exp_cor is not provided. 'sdate' by default. -#'@param memb_dim a character string indicating the dimension name where ensemble members are stored in the experimental arrays. 'member' by default. -#'@param window_dim a character string indicating the dimension name where samples has been stored. It can be NULL (default) in case all samples are used. -#'@param method a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used. -#'@param na.rm logical indicating ig missing values should be removed (FALSE by default). -#'@param ncores an integer that indicates the number of cores for parallel computations using multiApply function. The default value is NULL (1). -#'@param ... additional parameters to be used by the method choosen. See qmap package. -#' -#'@return an array containing the experimental data after applyingthe quantile mapping correction. -#') <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , -#'@import qmap -#'@import multiApply -#'@import s2dv +#'@param exp A multidimensional array with named dimensions. +#'@param obs A multidimensional array with named dimensions. +#'@param exp_cor A multidimensional array with named dimensions in which the +#' quantile mapping correction should be applied. If it is not specified, the +#' correction is applied in object 'exp'. +#'@param sdate_dim A character string indicating the dimension name in which +#' cross-validation would be applied when exp_cor is not provided. 'sdate' by +#' default. +#'@param memb_dim A character string indicating the dimension name where +#' ensemble members are stored in the experimental arrays. 'member' by +#' default. +#'@param window_dim A character string indicating the dimension name where +#' samples have been stored. It can be NULL (default) in case all samples are +#' used. +#'@param method A character string indicating the method to be used: 'PTF', +#' 'DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping +#' 'QUANT' is used. +#'@param na.rm A logical value indicating ig missing values should be removed +#' (FALSE by default). +#'@param ncores An integer indicating the number of cores for parallel +#' computation using multiApply function. The default value is NULL (1). +#'@param ... Additional parameters to be used by the method choosen. See qmap +#' package. #' -#'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} +#'@return An array containing the experimental data after applying the quantile +#'mapping correction. +#' +#'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} #'@examples #'exp <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) #'dim(exp) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , @@ -96,34 +127,56 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', #'exp <- lonlat_data$exp$data #'obs <- lonlat_data$obs$data #'res <- QuantileMapping(exp, obs) +#'@import qmap +#'@import multiApply +#'@import s2dv #'@export - QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', memb_dim = 'member', window_dim = NULL, method = 'QUANT', na.rm = FALSE, ncores = NULL, ...) { + # exp and obs obsdims <- names(dim(obs)) expdims <- names(dim(exp)) + if (!is.array(exp) || !is.numeric(exp)) { + stop("Parameter 'exp' must be a numeric array.") + } + if (!is.array(obs) || !is.numeric(obs)) { + stop("Parameter 'obs' must be a numeric array.") + } if (is.null(expdims)) { stop("Parameter 'exp' must have dimension names.") } if (is.null(obsdims)) { stop("Parameter 'obs' must have dimension names.") } + # sdate_dim + if (!is.character(sdate_dim) | length(sdate_dim) != 1) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% expdims | !sdate_dim %in% obsdims) { + stop("Parameter 'sdate_dim' is not found in 'exp' or 'obs' dimension.") + } + if (dim(exp)[sdate_dim] == 1 || dim(obs)[sdate_dim] == 1) { + stop("Parameter 'exp' and 'obs' must have dimension length of 'sdate_dim' bigger than 1.") + } + # exp_cor if (!is.null(exp_cor)) { - exp_cordims <- names(dim(exp_cor)) - if (is.null(exp_cordims)) { + if (is.null(names(dim(exp_cor)))) { stop("Parameter 'exp_cor' must have dimension names.") } } + # method if (!(method %in% c('PTF','DIST','RQUANT','QUANT','SSPLIN'))){ stop("Parameter 'method' must be one of the following methods: ", "'PTF','DIST','RQUANT','QUANT','SSPLIN'.") - } + } + # memb_dim if (!all(memb_dim %in% obsdims)) { obs <- InsertDim(obs, posdim = 1, lendim = 1, name = memb_dim[!(memb_dim %in% obsdims)]) } + # window_dim if (!is.null(window_dim)) { if (!(window_dim %in% obsdims)) { stop("Dimension 'window_dim' not found in 'obs'.") @@ -134,10 +187,12 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', warning("window_dim found in exp and it is merged to memb_dim.") } } + # sdate_dim and member_dim sample_dims <- c(memb_dim, sdate_dim) if (!all(sample_dims %in% expdims)) { stop("Parameter 'exp' requires 'sdate_dim' and 'member_dim' dimensions.") } + # method if (!is.character(method)) { warning("Parameter 'method' must be a character string indicating ", "one of the following methods: 'PTF', 'DIST', 'RQUANT', @@ -149,31 +204,52 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', " is used.") method <- method[1] } + if (method == 'DIST') { + stop("DIST method does not work. Use another method.") + } + # na.rm + if (!is.logical(na.rm) | length(na.rm) > 1) { + stop("Parameter 'na.rm' must be one logical value.") + } + # ncores + if (!is.null(ncores)) { + if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | + length(ncores) > 1) { + stop("Parameter 'ncores' must be either NULL or a positive integer.") + } + } + + ############################### + if (!is.null(exp_cor)) { qmaped <- Apply(list(exp, obs, exp_cor), target_dims = sample_dims, - fun = .qmapcor, method = method, na.rm = na.rm, ..., + fun = .qmapcor, method = method, sdate_dim = sdate_dim, + na.rm = na.rm, ..., ncores = ncores)$output1 } else { qmaped <- Apply(list(exp, obs), target_dims = sample_dims, fun = .qmapcor, exp_cor = NULL, method = method, - na.rm = na.rm, ..., + sdate_dim = sdate_dim, na.rm = na.rm, ..., ncores = ncores)$output1 } return(qmaped) } -.qmapcor <- function(exp, obs, exp_cor = NULL, method = 'QUANT', na.rm = FALSE, - ...) { + +.qmapcor <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', method = 'QUANT', + na.rm = FALSE, ...) { + # exp[memb, sdate] # obs[window, sdate] + if (is.null(exp_cor)) { applied <- exp * NA - for (sd in 1:dim(exp)['sdate']) { - if (na.rm) { + for (sd in 1:dim(exp)[sdate_dim]) { + if (na.rm) { # select start date for cross-val - nas_pos <- which(!is.na(exp[,sd])) - obs2 <- as.vector(obs[,-sd]) - exp2 <- as.vector(exp[,-sd]) - exp_cor2 <- as.vector(exp[,sd]) + nas_pos <- which(!is.na(exp[, sd])) + obs2 <- as.vector(obs[, -sd]) + exp2 <- as.vector(exp[, -sd]) + exp_cor2 <- as.vector(exp[, sd]) # remove NAs obs2 <- obs2[!is.na(obs2)] exp2 <- exp2[!is.na(exp2)] @@ -181,51 +257,47 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', tryCatch({ adjust <- fitQmap(obs2, exp2, method = method, ...) applied[nas_pos, sd] <- doQmap(exp_cor2, adjust, ...) - }, + }, error = function(error_message) { - return(applied[,sd]) + return(applied[, sd]) }) } else { # na.rm = FALSE shouldn't fail, just return NA - if (anyNA(obs[,-sd]) | anyNA(exp[,-sd])) { - applied[,sd] <- NA + if (anyNA(obs[, -sd]) | anyNA(exp[, -sd])) { + applied[, sd] <- NA } else { - adjust <- fitQmap(as.vector(obs[,-sd]), as.vector(exp[,-sd]), + adjust <- fitQmap(as.vector(obs[, -sd]), as.vector(exp[, -sd]), method = method, ...) - exp2 <- exp[,sd] + exp2 <- exp[, sd] if (sum(is.na(exp2)) >= 1) { app <- rep(NA, length(exp2)) nas_pos <- which(is.na(exp2)) exp2 <- exp2[!is.na(exp2)] - app[-nas_pos]<- doQmap(as.vector(exp2), adjust, ...) + app[-nas_pos] <- doQmap(as.vector(exp2), adjust, ...) } else { app <- doQmap(as.vector(exp2), adjust, ...) } - applied[,sd] <- app + applied[, sd] <- app } } } - names(dim(applied)) <- names(dim(exp)) } else { - applied <- exp_cor * NA - if (na.rm) { - tryCatch({ - adjust <- fitQmap(obs[!is.na(obs)], exp[!is.na(exp)], - method = method, ...) - applied[!is.na(exp_cor)] <- doQmap(exp_cor[!is.na(exp_cor)], - adjust, ...) - }, - error = function(error_message) { - return(applied) - }) - } else { - adjust <- fitQmap(as.vector(obs), as.vector(exp), method = method, ...) - applied <- doQmap(as.vector(exp_cor), adjust, ...) - } - dim(applied) <- length(applied) - names(dim(applied)) <- names(dim(exp_cor))[1] + applied <- exp_cor * NA + if (na.rm) { + tryCatch({ + adjust <- fitQmap(obs[!is.na(obs)], exp[!is.na(exp)], + method = method, ...) + applied[!is.na(exp_cor)] <- doQmap(exp_cor[!is.na(exp_cor)], + adjust, ...) + }, + error = function(error_message) { + return(applied) + }) + } else { + adjust <- fitQmap(as.vector(obs), as.vector(exp), method = method, ...) + applied <- doQmap(as.vector(exp_cor), adjust, ...) + } + dim(applied) <- dim(exp_cor) } return(applied) -} - - +} \ No newline at end of file diff --git a/man/CST_QuantileMapping.Rd b/man/CST_QuantileMapping.Rd index 736ea7d9..0fdd6e0c 100644 --- a/man/CST_QuantileMapping.Rd +++ b/man/CST_QuantileMapping.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/CST_QuantileMapping.R \name{CST_QuantileMapping} \alias{CST_QuantileMapping} -\title{Quaintiles Mapping for seasonal or decadal forecast data} +\title{Quantile Mapping for seasonal or decadal forecast data} \usage{ CST_QuantileMapping( exp, @@ -18,32 +18,46 @@ CST_QuantileMapping( ) } \arguments{ -\item{exp}{an object of class \code{s2dv_cube}} +\item{exp}{An object of class \code{s2dv_cube}.} -\item{obs}{an object of class \code{s2dv_cube}} +\item{obs}{An object of class \code{s2dv_cube}.} -\item{exp_cor}{an object of class \code{s2dv_cube} in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'.} +\item{exp_cor}{An object of class \code{s2dv_cube} in which the quantile +mapping correction should be applied. If it is not specified, the correction +is applied in object 'exp'.} -\item{sdate_dim}{a character string indicating the dimension name in which cross-validation would be applied when exp_cor is not provided. 'sdate' by default.} +\item{sdate_dim}{A character string indicating the dimension name in which +cross-validation would be applied when exp_cor is not provided. 'sdate' by +default.} -\item{memb_dim}{a character string indicating the dimension name where ensemble members are stored in the experimental arrays. 'member' by default.} +\item{memb_dim}{A character string indicating the dimension name where +ensemble members are stored in the experimental arrays. 'member' by default.} -\item{window_dim}{a character string indicating the dimension name where samples has been stored. It can be NULL (default) in case all samples are used.} +\item{window_dim}{A character string indicating the dimension name where +samples have been stored. It can be NULL (default) in case all samples are +used.} -\item{method}{a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used.} +\item{method}{A character string indicating the method to be used:'PTF', +'DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping +'QUANT' is used.} -\item{na.rm}{logical indicating ig missing values should be removed (FALSE by default).} +\item{na.rm}{A logical value indicating ig missing values should be removed +(FALSE by default).} -\item{ncores}{an integer that indicates the number of cores for parallel computations using multiApply function. The default value is NULL (1).} +\item{ncores}{An integer indicating the number of cores for parallel +computation using multiApply function. The default value is NULL (1).} -\item{...}{additional parameters to be used by the method choosen. See qmap package.} +\item{...}{Additional parameters to be used by the method choosen. See qmap +package.} } \value{ -an oject of class \code{s2dv_cube} containing the experimental data after applyingthe quantile mapping correction. -) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , +An oject of class \code{s2dv_cube} containing the experimental data +after applyingthe quantile mapping correction. ) <- c(dataset = 1, member = +10, sdate = 20, ftime = 60 , } \description{ -This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. +This function is a wrapper from fitQmap and doQmap from package +'qmap'to be applied in CSTools objects of class 's2dv_cube'. } \examples{ exp <- NULL diff --git a/man/QuantileMapping.Rd b/man/QuantileMapping.Rd index 1d797e5a..41a87bc3 100644 --- a/man/QuantileMapping.Rd +++ b/man/QuantileMapping.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/CST_QuantileMapping.R \name{QuantileMapping} \alias{QuantileMapping} -\title{Quaintiles Mapping for seasonal or decadal forecast data} +\title{Quantile Mapping for seasonal or decadal forecast data} \usage{ QuantileMapping( exp, @@ -18,32 +18,46 @@ QuantileMapping( ) } \arguments{ -\item{exp}{a multidimensional array with named dimensions.} +\item{exp}{A multidimensional array with named dimensions.} -\item{obs}{a multidimensional array with named dimensions.} +\item{obs}{A multidimensional array with named dimensions.} -\item{exp_cor}{a multidimensional array with named dimensions in which the quantile mapping correction should be applied. If it is not specified, the correction is applied in object 'exp'.} +\item{exp_cor}{A multidimensional array with named dimensions in which the +quantile mapping correction should be applied. If it is not specified, the +correction is applied in object 'exp'.} -\item{sdate_dim}{a character string indicating the dimension name in which cross-validation would be applied when exp_cor is not provided. 'sdate' by default.} +\item{sdate_dim}{A character string indicating the dimension name in which +cross-validation would be applied when exp_cor is not provided. 'sdate' by +default.} -\item{memb_dim}{a character string indicating the dimension name where ensemble members are stored in the experimental arrays. 'member' by default.} +\item{memb_dim}{A character string indicating the dimension name where +ensemble members are stored in the experimental arrays. 'member' by +default.} -\item{window_dim}{a character string indicating the dimension name where samples has been stored. It can be NULL (default) in case all samples are used.} +\item{window_dim}{A character string indicating the dimension name where +samples have been stored. It can be NULL (default) in case all samples are +used.} -\item{method}{a character string indicating the method to be used: 'PTF','DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used.} +\item{method}{A character string indicating the method to be used: 'PTF', +'DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping +'QUANT' is used.} -\item{na.rm}{logical indicating ig missing values should be removed (FALSE by default).} +\item{na.rm}{A logical value indicating ig missing values should be removed +(FALSE by default).} -\item{ncores}{an integer that indicates the number of cores for parallel computations using multiApply function. The default value is NULL (1).} +\item{ncores}{An integer indicating the number of cores for parallel +computation using multiApply function. The default value is NULL (1).} -\item{...}{additional parameters to be used by the method choosen. See qmap package.} +\item{...}{Additional parameters to be used by the method choosen. See qmap +package.} } \value{ -an array containing the experimental data after applyingthe quantile mapping correction. -) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , +An array containing the experimental data after applying the quantile +mapping correction. } \description{ -This function is a wrapper from fitQmap and doQmap from package 'qmap'to be applied in CSTools objects of class 's2dv_cube'. +This function is a wrapper from fitQmap and doQmap from package +'qmap'to be applied in CSTools objects of class 's2dv_cube'. } \examples{ exp <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) diff --git a/tests/testthat/test-CST_QuantileMapping.R b/tests/testthat/test-CST_QuantileMapping.R index 967f13d7..e8ad7e70 100644 --- a/tests/testthat/test-CST_QuantileMapping.R +++ b/tests/testthat/test-CST_QuantileMapping.R @@ -1,114 +1,213 @@ -context("Generic tests") -test_that("Sanity checks", { +context("CSTools::CST_QuantileMapping tests") + +############################################## + +# dat1 +exp1 <- list(data = array(1:20, dim = c(time = 20))) +class(exp1) <- 's2dv_cube' + +obs1 <- list(data = array(1:20, dim = c(time = 20))) +class(obs1) <- 's2dv_cube' + +exp1_2 <- list(data = array(1:20, dim = c(20))) +class(exp1_2) <- 's2dv_cube' + +obs1_2 <- list(data = array(1:20, dim = c(20))) +class(obs1_2) <- 's2dv_cube' + +exp_cor1 <- list(data = array(1:20, dim = c(20))) +class(exp_cor1) <- 's2dv_cube' + +# dat2 +exp2 <- list(data = array(1:20, dim = c(sdate = 5, member = 4))) +class(exp2) <- 's2dv_cube' + +obs2 <- list(data = array(1:20, dim = c(sdate = 5, time = 4))) +class(obs2) <- 's2dv_cube' + +res2 <- NULL +res2$data <- t(array(c(0, 2:20), c(sdate = 5, member = 4))) +names(dim(res2$data)) <- c('member', 'sdate') +class(res2) <- 's2dv_cube' + +# dat3 +exp3 <- exp2 +obs3 <- obs2 +res3 <- res2 +exp3$data[1] <- NA +obs3$data[1] <- NA +res3$data[c(1,2),1] <- c(NA,6) +res3$data[1,2] <- 0 + +res3_1 <- res3 +res3_1$data <- array(c(NA, 6, 11, 16, rep(NA, 16)), c(member = 4, sdate = 5)) + +exp3_2 <- exp3 +obs3_2 <- obs3 +res3_2 <- res3_1 +exp3_2$data[1] <- 1 +obs3_2$data[1] <- NA +res3_2$data[1] <- 0 + +# dat4 +exp4 <- lonlat_data$exp +exp4$data <- exp4$data[,,1:4,,,] +dim(exp4$data) <- c(dataset = 1, member = 15, sdate = 4, ftime = 3, + lat = 22, lon = 53) +obs4 <- lonlat_data$obs +obs4$data <- obs4$data[,,1:4, ,,] +dim(obs4$data) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3, + lat = 22, lon = 53) +exp_cor4 <- lonlat_data$exp +exp_cor4$data <- exp_cor4$data[,,5:6,,,] +dim(exp_cor4$data) <- c(dataset = 1, member = 15, sdate = 2, ftime = 3, + lat = 22, lon = 53) + +# dat5 +exp5 <- lonlat_data$exp +obs5 <- lonlat_data$obs +set.seed(1) +res5 <- NULL +res5$data <- array(rnorm(length(exp5)), dim = c(member = 15, sdate = 6, + dataset = 1, ftime = 3, lat = 22, lon = 53)) +class(res5) <- "s2dv_cube" + +res5_1 <- NULL +res5_1$data <- array(rnorm(length(res5_1)), dim = c(member = 15, ftime = 3, + dataset = 1, sdate = 6, lat = 22, lon = 53)) +class(res5_1) <- "s2dv_cube" + +# dat6 +exp6 <- lonlat_data$exp +obs6 <- lonlat_data$obs +obs6$data <- s2dv::InsertDim(obs6$data, pos = 1, len = 4, name = 'window') + +obs6_1 <- obs6 +obs6_1$data[2] <- NA + +exp6_1 <- exp6 + +exp6_1$data[1,,,1,1,1] <- NA +exp_cor6_1 <- exp6_1 +exp_cor6_1$data <- ClimProjDiags::Subset(exp_cor6_1$data, 'sdate', 1) + +############################################## + +test_that("1. Sanity checks", { + + # s2dv_cube expect_error( CST_QuantileMapping(exp = 1), paste0("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", - "as output by CSTools::CST_Load.")) - - exp <- 1 : 20 - dim(exp) <- c(time = 20) - exp <- list(data = exp) - class(exp) <- 's2dv_cube' + "as output by CSTools::CST_Load.") + ) expect_error( - CST_QuantileMapping(exp = exp), - 'argument "obs" is missing, with no default') + CST_QuantileMapping(exp = exp1), + 'argument "obs" is missing, with no default' + ) expect_error( - CST_QuantileMapping(exp = exp, obs = 1), + CST_QuantileMapping(exp = exp1, obs = 1), paste0("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", - "as output by CSTools::CST_Load.")) - - obs <- 1 : 20 - dim(obs) <- c(time = 20) - obs <- list(data = obs) - class(obs) <- 's2dv_cube' + "as output by CSTools::CST_Load.") + ) expect_error( - CST_QuantileMapping(exp = exp, obs = obs, exp_cor = 1), + CST_QuantileMapping(exp = exp1, obs = obs1, exp_cor = 1), paste0("Parameter 'exp_cor' must be of the class 's2dv_cube', as output ", - "by CSTools::CST_Load.")) - - exp <- 1 : 20 - dim(exp) <- 20 - exp <- list(data = exp) - class(exp) <- 's2dv_cube' + "by CSTools::CST_Load.") + ) + # exp and obs expect_error( - CST_QuantileMapping(exp = exp, obs = obs), - "Parameter 'exp' must have dimension names.") - - exp <- 1 : 20 - dim(exp) <- c(time = 20) - exp <- list(data = exp) - class(exp) <- 's2dv_cube' - obs <- 1 : 20 - dim(obs) <- 20 - obs <- list(data = obs) - class(obs) <- 's2dv_cube' + CST_QuantileMapping(exp = exp1_2, obs = obs1), + "Parameter 'exp' must have dimension names." + ) expect_error( - CST_QuantileMapping(exp = exp, obs = obs), - "Parameter 'obs' must have dimension names.") - - obs <- 1 : 20 - dim(obs) <- c(time = 20) - obs <- list(data = obs) - class(obs) <- 's2dv_cube' - exp_cor <- 1 : 20 - dim(exp_cor) <- 20 - exp_cor <- list(data = exp_cor) - class(exp_cor) <- 's2dv_cube' + CST_QuantileMapping(exp = exp1, obs = obs1_2), + "Parameter 'obs' must have dimension names." + ) + # exp_cor expect_error( - CST_QuantileMapping(exp = exp, obs = obs, exp_cor = exp_cor), - "Parameter 'exp_cor' must have dimension names.") + CST_QuantileMapping(exp = exp1, obs = obs1, exp_cor = exp_cor1, sdate_dim = 'time'), + "Parameter 'exp_cor' must have dimension names." + ) + # sdate_dim and member_dim expect_error( - CST_QuantileMapping(exp = exp, obs = obs), - paste0("Parameter 'exp' requires 'sdate_dim' and 'member_dim' dimensions.")) - - dim(exp$data) <- c(sdate = 5, member = 4) - dim(obs$data) <- c(sdate = 5, time = 4) + CST_QuantileMapping(exp = exp1, obs = obs1, sdate_dim = 'time'), + paste0("Parameter 'exp' requires 'sdate_dim' and 'member_dim' dimensions.") + ) + # window_dim expect_error( - CST_QuantileMapping(exp = exp, obs = obs, method = 'x'), - paste0("Parameter 'method' must be one of the following methods: ", - "'PTF','DIST','RQUANT','QUANT','SSPLIN'.")) + CST_QuantileMapping(exp = exp2, obs = obs2, window_dim = 'ftime'), + paste0("Dimension 'window_dim' not found in 'obs'.") + ) + # method expect_error( - CST_QuantileMapping(exp = exp, obs = obs, window_dim = 'ftime'), - paste0("Dimension 'window_dim' not found in 'obs'.")) - res <- NULL - res$data <- t(array(c(0, 2:20), c(sdate = 5, member = 4))) - names(dim(res$data)) <- c('member', 'sdate') - class(res) <- 's2dv_cube' - expect_equal(CST_QuantileMapping(exp = exp, obs = obs, window_dim = 'time'), - res) - - exp$data[1] <- NA - obs$data[1] <- NA - res$data[1,c(1,2)] <- c(NA, 0) - expect_equal( - CST_QuantileMapping(exp = exp, obs = obs, window_dim = 'time', na.rm = T), - res) - - res$data <- array(c(NA, 6, 11, 16, rep(NA, 16)), c(member = 4, sdate = 5)) + CST_QuantileMapping(exp = exp2, obs = obs2, method = 'x'), + paste0("Parameter 'method' must be one of the following methods: ", + "'PTF','DIST','RQUANT','QUANT','SSPLIN'.") + ) +}) + +############################################## + +test_that("2. dat2, dat3 and dat4", { + expect_equal( + CST_QuantileMapping(exp = exp2, obs = obs2, window_dim = 'time'), + res2 + ) + expect_equal( + CST_QuantileMapping(exp = exp3, obs = obs3, window_dim = 'time', na.rm = T), + res3 + ) + expect_equal( + CST_QuantileMapping(exp = exp3, obs = obs3, window_dim = 'time'), + res3_1 + ) + expect_equal( + CST_QuantileMapping(exp = exp3_2, obs = obs3_2, window_dim = 'time'), + res3_2 + ) + expect_equal( + length(CST_QuantileMapping(exp4, obs4, exp_cor4)), + 9 + ) +}) + +############################################## + +test_that("3. dat5", { expect_equal( - CST_QuantileMapping(exp = exp, obs = obs, window_dim = 'time'), - res) - - exp$data[1] <- 1 - obs$data[1] <- NA - res$data[1] <- 0 - expect_equal( - CST_QuantileMapping(exp = exp, obs = obs, window_dim = 'time'), - res) - obs$data[1] <- 1 - exp <- CSTools::lonlat_data$exp - exp$data <- exp$data[,,1:4,,,] - dim(exp$data) <- c(dataset = 1, member = 15, sdate = 4, ftime = 3, - lat = 22, lon = 53) - obs <- CSTools::lonlat_data$obs - obs$data <- obs$data[,,1:4, ,,] - dim(obs$data) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3, - lat = 22, lon = 53) - exp_cor <- CSTools::lonlat_data$exp - exp_cor$data <- exp_cor$data[,,5:6,,,] - dim(exp_cor$data) <- c(dataset = 1, member = 15, sdate = 2, ftime = 3, - lat = 22, lon = 53) - expect_equal(length(CST_QuantileMapping(exp, obs, exp_cor)), - 9) + dim(CST_QuantileMapping(exp5, obs5, ncores = 6)$data), + dim(res5$data) + ) + expect_equal( + dim(CST_QuantileMapping(exp5, obs5)$data), + dim(res5$data) + ) + expect_equal( + dim(CST_QuantileMapping(exp5, obs5, sdate_dim = "ftime", ncores = 6)$data), + dim(res5_1$data) + ) +}) +############################################## + +test_that("4. dat6", { + expect_equal( + CST_QuantileMapping(exp6, obs6, window_dim = 'window', ncores = 6), + CST_QuantileMapping(exp6, obs6, window_dim = 'window', ncores = 6, na.rm = TRUE) + ) + expect_equal( + dim(CST_QuantileMapping(exp6, obs6_1, window_dim = 'window', ncores = 6)), + dim(CST_QuantileMapping(exp6, obs6_1, window_dim = 'window', ncores = 6, na.rm = TRUE)) + ) + expect_equal( + sum(is.na(CST_QuantileMapping(exp6_1, obs6_1, exp_cor = exp_cor6_1, window_dim = 'window', ncores = 6, na.rm = TRUE)$data)), + sum(is.na(exp_cor6_1$data)) + ) + expect_error( + CST_QuantileMapping(exp6_1, obs6_1, exp_cor = exp_cor6_1, method = 'DIST', + window_dim = 'window', ncores = 6, + na.rm = TRUE, distr="berngamma"), + "DIST method does not work. Use another method." + ) }) -- GitLab From c08dc096a32de2c8b32b399b8a7c72406e239959 Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Mon, 10 Oct 2022 15:56:54 +0200 Subject: [PATCH 10/14] Correct initial checks for DIST method --- R/CST_QuantileMapping.R | 3 --- tests/testthat/test-CST_QuantileMapping.R | 10 +++++----- 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 9e042131..bdb183e4 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -204,9 +204,6 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', " is used.") method <- method[1] } - if (method == 'DIST') { - stop("DIST method does not work. Use another method.") - } # na.rm if (!is.logical(na.rm) | length(na.rm) > 1) { stop("Parameter 'na.rm' must be one logical value.") diff --git a/tests/testthat/test-CST_QuantileMapping.R b/tests/testthat/test-CST_QuantileMapping.R index e8ad7e70..7c4a2481 100644 --- a/tests/testthat/test-CST_QuantileMapping.R +++ b/tests/testthat/test-CST_QuantileMapping.R @@ -90,6 +90,8 @@ exp6_1 <- exp6 exp6_1$data[1,,,1,1,1] <- NA exp_cor6_1 <- exp6_1 exp_cor6_1$data <- ClimProjDiags::Subset(exp_cor6_1$data, 'sdate', 1) +exp_cor6_2 <- exp6 +exp_cor6_2$data <- ClimProjDiags::Subset(exp_cor6_2$data, 'member', 1:5) ############################################## @@ -204,10 +206,8 @@ test_that("4. dat6", { sum(is.na(CST_QuantileMapping(exp6_1, obs6_1, exp_cor = exp_cor6_1, window_dim = 'window', ncores = 6, na.rm = TRUE)$data)), sum(is.na(exp_cor6_1$data)) ) - expect_error( - CST_QuantileMapping(exp6_1, obs6_1, exp_cor = exp_cor6_1, method = 'DIST', - window_dim = 'window', ncores = 6, - na.rm = TRUE, distr="berngamma"), - "DIST method does not work. Use another method." + expect_equal( + dim(CST_QuantileMapping(exp6, obs6_1, exp_cor6_2, window_dim = 'window', ncores = 6)$data), + c(member = 5, sdate = 6, dataset = 1, ftime = 3, lat = 22, lon = 53) ) }) -- GitLab From 6d71d321d7504bd50924a5343b9a1c6365585e6e Mon Sep 17 00:00:00 2001 From: aho Date: Tue, 18 Oct 2022 17:15:28 +0200 Subject: [PATCH 11/14] Fix typos; improve sanity checks --- R/CST_QuantileMapping.R | 97 ++++++++++++++--------- man/CST_QuantileMapping.Rd | 15 ++-- man/QuantileMapping.Rd | 18 +++-- tests/testthat/test-CST_QuantileMapping.R | 30 ++++--- 4 files changed, 94 insertions(+), 66 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index bdb183e4..0a89692d 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -1,8 +1,10 @@ #'Quantile Mapping for seasonal or decadal forecast data #' -#'@description This function is a wrapper from fitQmap and doQmap from package -#''qmap'to be applied in CSTools objects of class 's2dv_cube'. -#' +#'@description This function is a wrapper of fitQmap and doQmap from package +#''qmap' to be applied on the object of class 's2dv_cube'. The quantile mapping +#'adjustment between an experiment, typically a hindcast, and observation is +#'applied to the experiment itself or to a provided forecast. + #'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} #'@param exp An object of class \code{s2dv_cube}. #'@param obs An object of class \code{s2dv_cube}. @@ -20,16 +22,16 @@ #'@param method A character string indicating the method to be used:'PTF', #' 'DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping #' 'QUANT' is used. -#'@param na.rm A logical value indicating ig missing values should be removed +#'@param na.rm A logical value indicating if missing values should be removed #' (FALSE by default). #'@param ncores An integer indicating the number of cores for parallel #' computation using multiApply function. The default value is NULL (1). #'@param ... Additional parameters to be used by the method choosen. See qmap -#' package. +#' package for details. +#' +#'@return An object of class \code{s2dv_cube} containing the experimental data +#'after applying the quantile mapping correction. #' -#'@return An oject of class \code{s2dv_cube} containing the experimental data -#'after applyingthe quantile mapping correction. ) <- c(dataset = 1, member = -#'10, sdate = 20, ftime = 60 , #'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} #'@examples #'exp <- NULL @@ -83,15 +85,19 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', #'Quantile Mapping for seasonal or decadal forecast data #' -#'@description This function is a wrapper from fitQmap and doQmap from package -#''qmap'to be applied in CSTools objects of class 's2dv_cube'. +#'@description This function is a wrapper of fitQmap and doQmap from package +#''qmap' to be applied on multi-dimensional arrays. The quantile mapping +#'adjustment between an experiment, typically a hindcast, and observation is +#'applied to the experiment itself or to a provided forecast. #' #'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} -#'@param exp A multidimensional array with named dimensions. -#'@param obs A multidimensional array with named dimensions. +#'@param exp A multidimensional array with named dimensions containing the +#' hindcast. +#'@param obs A multidimensional array with named dimensions containing the +#' reference dataset. #'@param exp_cor A multidimensional array with named dimensions in which the #' quantile mapping correction should be applied. If it is not specified, the -#' correction is applied in object 'exp'. +#' correction is applied on object 'exp'. #'@param sdate_dim A character string indicating the dimension name in which #' cross-validation would be applied when exp_cor is not provided. 'sdate' by #' default. @@ -104,12 +110,12 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', #'@param method A character string indicating the method to be used: 'PTF', #' 'DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping #' 'QUANT' is used. -#'@param na.rm A logical value indicating ig missing values should be removed +#'@param na.rm A logical value indicating if missing values should be removed #' (FALSE by default). #'@param ncores An integer indicating the number of cores for parallel #' computation using multiApply function. The default value is NULL (1). #'@param ... Additional parameters to be used by the method choosen. See qmap -#' package. +#' package for details. #' #'@return An array containing the experimental data after applying the quantile #'mapping correction. @@ -151,6 +157,8 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', stop("Parameter 'obs' must have dimension names.") } # sdate_dim +#-------NEW------------ + if (!is.null(exp_cor)) { if (!is.character(sdate_dim) | length(sdate_dim) != 1) { stop("Parameter 'sdate_dim' must be a character string.") } @@ -160,50 +168,60 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', if (dim(exp)[sdate_dim] == 1 || dim(obs)[sdate_dim] == 1) { stop("Parameter 'exp' and 'obs' must have dimension length of 'sdate_dim' bigger than 1.") } + } +#--------NEW_END---------- # exp_cor if (!is.null(exp_cor)) { if (is.null(names(dim(exp_cor)))) { stop("Parameter 'exp_cor' must have dimension names.") } } +#--------NEW----------- # method - if (!(method %in% c('PTF','DIST','RQUANT','QUANT','SSPLIN'))){ + if (!(method %in% c('PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN')) | length(method) != 1) { stop("Parameter 'method' must be one of the following methods: ", - "'PTF','DIST','RQUANT','QUANT','SSPLIN'.") + "'PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN'.") } +#--------NEW_END----------- # memb_dim if (!all(memb_dim %in% obsdims)) { obs <- InsertDim(obs, posdim = 1, lendim = 1, name = memb_dim[!(memb_dim %in% obsdims)]) } +#--------NEW----------- + if (any(!memb_dim %in% expdims)) { + stop("Parameter 'memb_dim' is not found in 'exp' dimensions.") + } + sample_dims <- c(memb_dim, sdate_dim) # window_dim if (!is.null(window_dim)) { if (!(window_dim %in% obsdims)) { - stop("Dimension 'window_dim' not found in 'obs'.") + stop("Parameter 'window_dim' is not found in 'obs'.") } obs <- CSTools::MergeDims(obs, c(memb_dim, window_dim)) if (window_dim %in% expdims) { exp <- CSTools::MergeDims(exp, c(memb_dim, window_dim)) - warning("window_dim found in exp and it is merged to memb_dim.") + warning("Parameter 'window_dim' is found in exp and is merged to 'memb_dim'.") } } - # sdate_dim and member_dim - sample_dims <- c(memb_dim, sdate_dim) - if (!all(sample_dims %in% expdims)) { - stop("Parameter 'exp' requires 'sdate_dim' and 'member_dim' dimensions.") - } - # method - if (!is.character(method)) { - warning("Parameter 'method' must be a character string indicating ", - "one of the following methods: 'PTF', 'DIST', 'RQUANT', - 'QUANT', 'SSPLIN'. Method 'QUANT' is being used.") - method = 'QUANT' - } - if (length(method) > 1) { - warning("Parameter 'method' has length > 1 and only the first element", - " is used.") - method <- method[1] - } +# # sdate_dim and member_dim +# sample_dims <- c(memb_dim, sdate_dim) +# if (!all(sample_dims %in% expdims)) { +# stop("Parameter 'exp' requires 'sdate_dim' and 'memb_dim' dimensions.") +# } +# # method +# if (!is.character(method)) { +# warning("Parameter 'method' must be a character string indicating ", +# "one of the following methods: 'PTF', 'DIST', 'RQUANT', +# 'QUANT', 'SSPLIN'. Method 'QUANT' is being used.") +# method <- 'QUANT' +# } +# if (length(method) > 1) { +# warning("Parameter 'method' has length > 1 and only the first element", +# " is used.") +# method <- method[1] +# } +#-------NEW_END------------- # na.rm if (!is.logical(na.rm) | length(na.rm) > 1) { stop("Parameter 'na.rm' must be one logical value.") @@ -235,8 +253,9 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', .qmapcor <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', method = 'QUANT', na.rm = FALSE, ...) { - # exp[memb, sdate] - # obs[window, sdate] + # exp: [memb (+ window), sdate] + # obs: [memb (+ window), sdate] + # exp_cor: NULL or [memb, sdate] if (is.null(exp_cor)) { applied <- exp * NA @@ -297,4 +316,4 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', dim(applied) <- dim(exp_cor) } return(applied) -} \ No newline at end of file +} diff --git a/man/CST_QuantileMapping.Rd b/man/CST_QuantileMapping.Rd index 0fdd6e0c..c0b7af16 100644 --- a/man/CST_QuantileMapping.Rd +++ b/man/CST_QuantileMapping.Rd @@ -41,23 +41,24 @@ used.} 'DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used.} -\item{na.rm}{A logical value indicating ig missing values should be removed +\item{na.rm}{A logical value indicating if missing values should be removed (FALSE by default).} \item{ncores}{An integer indicating the number of cores for parallel computation using multiApply function. The default value is NULL (1).} \item{...}{Additional parameters to be used by the method choosen. See qmap -package.} +package for details.} } \value{ -An oject of class \code{s2dv_cube} containing the experimental data -after applyingthe quantile mapping correction. ) <- c(dataset = 1, member = -10, sdate = 20, ftime = 60 , +An object of class \code{s2dv_cube} containing the experimental data +after applying the quantile mapping correction. } \description{ -This function is a wrapper from fitQmap and doQmap from package -'qmap'to be applied in CSTools objects of class 's2dv_cube'. +This function is a wrapper of fitQmap and doQmap from package +'qmap' to be applied on the object of class 's2dv_cube'. The quantile mapping +adjustment between an experiment, typically a hindcast, and observation is +applied to the experiment itself or to a provided forecast. } \examples{ exp <- NULL diff --git a/man/QuantileMapping.Rd b/man/QuantileMapping.Rd index 41a87bc3..197a5c4a 100644 --- a/man/QuantileMapping.Rd +++ b/man/QuantileMapping.Rd @@ -18,13 +18,15 @@ QuantileMapping( ) } \arguments{ -\item{exp}{A multidimensional array with named dimensions.} +\item{exp}{A multidimensional array with named dimensions containing the +hindcast.} -\item{obs}{A multidimensional array with named dimensions.} +\item{obs}{A multidimensional array with named dimensions containing the +reference dataset.} \item{exp_cor}{A multidimensional array with named dimensions in which the quantile mapping correction should be applied. If it is not specified, the -correction is applied in object 'exp'.} +correction is applied on object 'exp'.} \item{sdate_dim}{A character string indicating the dimension name in which cross-validation would be applied when exp_cor is not provided. 'sdate' by @@ -42,22 +44,24 @@ used.} 'DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping 'QUANT' is used.} -\item{na.rm}{A logical value indicating ig missing values should be removed +\item{na.rm}{A logical value indicating if missing values should be removed (FALSE by default).} \item{ncores}{An integer indicating the number of cores for parallel computation using multiApply function. The default value is NULL (1).} \item{...}{Additional parameters to be used by the method choosen. See qmap -package.} +package for details.} } \value{ An array containing the experimental data after applying the quantile mapping correction. } \description{ -This function is a wrapper from fitQmap and doQmap from package -'qmap'to be applied in CSTools objects of class 's2dv_cube'. +This function is a wrapper of fitQmap and doQmap from package +'qmap' to be applied on multi-dimensional arrays. The quantile mapping +adjustment between an experiment, typically a hindcast, and observation is +applied to the experiment itself or to a provided forecast. } \examples{ exp <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) diff --git a/tests/testthat/test-CST_QuantileMapping.R b/tests/testthat/test-CST_QuantileMapping.R index 7c4a2481..a596dd16 100644 --- a/tests/testthat/test-CST_QuantileMapping.R +++ b/tests/testthat/test-CST_QuantileMapping.R @@ -50,22 +50,22 @@ obs3_2$data[1] <- NA res3_2$data[1] <- 0 # dat4 -exp4 <- lonlat_data$exp +exp4 <- lonlat_temp$exp exp4$data <- exp4$data[,,1:4,,,] dim(exp4$data) <- c(dataset = 1, member = 15, sdate = 4, ftime = 3, lat = 22, lon = 53) -obs4 <- lonlat_data$obs +obs4 <- lonlat_temp$obs obs4$data <- obs4$data[,,1:4, ,,] dim(obs4$data) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3, lat = 22, lon = 53) -exp_cor4 <- lonlat_data$exp +exp_cor4 <- lonlat_temp$exp exp_cor4$data <- exp_cor4$data[,,5:6,,,] dim(exp_cor4$data) <- c(dataset = 1, member = 15, sdate = 2, ftime = 3, lat = 22, lon = 53) # dat5 -exp5 <- lonlat_data$exp -obs5 <- lonlat_data$obs +exp5 <- lonlat_temp$exp +obs5 <- lonlat_temp$obs set.seed(1) res5 <- NULL res5$data <- array(rnorm(length(exp5)), dim = c(member = 15, sdate = 6, @@ -78,8 +78,8 @@ res5_1$data <- array(rnorm(length(res5_1)), dim = c(member = 15, ftime = 3, class(res5_1) <- "s2dv_cube" # dat6 -exp6 <- lonlat_data$exp -obs6 <- lonlat_data$obs +exp6 <- lonlat_temp$exp +obs6 <- lonlat_temp$obs obs6$data <- s2dv::InsertDim(obs6$data, pos = 1, len = 4, name = 'window') obs6_1 <- obs6 @@ -131,21 +131,21 @@ test_that("1. Sanity checks", { CST_QuantileMapping(exp = exp1, obs = obs1, exp_cor = exp_cor1, sdate_dim = 'time'), "Parameter 'exp_cor' must have dimension names." ) - # sdate_dim and member_dim + # sdate_dim and memb_dim expect_error( CST_QuantileMapping(exp = exp1, obs = obs1, sdate_dim = 'time'), - paste0("Parameter 'exp' requires 'sdate_dim' and 'member_dim' dimensions.") + paste0("Parameter 'memb_dim' is not found in 'exp' dimensions.") ) # window_dim expect_error( CST_QuantileMapping(exp = exp2, obs = obs2, window_dim = 'ftime'), - paste0("Dimension 'window_dim' not found in 'obs'.") + paste0("Parameter 'window_dim' is not found in 'obs'.") ) # method expect_error( CST_QuantileMapping(exp = exp2, obs = obs2, method = 'x'), paste0("Parameter 'method' must be one of the following methods: ", - "'PTF','DIST','RQUANT','QUANT','SSPLIN'.") + "'PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN'.") ) }) @@ -199,13 +199,17 @@ test_that("4. dat6", { CST_QuantileMapping(exp6, obs6, window_dim = 'window', ncores = 6, na.rm = TRUE) ) expect_equal( - dim(CST_QuantileMapping(exp6, obs6_1, window_dim = 'window', ncores = 6)), - dim(CST_QuantileMapping(exp6, obs6_1, window_dim = 'window', ncores = 6, na.rm = TRUE)) + dim(CST_QuantileMapping(exp6, obs6_1, window_dim = 'window', ncores = 6)$data), + c(member = 15, sdate = 6, dataset = 1, ftime = 3, lat = 22, lon = 53) ) expect_equal( sum(is.na(CST_QuantileMapping(exp6_1, obs6_1, exp_cor = exp_cor6_1, window_dim = 'window', ncores = 6, na.rm = TRUE)$data)), sum(is.na(exp_cor6_1$data)) ) + expect_equal( + dim(CST_QuantileMapping(exp6, obs6_1, exp_cor6_1, window_dim = 'window', ncores = 6, na.rm = T)$data), + c(member = 15, sdate = 1, dataset = 1, ftime = 3, lat = 22, lon = 53) + ) expect_equal( dim(CST_QuantileMapping(exp6, obs6_1, exp_cor6_2, window_dim = 'window', ncores = 6)$data), c(member = 5, sdate = 6, dataset = 1, ftime = 3, lat = 22, lon = 53) -- GitLab From 89f5e6f46779fefc929336d3311cc36c02c9062b Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 19 Oct 2022 09:52:35 +0200 Subject: [PATCH 12/14] clean up code, correct if condition, change data name back to lonlat_data --- R/CST_QuantileMapping.R | 43 ++++++----------------- tests/testthat/test-CST_QuantileMapping.R | 14 ++++---- 2 files changed, 17 insertions(+), 40 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 0a89692d..65a3f1c7 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -157,38 +157,33 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', stop("Parameter 'obs' must have dimension names.") } # sdate_dim -#-------NEW------------ - if (!is.null(exp_cor)) { - if (!is.character(sdate_dim) | length(sdate_dim) != 1) { - stop("Parameter 'sdate_dim' must be a character string.") - } - if (!sdate_dim %in% expdims | !sdate_dim %in% obsdims) { - stop("Parameter 'sdate_dim' is not found in 'exp' or 'obs' dimension.") - } - if (dim(exp)[sdate_dim] == 1 || dim(obs)[sdate_dim] == 1) { - stop("Parameter 'exp' and 'obs' must have dimension length of 'sdate_dim' bigger than 1.") - } + if (is.null(exp_cor)) { + if (!is.character(sdate_dim) | length(sdate_dim) != 1) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% expdims | !sdate_dim %in% obsdims) { + stop("Parameter 'sdate_dim' is not found in 'exp' or 'obs' dimension.") + } + if (dim(exp)[sdate_dim] == 1 || dim(obs)[sdate_dim] == 1) { + stop("Parameter 'exp' and 'obs' must have dimension length of 'sdate_dim' bigger than 1.") + } } -#--------NEW_END---------- # exp_cor if (!is.null(exp_cor)) { if (is.null(names(dim(exp_cor)))) { stop("Parameter 'exp_cor' must have dimension names.") } } -#--------NEW----------- # method if (!(method %in% c('PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN')) | length(method) != 1) { stop("Parameter 'method' must be one of the following methods: ", "'PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN'.") } -#--------NEW_END----------- # memb_dim if (!all(memb_dim %in% obsdims)) { obs <- InsertDim(obs, posdim = 1, lendim = 1, name = memb_dim[!(memb_dim %in% obsdims)]) } -#--------NEW----------- if (any(!memb_dim %in% expdims)) { stop("Parameter 'memb_dim' is not found in 'exp' dimensions.") } @@ -204,24 +199,6 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', warning("Parameter 'window_dim' is found in exp and is merged to 'memb_dim'.") } } -# # sdate_dim and member_dim -# sample_dims <- c(memb_dim, sdate_dim) -# if (!all(sample_dims %in% expdims)) { -# stop("Parameter 'exp' requires 'sdate_dim' and 'memb_dim' dimensions.") -# } -# # method -# if (!is.character(method)) { -# warning("Parameter 'method' must be a character string indicating ", -# "one of the following methods: 'PTF', 'DIST', 'RQUANT', -# 'QUANT', 'SSPLIN'. Method 'QUANT' is being used.") -# method <- 'QUANT' -# } -# if (length(method) > 1) { -# warning("Parameter 'method' has length > 1 and only the first element", -# " is used.") -# method <- method[1] -# } -#-------NEW_END------------- # na.rm if (!is.logical(na.rm) | length(na.rm) > 1) { stop("Parameter 'na.rm' must be one logical value.") diff --git a/tests/testthat/test-CST_QuantileMapping.R b/tests/testthat/test-CST_QuantileMapping.R index a596dd16..23970cf1 100644 --- a/tests/testthat/test-CST_QuantileMapping.R +++ b/tests/testthat/test-CST_QuantileMapping.R @@ -50,22 +50,22 @@ obs3_2$data[1] <- NA res3_2$data[1] <- 0 # dat4 -exp4 <- lonlat_temp$exp +exp4 <- lonlat_data$exp exp4$data <- exp4$data[,,1:4,,,] dim(exp4$data) <- c(dataset = 1, member = 15, sdate = 4, ftime = 3, lat = 22, lon = 53) -obs4 <- lonlat_temp$obs +obs4 <- lonlat_data$obs obs4$data <- obs4$data[,,1:4, ,,] dim(obs4$data) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3, lat = 22, lon = 53) -exp_cor4 <- lonlat_temp$exp +exp_cor4 <- lonlat_data$exp exp_cor4$data <- exp_cor4$data[,,5:6,,,] dim(exp_cor4$data) <- c(dataset = 1, member = 15, sdate = 2, ftime = 3, lat = 22, lon = 53) # dat5 -exp5 <- lonlat_temp$exp -obs5 <- lonlat_temp$obs +exp5 <- lonlat_data$exp +obs5 <- lonlat_data$obs set.seed(1) res5 <- NULL res5$data <- array(rnorm(length(exp5)), dim = c(member = 15, sdate = 6, @@ -78,8 +78,8 @@ res5_1$data <- array(rnorm(length(res5_1)), dim = c(member = 15, ftime = 3, class(res5_1) <- "s2dv_cube" # dat6 -exp6 <- lonlat_temp$exp -obs6 <- lonlat_temp$obs +exp6 <- lonlat_data$exp +obs6 <- lonlat_data$obs obs6$data <- s2dv::InsertDim(obs6$data, pos = 1, len = 4, name = 'window') obs6_1 <- obs6 -- GitLab From 31d7f0cd85200da5c10028311ff8617a16054dc1 Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Wed, 19 Oct 2022 12:01:19 +0200 Subject: [PATCH 13/14] Add check and correct output metadata --- R/CST_QuantileMapping.R | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 65a3f1c7..17569860 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -73,10 +73,14 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', na.rm = na.rm, ncores = ncores, ...) if (is.null(exp_cor)) { exp$data <- QMapped + exp$Datasets <- c(exp$Datasets, obs$Datasets) + exp$source_files <- c(exp$source_files, obs$source_files) return(exp) } else { exp_cor$data <- QMapped + 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) } @@ -157,22 +161,23 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', stop("Parameter 'obs' must have dimension names.") } # sdate_dim - if (is.null(exp_cor)) { - if (!is.character(sdate_dim) | length(sdate_dim) != 1) { - stop("Parameter 'sdate_dim' must be a character string.") - } - if (!sdate_dim %in% expdims | !sdate_dim %in% obsdims) { - stop("Parameter 'sdate_dim' is not found in 'exp' or 'obs' dimension.") - } - if (dim(exp)[sdate_dim] == 1 || dim(obs)[sdate_dim] == 1) { - stop("Parameter 'exp' and 'obs' must have dimension length of 'sdate_dim' bigger than 1.") - } + if (!is.character(sdate_dim) | length(sdate_dim) != 1) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% expdims | !sdate_dim %in% obsdims) { + stop("Parameter 'sdate_dim' is not found in 'exp' or 'obs' dimension.") + } + if (dim(exp)[sdate_dim] == 1 || dim(obs)[sdate_dim] == 1) { + stop("Parameter 'exp' and 'obs' must have dimension length of 'sdate_dim' bigger than 1.") } # exp_cor if (!is.null(exp_cor)) { if (is.null(names(dim(exp_cor)))) { stop("Parameter 'exp_cor' must have dimension names.") } + if (!sdate_dim %in% names(dim(exp_cor))) { + stop("Parameter 'sdate_dim' is not found in 'exp_cor' dimension.") + } } # method if (!(method %in% c('PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN')) | length(method) != 1) { -- GitLab From 1aec981578e483b3f4a411ff9c0fe407fc46f62a Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Wed, 19 Oct 2022 12:35:55 +0200 Subject: [PATCH 14/14] tests/testthat/test-CST_QuantileMapping.R --- tests/testthat/test-CST_QuantileMapping.R | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/testthat/test-CST_QuantileMapping.R b/tests/testthat/test-CST_QuantileMapping.R index 23970cf1..a80babe9 100644 --- a/tests/testthat/test-CST_QuantileMapping.R +++ b/tests/testthat/test-CST_QuantileMapping.R @@ -178,7 +178,7 @@ test_that("2. dat2, dat3 and dat4", { test_that("3. dat5", { expect_equal( - dim(CST_QuantileMapping(exp5, obs5, ncores = 6)$data), + dim(CST_QuantileMapping(exp5, obs5)$data), dim(res5$data) ) expect_equal( @@ -186,7 +186,7 @@ test_that("3. dat5", { dim(res5$data) ) expect_equal( - dim(CST_QuantileMapping(exp5, obs5, sdate_dim = "ftime", ncores = 6)$data), + dim(CST_QuantileMapping(exp5, obs5, sdate_dim = "ftime")$data), dim(res5_1$data) ) }) @@ -195,23 +195,23 @@ test_that("3. dat5", { test_that("4. dat6", { expect_equal( - CST_QuantileMapping(exp6, obs6, window_dim = 'window', ncores = 6), - CST_QuantileMapping(exp6, obs6, window_dim = 'window', ncores = 6, na.rm = TRUE) + CST_QuantileMapping(exp6, obs6, window_dim = 'window'), + CST_QuantileMapping(exp6, obs6, window_dim = 'window', na.rm = TRUE) ) expect_equal( - dim(CST_QuantileMapping(exp6, obs6_1, window_dim = 'window', ncores = 6)$data), + dim(CST_QuantileMapping(exp6, obs6_1, window_dim = 'window')$data), c(member = 15, sdate = 6, dataset = 1, ftime = 3, lat = 22, lon = 53) ) expect_equal( - sum(is.na(CST_QuantileMapping(exp6_1, obs6_1, exp_cor = exp_cor6_1, window_dim = 'window', ncores = 6, na.rm = TRUE)$data)), + sum(is.na(CST_QuantileMapping(exp6_1, obs6_1, exp_cor = exp_cor6_1, window_dim = 'window', na.rm = TRUE)$data)), sum(is.na(exp_cor6_1$data)) ) expect_equal( - dim(CST_QuantileMapping(exp6, obs6_1, exp_cor6_1, window_dim = 'window', ncores = 6, na.rm = T)$data), + dim(CST_QuantileMapping(exp6, obs6_1, exp_cor6_1, window_dim = 'window', na.rm = T)$data), c(member = 15, sdate = 1, dataset = 1, ftime = 3, lat = 22, lon = 53) ) expect_equal( - dim(CST_QuantileMapping(exp6, obs6_1, exp_cor6_2, window_dim = 'window', ncores = 6)$data), + dim(CST_QuantileMapping(exp6, obs6_1, exp_cor6_2, window_dim = 'window')$data), c(member = 5, sdate = 6, dataset = 1, ftime = 3, lat = 22, lon = 53) ) }) -- GitLab