From 2bcbbc04d8a0fb2ed0b7e08c63cf064055fbc23e Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Tue, 13 Apr 2021 17:52:33 +0200
Subject: [PATCH 1/6] Transform Spread() from s2dverification

---
 NAMESPACE     |   4 +
 R/Spread.R    | 203 ++++++++++++++++++++++++++++++++++++++++++++++++++
 man/Spread.Rd | 100 +++++++++++++++++++++++++
 3 files changed, 307 insertions(+)
 create mode 100644 R/Spread.R
 create mode 100644 man/Spread.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 87937e2..ea090db 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -47,6 +47,7 @@ export(Reorder)
 export(SPOD)
 export(Season)
 export(Smoothing)
+export(Spread)
 export(TPI)
 export(ToyModel)
 export(Trend)
@@ -84,12 +85,14 @@ importFrom(grDevices,rainbow)
 importFrom(grDevices,rgb)
 importFrom(grDevices,svg)
 importFrom(grDevices,tiff)
+importFrom(stats,IQR)
 importFrom(stats,acf)
 importFrom(stats,anova)
 importFrom(stats,confint)
 importFrom(stats,cor)
 importFrom(stats,kmeans)
 importFrom(stats,lm)
+importFrom(stats,mad)
 importFrom(stats,median)
 importFrom(stats,na.fail)
 importFrom(stats,na.omit)
@@ -102,6 +105,7 @@ importFrom(stats,qnorm)
 importFrom(stats,qt)
 importFrom(stats,quantile)
 importFrom(stats,rnorm)
+importFrom(stats,runif)
 importFrom(stats,sd)
 importFrom(stats,setNames)
 importFrom(stats,ts)
diff --git a/R/Spread.R b/R/Spread.R
new file mode 100644
index 0000000..d808dc7
--- /dev/null
+++ b/R/Spread.R
@@ -0,0 +1,203 @@
+#'Compute interquartile range, maximum-minimum, standard deviation and median
+#'absolute deviation
+#'
+#'Compute interquartile range, maximum-minimum, standard deviation and median
+#'absolute deviation along the list of dimensions provided by the compute_dim 
+#'argument (typically along the ensemble member and start date dimension).
+#'The confidence interval is computed by bootstrapping. The input data can
+#'be the output of \code{Load()}, \code{Ano()}, or \code{Ano_CrossValid()}, for
+#'example.
+#'
+#'@param data A numeric vector or array with named dimensions to compute the 
+#'  statistics. The dimensions should at least include 'compute_dim'.
+#'@param compute_dim A vector of character strings of the dimension names along 
+#'  which to compute the statistics. The default value is 'member'.
+#'@param na.rm A logical value indicating if NAs should be removed (TRUE) or
+#'  kept (FALSE) for computation. The default value is TRUE.
+#'@param conf A logical value indicating whether to compute the confidence 
+#'  intervals or not. The default value is TRUE.
+#'@param conf.lev A numeric value of the confidence level for the computation. 
+#'  The default value is 0.95.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return 
+#'A list of numeric arrays with the same dimensions as 'data' but without
+#''compute_dim' and with the first dimension 'stats'. If 'conf' is TRUE, the 
+#'length of 'stats' is 3 corresponding to the lower limit of the confidence 
+#'interval, the spread, and the upper limit of the confidence interval. If 
+#''conf' is FALSE, the length of 'stats' is 1 corresponding to the spread.
+#'\item{$iqr}{
+#'  InterQuartile Range.
+#'}
+#'\item{$maxmin}{
+#'  Maximum - Minimum.
+#'}
+#'\item{$sd}{
+#'  Standard Deviation.
+#'}
+#'\item{$mad}{
+#'  Median Absolute Deviation.
+#'} 
+#'
+#'@examples
+#'# Load sample data as in Load() example:
+#'example(Load)
+#'clim <- Clim(sampleData$mod, sampleData$obs)
+#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+#'runmean_months <- 12
+#'smooth_ano_exp <- Smoothing(ano_exp, runmeanlen = runmean_months)
+#'smooth_ano_exp_m_sub <- smooth_ano_exp - InsertDim(MeanDims(smooth_ano_exp, 'member', 
+#'                                                            na.rm = TRUE), 
+#'                                                   posdim = 'member', 
+#'                                                   lendim = dim(smooth_ano_exp)['member'], 
+#'                                                   name = 'member')
+#'spread <- Spread(smooth_ano_exp_m_sub, compute_dim = c('member', 'sdate'))
+#'
+#'\donttest{
+#'PlotVsLTime(Reorder(spread$iqr, c('dataset', 'stats', 'ftime')), 
+#'            toptitle = "Inter-Quartile Range between ensemble members",
+#'            ytitle = "K", monini = 11, limits = NULL, 
+#'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
+#'            hlines = c(0), fileout = 'tos_iqr.png')
+#'PlotVsLTime(Reorder(spread$maxmin, c('dataset', 'stats', 'ftime')), 
+#'            toptitle = "Maximum minus minimum of the members", 
+#'            ytitle = "K", monini = 11, limits = NULL, 
+#'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
+#'            hlines = c(0), fileout = 'tos_maxmin.png')
+#'PlotVsLTime(Reorder(spread$sd, c('dataset', 'stats', 'ftime')), 
+#'            toptitle = "Standard deviation of the members", 
+#'            ytitle = "K", monini = 11, limits = NULL, 
+#'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
+#'            hlines = c(0), fileout = 'tos_sd.png')
+#'PlotVsLTime(Reorder(spread$mad, c('dataset', 'stats', 'ftime')), 
+#'            toptitle = "Median Absolute Deviation of the members",
+#'            ytitle = "K", monini = 11, limits = NULL, 
+#'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
+#'            hlines = c(0), fileout = 'tos_mad.png')
+#'}
+#'
+#'@import multiApply
+#'@importFrom stats IQR sd mad runif quantile
+#'@export
+Spread <- function(data, compute_dim = 'member', na.rm = TRUE, 
+                   conf = TRUE, conf.lev = 0.95, ncores = NULL) {
+
+  # Check inputs 
+  ## data
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
+  }
+  if (!is.numeric(data)) {
+    stop("Parameter 'data' must be a numeric array.")
+  }
+  if (is.null(dim(data))) {  #is vector
+    dim(data) <- c(length(data))
+    names(dim(data)) <- compute_dim[1]
+  }
+  if(any(is.null(names(dim(data))))| any(nchar(names(dim(data))) == 0)) {
+    stop("Parameter 'data' must have dimension names.")
+  }
+  ## compute_dim
+  if (!is.character(compute_dim)) {
+    stop("Parameter 'compute_dim' must be a character vector.")
+  }
+  if (any(!compute_dim %in% names(dim(data)))) {
+    stop("Parameter 'compute_dim' has some element not in 'data' dimension names.")
+  }
+  ## na.rm
+  if (!is.logical(na.rm) | length(na.rm) > 1) {
+    stop("Parameter 'na.rm' must be one logical value.")
+  }
+  ## conf
+  if (!is.logical(conf) | length(conf) > 1) {
+    stop("Parameter 'conf' must be one logical value.")
+  }
+  ## conf.lev
+  if (!is.numeric(conf.lev) | conf.lev < 0 | conf.lev > 1 | length(conf.lev) > 1) {
+    stop("Parameter 'conf.lev' must be a numeric number between 0 and 1.")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+      length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  } 
+
+  ###############################
+  # Calculate Spread
+
+  output <- Apply(list(data),
+                  target_dims = compute_dim,
+                  fun = .Spread,
+                  output_dims = list(iqr = 'stats', maxmin = 'stats', 
+                                     sd = 'stats', mad = 'stats'),
+                  na.rm = na.rm,
+                  conf = conf, conf.lev = conf.lev,
+                  ncores = ncores)
+
+  return(output)
+}
+
+.Spread <- function(data, compute_dim = 'member', na.rm = TRUE,
+                    conf = TRUE, conf.lev = 0.95) {
+
+  # data: compute_dim. [member] or [member, sdate] for example
+
+  # Compute spread
+  res_iqr <- IQR(data, na.rm = na.rm)
+  res_maxmin <- max(data, na.rm = na.rm) - min(data, na.rm = na.rm)
+  res_sd <- sd(data, na.rm = na.rm)
+  res_mad <- mad(data, na.rm = na.rm)
+
+  # Compute conf (bootstrapping)
+  if (conf) {
+    # The output length is 3, [conf.low, spread, conf.high]
+    res_iqr <- rep(res_iqr, 3)
+    res_maxmin <- rep(res_maxmin, 3)
+    res_sd <- rep(res_sd, 3)
+    res_mad <- rep(res_mad, 3)
+
+    conf_low <- (1 - conf.lev) / 2
+    conf_high <- 1 - conf_low
+  
+    # Create vector for saving bootstrap result
+    iqr_bs <- c()
+    maxmin_bs <- c()
+    sd_bs <- c()
+    mad_bs <- c()
+
+    # bootstrapping for 100 times
+    num <- length(data)
+    for (jmix in 1:100) {
+      drawings <- round(runif(num, 0.5, num + 0.5))
+      iqr_bs <- c(iqr_bs, IQR(data[drawings], na.rm = na.rm))
+      maxmin_bs <- c(maxmin_bs, max(data[drawings], na.rm = na.rm) - 
+                                min(data[drawings], na.rm = na.rm))
+      sd_bs <- c(sd_bs, sd(data[drawings], na.rm = na.rm))
+      mad_bs <- c(mad_bs, mad(data[drawings], na.rm = na.rm))
+    }
+    
+    # Calculate confidence interval with the bootstrapping results
+    res_iqr[1] <- quantile(iqr_bs, conf_low, na.rm = na.rm)
+    res_iqr[3] <- quantile(iqr_bs, conf_high, na.rm = na.rm)
+    res_maxmin[1] <- res_maxmin[2] + (quantile(maxmin_bs, conf_low, na.rm = na.rm) - 
+                                 quantile(maxmin_bs, conf_high, na.rm = na.rm)) / 2
+    res_maxmin[3] <- res_maxmin[2] - (quantile(maxmin_bs, conf_low, na.rm = na.rm) - 
+                                 quantile(maxmin_bs, conf_high, na.rm = na.rm)) / 2
+    res_sd[1] <- quantile(sd_bs, conf_low, na.rm = na.rm)
+    res_sd[3] <- quantile(sd_bs, conf_high, na.rm = na.rm)
+    res_mad[1] <- res_mad[2] + (quantile(mad_bs, conf_low, na.rm = na.rm) - 
+                                quantile(mad_bs, conf_high, na.rm = na.rm))
+    res_mad[3] <- res_mad[2] - (quantile(mad_bs, conf_low, na.rm = na.rm) - 
+                                quantile(mad_bs, conf_high, na.rm = na.rm))
+
+  }
+
+  # Turn infinite to NA
+  res_maxmin[which(is.infinite(res_maxmin))] <- NA
+
+
+  return(invisible(list(iqr = res_iqr, maxmin = res_maxmin, sd = res_sd, mad = res_mad)))
+}
diff --git a/man/Spread.Rd b/man/Spread.Rd
new file mode 100644
index 0000000..7eca295
--- /dev/null
+++ b/man/Spread.Rd
@@ -0,0 +1,100 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Spread.R
+\name{Spread}
+\alias{Spread}
+\title{Compute interquartile range, maximum-minimum, standard deviation and median
+absolute deviation}
+\usage{
+Spread(
+  data,
+  compute_dim = "member",
+  na.rm = TRUE,
+  conf = TRUE,
+  conf.lev = 0.95,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{A numeric vector or array with named dimensions to compute the 
+statistics. The dimensions should at least include 'compute_dim'.}
+
+\item{compute_dim}{A vector of character strings of the dimension names along 
+which to compute the statistics. The default value is 'member'.}
+
+\item{na.rm}{A logical value indicating if NAs should be removed (TRUE) or
+kept (FALSE) for computation. The default value is TRUE.}
+
+\item{conf}{A logical value indicating whether to compute the confidence 
+intervals or not. The default value is TRUE.}
+
+\item{conf.lev}{A numeric value of the confidence level for the computation. 
+The default value is 0.95.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+A list of numeric arrays with the same dimensions as 'data' but without
+'compute_dim' and with the first dimension 'stats'. If 'conf' is TRUE, the 
+length of 'stats' is 3 corresponding to the lower limit of the confidence 
+interval, the spread, and the upper limit of the confidence interval. If 
+'conf' is FALSE, the length of 'stats' is 1 corresponding to the spread.
+\item{$iqr}{
+ InterQuartile Range.
+}
+\item{$maxmin}{
+ Maximum - Minimum.
+}
+\item{$sd}{
+ Standard Deviation.
+}
+\item{$mad}{
+ Median Absolute Deviation.
+}
+}
+\description{
+Compute interquartile range, maximum-minimum, standard deviation and median
+absolute deviation along the list of dimensions provided by the compute_dim 
+argument (typically along the ensemble member and start date dimension).
+The confidence interval is computed by bootstrapping. The input data can
+be the output of \code{Load()}, \code{Ano()}, or \code{Ano_CrossValid()}, for
+example.
+}
+\examples{
+# Load sample data as in Load() example:
+example(Load)
+clim <- Clim(sampleData$mod, sampleData$obs)
+ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+runmean_months <- 12
+smooth_ano_exp <- Smoothing(ano_exp, runmeanlen = runmean_months)
+smooth_ano_exp_m_sub <- smooth_ano_exp - InsertDim(MeanDims(smooth_ano_exp, 'member', 
+                                                           na.rm = TRUE), 
+                                                  posdim = 'member', 
+                                                  lendim = dim(smooth_ano_exp)['member'], 
+                                                  name = 'member')
+spread <- Spread(smooth_ano_exp_m_sub, compute_dim = c('member', 'sdate'))
+
+\donttest{
+PlotVsLTime(Reorder(spread$iqr, c('dataset', 'stats', 'ftime')), 
+           toptitle = "Inter-Quartile Range between ensemble members",
+           ytitle = "K", monini = 11, limits = NULL, 
+           listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
+           hlines = c(0), fileout = 'tos_iqr.png')
+PlotVsLTime(Reorder(spread$maxmin, c('dataset', 'stats', 'ftime')), 
+           toptitle = "Maximum minus minimum of the members", 
+           ytitle = "K", monini = 11, limits = NULL, 
+           listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
+           hlines = c(0), fileout = 'tos_maxmin.png')
+PlotVsLTime(Reorder(spread$sd, c('dataset', 'stats', 'ftime')), 
+           toptitle = "Standard deviation of the members", 
+           ytitle = "K", monini = 11, limits = NULL, 
+           listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
+           hlines = c(0), fileout = 'tos_sd.png')
+PlotVsLTime(Reorder(spread$mad, c('dataset', 'stats', 'ftime')), 
+           toptitle = "Median Absolute Deviation of the members",
+           ytitle = "K", monini = 11, limits = NULL, 
+           listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
+           hlines = c(0), fileout = 'tos_mad.png')
+}
+
+}
-- 
GitLab


From 45758f3f50bc14e999729b681e35c76cb668b901 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Tue, 13 Apr 2021 17:52:47 +0200
Subject: [PATCH 2/6] Update indices .Rd

---
 man/AMV.Rd  | 29 ++++++++++++++---------------
 man/GMST.Rd | 51 ++++++++++++++++++++++++---------------------------
 man/GSAT.Rd | 30 +++++++++++++++---------------
 man/SPOD.Rd | 30 +++++++++++++++---------------
 man/TPI.Rd  | 30 +++++++++++++++---------------
 5 files changed, 83 insertions(+), 87 deletions(-)

diff --git a/man/AMV.Rd b/man/AMV.Rd
index 0f81652..3cc4113 100644
--- a/man/AMV.Rd
+++ b/man/AMV.Rd
@@ -18,16 +18,15 @@ AMV(
   indices_for_clim = NULL,
   year_dim = "year",
   month_dim = "month",
-  member_dim = "member",
+  na.rm = TRUE,
   ncores = NULL
 )
 }
 \arguments{
-\item{data}{A numerical array to be used for the index computation with the
-dimensions: 1) latitude, longitude, start date, forecast month, and member 
-(in case of decadal predictions), 2) latitude, longitude, year, month and 
-member (in case of historical simulations), or 3) latitude, longitude, year
-and month (in case of observations or reanalyses). This data has to be 
+\item{data}{A numerical array to be used for the index computation with, at least, the
+dimensions: 1) latitude, longitude, start date and forecast month 
+(in case of decadal predictions), 2) latitude, longitude, year and month 
+(in case of historical simulations or observations). This data has to be 
 provided, at least, over the whole region needed to compute the index.}
 
 \item{data_lats}{A numeric vector indicating the latitudes of the data.}
