From 07f3422a67006d9a9dbd64607fa76dec5f79fcf0 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 22 Sep 2020 15:41:18 +0200 Subject: [PATCH 01/13] testing qmap library on Depends --- DESCRIPTION | 4 +-- NAMESPACE | 1 + NEWS.md | 9 +++++++ R/CST_QuantileMapping.R | 33 +++++++++++++++++++++++ man/QuantileMapping.Rd | 60 +++++++++++++++++++++++++++++++++++++++++ 5 files changed, 105 insertions(+), 2 deletions(-) create mode 100644 man/QuantileMapping.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 07433f27..5050cd8c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -47,13 +47,13 @@ Description: Exploits dynamical seasonal forecasts in order to provide Yiou et al. (2013) . Depends: R (>= 3.4.0), - maps + maps, + qmap Imports: s2dverification, s2dv, rainfarmr, multiApply (>= 2.1.1), - qmap, ClimProjDiags, ncdf4, plyr, diff --git a/NAMESPACE b/NAMESPACE index 9f9a0e34..ba7186fe 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -33,6 +33,7 @@ export(PlotForecastPDF) export(PlotMostLikelyQuantileMap) export(PlotPDFsOLE) export(PlotTriangles4Categories) +export(QuantileMapping) export(RFSlope) export(RFTemp) export(RainFARM) diff --git a/NEWS.md b/NEWS.md index fab70c70..5f1c7ea1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,12 @@ +### CSTools 3.x.x +**Submission date to CRAN: XX-XX-2020** + +- New features: + + QuantileMapping is exposed to users + +- New fixes: + + qmap library moved from Imports to Depends + ### CSTools 3.1.0 **Submission date to CRAN: XX-06-2020** diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 3923569a..c0b9f8a6 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -91,6 +91,39 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, 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, ...) { obsdims <- names(dim(obs)) diff --git a/man/QuantileMapping.Rd b/man/QuantileMapping.Rd new file mode 100644 index 00000000..8771e549 --- /dev/null +++ b/man/QuantileMapping.Rd @@ -0,0 +1,60 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_QuantileMapping.R +\name{QuantileMapping} +\alias{QuantileMapping} +\title{Quantiles Mapping for seasonal or decadal forecast data} +\usage{ +QuantileMapping( + exp, + obs, + exp_cor = NULL, + sample_dims = "ftime", + sample_length = NULL, + method = "QUANT", + ncores = NULL, + ... +) +} +\arguments{ +\item{exp}{a multi-dimensional array with named dimensions containing the hindcast.} + +\item{obs}{a multi-dimensional array with named dimensions (the same as the provided in 'exp') containing the reference dataset.} + +\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{sample_dims}{a character vector indicating the dimensions that can be used as sample for the same distribution} + +\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{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{...}{additional arguments passed to the method specified by \code{method}.} +} +\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.}} +} +\seealso{ +\code{qmap::fitQmap} and \code{qmap::doQmap} +} +\author{ +Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +} -- GitLab From 6015933f083180c20ad24dc85d52dbc5aace9a94 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 13 Oct 2020 16:02:23 +0200 Subject: [PATCH 02/13] Check common exp_cor dimensions --- R/CST_QuantileMapping.R | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index c0b9f8a6..5d52d1f3 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -167,6 +167,24 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sample_dims = 'ftime', obs <- adrop(obs, drop = todrop) } } + 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, exp_cordims) + if (all(dim(exp_cor)[todroppos] != 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 = todroppos) + } + } + } + if (!all(sample_dims %in% obsdims)) { newobsdims <- sample_dims[!sample_dims %in% obsdims] -- GitLab From 4e0096226099681925d32266f1ef4272ffe5d6de Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 13 Oct 2020 18:24:28 +0200 Subject: [PATCH 03/13] exp_cor can have different length of sdates --- R/CST_QuantileMapping.R | 56 +++++++++++++++++++++++++---------------- 1 file changed, 35 insertions(+), 21 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 5d52d1f3..177574d2 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -79,13 +79,9 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, 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, ...) - pos <- match(dimnames, names(dim(QMapped))) - QMapped <- aperm(QMapped, pos) - names(dim(QMapped)) <- dimnames exp$data <- QMapped exp$Datasets <- c(exp$Datasets, obs$Datasets) exp$source_files <- c(exp$source_files, obs$source_files) @@ -167,30 +163,40 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sample_dims = 'ftime', 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))) { + 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", todrop, "wont be 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)] - todroppos <- match(todrop, exp_cordims) - if (all(dim(exp_cor)[todroppos] != 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 = todroppos) - } + 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 (!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(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.") @@ -225,8 +231,16 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sample_dims = 'ftime', 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) - dim(qmaped) <- dim(exp) + names(dim(qmaped)) <- out_names return(qmaped) } qmapcor <- function(exp, obs, exp_cor = NULL, sample_length = NULL, method = 'QUANT', -- GitLab From 9a4108bd9ced86e9e4d68d81975c078af435666f Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 13 Oct 2020 18:33:03 +0200 Subject: [PATCH 04/13] Fixing warning message --- R/CST_QuantileMapping.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 177574d2..9400fb80 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -177,7 +177,8 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sample_dims = 'ftime', todrop <- commondims[(commondims %in% sample_dims)] todroppos <- match(todrop, sample_dims) if (all(dim(exp_cor)[todrop] != 1)) { - warning(paste("The sample_dims", todrop, "wont be used when applying the", + 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 { -- GitLab From a8a96d6284c1db10b26600fdbc25d9b7eaa61227 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 13 Oct 2020 18:37:53 +0200 Subject: [PATCH 05/13] Add examples --- R/CST_QuantileMapping.R | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 9400fb80..fc0642e1 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -48,6 +48,14 @@ #'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] @@ -60,6 +68,27 @@ #'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')) #'} #'@export CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, -- GitLab From 750b8b6a47eba6764412813f5f67959a412a444e Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 13 Oct 2020 19:00:52 +0200 Subject: [PATCH 06/13] Fix tests and add details to documentation --- R/CST_QuantileMapping.R | 4 +++- tests/testthat/test-CST_QuantileMapping.R | 24 +++++++++++++++++++---- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index fc0642e1..64529460 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -23,6 +23,8 @@ #'\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 @@ -217,7 +219,7 @@ QuantileMapping <- function(exp, obs, exp_cor = NULL, sample_dims = 'ftime', 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", + 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.") diff --git a/tests/testthat/test-CST_QuantileMapping.R b/tests/testthat/test-CST_QuantileMapping.R index a03680f7..f8482967 100644 --- a/tests/testthat/test-CST_QuantileMapping.R +++ b/tests/testthat/test-CST_QuantileMapping.R @@ -1,6 +1,5 @@ context("Generic tests") test_that("Sanity checks", { -library(qmap) expect_error( CST_QuantileMapping(exp = 1), paste0("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", @@ -90,8 +89,25 @@ library(qmap) expect_equal(CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time', method = 'SSPLIN'), exp) library(CSTools) - expect_equal(CST_QuantileMapping(exp = lonlat_data$exp, obs = lonlat_data$obs), - CST_QuantileMapping(exp = lonlat_data$exp, obs = lonlat_data$obs, - exp_cor = lonlat_data$exp)) + 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$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[,,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'") }) -- GitLab From b00f743a04d4da61b8a59119cb9eda83734bf542 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 2 Nov 2020 11:46:08 +0100 Subject: [PATCH 07/13] If exp_cor, the output adds its information --- R/CST_QuantileMapping.R | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 64529460..3662afaa 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -113,9 +113,16 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, 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, ...) - exp$data <- QMapped - exp$Datasets <- c(exp$Datasets, obs$Datasets) - exp$source_files <- c(exp$source_files, obs$source_files) + 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) + } else { + exp_cor$data <- QMapped + exp_cor$Datsets <- c(exp$Datasets, obs$Datasets, exp_cor$Datasets) + exp_cor$source_files <- c(exp$source_files, obs$source_files, exp_cor$source_files) + exp <- exp_cor + } return(exp) } #'Quantiles Mapping for seasonal or decadal forecast data -- GitLab From 8477cba9970a36e741303e3870a149c1a729d8dc Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 2 Nov 2020 12:16:12 +0100 Subject: [PATCH 08/13] Fix rounding number in MultiMetric test --- tests/testthat/test-CST_MultiMetric.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/testthat/test-CST_MultiMetric.R b/tests/testthat/test-CST_MultiMetric.R index b08b68ba..cabb2a55 100644 --- a/tests/testthat/test-CST_MultiMetric.R +++ b/tests/testthat/test-CST_MultiMetric.R @@ -11,8 +11,7 @@ test_that("basic use case", { attr(exp, 'class') <- 's2dv_cube' attr(obs, 'class') <- 's2dv_cube' - result <- list(data = array(rep(c(rep(1, 9), -rep(0.89999999999999991118215802998747676610950, 3)), 48), + result <- list(data = array(rep(c(rep(1, 9), rep(0.9, 3)), 48), dim = c(dataset = 3, dataset = 1, statistics = 4, lat = 6, lon = 8)), lat = lat, lon = lon) -- GitLab From e9efbf030f0352f53f515d0cddcc6adc37da47e0 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 2 Nov 2020 12:55:39 +0100 Subject: [PATCH 09/13] p-value computation has changed in s2dverification --- tests/testthat/test-CST_MultiMetric.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-CST_MultiMetric.R b/tests/testthat/test-CST_MultiMetric.R index cabb2a55..1a41a83a 100644 --- a/tests/testthat/test-CST_MultiMetric.R +++ b/tests/testthat/test-CST_MultiMetric.R @@ -11,7 +11,7 @@ test_that("basic use case", { attr(exp, 'class') <- 's2dv_cube' attr(obs, 'class') <- 's2dv_cube' - result <- list(data = array(rep(c(rep(1, 9), rep(0.9, 3)), 48), + result <- list(data = array(rep(c(rep(1, 9), rep(0, 3)), 48), dim = c(dataset = 3, dataset = 1, statistics = 4, lat = 6, lon = 8)), lat = lat, lon = lon) -- GitLab From 89b3edfc5904d9748ca637f9b9f951f298bd09e7 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 2 Nov 2020 20:37:06 +0100 Subject: [PATCH 10/13] only 1 dataset to correct save output with SaveExp --- R/CST_QuantileMapping.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/R/CST_QuantileMapping.R b/R/CST_QuantileMapping.R index 3662afaa..b27cbd4c 100644 --- a/R/CST_QuantileMapping.R +++ b/R/CST_QuantileMapping.R @@ -115,11 +115,9 @@ CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, method = method, 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) } else { exp_cor$data <- QMapped - exp_cor$Datsets <- c(exp$Datasets, obs$Datasets, exp_cor$Datasets) exp_cor$source_files <- c(exp$source_files, obs$source_files, exp_cor$source_files) exp <- exp_cor } -- GitLab From 61ba7f405890ea5c72452e488bf9dbd6e4bfdd91 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 4 Nov 2020 10:40:46 +0100 Subject: [PATCH 11/13] Fix initialization dates --- R/CST_SaveExp.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/CST_SaveExp.R b/R/CST_SaveExp.R index 1bdd4779..1b1c0202 100644 --- a/R/CST_SaveExp.R +++ b/R/CST_SaveExp.R @@ -40,7 +40,7 @@ CST_SaveExp <- function(data, destination = "./CST_Data") { "as output by CSTools::CST_Load.") } sdates <- lapply(1:length(data$Datasets), function(x) { - data$Datasets[[x]]$InitializationDates[[1]]})[[1]] + unique(data$Datasets[[x]]$InitializationDates)})[[1]] if (!is.character(attributes(data$Variable)$units)) { units <- attributes(data$Variable)$variable$units } else { -- GitLab From f8870f8300fdf23d38c943d8796e347fa4f6142c Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 19 Nov 2020 21:56:18 +0100 Subject: [PATCH 12/13] replace NAs without lapply --- R/CST_RainFARM.R | 47 ++++++++++---------------------------- man/CST_QuantileMapping.Rd | 30 ++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 35 deletions(-) diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R index 73b075c4..0c94650f 100644 --- a/R/CST_RainFARM.R +++ b/R/CST_RainFARM.R @@ -383,42 +383,19 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1, # restoring NA values in their position: if (!is.null(posna)) { pos <- which(posna == FALSE) - if (pos[1] == 1) { - result_dims['rainfarm_samples'] <- 1 - x <- array(rep(NA,prod(result_dims)), result_dims) - r <- abind(x, r, along = 3) - names(dim(r)) <- names(result_dims) - pos <- pos[-1] - } - result_dims['rainfarm_samples'] <- 1 - if ((pos[1] - 1) < dim(r)['rainfarm_samples']) { - x <- max(which((pos - 1) < dim(r)['rainfarm_samples'])) - pos_internal <- pos[1:x] - if (pos[length(pos)] > pos[x]) { - pos <- pos[(x + 1) : length(pos)] - } else { - pos <- NULL - } - } else { - pos_internal <- NULL - } - if (!is.null(pos_internal)) { - r <- lapply(pos_internal, - function(x) { - rrr <- Subset(r, along = 'rainfarm_samples', indices = 1:(x-1)) - rrrr <- Subset(r, along = 'rainfarm_samples', indices = x:dim(r)['rainfarm_samples']) - r <- abind(rrr, array(NA, result_dims), rrrr, along = 3) - names(dim(r)) <- names(result_dims) - return(r) - })[[length(pos)]] - } - if (!is.null(pos)) { - result_dims['rainfarm_samples'] <- length(pos) - x <- array(rep(NA,prod(result_dims)), result_dims) - r <- abind(r, x, along = 3) - } - names(dim(r)) <- names(result_dims) + dimdata <- dim(r) + xdim <- which(names(dimdata) == 'rainfarm_samples') + dimdata[xdim] <- dimdata[xdim] + length(pos) + new <- array(NA, dimdata) + posT <- which(posna == TRUE) + i = 1 + invisible(lapply(posT, function(x) { + new[,,x,] <<- r[,,i,] + i <<- i + 1 + })) + #names(dim(r)) <- names(result_dims) warning("Missing values found in the samples.") + r <- new } return(r) } diff --git a/man/CST_QuantileMapping.Rd b/man/CST_QuantileMapping.Rd index e78c8d56..ec5fc8a3 100644 --- a/man/CST_QuantileMapping.Rd +++ b/man/CST_QuantileMapping.Rd @@ -51,6 +51,7 @@ 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. } \examples{ library(qmap) @@ -70,6 +71,14 @@ 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] @@ -82,6 +91,27 @@ 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')) } } \seealso{ -- GitLab From 5271b755acf6509b03abb3757aa9324687341dcd Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 20 Nov 2020 18:41:12 +0100 Subject: [PATCH 13/13] CST_SplitDim includes parameter insert_ftime --- R/CST_SplitDim.R | 39 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/R/CST_SplitDim.R b/R/CST_SplitDim.R index 3326d435..dc1718ea 100644 --- a/R/CST_SplitDim.R +++ b/R/CST_SplitDim.R @@ -35,11 +35,44 @@ #'dim(new_data$data) #'@export CST_SplitDim <- function(data, split_dim = 'time', indices = NULL, - freq = 'monthly', new_dim_name = NULL) { + freq = 'monthly', new_dim_name = NULL, insert_ftime = NULL) { if (!inherits(data, 's2dv_cube')) { - stop("Parameter 'data' must be of the class 's2dv_cube', ", - "as output by CSTools::CST_Load.") + stop("Parameter 'data' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") } + if (!is.null(insert_ftime)) { + if (!is.numeric(insert_ftime)) { + stop("Parameter 'insert_ftime' should be an integer.") + } else { + if (length(insert_ftime) > 1) { + warning("Parameter 'insert_ftime' must be of length 1, and only the", + " first element will be used.") + insert_ftime <- insert_ftime[1] + } + # adding NAs at the begining of the data in ftime dim + ftimedim <- which(names(dim(data$data)) == 'ftime') + dims <- dim(data$data) + dims[ftimedim] <- insert_ftime + empty_array <- array(NA, dims) + data$data <- abind(empty_array, data$data, along = ftimedim) + names(dim(data$data)) <- names(dims) + # adding dates to Dates for the new NAs introduced + if ((data$Dates[[1]][2] - data$Dates[[1]][1]) == 1) { + timefreq <- 'days' + } else { + timefreq <- 'months' + warning("Time frequency of forecast time is considered monthly.") + } + start <- data$Dates[[1]] + dim(start) <- c(ftime = length(start)/dims['sdate'], sdate = dims['sdate']) + #new <- array(NA, prod(dim(data$data)[c('ftime', 'sdate')])) + # Pending fix transform to UTC when concatenaiting + data$Dates$start <- do.call(c, lapply(1:dim(start)[2], function(x) { + seq(start[1,x] - as.difftime(insert_ftime, + unit = timefreq), + start[dim(start)[1],x], by = timefreq, tz = "UTC")})) + } + } if (is.null(indices)) { if (any(split_dim %in% c('ftime', 'time', 'sdate'))) { if (is.list(data$Dates)) { -- GitLab