From 7df48b25d46c9a839bf8921b687861e5a26b81b5 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 10 Nov 2020 09:56:25 +0100 Subject: [PATCH 1/5] adding parameter cex_axis --- R/PlotTriangles4Categories.R | 6 ++++-- man/PlotTriangles4Categories.Rd | 3 +++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/R/PlotTriangles4Categories.R b/R/PlotTriangles4Categories.R index fcfb36bb..341d3a1d 100644 --- a/R/PlotTriangles4Categories.R +++ b/R/PlotTriangles4Categories.R @@ -40,6 +40,7 @@ #'@param cex_leg a number to indicate the increase/reductuion of the lab_legend used #' to represent sig_data. #'@param col_leg color of the legend (triangles). +#'@param cex_axis a number to indicate the increase/reduction of the axis labels. #'@param fileout A string of full directory path and file name indicating where #' to save the plot. If not specified (default), a graphics device will pop up. #'@param size_units A string indicating the units of the size of the device @@ -81,6 +82,7 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, ylabels = NULL, ytitle = NULL, legend = TRUE, lab_legend = NULL, cex_leg = 1, col_leg = 'black', + cex_axis = 1.5, fileout = NULL, size_units = 'px', res = 100, figure.width = 1, ...) { # Checking the dimensions @@ -202,10 +204,10 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, } if (xlab){ - axis(1, at =(1:ncol) - 0.5, las = 2, labels = xlabels, cex.axis = 1.5) + axis(1, at =(1:ncol) - 0.5, las = 2, labels = xlabels, cex.axis = cex_axis) } if (ylab){ - axis(2, at = (1:nrow) - 0.5, las = 2, labels = ylabels, cex.axis = 1.5) + axis(2, at = (1:nrow) - 0.5, las = 2, labels = ylabels, cex.axis = cex_axis) } diff --git a/man/PlotTriangles4Categories.Rd b/man/PlotTriangles4Categories.Rd index 6e95df38..9c47fc62 100644 --- a/man/PlotTriangles4Categories.Rd +++ b/man/PlotTriangles4Categories.Rd @@ -23,6 +23,7 @@ PlotTriangles4Categories( lab_legend = NULL, cex_leg = 1, col_leg = "black", + cex_axis = 1.5, fileout = NULL, size_units = "px", res = 100, @@ -84,6 +85,8 @@ to represent sig_data.} \item{col_leg}{color of the legend (triangles).} +\item{cex_axis}{a number to indicate the increase/reduction of the axis labels.} + \item{fileout}{A string of full directory path and file name indicating where to save the plot. If not specified (default), a graphics device will pop up.} -- GitLab From 5b026db737eeffa9665067d0c22cb600bee983f0 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 11 Nov 2020 13:15:18 +0100 Subject: [PATCH 2/5] PlotTriangles has parameter mar --- R/PlotTriangles4Categories.R | 9 ++++++--- man/PlotTriangles4Categories.Rd | 3 +++ 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/R/PlotTriangles4Categories.R b/R/PlotTriangles4Categories.R index 341d3a1d..f56e8bd6 100644 --- a/R/PlotTriangles4Categories.R +++ b/R/PlotTriangles4Categories.R @@ -43,6 +43,7 @@ #'@param cex_axis a number to indicate the increase/reduction of the axis labels. #'@param fileout A string of full directory path and file name indicating where #' to save the plot. If not specified (default), a graphics device will pop up. +#'@param mar A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot. #'@param size_units A string indicating the units of the size of the device #' (file or window) to plot in. Set 'px' as default. See ?Devices and the #' creator function of the corresponding device. @@ -82,7 +83,7 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, ylabels = NULL, ytitle = NULL, legend = TRUE, lab_legend = NULL, cex_leg = 1, col_leg = 'black', - cex_axis = 1.5, + cex_axis = 1.5, mar = c(5, 4, 0, 0), fileout = NULL, size_units = 'px', res = 100, figure.width = 1, ...) { # Checking the dimensions @@ -102,7 +103,9 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, stop("Parameter 'data' should contain 'dimx', 'dimy' and 'dimcat' dimension names. ") } } - + if (!is.vector(mar) & length(mar) != 4) { + stop("Parameter 'mar' must be a vector of length 4.") + } if (!is.null(sig_data)) { if (!is.logical(sig_data)) { stop("Parameter 'sig_data' array must be logical.")} @@ -181,7 +184,7 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, if(legend){ layout(matrix(c(1, 2, 1, 3), 2, 2, byrow = T), widths = c(10, 3.4), heights = c(10, 3.5)) - par(oma = c(1, 1, 1, 1), mar = c(5, 4, 0, 0)) + par(oma = c(1, 1, 1, 1), mar = mar) if(is.null(lab_legend)) { lab_legend = 1:ncat } diff --git a/man/PlotTriangles4Categories.Rd b/man/PlotTriangles4Categories.Rd index 9c47fc62..abd3085e 100644 --- a/man/PlotTriangles4Categories.Rd +++ b/man/PlotTriangles4Categories.Rd @@ -24,6 +24,7 @@ PlotTriangles4Categories( cex_leg = 1, col_leg = "black", cex_axis = 1.5, + mar = c(5, 4, 0, 0), fileout = NULL, size_units = "px", res = 100, @@ -87,6 +88,8 @@ to represent sig_data.} \item{cex_axis}{a number to indicate the increase/reduction of the axis labels.} +\item{mar}{A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.} + \item{fileout}{A string of full directory path and file name indicating where to save the plot. If not specified (default), a graphics device will pop up.} -- GitLab From 08ce8fa3a66436012f09e8245cfe67e0248319cc Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 1 Dec 2020 12:16:33 +0100 Subject: [PATCH 3/5] Update news --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 3a8dedb3..d21f29e6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,7 @@ + CST_RainFARM accepts multi-dimensional slopes and weights and handless missing values in sample dimensions. + CST_MultiMetric includes 'rpss' metric and it is addapted to s2dv. + PlotMostLikelyQuantileMap vignette + + PlotTriangles4Categories includes two parameters to adjust axis and margins. - Fixes: + PlotForecastPDF correctly displays terciles labels -- GitLab From 01bda57992f889b973bbf4da851127323421eb6b Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 1 Dec 2020 13:25:59 +0100 Subject: [PATCH 4/5] InsertDim only from s2dv --- NAMESPACE | 2 +- NEWS.md | 1 + R/CST_Anomaly.R | 11 ++++---- R/CST_Calibration.R | 14 +++++----- R/CST_CategoricalEnsCombination.R | 44 ++++++++++++++++++++++++------- R/CST_RegimesAssign.R | 20 +++++++------- R/CST_SaveExp.R | 8 +++--- R/PlotForecastPDF.R | 6 ++--- man/CST_MultiMetric.Rd | 2 +- man/CategoricalEnsCombination.Rd | 43 ++++++++++++++++++++++++++++++ 10 files changed, 109 insertions(+), 42 deletions(-) create mode 100644 man/CategoricalEnsCombination.Rd diff --git a/NAMESPACE b/NAMESPACE index 29dc0811..91e34471 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -25,6 +25,7 @@ export(CST_SaveExp) export(CST_SplitDim) export(CST_WeatherRegimes) export(Calibration) +export(CategoricalEnsCombination) export(EnsClustering) export(MergeDims) export(MultiEOF) @@ -101,7 +102,6 @@ importFrom(s2dverification,Ano_CrossValid) importFrom(s2dverification,Clim) importFrom(s2dverification,EOF) importFrom(s2dverification,Eno) -importFrom(s2dverification,InsertDim) importFrom(s2dverification,Load) importFrom(s2dverification,Mean1Dim) importFrom(s2dverification,Subset) diff --git a/NEWS.md b/NEWS.md index d21f29e6..7e06bfcd 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,7 @@ + CST_MultiMetric includes 'rpss' metric and it is addapted to s2dv. + PlotMostLikelyQuantileMap vignette + PlotTriangles4Categories includes two parameters to adjust axis and margins. + + CategoricalEnsCombination is exposed to users - Fixes: + PlotForecastPDF correctly displays terciles labels diff --git a/R/CST_Anomaly.R b/R/CST_Anomaly.R index ef616225..52a786fb 100644 --- a/R/CST_Anomaly.R +++ b/R/CST_Anomaly.R @@ -18,7 +18,8 @@ #' #' @return A list with two S3 objects, 'exp' and 'obs', of the class 's2dv_cube', containing experimental and date-corresponding observational anomalies, respectively. These 's2dv_cube's can be ingested by other functions in CSTools. #' -#'@importFrom s2dverification Clim Ano_CrossValid InsertDim +#'@importFrom s2dverification Clim Ano_CrossValid +#'@importFrom s2dv InsertDim #' #'@examples #'# Example 1: @@ -165,12 +166,12 @@ CST_Anomaly <- function(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE, clim_exp <- tmp$clim_exp clim_obs <- tmp$clim_obs } else { - clim_exp <- s2dverification::InsertDim(tmp$clim_exp, 2, dim_exp[2]) - clim_obs <- s2dverification::InsertDim(tmp$clim_obs, 2, dim_obs[2]) + clim_exp <- InsertDim(tmp$clim_exp, 2, dim_exp[2]) + clim_obs <- InsertDim(tmp$clim_obs, 2, dim_obs[2]) } - clim_exp <- s2dverification::InsertDim(clim_exp, 3, dim_exp[3]) - clim_obs <- s2dverification::InsertDim(clim_obs, 3, dim_obs[3]) + clim_exp <- InsertDim(clim_exp, 3, dim_exp[3]) + clim_obs <- InsertDim(clim_obs, 3, dim_obs[3]) ano <- NULL ano$ano_exp <- exp$data - clim_exp ano$ano_obs <- obs$data - clim_obs diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R index fd381790..874c351e 100644 --- a/R/CST_Calibration.R +++ b/R/CST_Calibration.R @@ -13,7 +13,7 @@ #'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one. #'@return an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions as the one in the exp object. #' -#'@importFrom s2dverification InsertDim +#'@importFrom s2dv InsertDim #'@import abind #' #'@seealso \code{\link{CST_Load}} @@ -84,7 +84,7 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", #'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one. #'@return an array containing the calibrated forecasts with the same dimensions as the \code{exp} array. #' -#'@importFrom s2dverification InsertDim +#'@importFrom s2dv InsertDim #'@import abind #' #'@seealso \code{\link{CST_Load}} @@ -251,7 +251,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o .calc.obs.fc.quant <- function(obs, fc){ #function to calculate different quantities of a series of ensemble forecasts and corresponding observations amt.mbr <- dim(fc)[1] - obs.per.ens <- s2dverification::InsertDim(obs, posdim = 1, lendim = amt.mbr) + obs.per.ens <- InsertDim(obs, posdim = 1, lendim = amt.mbr) fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs") obs.av <- mean(obs, na.rm = TRUE) @@ -271,7 +271,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o .calc.obs.fc.quant.ext <- function(obs, fc){ #extended function to calculate different quantities of a series of ensemble forecasts and corresponding observations amt.mbr <- dim(fc)[1] - obs.per.ens <- s2dverification::InsertDim(obs, posdim = 1, lendim = amt.mbr) + obs.per.ens <- InsertDim(obs, posdim = 1, lendim = amt.mbr) fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs") obs.av <- mean(obs, na.rm = TRUE) @@ -296,7 +296,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) fc.ens.av.av <- mean(fc.ens.av, na.rm = TRUE) fc.ens.av.sd <- sd(fc.ens.av, na.rm = TRUE) - fc.ens.av.per.ens <- s2dverification::InsertDim(fc.ens.av, posdim = 1, lendim = amt.mbr) + fc.ens.av.per.ens <- InsertDim(fc.ens.av, posdim = 1, lendim = amt.mbr) fc.ens.sd <- apply(fc, c(2), sd, na.rm = TRUE) fc.ens.var.av.sqrt <- sqrt(mean(fc.ens.sd^2, na.rm = TRUE)) fc.dev <- fc - fc.ens.av.per.ens @@ -322,10 +322,10 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o .calc.fc.quant.ext <- function(fc){ #extended function to calculate different quantities of a series of ensemble forecasts amt.mbr <- dim(fc)[1] - repmat1.tmp <- s2dverification::InsertDim(fc, posdim = 1, lendim = amt.mbr) + repmat1.tmp <- InsertDim(fc, posdim = 1, lendim = amt.mbr) repmat2.tmp <- aperm(repmat1.tmp, c(2, 1, 3)) spr.abs <- apply(abs(repmat1.tmp - repmat2.tmp), c(3), mean, na.rm = TRUE) - spr.abs.per.ens <- s2dverification::InsertDim(spr.abs, posdim = 1, lendim = amt.mbr) + spr.abs.per.ens <- InsertDim(spr.abs, posdim = 1, lendim = amt.mbr) return( append(.calc.fc.quant(fc), diff --git a/R/CST_CategoricalEnsCombination.R b/R/CST_CategoricalEnsCombination.R index 26ac5449..6664a2f1 100644 --- a/R/CST_CategoricalEnsCombination.R +++ b/R/CST_CategoricalEnsCombination.R @@ -55,7 +55,7 @@ #'@references Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744. #'@references Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290). #' -#'@importFrom s2dverification InsertDim +#'@importFrom s2dv InsertDim #'@import abind #'@examples #' @@ -86,19 +86,45 @@ CST_CategoricalEnsCombination <- function(exp, obs, cat.method = "pool", eval.me "of the parameter 'obs' must be equal to 1.") } names.dim.tmp <- names(dim(exp$data)) - exp$data <- .CategoricalEnsCombination.wrap(fc = exp$data, obs = obs$data, cat.method = cat.method, eval.method = eval.method, amt.cat = amt.cat, ...) + exp$data <- CategoricalEnsCombination(fc = exp$data, obs = obs$data, cat.method = cat.method, + eval.method = eval.method, amt.cat = amt.cat, ...) names.dim.tmp[which(names.dim.tmp == "member")] <- "category" names(dim(exp$data)) <- names.dim.tmp - exp$data <- s2dverification::InsertDim(exp$data, lendim = 1, posdim = 2) + exp$data <- InsertDim(exp$data, lendim = 1, posdim = 2) names(dim(exp$data))[2] <- "member" exp$Datasets <- c(exp$Datasets, obs$Datasets) exp$source_files <- c(exp$source_files, obs$source_files) return(exp) } +#' Make categorical forecast based on a multi-model forecast with potential for calibrate +#' +#'@author Bert Van Schaeybroeck, \email{bertvs@meteo.be} +#'@description This function converts a multi-model ensemble forecast +#' into a categorical forecast by giving the probability +#' for each category. Different methods are available to combine +#' the different ensemble forecasting models into +#' probabilistic categorical forecasts. +#' +#' See details in ?CST_CategoricalEnsCombination +#'@param fc a multi-dimensional array with named dimensions containing the seasonal forecast experiment data in the element named \code{$data}. The amount of forecasting models is equal to the size of the \code{dataset} dimension of the data array. The amount of members per model may be different. The size of the \code{member} dimension of the data array is equal to the maximum of the ensemble members among the models. Models with smaller ensemble sizes have residual indices of \code{member} dimension in the data array filled with NA values. +#'@param obs a multidimensional array with named dimensions containing the observed data in the element named \code{$data}. +#'@param amt.cat is the amount of categories. Equally-sized quantiles will be calculated based on the amount of categories. +#'@param cat.method method used to produce the categorical forecast, can be either \code{pool}, \code{comb}, \code{mmw} or \code{obs}. The method pool assumes equal weight for all ensemble members while the method comb assumes equal weight for each model. The weighting method is descirbed in Rajagopalan et al. (2002), Robertson et al. (2004) and Van Schaeybroeck and Vannitsem (2019). Finally, the \code{obs} method classifies the observations into the different categories and therefore contains only 0 and 1 values. +#'@param eval.method is the sampling method used, can be either \code{"in-sample"} or \code{"leave-one-out"}. Default value is the \code{"leave-one-out"} cross validation. +#'@param ... other parameters to be passed on to the calibration procedure. +#' +#'@return an array containing the categorical forecasts in the element called \code{$data}. The first two dimensions of the returned object are named dataset and member and are both of size one. An additional dimension named category is introduced and is of size amt.cat. +#' +#'@references Rajagopalan, B., Lall, U., & Zebiak, S. E. (2002). Categorical climate forecasts through regularization and optimal combination of multiple GCM ensembles. Monthly Weather Review, 130(7), 1792-1811. +#'@references Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744. +#'@references Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290). +#' +#'@importFrom s2dv InsertDim +#'@import abind +#'@export - -.CategoricalEnsCombination.wrap <- function (fc, obs, cat.method, eval.method, amt.cat, ...) { +CategoricalEnsCombination <- function (fc, obs, cat.method, eval.method, amt.cat, ...) { if (!all(c("member", "sdate") %in% names(dim(fc)))) { stop("Parameter 'exp' must have the dimensions 'member' and 'sdate'.") @@ -277,8 +303,8 @@ comb.dims <- function(arr.in, dims.to.combine){ ui = constr.mtrx, ci = constr.vec, freq.per.mdl.at.obs = freq.per.mdl.at.obs) init.par <- optim.tmp$par * (1 - abs(rnorm(amt.coeff, 0, 0.01))) - var.cat.fc[ , eval.dexes] <- apply(s2dverification::InsertDim( - s2dverification::InsertDim(optim.tmp$par, lendim = amt.cat, posdim = 2), + var.cat.fc[ , eval.dexes] <- apply(InsertDim( + InsertDim(optim.tmp$par, lendim = amt.cat, posdim = 2), lendim = amt.sdate.ev, posdim = 3) * freq.per.mdl.ev[ , , , drop = FALSE], c(2,3), sum, na.rm = TRUE) } else if (cat.method == "comb") { @@ -363,7 +389,7 @@ comb.dims <- function(arr.in, dims.to.combine){ .funct.optim.grad <- function(par, freq.per.mdl.at.obs){ amt.model <- dim(freq.per.mdl.at.obs)[1] - return(-apply(freq.per.mdl.at.obs/s2dverification::InsertDim(drop(par %*% + return(-apply(freq.per.mdl.at.obs/InsertDim(drop(par %*% freq.per.mdl.at.obs), lendim = amt.model, posdim = 1), c(1), mean, na.rm = TRUE)) } @@ -374,7 +400,7 @@ freq.per.mdl.at.obs), amt.mdl <- mdl.feat$amt.mdl mdl.msk.tmp <- mdl.feat$mdl.msk amt.coeff <- amt.mdl + 1 - msk.fc.obs <- (cat.fc == s2dverification::InsertDim(cat.obs, posdim = 1, lendim = amt.mbr)) + msk.fc.obs <- (cat.fc == InsertDim(cat.obs, posdim = 1, lendim = amt.mbr)) freq.per.mdl.at.obs <- array(NA, c(amt.coeff, amt.sdate)) for (i.mdl in seq(1, amt.mdl)){ freq.per.mdl.at.obs[i.mdl, ] <- apply(msk.fc.obs[mdl.msk.tmp[i.mdl, ], , drop = FALSE], diff --git a/R/CST_RegimesAssign.R b/R/CST_RegimesAssign.R index 72ae69f1..8985d335 100644 --- a/R/CST_RegimesAssign.R +++ b/R/CST_RegimesAssign.R @@ -27,7 +27,8 @@ #' that accounts for the serial dependence of the data with the same structure as Composite.)(only when composite = 'TRUE'), #' \code{$cluster} (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.) , #' \code{$frequency} (A vector of integers (from k=1,...k n reference maps) indicating the percentage of assignations corresponding to each map.), -#'@importFrom s2dverification ACC Mean1Dim InsertDim +#'@importFrom s2dverification ACC Mean1Dim +#'@importFrom s2dv InsertDim #'@import multiApply #'@examples #'\dontrun{ @@ -102,7 +103,8 @@ CST_RegimesAssign <- function(data, ref_maps, #' \code{$cluster} (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.) , #' \code{$frequency} (A vector of integers (from k = 1, ... k n reference maps) indicating the percentage of assignations corresponding to each map.), #' -#'@importFrom s2dverification ACC Mean1Dim Eno InsertDim +#'@importFrom s2dverification ACC Mean1Dim Eno +#'@importFrom s2dv InsertDim #'@import multiApply #'@examples #'\dontrun{ @@ -190,9 +192,9 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite = if (any(is.na(index))) { recon <-list( - composite = s2dverification::InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), + composite = InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), postime, dim(ref_maps)['composite.cluster']), - pvalue = s2dverification::InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), + pvalue = InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), postime, dim(ref_maps)['composite.cluster'])) } else { if (memb) { @@ -255,7 +257,7 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite = which(names(dim(target)) == 'lon'))) # weights are defined - latWeights <- s2dverification::InsertDim(sqrt(cos(lat * pi / 180)), 2, dim(ref)[3]) + latWeights <- InsertDim(sqrt(cos(lat * pi / 180)), 2, dim(ref)[3]) rmsdiff <- function(x, y) { @@ -278,12 +280,8 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite = corr <- rep(NA, nclust) for (i in 1:nclust) { corr[i] <- - ACC(s2dverification::InsertDim(InsertDim( - s2dverification::InsertDim(ref[i, , ] * latWeights, 1, 1), 2, 1 - ), 3, 1), - s2dverification::InsertDim(InsertDim( - s2dverification::InsertDim(target * latWeights, 1, 1), 2, 1 - ), 3, 1))$ACC[2] + ACC(InsertDim(InsertDim(InsertDim(ref[i, , ] * latWeights, 1, 1), 2, 1), 3, 1), + InsertDim(InsertDim(InsertDim(target * latWeights, 1, 1), 2, 1), 3, 1))$ACC[2] } assign <- which(corr == max(corr)) } diff --git a/R/CST_SaveExp.R b/R/CST_SaveExp.R index 4ae560bf..ac377101 100644 --- a/R/CST_SaveExp.R +++ b/R/CST_SaveExp.R @@ -17,8 +17,7 @@ #'@seealso \code{\link{CST_Load}}, \code{\link{as.s2dv_cube}} and \code{\link{s2dv_cube}} #' #'@import ncdf4 -#'@importFrom s2dv Reorder -#'@importFrom s2dverification InsertDim +#'@importFrom s2dv Reorder InsertDim #'@import multiApply #' #'@examples @@ -79,8 +78,7 @@ CST_SaveExp <- function(data, destination = "./CST_Data") { #' The path will be created with the name of the variable and each Datasets. #' #'@import ncdf4 -#'@importFrom s2dv Reorder -#'@importFrom s2dverification InsertDim +#'@importFrom s2dv Reorder InsertDim #'@import multiApply #' #'@examples @@ -144,7 +142,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, if (length(dataset_pos) == 0) { warning("Element 'data' in parameter 'data' hasn't 'dataset' dimension. ", "All data is stored in the same 'dataset' folder.") - data$data <- s2dverification::InsertDim(data, posdim = 1, lendim = 1) + data$data <- InsertDim(data, posdim = 1, lendim = 1) names(dim(data))[1] <- "dataset" dimname <- c("dataset", dimname) dataset_pos = 1 diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 8c4b1068..4cf6de6c 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -23,7 +23,7 @@ #'@importFrom reshape2 melt #'@importFrom plyr . #'@importFrom plyr dlply -#'@importFrom s2dverification InsertDim +#'@importFrom s2dv InsertDim #'@examples #'fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2), #' fcst3 = rnorm(10, -0.5, 0.9)) @@ -108,7 +108,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N if (length(tercile.limits) != 2) { stop("Provide two tercile limits") } - tercile.limits <- s2dverification::InsertDim(tercile.limits, 1, npanels) + tercile.limits <- InsertDim(tercile.limits, 1, npanels) } else if (is.array(tercile.limits)) { if (length(dim(tercile.limits)) == 2) { if (dim(tercile.limits)[2] != 2) { @@ -135,7 +135,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N if (length(extreme.limits) != 2) { stop("Provide two extreme limits") } - extreme.limits <- s2dverification::InsertDim(extreme.limits, 1, npanels) + extreme.limits <- InsertDim(extreme.limits, 1, npanels) } else if (is.array(extreme.limits)) { if (length(dim(extreme.limits)) == 2) { if (dim(extreme.limits)[2] != 2) { diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd index e8e047dc..dc2f7566 100644 --- a/man/CST_MultiMetric.Rd +++ b/man/CST_MultiMetric.Rd @@ -30,7 +30,7 @@ CST_MultiMetric( \item{sdate_dim}{name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'.} } \value{ -an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean. +an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the first position in the first 'nexp' dimension correspons to the Multi-Model Mean. } \description{ This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations. diff --git a/man/CategoricalEnsCombination.Rd b/man/CategoricalEnsCombination.Rd new file mode 100644 index 00000000..2f5ad14d --- /dev/null +++ b/man/CategoricalEnsCombination.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_CategoricalEnsCombination.R +\name{CategoricalEnsCombination} +\alias{CategoricalEnsCombination} +\title{Make categorical forecast based on a multi-model forecast with potential for calibrate} +\usage{ +CategoricalEnsCombination(fc, obs, cat.method, eval.method, amt.cat, ...) +} +\arguments{ +\item{fc}{a multi-dimensional array with named dimensions containing the seasonal forecast experiment data in the element named \code{$data}. The amount of forecasting models is equal to the size of the \code{dataset} dimension of the data array. The amount of members per model may be different. The size of the \code{member} dimension of the data array is equal to the maximum of the ensemble members among the models. Models with smaller ensemble sizes have residual indices of \code{member} dimension in the data array filled with NA values.} + +\item{obs}{a multidimensional array with named dimensions containing the observed data in the element named \code{$data}.} + +\item{cat.method}{method used to produce the categorical forecast, can be either \code{pool}, \code{comb}, \code{mmw} or \code{obs}. The method pool assumes equal weight for all ensemble members while the method comb assumes equal weight for each model. The weighting method is descirbed in Rajagopalan et al. (2002), Robertson et al. (2004) and Van Schaeybroeck and Vannitsem (2019). Finally, the \code{obs} method classifies the observations into the different categories and therefore contains only 0 and 1 values.} + +\item{eval.method}{is the sampling method used, can be either \code{"in-sample"} or \code{"leave-one-out"}. Default value is the \code{"leave-one-out"} cross validation.} + +\item{amt.cat}{is the amount of categories. Equally-sized quantiles will be calculated based on the amount of categories.} + +\item{...}{other parameters to be passed on to the calibration procedure.} +} +\value{ +an array containing the categorical forecasts in the element called \code{$data}. The first two dimensions of the returned object are named dataset and member and are both of size one. An additional dimension named category is introduced and is of size amt.cat. +} +\description{ +This function converts a multi-model ensemble forecast +into a categorical forecast by giving the probability +for each category. Different methods are available to combine +the different ensemble forecasting models into +probabilistic categorical forecasts. + +See details in ?CST_CategoricalEnsCombination +} +\references{ +Rajagopalan, B., Lall, U., & Zebiak, S. E. (2002). Categorical climate forecasts through regularization and optimal combination of multiple GCM ensembles. Monthly Weather Review, 130(7), 1792-1811. + +Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744. + +Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290). +} +\author{ +Bert Van Schaeybroeck, \email{bertvs@meteo.be} +} -- GitLab From 39f5a6b86d3ce857065af9f4e27ea0125bc193f5 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 1 Dec 2020 15:24:06 +0100 Subject: [PATCH 5/5] Fix CategoricalEnsCombination to use InsertDIm from s2dv --- R/CST_CategoricalEnsCombination.R | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/R/CST_CategoricalEnsCombination.R b/R/CST_CategoricalEnsCombination.R index 6664a2f1..909792ee 100644 --- a/R/CST_CategoricalEnsCombination.R +++ b/R/CST_CategoricalEnsCombination.R @@ -90,7 +90,7 @@ CST_CategoricalEnsCombination <- function(exp, obs, cat.method = "pool", eval.me eval.method = eval.method, amt.cat = amt.cat, ...) names.dim.tmp[which(names.dim.tmp == "member")] <- "category" names(dim(exp$data)) <- names.dim.tmp - exp$data <- InsertDim(exp$data, lendim = 1, posdim = 2) + exp$data <- suppressWarnings(InsertDim(exp$data, lendim = 1, posdim = 2)) names(dim(exp$data))[2] <- "member" exp$Datasets <- c(exp$Datasets, obs$Datasets) exp$source_files <- c(exp$source_files, obs$source_files) @@ -302,10 +302,10 @@ comb.dims <- function(arr.in, dims.to.combine){ optim.tmp <- constrOptim(theta = init.par, f = .funct.optim, grad = .funct.optim.grad, ui = constr.mtrx, ci = constr.vec, freq.per.mdl.at.obs = freq.per.mdl.at.obs) - init.par <- optim.tmp$par * (1 - abs(rnorm(amt.coeff, 0, 0.01))) - var.cat.fc[ , eval.dexes] <- apply(InsertDim( + init.par <- optim.tmp$par * (1 - abs(rnorm(amt.coeff, 0, 0.01))) + var.cat.fc[ , eval.dexes] <- apply(suppressWarnings(InsertDim( InsertDim(optim.tmp$par, lendim = amt.cat, posdim = 2), - lendim = amt.sdate.ev, posdim = 3) * + lendim = amt.sdate.ev, posdim = 3)) * freq.per.mdl.ev[ , , , drop = FALSE], c(2,3), sum, na.rm = TRUE) } else if (cat.method == "comb") { freq.per.mdl.ev <- .calc.freq.per.mdl(cat.fc = cat.fc.ev, mdl.feat = mdl.feat, amt.cat = amt.cat) @@ -389,9 +389,12 @@ comb.dims <- function(arr.in, dims.to.combine){ .funct.optim.grad <- function(par, freq.per.mdl.at.obs){ amt.model <- dim(freq.per.mdl.at.obs)[1] - return(-apply(freq.per.mdl.at.obs/InsertDim(drop(par %*% -freq.per.mdl.at.obs), - lendim = amt.model, posdim = 1), c(1), mean, na.rm = TRUE)) + preprocess <- drop(par %*% freq.per.mdl.at.obs) + if (is.null(dim(preprocess))) { + dim(preprocess) <- c(dim = length(preprocess)) + } + return(-apply(freq.per.mdl.at.obs/suppressWarnings(InsertDim(preprocess, + lendim = as.numeric(amt.model), posdim = 1)), c(1), mean, na.rm = TRUE)) } .calc.freq.per.mdl.at.obs <- function(cat.obs, cat.fc, amt.cat, mdl.feat){ -- GitLab