@@ -70,7 +69,7 @@ climatology is calculated over the whole period. If the data are already
 anomalies, set it to FALSE. The default value is NULL.\cr
 In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative 
 to the first forecast year, and the climatology is automatically computed 
-over the actual common period for the different forecast years.}
+over the common calendar period for the different forecast years.}
 
 \item{year_dim}{A character string indicating the name of the year dimension
 The default value is 'year'. Only used if parameter 'type' is 'hist' or 
@@ -80,18 +79,17 @@ The default value is 'year'. Only used if parameter 'type' is 'hist' or
 dimension. The default value is 'month'. Only used if parameter 'type' is 
 'hist' or 'obs'.}
 
-\item{member_dim}{A character string indicating the name of the member 
-dimension. The default value is 'member'. Only used if parameter 'type' is
-'dcpp' or 'hist'.}
+\item{na.rm}{A logical value indicanting whether to remove NA values. The default 
+value is TRUE.}
 
 \item{ncores}{An integer indicating the number of cores to use for parallel 
 computation. The default value is NULL.}
 }
 \value{
-A numerical array of the AMV index with the dimensions of:
- 1)  sdate, forecast year, and member (in case of decadal predictions);
- 2)  year and member (in case of historical simulations); or 
- 3)  year (in case of observations or reanalyses).
+A numerical array with the AMV index with the same dimensions as data except 
+ the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions 
+ (historical simulations or observations). In case of decadal predictions, a new dimension
+ 'fyear' is added.
 }
 \description{
 The Atlantic Multidecadal Variability (AMV), also known as Atlantic
@@ -100,7 +98,8 @@ surface temperatures (SST) over the North Atlantic Ocean on multi-decadal
 time scales. The AMV index is computed as the difference of weighted-averaged
 SST anomalies over the North Atlantic region (0ºN-60ºN, 280ºE-360ºE) and the 
 weighted-averaged SST anomalies over 60ºS-60ºN, 0ºE-360ºE (Trenberth & 
-Dennis, 2005; Doblas-Reyes et al., 2013).
+Dennis, 2005; Doblas-Reyes et al., 2013). If different members and/or datasets are provided, 
+the climatology (used to calculate the anomalies) is computed individually for all of them.
 }
 \examples{
 ## Observations or reanalyses
diff --git a/man/GMST.Rd b/man/GMST.Rd
index 8021943..f23e60d 100644
--- a/man/GMST.Rd
+++ b/man/GMST.Rd
@@ -21,28 +21,25 @@ GMST(
   indices_for_clim = NULL,
   year_dim = "year",
   month_dim = "month",
-  member_dim = "member",
+  na.rm = TRUE,
   ncores = NULL
 )
 }
 \arguments{
-\item{data_tas}{A numerical array indicating the surface air temperature data
-to be used for the index computation with the dimensions: 1) latitude, 
-longitude, start date, forecast month, and member (in case of decadal 
-predictions), 2) latitude, longitude, year, month and member (in case of 
-historical simulations), or 3) latitude, longitude, year and month (in case
-of observations or reanalyses). This data has to be provided, at least, 
-over the whole region needed to compute the index. The dimensions must be 
-identical to those of data_tos.}
-
-\item{data_tos}{A numerical array indicating the sea surface temperature data
-to be used for the index computation with the dimensions: 1) latitude, 
-longitude, start date, forecast month, and member (in case of decadal 
-predictions), 2) latitude, longitude, year, month and member (in case of 
-historical simulations), or 3) latitude, longitude, year and month (in case
-of observations or reanalyses). This data has to be provided, at least, 
-over the whole region needed to compute the index. The dimensions must be
-identical to those of data_tas.}
+\item{data_tas}{A numerical array  with the surface air temperature data
+to be used for the index computation with, at least, the
+dimensions: 1) latitude, longitude, start date and forecast month 
+(in case of decadal predictions), 2) latitude, longitude, year and month 
+(in case of historical simulations or observations). This data has to be 
+provided, at least, over the whole region needed to compute the index.
+The dimensions must be identical to thos of data_tos.
+#'@param data_tos A numerical array  with the sea surface temperature data
+to be used for the index computation with, at least, the
+dimensions: 1) latitude, longitude, start date and forecast month 
+(in case of decadal predictions), 2) latitude, longitude, year and month 
+(in case of historical simulations or observations). This data has to be 
+provided, at least, over the whole region needed to compute the index.
+The dimensions must be identical to thos of data_tas.}
 
 \item{data_lats}{A numeric vector indicating the latitudes of the data.}
 
@@ -90,7 +87,7 @@ climatology is calculated over the whole period. If the data are already
 anomalies, set it to FALSE. The default value is NULL.\cr
 In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative 
 to the first forecast year, and the climatology is automatically computed 
-over the actual common period for the different forecast years.}
+over the common calendar period for the different forecast years.}
 
 \item{year_dim}{A character string indicating the name of the year dimension
 The default value is 'year'. Only used if parameter 'type' is 'hist' or 
@@ -100,23 +97,23 @@ The default value is 'year'. Only used if parameter 'type' is 'hist' or
 dimension. The default value is 'month'. Only used if parameter 'type' is 
 'hist' or 'obs'.}
 
-\item{member_dim}{A character string indicating the name of the member 
-dimension. The default value is 'member'. Only used if parameter 'type' is
-'dcpp' or 'hist'.}
+\item{na.rm}{A logical value indicanting whether to remove NA values. The default 
+value is TRUE.}
 
 \item{ncores}{An integer indicating the number of cores to use for parallel 
 computation. The default value is NULL.}
 }
 \value{
-A numerical array of the GMST anomalies with the dimensions of:
- 1) sdate, forecast year, and member (in case of decadal predictions);
- 2) year and member (in case of historical simulations); or
- 3) year (in case of observations or reanalyses).
+A numerical array with the GMST anomalies with the same dimensions as data_tas except 
+ the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions 
+ (historical simulations or observations). In case of decadal predictions, a new dimension
+ 'fyear' is added.
 }
 \description{
 The Global Mean Surface Temperature (GMST) anomalies are computed as the
 weighted-averaged surface air temperature anomalies over land and sea surface
-temperature anomalies over the ocean.
+temperature anomalies over the ocean. If different members and/or datasets are provided, 
+the climatology (used to calculate the anomalies) is computed individually for all of them.
 }
 \examples{
 ## Observations or reanalyses
diff --git a/man/GSAT.Rd b/man/GSAT.Rd
index d7fe3cf..9f03ff2 100644
--- a/man/GSAT.Rd
+++ b/man/GSAT.Rd
@@ -18,16 +18,15 @@ GSAT(
   indices_for_clim = NULL,
   year_dim = "year",
   month_dim = "month",
-  member_dim = "member",
+  na.rm = TRUE,
   ncores = NULL
 )
 }
 \arguments{
-\item{data}{A numerical array to be used for the index computation with the 
-dimensions: 1) latitude, longitude, start date, forecast month, and member 
-(in case of decadal predictions), 2) latitude, longitude, year, month and 
-member (in case of historical simulations), or 3) latitude, longitude, year
-and month (in case of observations or reanalyses). This data has to be 
+\item{data}{A numerical array to be used for the index computation with, at least, the
+dimensions: 1) latitude, longitude, start date and forecast month 
+(in case of decadal predictions), 2) latitude, longitude, year and month 
+(in case of historical simulations or observations). This data has to be 
 provided, at least, over the whole region needed to compute the index.}
 
 \item{data_lats}{A numeric vector indicating the latitudes of the data.}
@@ -70,7 +69,7 @@ climatology is calculated over the whole period. If the data are already
 anomalies, set it to FALSE. The default value is NULL.\cr
 In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative 
 to the first forecast year, and the climatology is automatically computed 
-over the actual common period for the different forecast years.}
+over the common calendar period for the different forecast years.}
 
 \item{year_dim}{A character string indicating the name of the year dimension
 The default value is 'year'. Only used if parameter 'type' is 'hist' or 
@@ -80,22 +79,23 @@ The default value is 'year'. Only used if parameter 'type' is 'hist' or
 dimension. The default value is 'month'. Only used if parameter 'type' is 
 'hist' or 'obs'.}
 
-\item{member_dim}{A character string indicating the name of the member 
-dimension. The default value is 'member'. Only used if parameter 'type' is
-'dcpp' or 'hist'.}
+\item{na.rm}{A logical value indicanting whether to remove NA values. The default 
+value is TRUE.}
 
 \item{ncores}{An integer indicating the number of cores to use for parallel 
 computation. The default value is NULL.}
 }
 \value{
-A numerical array of the GSAT anomalies with the dimensions of:
- 1) sdate, forecast year, and member (in case of decadal predictions);
- 2) year and member (in case of historical simulations); or 
- 3) year (in case of observations or reanalyses).
+A numerical array with the GSAT anomalies with the same dimensions as data except 
+ the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions 
+ (historical simulations or observations). In case of decadal predictions, a new dimension
+ 'fyear' is added.
 }
 \description{
 The Global Surface Air Temperature (GSAT) anomalies are computed as the 
-weighted-averaged surface air temperature anomalies over the global region.
+weighted-averaged surface air temperature anomalies over the global region. 
+If different members and/or datasets are provided, the climatology (used to 
+calculate the anomalies) is computed individually for all of them.
 }
 \examples{
 ## Observations or reanalyses
diff --git a/man/SPOD.Rd b/man/SPOD.Rd
index 0491739..4c4ed29 100644
--- a/man/SPOD.Rd
+++ b/man/SPOD.Rd
@@ -18,16 +18,15 @@ SPOD(
   indices_for_clim = NULL,
   year_dim = "year",
   month_dim = "month",
-  member_dim = "member",
+  na.rm = TRUE,
   ncores = NULL
 )
 }
 \arguments{
-\item{data}{A numerical array to be used for the index computation with the
-dimensions: 1) latitude, longitude, start date, forecast month, and member 
-(in case of decadal predictions), 2) latitude, longitude, year, month and 
-member (in case of historical simulations), or 3) latitude, longitude, year
-and month (in case of observations or reanalyses). This data has to be 
+\item{data}{A numerical array to be used for the index computation with, at least, the
+dimensions: 1) latitude, longitude, start date and forecast month 
+(in case of decadal predictions), 2) latitude, longitude, year and month 
+(in case of historical simulations or observations). This data has to be 
 provided, at least, over the whole region needed to compute the index.}
 
 \item{data_lats}{A numeric vector indicating the latitudes of the data.}
@@ -70,7 +69,7 @@ climatology is calculated over the whole period. If the data are already
 anomalies, set it to FALSE. The default value is NULL.\cr
 In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative 
 to the first forecast year, and the climatology is automatically computed 
-over the actual common period for the different forecast years.}
+over the common calendar period for the different forecast years.}
 
 \item{year_dim}{A character string indicating the name of the year dimension
 The default value is 'year'. Only used if parameter 'type' is 'hist' or 
@@ -80,25 +79,26 @@ The default value is 'year'. Only used if parameter 'type' is 'hist' or
 dimension. The default value is 'month'. Only used if parameter 'type' is 
 'hist' or 'obs'.}
 
-\item{member_dim}{A character string indicating the name of the member 
-dimension. The default value is 'member'. Only used if parameter 'type' is
-'dcpp' or 'hist'.}
+\item{na.rm}{A logical value indicanting whether to remove NA values. The default 
+value is TRUE.}
 
 \item{ncores}{An integer indicating the number of cores to use for parallel 
 computation. The default value is NULL.}
 }
 \value{
-A numerical array of the SPOD index with the dimensions of:
- 1) sdate, forecast year, and member (in case of decadal predictions); 
- 2) year and member (in case of historical simulations); or
- 3) year (in case of observations or reanalyses).
+A numerical array with the SPOD index with the same dimensions as data except 
+ the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions 
+ (historical simulations or observations). In case of decadal predictions, a new dimension
+ 'fyear' is added.
 }
 \description{
 The South Pacific Ocean Dipole (SPOD) index is related to the El 
 Nino-Southern Oscillation (ENSO) and the Inderdecadal Pacific Oscillation 
 (IPO). The SPOD index is computed as the difference of weighted-averaged SST
 anomalies over 20ºS-48ºS, 165ºE-190ºE (NW pole) and the weighted-averaged SST
-anomalies over 44ºS-65ºS, 220ºE-260ºE (SE pole) (Saurral et al., 2020).
+anomalies over 44ºS-65ºS, 220ºE-260ºE (SE pole) (Saurral et al., 2020). 
+If different members and/or datasets are provided, the climatology (used to 
+calculate the anomalies) is computed individually for all of them.
 }
 \examples{
 ## Observations or reanalyses
diff --git a/man/TPI.Rd b/man/TPI.Rd
index 3bdc17c..5e8f716 100644
--- a/man/TPI.Rd
+++ b/man/TPI.Rd
@@ -18,16 +18,15 @@ TPI(
   indices_for_clim = NULL,
   year_dim = "year",
   month_dim = "month",
-  member_dim = "member",
+  na.rm = TRUE,
   ncores = NULL
 )
 }
 \arguments{
-\item{data}{A numerical array to be used for the index computation with the
-dimensions: 1) latitude, longitude, start date, forecast month, and member 
-(in case of decadal predictions), 2) latitude, longitude, year, month and 
-member (in case of historical simulations), or 3) latitude, longitude, year
-and month (in case of observations or reanalyses). This data has to be 
+\item{data}{A numerical array to be used for the index computation with, at least, the
+dimensions: 1) latitude, longitude, start date and forecast month 
+(in case of decadal predictions), 2) latitude, longitude, year and month 
+(in case of historical simulations or observations). This data has to be 
 provided, at least, over the whole region needed to compute the index.}
 
 \item{data_lats}{A numeric vector indicating the latitudes of the data.}
@@ -70,7 +69,7 @@ climatology is calculated over the whole period. If the data are already
 anomalies, set it to FALSE. The default value is NULL.\cr
 In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative 
 to the first forecast year, and the climatology is automatically computed 
-over the actual common period for the different forecast years.}
+over the common calendar period for the different forecast years.}
 
 \item{year_dim}{A character string indicating the name of the year dimension
 The default value is 'year'. Only used if parameter 'type' is 'hist' or 
@@ -80,24 +79,25 @@ The default value is 'year'. Only used if parameter 'type' is 'hist' or
 dimension. The default value is 'month'. Only used if parameter 'type' is 
 'hist' or 'obs'.}
 
-\item{member_dim}{A character string indicating the name of the member 
-dimension. The default value is 'member'. Only used if parameter 'type' is
-'dcpp' or 'hist'.}
+\item{na.rm}{A logical value indicanting whether to remove NA values. The default 
+value is TRUE.}
 
 \item{ncores}{An integer indicating the number of cores to use for parallel 
 computation. The default value is NULL.}
 }
 \value{
-A numerical array of the TPI index with the dimensions of:
- 1)  sdate, forecast year, and member (in case of decadal predictions);
- 2)  year and member (in case of historical simulations); or 
- 3)  year (in case of observations or reanalyses).
+A numerical array with the TPI index with the same dimensions as data except 
+ the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions 
+ (historical simulations or observations). In case of decadal predictions, a new dimension
+ 'fyear' is added.
 }
 \description{
 The Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO) is 
 computed as the difference of weighted-averaged SST anomalies over 10ºS-10ºN, 
 170ºE-270ºE minus the mean of the weighted-averaged SST anomalies over 
-25ºN-45ºN, 140ºE-215ºE and 50ºS-15ºS, 150ºE-200ºE (Henley et al., 2015).
+25ºN-45ºN, 140ºE-215ºE and 50ºS-15ºS, 150ºE-200ºE (Henley et al., 2015). 
+If different members and/or datasets are provided, the climatology (used to 
+calculate the anomalies) is computed individually for all of them.
 }
 \examples{
 ## Observations or reanalyses
-- 
GitLab


From 52529d35be4d6830435cf9b34ce844bd9a56d328 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Tue, 13 Apr 2021 20:41:32 +0200
Subject: [PATCH 3/6] Create unit test and improve document

---
 R/Spread.R                   |  10 +--
 man/Spread.Rd                |   6 +-
 tests/testthat/test-Spread.R | 164 +++++++++++++++++++++++++++++++++++
 3 files changed, 172 insertions(+), 8 deletions(-)
 create mode 100644 tests/testthat/test-Spread.R

diff --git a/R/Spread.R b/R/Spread.R
index d808dc7..9589fae 100644
--- a/R/Spread.R
+++ b/R/Spread.R
@@ -4,9 +4,9 @@
 #'Compute interquartile range, maximum-minimum, standard deviation and median
 #'absolute deviation along the list of dimensions provided by the compute_dim 
 #'argument (typically along the ensemble member and start date dimension).
-#'The confidence interval is computed by bootstrapping. The input data can
-#'be the output of \code{Load()}, \code{Ano()}, or \code{Ano_CrossValid()}, for
-#'example.
+#'The confidence interval is computed by bootstrapping by 100 times. The input
+#'data can be the output of \code{Load()}, \code{Ano()}, or 
+#'\code{Ano_CrossValid()}, for example.
 #'
 #'@param data A numeric vector or array with named dimensions to compute the 
 #'  statistics. The dimensions should at least include 'compute_dim'.
@@ -198,6 +198,6 @@ Spread <- function(data, compute_dim = 'member', na.rm = TRUE,
   # Turn infinite to NA
   res_maxmin[which(is.infinite(res_maxmin))] <- NA
 
-
-  return(invisible(list(iqr = res_iqr, maxmin = res_maxmin, sd = res_sd, mad = res_mad)))
+  return(invisible(list(iqr = as.array(res_iqr), maxmin = as.array(res_maxmin), 
+                        sd = as.array(res_sd), mad = as.array(res_mad))))
 }
diff --git a/man/Spread.Rd b/man/Spread.Rd
index 7eca295..e5f3edd 100644
--- a/man/Spread.Rd
+++ b/man/Spread.Rd
@@ -56,9 +56,9 @@ interval, the spread, and the upper limit of the confidence interval. If
 Compute interquartile range, maximum-minimum, standard deviation and median
 absolute deviation along the list of dimensions provided by the compute_dim 
 argument (typically along the ensemble member and start date dimension).
-The confidence interval is computed by bootstrapping. The input data can
-be the output of \code{Load()}, \code{Ano()}, or \code{Ano_CrossValid()}, for
-example.
+The confidence interval is computed by bootstrapping by 100 times. The input
+data can be the output of \code{Load()}, \code{Ano()}, or 
+\code{Ano_CrossValid()}, for example.
 }
 \examples{
 # Load sample data as in Load() example:
diff --git a/tests/testthat/test-Spread.R b/tests/testthat/test-Spread.R
new file mode 100644
index 0000000..1d299a6
--- /dev/null
+++ b/tests/testthat/test-Spread.R
@@ -0,0 +1,164 @@
+context("s2dv::Spread test")
+
+##############################################
+  # dat1
+  set.seed(1)
+  dat1 <- array(rnorm(240), c(member = 3, sdate = 4, ftime = 20))
+
+  # dat2
+  set.seed(2)
+  dat2 <- rnorm(20)
+
+##############################################
+
+test_that("1. Input checks", {
+
+  # data
+  expect_error(
+  Spread(c()),
+  "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    Spread(c(NA, NA)),
+    "Parameter 'data' must be a numeric array."
+  )
+  expect_error(
+    Spread(array(1:10, dim = c(2, 5))),
+    "Parameter 'data' must have dimension names."
+  )
+  # compute_dim
+  expect_error(
+  Spread(dat1, compute_dim = 1),
+  "Parameter 'compute_dim' must be a character vector."
+  )
+  expect_error(
+  Spread(dat1, compute_dim = 'memb'),
+  "Parameter 'compute_dim' has some element not in 'data' dimension names."
+  )
+  # na.rm
+  expect_error(
+  Spread(dat1, na.rm = 1),
+  "Parameter 'na.rm' must be one logical value."
+  )
+  # conf
+  expect_error(
+  Spread(dat1, conf = 0.1),
+  "Parameter 'conf' must be one logical value."
+  )
+  # conf.lev
+  expect_error(
+  Spread(dat1, conf.lev = c(0.05, 0.95)),
+  "Parameter 'conf.lev' must be a numeric number between 0 and 1."
+  )
+  # ncores
+  expect_error(
+  Spread(dat1, ncores = 0),
+  "Parameter 'ncores' must be a positive integer."
+  )
+
+})
+
+##############################################
+test_that("2. dat1", {
+
+res1_1 <- Spread(dat1)
+
+expect_equal(
+names(res1_1),
+c('iqr', 'maxmin', 'sd', 'mad')
+)
+expect_equal(
+dim(res1_1$iqr),
+c(stats = 3, sdate = 4, ftime = 20)
+)
+expect_equal(
+dim(res1_1$maxmin),
+c(stats = 3, sdate = 4, ftime = 20)
+)
+expect_equal(
+dim(res1_1$sd),
+c(stats = 3, sdate = 4, ftime = 20)
+)
+expect_equal(
+dim(res1_1$mad),
+c(stats = 3, sdate = 4, ftime = 20)
+)
+
+res1_2 <- Spread(dat1, conf = F)
+
+expect_equal(
+names(res1_2),
+c('iqr', 'maxmin', 'sd', 'mad')
+)
+expect_equal(
+dim(res1_2$iqr),
+c(stats = 1, sdate = 4, ftime = 20)
+)
+expect_equal(
+dim(res1_2$iqr),
+dim(res1_2$maxmin)
+)
+expect_equal(
+dim(res1_2$iqr),
+dim(res1_2$sd)
+)
+expect_equal(
+dim(res1_2$iqr),
+dim(res1_2$mad)
+)
+
+})
+
+##############################################
+
+test_that("3. dat2", {
+
+res2_1 <- Spread(dat2)
+
+expect_equal(
+names(res2_1),
+c('iqr', 'maxmin', 'sd', 'mad')
+)
+expect_equal(
+dim(res2_1$iqr),
+c(stats = 3)
+)
+expect_equal(
+dim(res2_1$maxmin),
+c(stats = 3)
+)
+expect_equal(
+dim(res2_1$sd),
+c(stats = 3)
+)
+expect_equal(
+dim(res2_1$mad),
+c(stats = 3)
+)
+
+res2_2 <- Spread(dat2, conf = F)
+
+expect_equal(
+names(res2_2),
+c('iqr', 'maxmin', 'sd', 'mad')
+)
+expect_equal(
+dim(res2_2$iqr),
+c(stats = 1)
+)
+expect_equal(
+dim(res2_2$maxmin),
+c(stats = 1)
+)
+expect_equal(
+dim(res2_2$sd),
+c(stats = 1)
+)
+expect_equal(
+dim(res2_2$mad),
+c(stats = 1)
+)
+
+})
+
+
-- 
GitLab


From b5452b66cdb10b1bdaffc51e1516e7bbf448f5b7 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 15 Apr 2021 16:22:16 +0200
Subject: [PATCH 4/6] Correct example

---
 R/Spread.R    | 2 +-
 man/Spread.Rd | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/Spread.R b/R/Spread.R
index 9589fae..4b3bc6b 100644
--- a/R/Spread.R
+++ b/R/Spread.R
@@ -49,7 +49,7 @@
 #'smooth_ano_exp <- Smoothing(ano_exp, runmeanlen = runmean_months)
 #'smooth_ano_exp_m_sub <- smooth_ano_exp - InsertDim(MeanDims(smooth_ano_exp, 'member', 
 #'                                                            na.rm = TRUE), 
-#'                                                   posdim = 'member', 
+#'                                                   posdim = 3, 
 #'                                                   lendim = dim(smooth_ano_exp)['member'], 
 #'                                                   name = 'member')
 #'spread <- Spread(smooth_ano_exp_m_sub, compute_dim = c('member', 'sdate'))
diff --git a/man/Spread.Rd b/man/Spread.Rd
index e5f3edd..26e289e 100644
--- a/man/Spread.Rd
+++ b/man/Spread.Rd
@@ -69,7 +69,7 @@ runmean_months <- 12
 smooth_ano_exp <- Smoothing(ano_exp, runmeanlen = runmean_months)
 smooth_ano_exp_m_sub <- smooth_ano_exp - InsertDim(MeanDims(smooth_ano_exp, 'member', 
                                                            na.rm = TRUE), 
-                                                  posdim = 'member', 
+                                                  posdim = 3, 
                                                   lendim = dim(smooth_ano_exp)['member'], 
                                                   name = 'member')
 spread <- Spread(smooth_ano_exp_m_sub, compute_dim = c('member', 'sdate'))
-- 
GitLab


From 9299abbbad2930a764a99bf0ef340a32f1b7aad5 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 15 Apr 2021 17:02:33 +0200
Subject: [PATCH 5/6] Include PlotVsLTime.R

---
 NAMESPACE          |   1 +
 R/PlotVsLTime.R    | 265 +++++++++++++++++++++++++++++++++++++++++++++
 man/PlotVsLTime.Rd | 144 ++++++++++++++++++++++++
 3 files changed, 410 insertions(+)
 create mode 100644 R/PlotVsLTime.R
 create mode 100644 man/PlotVsLTime.Rd

diff --git a/NAMESPACE b/NAMESPACE
index ea090db..c9ac3a1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -39,6 +39,7 @@ export(PlotLayout)
 export(PlotMatrix)
 export(PlotSection)
 export(PlotStereoMap)
+export(PlotVsLTime)
 export(RMS)
 export(RMSSS)
 export(RandomWalkTest)
diff --git a/R/PlotVsLTime.R b/R/PlotVsLTime.R
new file mode 100644
index 0000000..94c82e0
--- /dev/null
+++ b/R/PlotVsLTime.R
@@ -0,0 +1,265 @@
+#'Plot a score along the forecast time with its confidence interval
+#'
+#'Plot the correlation (\code{Corr()}), the root mean square error 
+#'(\code{RMS()}) between the forecast values and their observational 
+#'counterpart, the slope of their trend (\code{Trend()}), the 
+#'InterQuartile range, maximum-mininum, standard deviation or median absolute 
+#'Deviation of the ensemble members (\code{Spread()}), or the ratio between 
+#'the ensemble spread and the RMSE of the ensemble mean (\code{RatioSDRMS()}) 
+#'along the forecast time for all the input experiments on the same figure 
+#'with their confidence intervals.
+#'
+#'@param var Matrix containing any Prediction Score with dimensions:\cr
+#'  (nexp/nmod, 3/4 ,nltime)\cr
+#'  or (nexp/nmod, nobs, 3/4 ,nltime).
+#'@param toptitle Main title, optional.
+#'@param ytitle Title of Y-axis, optional.
+#'@param monini Starting month between 1 and 12. Default = 1.
+#'@param freq 1 = yearly, 12 = monthly, 4 = seasonal, ... Default = 12.
+#'@param nticks Number of ticks and labels on the x-axis, optional.
+#'@param limits c(lower limit, upper limit): limits of the Y-axis, optional.
+#'@param listexp List of experiment names, optional.
+#'@param listobs List of observation names, optional.
+#'@param biglab TRUE/FALSE for presentation/paper plot. Default = FALSE.
+#'@param hlines c(a,b, ..) Add horizontal black lines at Y-positions a,b, ...\cr
+#'  Default = NULL.
+#'@param leg TRUE/FALSE if legend should be added or not to the plot. 
+#'  Default = TRUE.
+#'@param siglev TRUE/FALSE if significance level should replace confidence 
+#'  interval.\cr
+#'  Default = FALSE.
+#'@param sizetit Multiplicative factor to change title size, optional.
+#'@param show_conf TRUE/FALSE to show/not confidence intervals for input 
+#'  variables.
+#'@param fileout Name of output file. Extensions allowed: eps/ps, jpeg, png, 
+#'  pdf, bmp and tiff. The default value is NULL.
+#'@param width File width, in the units specified in the parameter size_units 
+#'  (inches by default). Takes 8 by default.
+#'@param height File height, in the units specified in the parameter 
+#'  size_units (inches by default). Takes 5 by default.
+#'@param size_units Units of the size of the device (file or window) to plot 
+#'  in. Inches ('in') by default. See ?Devices and the creator function of the 
+#'  corresponding device.
+#'@param res Resolution of the device (file or window) to plot in. See 
+#'  ?Devices and the creator function of the corresponding device.
+#'@param ... Arguments to be passed to the method. Only accepts the following
+#'  graphical parameters:\cr
+#'  adj ann ask bg bty cex.sub cin col.axis col.lab col.main col.sub cra crt 
+#'  csi cxy err family fg fig font font.axis font.lab font.main font.sub 
+#'  lheight ljoin lmitre mar mex mfcol mfrow mfg mkh oma omd omi page pch plt 
+#'  smo srt tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog \cr
+#'  For more information about the parameters see `par`.
+#'
+#'@details
+#'Examples of input:\cr
+#'Model and observed output from \code{Load()} then \code{Clim()} then 
+#'\code{Ano()} then \code{Smoothing()}:\cr
+#'(nmod, nmemb, nsdate, nltime) and (nobs, nmemb, nsdate, nltime)\cr
+#'then averaged over the members\cr
+#'\code{Mean1Dim(var_exp/var_obs, posdim = 2)}:\cr
+#'(nmod, nsdate, nltime) and (nobs, nsdate, nltime)\cr
+#'then passed through\cr
+#'  \code{Corr(exp, obs, posloop = 1, poscor = 2)} or\cr
+#'  \code{RMS(exp, obs, posloop = 1, posRMS = 2)}:\cr
+#'  (nmod, nobs, 3, nltime)\cr
+#'would plot the correlations or RMS between each exp & each obs as a function 
+#'of the forecast time.
+#'
+#'@examples
+#'# Load sample data as in Load() example:
+#'example(Load)
+#'clim <- Clim(sampleData$mod, sampleData$obs)
+#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+#'runmean_months <- 12
+#'smooth_ano_exp <- Smoothing(data = ano_exp, runmeanlen = runmean_months)
+#'smooth_ano_obs <- Smoothing(data = ano_obs, runmeanlen = runmean_months)
+#'dim_to_mean <- 'member'  # mean along members
+#'required_complete_row <- 'ftime'  # discard startdates for which there are NA leadtimes
+#'leadtimes_per_startdate <- 60
+#'corr <- Corr(MeanDims(smooth_ano_exp, dim_to_mean), 
+#'             MeanDims(smooth_ano_obs, dim_to_mean), 
+#'             comp_dim = required_complete_row, 
+#'             limits = c(ceiling((runmean_months + 1) / 2), 
+#'                        leadtimes_per_startdate - floor(runmean_months / 2)))
+#'# Combine corr results for plotting
+#'corr_combine <- abind::abind(corr$conf.lower, corr$corr, corr$conf.upper, corr$p.val, along = 0)
+#'corr_combine <- Reorder(corr_combine, c(2, 3, 1, 4))
+#'\donttest{
+#'PlotVsLTime(corr_combine, toptitle = "correlations", ytitle = "correlation", 
+#'            monini = 11, limits = c(-1, 2), listexp = c('CMIP5 IC3'), 
+#'            listobs = c('ERSST'), biglab = FALSE, hlines = c(-1, 0, 1))
+#'  }
+#'
+#'@importFrom grDevices dev.cur dev.new dev.off
+#'@importFrom stats ts
+#'@export
+PlotVsLTime <- function(var, toptitle = '', ytitle = '', monini = 1, freq = 12, 
+               nticks = NULL, limits = NULL, 
+               listexp = c('exp1', 'exp2', 'exp3'), 
+               listobs = c('obs1', 'obs2', 'obs3'), biglab = FALSE, hlines = NULL, 
+               leg = TRUE, siglev = FALSE, sizetit = 1, show_conf = TRUE,
+               fileout = NULL, 
+               width = 8, height = 5, size_units = 'in', res = 100, ...) {
+  # Process the user graphical parameters that may be passed in the call
+  ## Graphical parameters to exclude
+  excludedArgs <- c("cex", "cex.axis", "cex.lab", "cex.main", "col", "fin", "lab", "las", "lend", "lty", "lwd", "mai", "mgp", "new", "pin", "ps", "pty")
+  userArgs <- .FilterUserGraphicArgs(excludedArgs, ...)
+
+  # If there is any filenames to store the graphics, process them
+  # to select the right device 
+  if (!is.null(fileout)) {
+    deviceInfo <- .SelectDevice(fileout = fileout, width = width, height = height, units = size_units, res = res)
+    saveToFile <- deviceInfo$fun
+    fileout <- deviceInfo$files
+  }
+
+  #
+  #  Get some arguments
+  # ~~~~~~~~~~~~~~~~~~~~
+  #
+  if (length(dim(var)) == 3) {
+    var <- InsertDim(var, posdim = 2, lendim = 1)
+  } else if (length(dim(var)) != 4) {
+    stop("Parameter 'var' should have 3 or 4 dimensions: c(n. exp[, n. obs], 3/4, n. lead-times)")
+  }
+  nleadtime <- dim(var)[4]
+  nexp <- dim(var)[1]
+  nobs <- dim(var)[2]
+  if (is.null(limits) == TRUE) {
+	if (all(is.na(var > 0))) {
+      ll <- ul <- 0
+	} else {
+      ll <- min(var, na.rm = TRUE)
+      ul <- max(var, na.rm = TRUE)
+    }
+    if (biglab) {
+      ul <- ul + 0.4 * (ul - ll)
+    } else {
+      ul <- ul + 0.3 * (ul - ll)
+    }
+  } else {
+    ll <- limits[1]
+    ul <- limits[2]
+  }
+  lastyear <- (monini + (nleadtime - 1) * 12 / freq - 1) %/% 12
+  lastmonth <- (monini + (nleadtime - 1) * 12 / freq - 1) %% 12 + 1
+  empty_ts <- ts(start = c(0000, (monini - 1) %/% (12 / freq) + 1), 
+                 end = c(lastyear, (lastmonth - 1) %/% (12 / freq) + 1),
+                 frequency = freq)
+  empty <- array(dim = length(empty_ts))
+  #
+  #  Define some plot parameters
+  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  #
+  if (is.null(nticks)) {
+    if (biglab) {
+      nticks <- 5
+    } else {
+      nticks <- 10
+    }
+  }
+  labind <- seq(1, nleadtime, max(nleadtime %/% nticks, 1))
+  months <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep",
+              "Oct", "Nov", "Dec")
+  labyear <- ((labind - 1) * 12 / freq + monini - 1) %/% 12
+  labmonth <- months[((labind - 1) * 12 / freq + monini -1 ) %% 12 + 1]
+  for (jx in 1:length(labmonth)) {
+    y2o3dig <- paste("0", as.character(labyear[jx]), sep = "")
+    labmonth[jx] <- paste(labmonth[jx], "\nYr ", substr(y2o3dig, nchar(y2o3dig)
+                    - 1, nchar(y2o3dig)), sep = "")
+  }
+  color <- c("red1", "dodgerblue1", "green1", "orange1", "lightblue1", 
+             "deeppink1", "mediumpurple1", "lightgoldenrod1", "olivedrab1", 
+             "mediumorchid1")
+  type <- c(1, 3, 2, 4)
+  thickness <- array(dim = c(4, 4))
+  thickness[, 1] <- c(1, 2, 1, 1.5)
+  thickness[, 2] <- c(8, 12, 8, 10)
+  thickness[, 3] <- thickness[, 1]
+  thickness[, 4] <- c(4, 6, 4, 5)
+  if (siglev == TRUE) {
+    lines <- c("n", "l", "n", "l")
+  } else {
+    lines <- c("l", "l", "l", "n")
+  }
+  #
+  #  Define plot layout
+  # ~~~~~~~~~~~~~~~~~~~~
+  #
+
+  # Open connection to graphical device
+  if (!is.null(fileout)) {
+    saveToFile(fileout)
+  } else if (names(dev.cur()) == 'null device') {
+    dev.new(units = size_units, res = res, width = width, height = height)
+  }
+
+  # Load the user parameters
+  par(userArgs)
+  
+  if (biglab) {
+    par(mai = c(1.25, 1.4, 0.5, 1), mgp = c(4, 2.5, 0))
+    par(cex = 1.3, cex.lab = 2, cex.axis = 1.8)
+    cexmain <- 2.2
+    legsize <- 1.5
+  } else {
+    par(mai = c(1, 1.1, 0.5, 0), mgp = c(3, 1.8, 0))
+    par(cex = 1.3, cex.lab = 1.5, cex.axis = 1.1)
+    cexmain <- 1.5
+    legsize <- 1
+  }
+  plot(empty, ylim = c(ll, ul), xlab = "Time (months)", ylab = ytitle, 
+       main = toptitle, cex.main = cexmain*sizetit, axes = FALSE)
+  axis(1, at = labind, labels = labmonth)
+  axis(2)
+  box()
+  if (is.null(hlines) != TRUE) { 
+    for (jy in 1:length(hlines)) {
+      par(new = TRUE) 
+      abline(h = hlines[jy])
+    }
+  }
+  #
+  #  Loop on experimental & observational data
+  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  #                
+  legendnames <- array(dim = nobs * nexp)
+  legendthick <- array(dim = nobs * nexp)
+  legendsty <- array(dim = nobs * nexp)
+  legendcol <- array(dim = nobs * nexp)
+  ind <- 1
+  if (show_conf == TRUE) {
+    start_line <- dim(var)[3]
+    end_line <- 1
+  } else {
+    start_line <- 2
+    end_line <- 2
+  }
+  for (jt in seq(start_line, end_line, -1)) {
+    ind <- 1
+    for (jexp in 1:nexp) {
+      for (jobs in 1:nobs) {
+        par(new = TRUE)
+        plot(var[jexp, jobs, jt, ], type = lines[jt], ylim = c(ll, ul), 
+             col = color[jexp], lty = type[jobs], lwd = thickness[jobs, jt], 
+             ylab = "", xlab = "", axes = FALSE)
+        legendnames[ind] <- paste(listexp[jexp], 'vs', listobs[jobs])
+        legendthick[ind] <- thickness[jobs, 1] * 3
+        legendsty[ind] <- type[jobs]
+        legendcol[ind] <- color[jexp]
+        ind <- ind + 1
+      }
+    }
+  }
+  if (leg) {
+    if (nobs == 1) {
+      legendnames <- listexp[1:nexp]
+    }  
+    legend(1, ul, legendnames, lty = legendsty, lwd = legendthick,
+           col = legendcol, cex = legsize)
+  }
+  
+  # If the graphic was saved to file, close the connection with the device
+  if(!is.null(fileout)) dev.off()
+}
diff --git a/man/PlotVsLTime.Rd b/man/PlotVsLTime.Rd
new file mode 100644
index 0000000..05e2b42
--- /dev/null
+++ b/man/PlotVsLTime.Rd
@@ -0,0 +1,144 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PlotVsLTime.R
+\name{PlotVsLTime}
+\alias{PlotVsLTime}
+\title{Plot a score along the forecast time with its confidence interval}
+\usage{
+PlotVsLTime(
+  var,
+  toptitle = "",
+  ytitle = "",
+  monini = 1,
+  freq = 12,
+  nticks = NULL,
+  limits = NULL,
+  listexp = c("exp1", "exp2", "exp3"),
+  listobs = c("obs1", "obs2", "obs3"),
+  biglab = FALSE,
+  hlines = NULL,
+  leg = TRUE,
+  siglev = FALSE,
+  sizetit = 1,
+  show_conf = TRUE,
+  fileout = NULL,
+  width = 8,
+  height = 5,
+  size_units = "in",
+  res = 100,
+  ...
+)
+}
+\arguments{
+\item{var}{Matrix containing any Prediction Score with dimensions:\cr
+(nexp/nmod, 3/4 ,nltime)\cr
+or (nexp/nmod, nobs, 3/4 ,nltime).}
+
+\item{toptitle}{Main title, optional.}
+
+\item{ytitle}{Title of Y-axis, optional.}
+
+\item{monini}{Starting month between 1 and 12. Default = 1.}
+
+\item{freq}{1 = yearly, 12 = monthly, 4 = seasonal, ... Default = 12.}
+
+\item{nticks}{Number of ticks and labels on the x-axis, optional.}
+
+\item{limits}{c(lower limit, upper limit): limits of the Y-axis, optional.}
+
+\item{listexp}{List of experiment names, optional.}
+
+\item{listobs}{List of observation names, optional.}
+
+\item{biglab}{TRUE/FALSE for presentation/paper plot. Default = FALSE.}
+
+\item{hlines}{c(a,b, ..) Add horizontal black lines at Y-positions a,b, ...\cr
+Default = NULL.}
+
+\item{leg}{TRUE/FALSE if legend should be added or not to the plot. 
+Default = TRUE.}
+
+\item{siglev}{TRUE/FALSE if significance level should replace confidence 
+interval.\cr
+Default = FALSE.}
+
+\item{sizetit}{Multiplicative factor to change title size, optional.}
+
+\item{show_conf}{TRUE/FALSE to show/not confidence intervals for input 
+variables.}
+
+\item{fileout}{Name of output file. Extensions allowed: eps/ps, jpeg, png, 
+pdf, bmp and tiff. The default value is NULL.}
+
+\item{width}{File width, in the units specified in the parameter size_units 
+(inches by default). Takes 8 by default.}
+
+\item{height}{File height, in the units specified in the parameter 
+size_units (inches by default). Takes 5 by default.}
+
+\item{size_units}{Units of the size of the device (file or window) to plot 
+in. Inches ('in') by default. See ?Devices and the creator function of the 
+corresponding device.}
+
+\item{res}{Resolution of the device (file or window) to plot in. See 
+?Devices and the creator function of the corresponding device.}
+
+\item{...}{Arguments to be passed to the method. Only accepts the following
+graphical parameters:\cr
+adj ann ask bg bty cex.sub cin col.axis col.lab col.main col.sub cra crt 
+csi cxy err family fg fig font font.axis font.lab font.main font.sub 
+lheight ljoin lmitre mar mex mfcol mfrow mfg mkh oma omd omi page pch plt 
+smo srt tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog \cr
+For more information about the parameters see `par`.}
+}
+\description{
+Plot the correlation (\code{Corr()}), the root mean square error 
+(\code{RMS()}) between the forecast values and their observational 
+counterpart, the slope of their trend (\code{Trend()}), the 
+InterQuartile range, maximum-mininum, standard deviation or median absolute 
+Deviation of the ensemble members (\code{Spread()}), or the ratio between 
+the ensemble spread and the RMSE of the ensemble mean (\code{RatioSDRMS()}) 
+along the forecast time for all the input experiments on the same figure 
+with their confidence intervals.
+}
+\details{
+Examples of input:\cr
+Model and observed output from \code{Load()} then \code{Clim()} then 
+\code{Ano()} then \code{Smoothing()}:\cr
+(nmod, nmemb, nsdate, nltime) and (nobs, nmemb, nsdate, nltime)\cr
+then averaged over the members\cr
+\code{Mean1Dim(var_exp/var_obs, posdim = 2)}:\cr
+(nmod, nsdate, nltime) and (nobs, nsdate, nltime)\cr
+then passed through\cr
+ \code{Corr(exp, obs, posloop = 1, poscor = 2)} or\cr
+ \code{RMS(exp, obs, posloop = 1, posRMS = 2)}:\cr
+ (nmod, nobs, 3, nltime)\cr
+would plot the correlations or RMS between each exp & each obs as a function 
+of the forecast time.
+}
+\examples{
+# Load sample data as in Load() example:
+example(Load)
+clim <- Clim(sampleData$mod, sampleData$obs)
+ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+runmean_months <- 12
+smooth_ano_exp <- Smoothing(data = ano_exp, runmeanlen = runmean_months)
+smooth_ano_obs <- Smoothing(data = ano_obs, runmeanlen = runmean_months)
+dim_to_mean <- 'member'  # mean along members
+required_complete_row <- 'ftime'  # discard startdates for which there are NA leadtimes
+leadtimes_per_startdate <- 60
+corr <- Corr(MeanDims(smooth_ano_exp, dim_to_mean), 
+            MeanDims(smooth_ano_obs, dim_to_mean), 
+            comp_dim = required_complete_row, 
+            limits = c(ceiling((runmean_months + 1) / 2), 
+                       leadtimes_per_startdate - floor(runmean_months / 2)))
+# Combine corr results for plotting
+corr_combine <- abind::abind(corr$conf.lower, corr$corr, corr$conf.upper, corr$p.val, along = 0)
+corr_combine <- Reorder(corr_combine, c(2, 3, 1, 4))
+\donttest{
+PlotVsLTime(corr_combine, toptitle = "correlations", ytitle = "correlation", 
+           monini = 11, limits = c(-1, 2), listexp = c('CMIP5 IC3'), 
+           listobs = c('ERSST'), biglab = FALSE, hlines = c(-1, 0, 1))
+ }
+
+}
-- 
GitLab


From 654aae7b7e12074cd91a02fca897376be5812ac5 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 15 Apr 2021 17:41:22 +0200
Subject: [PATCH 6/6] Include Plot2VarsVsLTime.R

---
 NAMESPACE               |   1 +
 R/Plot2VarsVsLTime.R    | 256 ++++++++++++++++++++++++++++++++++++++++
 man/Plot2VarsVsLTime.Rd | 136 +++++++++++++++++++++
 3 files changed, 393 insertions(+)
 create mode 100644 R/Plot2VarsVsLTime.R
 create mode 100644 man/Plot2VarsVsLTime.Rd

diff --git a/NAMESPACE b/NAMESPACE
index c9ac3a1..19648b8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -31,6 +31,7 @@ export(LeapYear)
 export(Load)
 export(MeanDims)
 export(Persistence)
+export(Plot2VarsVsLTime)
 export(PlotACC)
 export(PlotAno)
 export(PlotClim)
diff --git a/R/Plot2VarsVsLTime.R b/R/Plot2VarsVsLTime.R
new file mode 100644
index 0000000..1c784dd
--- /dev/null
+++ b/R/Plot2VarsVsLTime.R
@@ -0,0 +1,256 @@
+#'Plot two scores with confidence intervals in a common plot
+#'
+#'Plot two input variables that have the same dimensions in a common plot.
+#'One plot for all experiments.
+#'The input variables should have dimensions (nexp/nmod, nltime).
+#'
+#'@param var1 Matrix of dimensions (nexp/nmod, nltime).
+#'@param var2 Matrix of dimensions (nexp/nmod, nltime).
+#'@param toptitle Main title, optional.
+#'@param ytitle Title of Y-axis, optional.
+#'@param monini Starting month between 1 and 12. Default = 1.
+#'@param freq 1 = yearly, 12 = monthly, 4 = seasonal, ... Default = 12.
+#'@param nticks Number of ticks and labels on the x-axis, optional.
+#'@param limits c(lower limit, upper limit): limits of the Y-axis, optional.
+#'@param listexp List of experiment names, up to three, optional.
+#'@param listvars List of names of input variables, optional.
+#'@param biglab TRUE/FALSE for presentation/paper plot. Default = FALSE.
+#'@param hlines c(a, b, ...) Add horizontal black lines at Y-positions a, b, 
+#'  ... The default value is NULL.
+#'@param leg TRUE/FALSE if legend should be added or not to the plot. 
+#'  Default = TRUE.
+#'@param siglev TRUE/FALSE if significance level should replace confidence 
+#'  interval.\cr 
+#'  Default = FALSE.
+#'@param sizetit Multiplicative factor to change title size, optional.
+#'@param show_conf TRUE/FALSE to show/not confidence intervals for input 
+#'  variables.
+#'@param fileout Name of output file. Extensions allowed: eps/ps, jpeg, png, 
+#'  pdf, bmp and tiff. The default value is NULL.
+#'@param width File width, in the units specified in the parameter size_units 
+#'  (inches by default). Takes 8 by default.
+#'@param height File height, in the units specified in the parameter 
+#'  size_units (inches by default). Takes 5 by default.
+#'@param size_units Units of the size of the device (file or window) to plot 
+#'  in. Inches ('in') by default. See ?Devices and the creator function of the 
+#'  corresponding device.
+#'@param res Resolution of the device (file or window) to plot in. See 
+#'  ?Devices and the creator function of the corresponding device.
+#'@param ... Arguments to be passed to the method. Only accepts the following
+#'  graphical parameters:\cr  
+#'  adj ann ask bg bty cex.sub cin col.axis col.lab col.main col.sub cra crt 
+#'  csi cxy err family fg fig font font.axis font.lab font.main font.sub lend 
+#'  lheight ljoin lmitre mar mex mfcol mfrow mfg mkh oma omd omi page pch plt 
+#'  smo srt tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog \cr
+#'  For more information about the parameters see `par`.
+#'
+#'@details
+#'Examples of input:\cr
+#'------------------\cr\cr  
+#'RMSE error for a number of experiments and along lead-time: (nexp, nltime)
+#'
+#'@examples
+#'# Load sample data as in Load() example:
+#'example(Load)
+#'clim <- Clim(sampleData$mod, sampleData$obs)
+#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+#'runmean_months <- 12
+#'smooth_ano_exp <- Smoothing(data = ano_exp, runmeanlen = runmean_months)
+#'smooth_ano_obs <- Smoothing(data = ano_obs, runmeanlen = runmean_months)
+#'dim_to_mean <- 'member'  # mean along members
+#'required_complete_row <- 'ftime'  # discard startdates for which there are NA leadtimes
+#'leadtimes_per_startdate <- 60
+#'rms <- RMS(MeanDims(smooth_ano_exp, dim_to_mean), 
+#'           MeanDims(smooth_ano_obs, dim_to_mean), 
+#'           comp_dim = required_complete_row, 
+#'           limits = c(ceiling((runmean_months + 1) / 2), 
+#'                      leadtimes_per_startdate - floor(runmean_months / 2)))
+#'smooth_ano_exp_m_sub <- smooth_ano_exp - InsertDim(MeanDims(smooth_ano_exp, 'member', 
+#'                                                            na.rm = TRUE), 
+#'                                                   posdim = 3, 
+#'                                                   lendim = dim(smooth_ano_exp)['member'], 
+#'                                                   name = 'member')
+#'spread <- Spread(smooth_ano_exp_m_sub, compute_dim = c('member', 'sdate'))
+#'#Combine rms outputs into one array 
+#'rms_combine <- abind::abind(rms$conf.lower, rms$rms, rms$conf.upper, along = 0)
+#'rms_combine <- Reorder(rms_combine, c(2, 3, 1, 4))
+#'  \donttest{
+#'Plot2VarsVsLTime(InsertDim(rms_combine[, , , ], 1, 1), Reorder(spread$sd, c(1, 3, 2)), 
+#'                 toptitle = 'RMSE and spread', monini = 11, freq = 12, 
+#'                 listexp = c('CMIP5 IC3'), listvar = c('RMSE', 'spread'))
+#'  }
+#'
+#'@importFrom grDevices png jpeg postscript pdf svg bmp tiff postscript dev.cur dev.new dev.off
+#'@importFrom stats ts 
+#'@export
+Plot2VarsVsLTime <- function(var1, var2, toptitle = '', ytitle = '', monini = 1,
+                             freq = 12, nticks = NULL, limits = NULL, listexp =
+                             c('exp1', 'exp2', 'exp3'), listvars = c('var1',
+                             'var2'), biglab = FALSE, hlines = NULL, leg = TRUE,
+                             siglev = FALSE, sizetit = 1, show_conf = TRUE,
+                             fileout = NULL, 
+                             width = 8, height = 5, size_units = 'in', res = 100, ...) {
+  # Process the user graphical parameters that may be passed in the call
+  ## Graphical parameters to exclude
+  excludedArgs <- c("cex", "cex.axis", "cex.lab", "cex.main", "col", "fin", "lab", "las", "lty", "lwd", "mai", "mgp", "new", "pin", "ps", "pty")
+  userArgs <- .FilterUserGraphicArgs(excludedArgs, ...)
+
+  # If there is any filenames to store the graphics, process them
+  # to select the right device 
+  if (!is.null(fileout)) {
+    deviceInfo <- .SelectDevice(fileout = fileout, width = width, height = height, units = size_units, res = res)
+    saveToFile <- deviceInfo$fun
+    fileout <- deviceInfo$files
+  }
+
+  #
+  nvars <- 2
+
+  if (length(dim(var1)) != length(dim(var2))) { 
+    print("the two input variables should have the same dimensions")
+    stop()
+  }
+  if (length(dim(var1)) >= 4) { 
+    print("dimensions of input variables should be 3")
+    stop()
+  }
+  nleadtime <- dim(var1)[3]
+  nexp <- dim(var1)[1]
+  var <- array(dim = c(nvars, nexp, 3, nleadtime))
+  for (jvar in 1:nvars) {
+    varname <- paste("var", as.character(jvar), sep = "")
+    var[jvar, , , ] <- get(varname)
+    rm(varname)
+  }
+
+  if (is.null(limits) == TRUE) {
+    ll <- min(var1, na.rm = TRUE)
+    ul <- max(var1, na.rm = TRUE)
+    if (biglab) {
+      ul <- ul + 0.4 * (ul - ll)
+    } else {
+      ul <- ul + 0.3 * (ul - ll)
+    }
+  } else {
+    ll <- limits[1]
+    ul <- limits[2]
+  }
+  lastyear <- (monini + (nleadtime - 1) * 12 / freq - 1) %/% 12
+  lastmonth <- (monini + (nleadtime - 1) * 12 / freq - 1) %% 12 + 1
+  empty_ts <- ts(start = c(0000, (monini - 1) %/% (12 / freq) + 1), 
+                 end = c(lastyear, (lastmonth - 1) %/% (12 / freq) + 1), 
+                 frequency = freq)
+  empty <- array(dim = length(empty_ts))
+  #
+  #  Define some plot parameters
+  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  #
+  if (is.null(nticks)) {
+    if (biglab) {
+      nticks <- 5
+    } else {
+      nticks <- 10
+    }
+  }
+  labind <- seq(1, nleadtime, max(nleadtime %/% nticks, 1))
+  months <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep",
+              "Oct", "Nov", "Dec")
+  labyear <- ((labind - 1) * 12 / freq + monini - 1) %/% 12
+  labmonth <- months[((labind - 1) * 12 / freq + monini - 1) %% 12 + 1]
+  for (jx in 1:length(labmonth)) {
+    y2o3dig <- paste("0", as.character(labyear[jx]), sep = "")
+    labmonth[jx] <- paste(labmonth[jx], "\nYr ", substr(y2o3dig, nchar(y2o3dig)
+                    - 1, nchar(y2o3dig)), sep = "")
+  }
+  color <- c("red1", "dodgerblue1", "green1", "orange1", "lightblue1", 
+             "deeppink1", "mediumpurple1", "lightgoldenrod1", "olivedrab1",
+             "mediumorchid1")
+  type <- c(1, 3)
+  if (siglev == TRUE) {
+    lines <- c("n", "l", "n")
+  }
+  else{
+    lines <- c("l", "l", "l")
+  }
+  thickness <- array(dim = c(3))
+  thickness[1] <- c(1)
+  thickness[2] <- c(8)
+  thickness[3] <- thickness[1]
+  
+  #
+  #  Define plot layout
+  # ~~~~~~~~~~~~~~~~~~~~
+  #
+
+  # Open connection to graphical device
+  if (!is.null(fileout)) {
+    saveToFile(fileout)
+  } else if (names(dev.cur()) == 'null device') {
+    dev.new(units = size_units, res = res, width = width, height = height)
+  }
+
+  # Load the user parameters
+  par(userArgs)
+
+  if (biglab) {
+    par(mai = c(1.25, 1.4, 0.5, 1), mgp = c(4, 2.5, 0))
+    par(cex = 1.3, cex.lab = 2, cex.axis = 1.8)
+    cexmain <- 2.2
+    legsize <- 1.5
+  } else {
+    par(mai = c(1, 1.1, 0.5, 0), mgp = c(3, 1.8, 0))
+    par(cex = 1.3, cex.lab = 1.5, cex.axis = 1.1)
+    cexmain <- 1.5
+    legsize <- 1
+  }
+  plot(empty, ylim = c(ll, ul), xlab = "Time (months)", ylab = ytitle, 
+       main = toptitle, cex.main = cexmain * sizetit, axes = FALSE)
+  axis(1, at = labind, labels = labmonth)
+  axis(2)
+  box()
+  if (is.null(hlines) != TRUE) { 
+    for (jy in 1:length(hlines)) {
+      par(new = TRUE) 
+      abline(h = hlines[jy])
+    }
+  }
+  #
+  #  Loop on experimental data
+  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  #                 
+  legendnames <- array(dim = nexp * nvars)
+  legendthick <- array(dim = nexp * nvars)
+  legendsty <- array(dim = nexp * nvars)
+  legendcol <- array(dim = nexp * nvars)
+  if (show_conf == TRUE) {
+    start_line <- 3
+    end_line <- 1
+  } else {
+    start_line <- 2
+    end_line <- 2
+  }
+  for (jint in seq(start_line, end_line, -1)) {
+    ind <- 1
+    for (jexp in 1:nexp) {
+      for (jvar in 1:nvars) {
+        par(new = TRUE)
+        plot(var[jvar, jexp, jint, ], type = lines[jint], ylim = c(ll, ul), 
+             col = color[jexp], lty = type[jvar], lwd = thickness[jint],
+             ylab = "", xlab = "", axes = FALSE)
+        legendnames[ind] <- paste(listexp[jexp], listvars[jvar])
+        legendthick[ind] <- 2
+        legendsty[ind] <- type[jvar]
+        legendcol[ind] <- color[jexp]
+        ind <- ind + 1
+      }
+    }
+  }
+  if (leg) {
+    legend(1, ul, legendnames, lty = legendsty, lwd = legendthick,
+           col = legendcol, cex = legsize)
+  }
+
+  # If the graphic was saved to file, close the connection with the device
+  if(!is.null(fileout)) dev.off()
+}
diff --git a/man/Plot2VarsVsLTime.Rd b/man/Plot2VarsVsLTime.Rd
new file mode 100644
index 0000000..46b9cd5
--- /dev/null
+++ b/man/Plot2VarsVsLTime.Rd
@@ -0,0 +1,136 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Plot2VarsVsLTime.R
+\name{Plot2VarsVsLTime}
+\alias{Plot2VarsVsLTime}
+\title{Plot two scores with confidence intervals in a common plot}
+\usage{
+Plot2VarsVsLTime(
+  var1,
+  var2,
+  toptitle = "",
+  ytitle = "",
+  monini = 1,
+  freq = 12,
+  nticks = NULL,
+  limits = NULL,
+  listexp = c("exp1", "exp2", "exp3"),
+  listvars = c("var1", "var2"),
+  biglab = FALSE,
+  hlines = NULL,
+  leg = TRUE,
+  siglev = FALSE,
+  sizetit = 1,
+  show_conf = TRUE,
+  fileout = NULL,
+  width = 8,
+  height = 5,
+  size_units = "in",
+  res = 100,
+  ...
+)
+}
+\arguments{
+\item{var1}{Matrix of dimensions (nexp/nmod, nltime).}
+
+\item{var2}{Matrix of dimensions (nexp/nmod, nltime).}
+
+\item{toptitle}{Main title, optional.}
+
+\item{ytitle}{Title of Y-axis, optional.}
+
+\item{monini}{Starting month between 1 and 12. Default = 1.}
+
+\item{freq}{1 = yearly, 12 = monthly, 4 = seasonal, ... Default = 12.}
+
+\item{nticks}{Number of ticks and labels on the x-axis, optional.}
+
+\item{limits}{c(lower limit, upper limit): limits of the Y-axis, optional.}
+
+\item{listexp}{List of experiment names, up to three, optional.}
+
+\item{listvars}{List of names of input variables, optional.}
+
+\item{biglab}{TRUE/FALSE for presentation/paper plot. Default = FALSE.}
+
+\item{hlines}{c(a, b, ...) Add horizontal black lines at Y-positions a, b, 
+... The default value is NULL.}
+
+\item{leg}{TRUE/FALSE if legend should be added or not to the plot. 
+Default = TRUE.}
+
+\item{siglev}{TRUE/FALSE if significance level should replace confidence 
+interval.\cr 
+Default = FALSE.}
+
+\item{sizetit}{Multiplicative factor to change title size, optional.}
+
+\item{show_conf}{TRUE/FALSE to show/not confidence intervals for input 
+variables.}
+
+\item{fileout}{Name of output file. Extensions allowed: eps/ps, jpeg, png, 
+pdf, bmp and tiff. The default value is NULL.}
+
+\item{width}{File width, in the units specified in the parameter size_units 
+(inches by default). Takes 8 by default.}
+
+\item{height}{File height, in the units specified in the parameter 
+size_units (inches by default). Takes 5 by default.}
+
+\item{size_units}{Units of the size of the device (file or window) to plot 
+in. Inches ('in') by default. See ?Devices and the creator function of the 
+corresponding device.}
+
+\item{res}{Resolution of the device (file or window) to plot in. See 
+?Devices and the creator function of the corresponding device.}
+
+\item{...}{Arguments to be passed to the method. Only accepts the following
+graphical parameters:\cr  
+adj ann ask bg bty cex.sub cin col.axis col.lab col.main col.sub cra crt 
+csi cxy err family fg fig font font.axis font.lab font.main font.sub lend 
+lheight ljoin lmitre mar mex mfcol mfrow mfg mkh oma omd omi page pch plt 
+smo srt tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog \cr
+For more information about the parameters see `par`.}
+}
+\description{
+Plot two input variables that have the same dimensions in a common plot.
+One plot for all experiments.
+The input variables should have dimensions (nexp/nmod, nltime).
+}
+\details{
+Examples of input:\cr
+------------------\cr\cr  
+RMSE error for a number of experiments and along lead-time: (nexp, nltime)
+}
+\examples{
+# Load sample data as in Load() example:
+example(Load)
+clim <- Clim(sampleData$mod, sampleData$obs)
+ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+runmean_months <- 12
+smooth_ano_exp <- Smoothing(data = ano_exp, runmeanlen = runmean_months)
+smooth_ano_obs <- Smoothing(data = ano_obs, runmeanlen = runmean_months)
+dim_to_mean <- 'member'  # mean along members
+required_complete_row <- 'ftime'  # discard startdates for which there are NA leadtimes
+leadtimes_per_startdate <- 60
+rms <- RMS(MeanDims(smooth_ano_exp, dim_to_mean), 
+          MeanDims(smooth_ano_obs, dim_to_mean), 
+          comp_dim = required_complete_row, 
+          limits = c(ceiling((runmean_months + 1) / 2), 
+                     leadtimes_per_startdate - floor(runmean_months / 2)))
+smooth_ano_exp_m_sub <- smooth_ano_exp - InsertDim(MeanDims(smooth_ano_exp, 'member', 
+                                                           na.rm = TRUE), 
+                                                  posdim = 3, 
+                                                  lendim = dim(smooth_ano_exp)['member'], 
+                                                  name = 'member')
+spread <- Spread(smooth_ano_exp_m_sub, compute_dim = c('member', 'sdate'))
+#Combine rms outputs into one array 
+rms_combine <- abind::abind(rms$conf.lower, rms$rms, rms$conf.upper, along = 0)
+rms_combine <- Reorder(rms_combine, c(2, 3, 1, 4))
+ \donttest{
+Plot2VarsVsLTime(InsertDim(rms_combine[, , , ], 1, 1), Reorder(spread$sd, c(1, 3, 2)), 
+                toptitle = 'RMSE and spread', monini = 11, freq = 12, 
+                listexp = c('CMIP5 IC3'), listvar = c('RMSE', 'spread'))
+ }
+
+}
-- 
GitLab