diff --git a/.Rbuildignore b/.Rbuildignore
index f6430a4587524c56d276b040b1318e3561853bcf..fef43e49c022fe079bca1257de840ce480dbe2f8 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -2,6 +2,8 @@
 .*\.gitignore$
 .*\.tar.gz$
 .*\.pdf$
+.*\.png$
+.*\.gif$
 .*^(?!inst)\.nc$
 sample_data
 README\.Rmd$
diff --git a/DESCRIPTION b/DESCRIPTION
index 86859f21f36e58d19e966fc6f68d86c7197f7a42..5327ec3f7a136ba315e6d41022457f286285f7e4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -16,14 +16,12 @@ Imports:
     grDevices,
     maps,
     mapproj,
-    methods,
     ClimProjDiags,
     sf,
     ggplot2,
     rnaturalearth,
     cowplot,
     s2dv,
-    multiApply,
     RColorBrewer,
     CSTools,
     easyNCDF,
@@ -35,7 +33,8 @@ Imports:
     stats,
     utils,
     foreach,
-    doParallel
+    doParallel,
+    rlang
 Suggests:
     testthat
 License: GPL-3
diff --git a/NAMESPACE b/NAMESPACE
index 9139e77251726d2938431ba60c824aff6f61d8e5..6601969aa0e3768771c9f42eb7a5a54969b33828 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,24 +5,15 @@ export(ClimPalette)
 export(ColorBarContinuous)
 export(ColorBarDiscrete)
 export(ShapeToMask)
-export(Viz2VarsVsLTime)
-export(VizACC)
-export(VizAnimateMap)
-export(VizAno)
-export(VizBoxWhisker)
-export(VizClim)
 export(VizCombinedMap)
 export(VizEquiMap)
 export(VizForecastPDF)
 export(VizLayout)
-export(VizMatrix)
 export(VizMostLikelyQuantileMap)
 export(VizPDFsOLE)
 export(VizRobinson)
-export(VizSection)
 export(VizStereoMap)
 export(VizTriangles4Categories)
-export(VizVsLTime)
 export(VizWeeklyClim)
 import(RColorBrewer)
 import(cowplot)
@@ -79,6 +70,7 @@ importFrom(maps,map)
 importFrom(plyr,.)
 importFrom(plyr,dlply)
 importFrom(reshape2,melt)
+importFrom(rlang,.data)
 importFrom(s2dv,InsertDim)
 importFrom(s2dv,MeanDims)
 importFrom(s2dv,Reorder)
diff --git a/R/Viz2VarsVsLTime.R b/R/Viz2VarsVsLTime.R
index 47aae43ce632c5e7c6f9b9da17e182363ee3b5d6..f4a634585ae0fee15c6b599d3a9543f612884320 100644
--- a/R/Viz2VarsVsLTime.R
+++ b/R/Viz2VarsVsLTime.R
@@ -63,7 +63,7 @@
 #'input_rms[, 1, ] <- rms_ano$conf.lower[, 1, 1, ]
 #'input_rms[, 2, ] <- rms_ano$rms[, 1, 1, ]
 #'input_rms[, 3, ] <- rms_ano$conf.upper[, 1, 1, ]
-#'Viz2VarsVsLTime(input_cor, input_rms,
+#'esviz:::Viz2VarsVsLTime(input_cor, input_rms,
 #'                toptitle = "Time correlation and RMSE with ERA5",
 #'                ytitle = "K", sizetit = 0.7,
 #'                monini = 11, freq = 1, limits = c(-1, 5),
@@ -72,7 +72,6 @@
 #'
 #'@importFrom grDevices png jpeg postscript pdf svg bmp tiff postscript dev.cur dev.new dev.off
 #'@importFrom stats ts 
-#'@export
 Viz2VarsVsLTime <- function(var1, var2, toptitle = '', ytitle = '', monini = 1,
                              freq = 12, nticks = NULL, limits = NULL, listexp =
                              c('exp1', 'exp2', 'exp3'), listvars = c('var1',
diff --git a/R/VizACC.R b/R/VizACC.R
index a70a980a60babd5846285a417d101711124c0937..2787383b4f5a65b5cca08bf2a04669923195e2fb 100644
--- a/R/VizACC.R
+++ b/R/VizACC.R
@@ -49,6 +49,7 @@
 #'  For more information about the parameters see `par`.
 #'
 #'@examples
+#'  \dontrun{
 #'ano <- s2dv::Ano_CrossValid(map_temp$exp, map_temp$obs, memb = FALSE, 
 #'                            dat_dim = c('dat', 'member'), memb_dim = 'member')
 #'lats <- attr(map_temp$exp, "Variables")$common$lat
@@ -61,14 +62,13 @@
 #'input_acc[, , , , 3] <- acc$conf.upper
 #'input_acc[, , , , 4] <- acc$p.val
 #'sdates <- paste0(2000:2005, '1101')
-#'VizACC(input_acc, sdates, 
+#'esviz:::VizACC(input_acc, sdates, 
 #'       toptitle = "Spatial anomaly corr. coeff. with ERA5",
 #'       ytitle = "K", sizetit = 0.7, freq = 12,
 #'       legends = 'SEAS5', fileout = NULL)
-#'
+#'  }
 #'@importFrom grDevices dev.cur dev.new dev.off 
 #'@importFrom stats ts
-#'@export
 VizACC <- function(ACC, sdates, toptitle = "", sizetit = 1, ytitle = "", 
                     limits = NULL, legends = NULL, freq = 12, biglab = FALSE, 
                     fill = FALSE, linezero = FALSE, points = TRUE, vlines = NULL, 
diff --git a/R/VizAnimateMap.R b/R/VizAnimateMap.R
index 748c1057f2297afd3eea5c1e932e94874284c677..061ca1e9604767c56e458f096947f6e82e36c172 100644
--- a/R/VizAnimateMap.R
+++ b/R/VizAnimateMap.R
@@ -94,15 +94,14 @@
 #'                   dat_dim = c('dat', 'member'), memb_dim = 'member')
 #'lats <- attr(map_temp$exp, "Variables")$common$lat
 #'lons <- attr(map_temp$exp, "Variables")$common$lon
-#'VizAnimateMap(clim$clim_exp[1, 1, , , ], lon = lons, lat = lats,
+#'  \dontrun{
+#'esviz:::VizAnimateMap(clim$clim_exp[1, 1, , , ], lon = lons, lat = lats,
 #'  toptitle = "climatology of decadal prediction", sizetit = 1, 
 #'  units = "K", brks = seq(270, 300, 3), monini = 11, freq = 12, 
-#'  msk95lev = FALSE, filled.continents = FALSE, intlon = 10, intlat = 10,
-#'  fileout = NULL)
-#'
+#'  msk95lev = FALSE, filled.continents = FALSE, intlon = 10, intlat = 10)
+#'  }
 #'@importFrom grDevices postscript dev.off
-#' @importFrom s2dv InsertDim
-#'@export
+#'@importFrom s2dv InsertDim
 VizAnimateMap <- function(var, lon, lat, toptitle = rep("", 11), sizetit = 1, 
                        units = "", monini = 1, freq = 12, msk95lev = FALSE, 
                        brks = NULL, cols = NULL, filled.continents = FALSE, 
diff --git a/R/VizAno.R b/R/VizAno.R
index 86d6c33a848fdad567ed8e262443aaa25efc9a0a..7a034209d4d5cb4c6926484286b5519ca9a037cd 100644
--- a/R/VizAno.R
+++ b/R/VizAno.R
@@ -55,13 +55,12 @@
 #'dim(ts_temp$exp) <- c(dat = 1, member = 3, sdate = 6, time = 5)
 #'dim(ts_temp$obs) <- c(dat = 1, member = 1, sdate = 6, time = 5)
 #'sdates <- paste0(2000:2005, '1101')
-#'VizAno(ts_temp$exp, ts_temp$obs, sdates,
+#'esviz:::VizAno(ts_temp$exp, ts_temp$obs, sdates,
 #'       toptitle = "Raw 'tas'", ytitle = 'K', legends = 'ERA5', biglab = FALSE)
 #'
 #'@importFrom grDevices dev.cur dev.new dev.off 
 #'@importFrom stats ts
 #'@importFrom s2dv MeanDims
-#'@export
 VizAno <- function(exp_ano, obs_ano = NULL, sdates, toptitle = rep('', 15),
                    ytitle = rep('', 15), limits = NULL, legends = NULL, 
                    freq = 12, biglab = FALSE, fill = TRUE, memb = TRUE, 
diff --git a/R/VizBoxWhisker.R b/R/VizBoxWhisker.R
index 09dd252a03f2c7076b386805194b4d5bafcb578a..ec1d8d21be5eef12dfe32ee5344b4ac759643dfd 100644
--- a/R/VizBoxWhisker.R
+++ b/R/VizBoxWhisker.R
@@ -50,6 +50,7 @@
 #'@examples
 #'# No example data is available over NAO region, so in this example we will
 #'# tweak the longitude and latitude.
+#'  \dontrun{
 #'ano <- s2dv::Ano_CrossValid(map_temp$exp, map_temp$obs, memb = FALSE,
 #'                            dat_dim = c('dat', 'member'), memb_dim = 'member')
 #'nao <- s2dv::NAO(ano$exp, ano$obs, lat = seq(20, 80, length.out = 11),
@@ -57,13 +58,12 @@
 #'                 ftime_dim = "time")
 #'nao$exp <- drop(aperm(nao$exp, c(2, 1, 3, 4)))
 #'nao$obs <- drop(nao$obs)
-#'VizBoxWhisker(nao$exp, nao$obs, toptitle = "NAO index", 
+#'esviz:::VizBoxWhisker(nao$exp, nao$obs, toptitle = "NAO index", 
 #'              ytitle = "NAO index (PC1) TOS", monini = 11, freq = 1,
 #'               yearini = 2000, expname = "SEAS5", obsname = "ERA5")
-#'
+#'  }
 #'@importFrom grDevices dev.cur dev.new dev.off 
 #'@importFrom stats cor
-#'@export
 VizBoxWhisker <- function(exp, obs, toptitle = '', ytitle = '', monini = 1, 
                            yearini = 0, freq = 1, expname = "exp 1", 
                            obsname = "obs 1", drawleg = TRUE,
diff --git a/R/VizClim.R b/R/VizClim.R
index 412122bf3b5e25c3003a69f2b469a484015f18a0..84dcd4616a98241f2131e60a8cc731ebba8e88b9 100644
--- a/R/VizClim.R
+++ b/R/VizClim.R
@@ -46,14 +46,13 @@
 #'                   dat_dim = c("dat", "member"))
 #'dim(clim$clim_exp) <- dim(clim$clim_exp)[-3]
 #'dim(clim$clim_obs) <- dim(clim$clim_obs)[-2]
-#'VizClim(clim$clim_exp, clim$clim_obs, toptitle = 'Climatologies',
+#'esviz:::VizClim(clim$clim_exp, clim$clim_obs, toptitle = 'Climatologies',
 #'        ytitle = 'K', monini = 11, listexp = 'SEAS5',
 #'        listobs = 'ERA5', biglab = FALSE, fileout = NULL)
 #'
 #'@importFrom grDevices dev.cur dev.new dev.off
 #'@importFrom stats ts
 #'@importFrom s2dv InsertDim
-#'@export
 VizClim <- function(exp_clim, obs_clim = NULL, toptitle = '', ytitle = '',
                     monini = 1, freq = 12, limits = NULL,
                     listexp = c('exp1', 'exp2', 'exp3'),
diff --git a/R/VizEquiMap.R b/R/VizEquiMap.R
index deafcfbf7060f7cf28f9db6c485b17ead1d44098..7449f7ce532ca09b1d5fc491eff99b3e7a7d571b 100644
--- a/R/VizEquiMap.R
+++ b/R/VizEquiMap.R
@@ -225,6 +225,7 @@
 #'}
 #'
 #'@examples
+#'  \dontrun{
 #'ano <- s2dv::Ano_CrossValid(map_temp$exp, map_temp$obs, memb = FALSE,
 #'                            dat_dim = c('dat', 'member'), memb_dim = 'member')
 #'var <- s2dv::MeanDims(ano$exp, "member")
@@ -234,7 +235,7 @@
 #'VizEquiMap(var[1, 1, 1, 1, , ], lon = lons, lat = lats,
 #'           toptitle = 'Near-surface temperature anomaly, Nov. 2000',
 #'           filled.continents = FALSE, title_scale = 0.7)
-#'
+#'  }
 #'@import graphics maps utils
 #'@importFrom grDevices dev.cur dev.new dev.off gray
 #'@importFrom stats cor
diff --git a/R/VizLayout.R b/R/VizLayout.R
index c9984f76989871eac0e17adcd7a10710ee229f36..0e4e64128b17b84e12e34fc6aea46f6d3aa909bc 100644
--- a/R/VizLayout.R
+++ b/R/VizLayout.R
@@ -179,6 +179,7 @@
 #'}
 #'
 #'@examples
+#'  \dontrun{
 #'ano <- s2dv::Ano_CrossValid(map_temp$exp, map_temp$obs, memb = FALSE,
 #'                            dat_dim = c('dat', 'member'), memb_dim = 'member')
 #'var <- s2dv::MeanDims(ano$exp, "member")
@@ -190,6 +191,7 @@
 #'          filled.continents = FALSE, 
 #'          toptitle = 'Near-surface temperature Nov.',
 #'          titles = paste(2000:2005))
+#'  }
 #'@import utils
 #'@importFrom grDevices dev.cur dev.new dev.off
 #'@importFrom s2dv Reorder
diff --git a/R/VizMatrix.R b/R/VizMatrix.R
index 361a255376254e9a44738e5d25409f2a25942585..1fd9ac0201ea5a2ce53f80982851caf0c626d0e3 100644
--- a/R/VizMatrix.R
+++ b/R/VizMatrix.R
@@ -51,7 +51,8 @@
 #'
 #'@examples
 #'#Example with random data
-#' VizMatrix(var = matrix(rnorm(n = 120, mean = 0.3), 10, 12),
+#'  \dontrun{
+#' esviz:::VizMatrix(var = matrix(rnorm(n = 120, mean = 0.3), 10, 12),
 #'            cols = c('white','#fef0d9','#fdd49e','#fdbb84','#fc8d59',
 #'                      '#e34a33','#b30000', '#7f0000'),
 #'            brks = c(-1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1),
@@ -59,9 +60,9 @@
 #'            xtitle = "Forecast time (month)", ytitle = "Start date",
 #'            xlabels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", 
 #'                        "Aug", "Sep", "Oct", "Nov", "Dec"))
+#'  }
 #'
 #'@importFrom grDevices dev.new dev.off dev.cur 
-#'@export
 VizMatrix <- function(var, brks = NULL, cols = NULL, 
                        toptitle = NULL, title.color = "royalblue4", 
                        xtitle = NULL, ytitle = NULL, xlabels = NULL, xvert = FALSE, 
diff --git a/R/VizMostLikelyQuantileMap.R b/R/VizMostLikelyQuantileMap.R
index 2de037c050619d9c427dece34be4f023ccb97fd2..1796e392521451f51c9b438c6ba5ca8aee6f5297 100644
--- a/R/VizMostLikelyQuantileMap.R
+++ b/R/VizMostLikelyQuantileMap.R
@@ -57,47 +57,17 @@
 #'}
 #'
 #'# More complex example
-#'n_lons <- 40
-#'n_lats <- 20
-#'n_timesteps <- 100
-#'n_bins <- 4
 #'
 #'# 1. Generation of sample data
-#'lons <- seq(0, 359.5, length = n_lons)
-#'lats <- seq(-89.5, 89.5, length = n_lats)
-#'
-#'# This function builds a 3-D gaussian at a specified point in the map.
-#'make_gaussian <- function(lon, sd_lon, lat, sd_lat) {
-#'  w <- outer(lons, lats, function(x, y) dnorm(x, lon, sd_lon) * dnorm(y, lat, sd_lat))
-#'  min_w <- min(w)
-#'  w <- w - min_w
-#'  w <- w / max(w)
-#'  w <- t(w)
-#'  names(dim(w)) <- c('lat', 'lon')
-#'  w
-#'}
-#'
-#'# This function generates random time series (with values ranging 1 to 5) 
-#'# according to 2 input weights.
-#'gen_data <- function(w1, w2, n) {
-#'  r <- sample(1:5, n, 
-#'              prob = c(.05, .9 * w1, .05, .05, .9 * w2), 
-#'              replace = TRUE)
-#'  r <- r + runif(n, -0.5, 0.5)
-#'  dim(r) <- c(time = n)
-#'  r
-#'}
-#'
-#'# We build two 3-D gaussians.
-#'w1 <- make_gaussian(120, 80, 20, 30)
-#'w2 <- make_gaussian(260, 60, -10, 40)
+#'lons <- seq(0, 359.5, length = 40)
+#'lats <- seq(-89.5, 89.5, length = 20)
 #'
-#'# We generate sample data (with dimensions time, lat, lon) according 
-#'# to the generated gaussians
-#'sample_data <- multiApply::Apply(list(w1, w2), NULL, 
-#'                                 gen_data, n = n_timesteps)$output1
+#'# Generate sample data (with dimensions time, lat, lon)
+#'sample_data <- sample(1:5, 50 * 20 * 40, replace = TRUE)
+#'dim(sample_data) <- c(time = 50, lat = 20, lon = 40)
 #'
 #'# 2. Binning sample data
+#'n_bins <- 4
 #'prob_thresholds <- 1:n_bins / n_bins
 #'prob_thresholds <- prob_thresholds[1:(n_bins - 1)]
 #'thresholds <- quantile(sample_data, prob_thresholds)
@@ -105,7 +75,6 @@
 #'binning <- function(x, thresholds) {
 #'  n_samples <- length(x)
 #'  n_bins <- length(thresholds) + 1
-#'
 #'  thresholds <- c(thresholds, max(x))
 #'  result <- 1:n_bins
 #'  lower_threshold <- min(x) - 1
@@ -113,20 +82,17 @@
 #'    result[i] <- sum(x > lower_threshold & x <= thresholds[i]) / n_samples
 #'    lower_threshold <- thresholds[i]
 #'  }
-#'
 #'  dim(result) <- c(bin = n_bins)
 #'  result
 #'}
+#'bins <- apply(sample_data, 2:3, binning, thresholds)
+#'names(dim(bins))[1] <- "bin"
 #'
-#'bins <- multiApply::Apply(sample_data, 'time', binning, thresholds)$output1
-#'
-#'# 3. Plotting most likely quantile/bin
 #'\dontrun{
 #'VizMostLikelyQuantileMap(bins, lons, lats, 
 #'                         toptitle = 'Most likely quantile map',
 #'                         bar_titles = paste('% of belonging to', letters[1:n_bins]),
-#'                         mask = 1 - (w1 + w2 / max(c(w1, w2))), 
-#'                         brks = 20, width = 10, height = 8)
+#'                         brks = 20, width = 10, height = 10)
 #'}
 #'@importFrom maps map 
 #'@importFrom graphics box image layout mtext par plot.new
diff --git a/R/VizRobinson.R b/R/VizRobinson.R
index 7130f32fd44d1fd1e6768fd00569f5066d46b01d..c1a0cbb1941f8be6e8f4e2dd731d781eb46ff855 100644
--- a/R/VizRobinson.R
+++ b/R/VizRobinson.R
@@ -101,6 +101,7 @@
 #'dots <- data
 #'dots[which(dots < 4 & dots > -4)] <- 0
 #'dots[which(dots != 0)] <- 1
+#'  \dontrun{
 #'VizRobinson(data, lon = 0:359, lat = -90:90, dots = dots,
 #'            brks = seq(-10, 10, length.out = 11),
 #'            toptitle = 'synthetic example', vertical = FALSE,
@@ -111,13 +112,14 @@
 #'           color_fun = ClimPalette("purpleorange"), colNA = 'green',
 #'           toptitle = 'synthetic example', caption = 'Mollweide Projection',
 #'           width = 8, height = 6)
-#'
+#'  }
 #'@import sf ggplot2 rnaturalearth cowplot utils
 #'@importFrom dplyr mutate group_by summarise
 #'@importFrom ClimProjDiags Subset
+#' @importFrom rlang .data
 #'@export
 VizRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL,
-                        target_proj = 54030, legend = 'bar', style = 'point', 
+                        target_proj = "ESRI:54030", legend = 'bar', style = 'point', 
                         dots = NULL, mask = NULL, brks = NULL, cols = NULL, bar_limits = NULL,
                         triangle_ends = NULL, col_inf = NULL, col_sup = NULL, colNA = NULL, 
                         color_fun = ClimPalette(), bar_extra_margin = rep(0, 4), vertical = TRUE,
@@ -384,7 +386,7 @@ VizRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL,
     dots_df <- dots_df %>%
       dplyr::mutate(long = st_coordinates(dots_df)[, 1],
                     lat = st_coordinates(dots_df)[, 2])
-    dots_df <- subset(dots_df, dot == FALSE)
+    dots_df <- dplyr::filter(dots_df, .data$dot == FALSE)
   }
 
   # coastlines
@@ -437,9 +439,9 @@ VizRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL,
     datapoly <- merge(values, positions, by = "id")
     datapoly <- st_as_sf(datapoly, coords = c("x", "y"), crs = target_proj)
     datapoly <- datapoly %>%
-                dplyr::group_by(id) %>%
+                dplyr::group_by(.data$id) %>%
                 dplyr::summarise() %>%  #NOTE: VERY SLOW if plot global
-                dplyr::mutate(value = values[order(values$id), ]$value) %>%
+                dplyr::mutate(value = .data[[values[order(values$id), ]$value]]) %>%
                 st_cast("POLYGON") %>%
                 st_convex_hull()  # maintain outer polygen (no bowtie shape)
   }
@@ -448,12 +450,12 @@ VizRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL,
   if (style == 'polygon') {
     res_p <- ggplot(data = data_df) +  #NOTE: must be data_df?
              geom_sf(data = datapoly,
-                     aes(col = cut(value, breaks = brks_ggplot, include.lowest = T),
-                         fill = cut(value, breaks = brks_ggplot, include.lowest = T)))
+                     aes(col = cut(.data$value, breaks = brks_ggplot, include.lowest = T),
+                         fill = cut(.data$value, breaks = brks_ggplot, include.lowest = T)))
   } else if (style == 'point') {
     res_p <- ggplot(data = data_df) +
-             geom_point(aes(x = long, y = lat,
-                            col = cut(dat, breaks = brks_ggplot, include.lowest = T)),
+             geom_point(aes(x = .data$long, y = .data$lat,
+                            col = cut(.data$dat, breaks = brks_ggplot, include.lowest = T)),
                         #NOTE: These two lines make point size vary with lat
                         #size = point_size /  (data_df$lat / min(data_df$lat))) +
                         #size = (sort(rep(as.vector(lat), length(lon))) / max(lat)) * point_size) +
@@ -474,7 +476,7 @@ VizRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL,
                     expand = TRUE, datum = target_proj)
 
   if (!is.null(dots)) {
-    res_p <- res_p + geom_point(data = dots_df, aes(x = long, y = lat),
+    res_p <- res_p + geom_point(data = dots_df, aes(x = .data$long, y = .data$lat),
                                 shape = dots_shape, size = dots_size)
                                 #NOTE: This line makes point size vary with lat
                                 #size = dots_size / (dots_df$lat / min(dots_df$lat)))
diff --git a/R/VizSection.R b/R/VizSection.R
index b375ea78bbd1d702e8b8f240c54c985a8d31c7f9..f77233365cbb1260d151078f282893a78a7a2e58 100644
--- a/R/VizSection.R
+++ b/R/VizSection.R
@@ -40,10 +40,9 @@
 #'# Synthetic data
 #'data <- array(rep(seq(25, 10, length.out = 7), each = 21) - rnorm(147),
 #'              dim = c(lat = 21, depth = 7))
-#'VizSection(data, horiz = 0:20, depth = seq(0, 300, length.out = 7),
+#'esviz:::VizSection(data, horiz = 0:20, depth = seq(0, 300, length.out = 7),
 #'           toptitle = 'Temperature cross-section', units = "degC")
 #'@importFrom grDevices dev.cur dev.new dev.off rainbow
-#'@export
 VizSection <- function(var, horiz, depth, toptitle = '', sizetit = 1,
                         units = '', brks = NULL, cols = NULL, axelab = TRUE,
                         intydep = 200, intxhoriz = 20, drawleg = TRUE,
diff --git a/R/VizVsLTime.R b/R/VizVsLTime.R
index cff8550fbe7e268c39687879318ef01155da96f8..7fde22699f45bebcd3b501dba3fc4ef70bea98cd 100644
--- a/R/VizVsLTime.R
+++ b/R/VizVsLTime.R
@@ -77,14 +77,13 @@
 #'input_cor[, 2, ] <- corr_ano$corr[, 1, 1, ]
 #'input_cor[, 3, ] <- corr_ano$conf.upper[, 1, 1, ]
 #'input_cor[, 4, ] <- corr_ano$p.val[, 1, 1, ]
-#'VizVsLTime(input_cor, toptitle = "Correlation",
+#'esviz:::VizVsLTime(input_cor, toptitle = "Correlation",
 #'           monini = 11, limits = c(-1, 2), listexp = 'SEAS5', 
 #'           listobs = 'ERA5', biglab = FALSE, hlines = c(-1, 0, 1))
 #'
 #'@importFrom grDevices dev.cur dev.new dev.off
 #'@importFrom stats ts
 #'@importFrom s2dv InsertDim
-#'@export
 VizVsLTime <- function(var, toptitle = '', ytitle = '', monini = 1, freq = 12,
                nticks = NULL, limits = NULL,
                listexp = c('exp1', 'exp2', 'exp3'),
diff --git a/R/VizWeeklyClim.R b/R/VizWeeklyClim.R
index 01f71d84dbfdc1592d22558aa9d49bce6eb30125..663e4c5da2185e694c0ead1e778ce97ac940ac67 100644
--- a/R/VizWeeklyClim.R
+++ b/R/VizWeeklyClim.R
@@ -75,6 +75,7 @@
 #'@importFrom ClimProjDiags Subset
 #'@importFrom s2dv MeanDims
 #'@importFrom CSTools SplitDim
+#' @importFrom rlang .data
 #'@export
 VizWeeklyClim <- function(data, first_date, ref_period, last_date = NULL, 
                           data_years = NULL, time_dim = 'time', 
@@ -267,18 +268,18 @@ VizWeeklyClim <- function(data, first_date, ref_period, last_date = NULL,
   ## Create a ggplot object
   cols <- colorRampPalette(brewer.pal(9, palette))(6)
 
-  p <- ggplot(all, aes(x = day)) +
-       geom_ribbon(aes(ymin = p10, ymax = p90, group = week, fill = "p10-p90"), 
+  p <- ggplot(all, aes(x = .data$day)) +
+       geom_ribbon(aes(ymin = p10, ymax = p90, group = .data$week, fill = "p10-p90"), 
                    alpha = 0.7, show.legend = legend) + # extremes clim
-       geom_ribbon(aes(ymin = p33, ymax = p66, group = week, fill = "p33-p66"),
+       geom_ribbon(aes(ymin = p33, ymax = p66, group = .data$week, fill = "p33-p66"),
                    alpha = 0.7, show.legend = legend) +  # terciles clim
-       geom_line(aes(y = clim, group = week, color = "climatological mean",
+       geom_line(aes(y = clim, group = .data$week, color = "climatological mean",
                  linetype = "climatological mean"),
                  alpha = 1.0, linewidth = 0.7, show.legend = legend) + # mean clim
        geom_line(aes(y = data, color = "observed daily mean",
                  linetype = "observed daily mean"),
                  alpha = 1, linewidth = 0.2, show.legend = legend) + # daily evolution
-       geom_line(aes(y = week_mean, group = week, color = "observed weekly mean",
+       geom_line(aes(y = week_mean, group = .data$week, color = "observed weekly mean",
                  linetype = "observed weekly mean"),
                  alpha = 1, linewidth = 0.7, show.legend = legend) + 		 # weekly evolution
        theme_bw() + ylab(ytitle) + xlab(NULL) + 
diff --git a/man/Viz2VarsVsLTime.Rd b/man/Viz2VarsVsLTime.Rd
index c74534b87628b2aa842a27f8a4abc219a5099e85..8b476183b23da0b56d04ce9213da7c2d6ce4b5f2 100644
--- a/man/Viz2VarsVsLTime.Rd
+++ b/man/Viz2VarsVsLTime.Rd
@@ -115,7 +115,7 @@ input_rms <- array(dim = c(dat = 1, 3, time = 5))
 input_rms[, 1, ] <- rms_ano$conf.lower[, 1, 1, ]
 input_rms[, 2, ] <- rms_ano$rms[, 1, 1, ]
 input_rms[, 3, ] <- rms_ano$conf.upper[, 1, 1, ]
-Viz2VarsVsLTime(input_cor, input_rms,
+esviz:::Viz2VarsVsLTime(input_cor, input_rms,
                toptitle = "Time correlation and RMSE with ERA5",
                ytitle = "K", sizetit = 0.7,
                monini = 11, freq = 1, limits = c(-1, 5),
diff --git a/man/VizACC.Rd b/man/VizACC.Rd
index 81a4b3e5ffa5715c71e7ddaf7479e380fc3dd2f4..51ff6e65a38ff16085cb824af4c978de6085b137 100644
--- a/man/VizACC.Rd
+++ b/man/VizACC.Rd
@@ -95,6 +95,7 @@ confidence interval and the 95\% significance level given by a one-sided
 T-test.
 }
 \examples{
+ \dontrun{
 ano <- s2dv::Ano_CrossValid(map_temp$exp, map_temp$obs, memb = FALSE, 
                            dat_dim = c('dat', 'member'), memb_dim = 'member')
 lats <- attr(map_temp$exp, "Variables")$common$lat
@@ -107,9 +108,9 @@ input_acc[, , , , 2] <- acc$acc
 input_acc[, , , , 3] <- acc$conf.upper
 input_acc[, , , , 4] <- acc$p.val
 sdates <- paste0(2000:2005, '1101')
-VizACC(input_acc, sdates, 
+esviz:::VizACC(input_acc, sdates, 
       toptitle = "Spatial anomaly corr. coeff. with ERA5",
       ytitle = "K", sizetit = 0.7, freq = 12,
       legends = 'SEAS5', fileout = NULL)
-
+ }
 }
diff --git a/man/VizAnimateMap.Rd b/man/VizAnimateMap.Rd
index 23c0d4fce9aed6f9663bbf17d817e206b3b35e67..29b8effe6d93e0cc2dda169263e82659febdcfdd 100644
--- a/man/VizAnimateMap.Rd
+++ b/man/VizAnimateMap.Rd
@@ -151,10 +151,10 @@ clim <- s2dv::Clim(map_temp$exp, map_temp$obs, memb = FALSE,
                   dat_dim = c('dat', 'member'), memb_dim = 'member')
 lats <- attr(map_temp$exp, "Variables")$common$lat
 lons <- attr(map_temp$exp, "Variables")$common$lon
-VizAnimateMap(clim$clim_exp[1, 1, , , ], lon = lons, lat = lats,
+ \dontrun{
+esviz:::VizAnimateMap(clim$clim_exp[1, 1, , , ], lon = lons, lat = lats,
  toptitle = "climatology of decadal prediction", sizetit = 1, 
  units = "K", brks = seq(270, 300, 3), monini = 11, freq = 12, 
- msk95lev = FALSE, filled.continents = FALSE, intlon = 10, intlat = 10,
- fileout = NULL)
-
+ msk95lev = FALSE, filled.continents = FALSE, intlon = 10, intlat = 10)
+ }
 }
diff --git a/man/VizAno.Rd b/man/VizAno.Rd
index 3db84ec96d4e6e9097d159330be8e63324a2b6ae..36f27737a8764688b9057a0374cfd350ca877bed 100644
--- a/man/VizAno.Rd
+++ b/man/VizAno.Rd
@@ -107,7 +107,7 @@ Plots time series of raw or smoothed anomalies of any variable output from
 dim(ts_temp$exp) <- c(dat = 1, member = 3, sdate = 6, time = 5)
 dim(ts_temp$obs) <- c(dat = 1, member = 1, sdate = 6, time = 5)
 sdates <- paste0(2000:2005, '1101')
-VizAno(ts_temp$exp, ts_temp$obs, sdates,
+esviz:::VizAno(ts_temp$exp, ts_temp$obs, sdates,
       toptitle = "Raw 'tas'", ytitle = 'K', legends = 'ERA5', biglab = FALSE)
 
 }
diff --git a/man/VizBoxWhisker.Rd b/man/VizBoxWhisker.Rd
index 21f1a87e6222d754baa3f84b8a5e3941a8523825..3765c118146d4848edec9a02e809b04a893bdd86 100644
--- a/man/VizBoxWhisker.Rd
+++ b/man/VizBoxWhisker.Rd
@@ -89,6 +89,7 @@ n-monthly to n-yearly time series.
 \examples{
 # No example data is available over NAO region, so in this example we will
 # tweak the longitude and latitude.
+ \dontrun{
 ano <- s2dv::Ano_CrossValid(map_temp$exp, map_temp$obs, memb = FALSE,
                            dat_dim = c('dat', 'member'), memb_dim = 'member')
 nao <- s2dv::NAO(ano$exp, ano$obs, lat = seq(20, 80, length.out = 11),
@@ -96,10 +97,10 @@ nao <- s2dv::NAO(ano$exp, ano$obs, lat = seq(20, 80, length.out = 11),
                 ftime_dim = "time")
 nao$exp <- drop(aperm(nao$exp, c(2, 1, 3, 4)))
 nao$obs <- drop(nao$obs)
-VizBoxWhisker(nao$exp, nao$obs, toptitle = "NAO index", 
+esviz:::VizBoxWhisker(nao$exp, nao$obs, toptitle = "NAO index", 
              ytitle = "NAO index (PC1) TOS", monini = 11, freq = 1,
               yearini = 2000, expname = "SEAS5", obsname = "ERA5")
-
+ }
 }
 \seealso{
 EOF, ProjectField, NAO
diff --git a/man/VizClim.Rd b/man/VizClim.Rd
index e2d0c5274d559af6a9858429ca5f748f7450afe0..b207aaab2b3775bbaf3d348ccc297ca7b6f58d94 100644
--- a/man/VizClim.Rd
+++ b/man/VizClim.Rd
@@ -90,7 +90,7 @@ clim <- s2dv::Clim(ts_temp$exp, ts_temp$obs, time_dim = "sdate",
                   dat_dim = c("dat", "member"))
 dim(clim$clim_exp) <- dim(clim$clim_exp)[-3]
 dim(clim$clim_obs) <- dim(clim$clim_obs)[-2]
-VizClim(clim$clim_exp, clim$clim_obs, toptitle = 'Climatologies',
+esviz:::VizClim(clim$clim_exp, clim$clim_obs, toptitle = 'Climatologies',
        ytitle = 'K', monini = 11, listexp = 'SEAS5',
        listobs = 'ERA5', biglab = FALSE, fileout = NULL)
 
diff --git a/man/VizEquiMap.Rd b/man/VizEquiMap.Rd
index 57ba1c60f66afc8998fd6fe4ee05bfb66ba295e2..b342fc397f58f6230fd55f82837bfc30e54a45cd 100644
--- a/man/VizEquiMap.Rd
+++ b/man/VizEquiMap.Rd
@@ -376,6 +376,7 @@ include continents, oceans, and lakes. This plot function is compatible with
 figure layouts if colour bar is disabled.
 }
 \examples{
+ \dontrun{
 ano <- s2dv::Ano_CrossValid(map_temp$exp, map_temp$obs, memb = FALSE,
                            dat_dim = c('dat', 'member'), memb_dim = 'member')
 var <- s2dv::MeanDims(ano$exp, "member")
@@ -385,5 +386,5 @@ lons <- attr(map_temp$exp, "Variables")$common$lon
 VizEquiMap(var[1, 1, 1, 1, , ], lon = lons, lat = lats,
           toptitle = 'Near-surface temperature anomaly, Nov. 2000',
           filled.continents = FALSE, title_scale = 0.7)
-
+ }
 }
diff --git a/man/VizLayout.Rd b/man/VizLayout.Rd
index 525b11e20611eb86301dbc0417cbf06e0e80eaf2..0fd373c4165ef4128f4425a33d04894408a48db0 100644
--- a/man/VizLayout.Rd
+++ b/man/VizLayout.Rd
@@ -265,6 +265,7 @@ This function pops up a blank new device and fills it in, so it cannot be
 nested in complex layouts.
 }
 \examples{
+ \dontrun{
 ano <- s2dv::Ano_CrossValid(map_temp$exp, map_temp$obs, memb = FALSE,
                            dat_dim = c('dat', 'member'), memb_dim = 'member')
 var <- s2dv::MeanDims(ano$exp, "member")
@@ -276,4 +277,5 @@ VizLayout(fun = VizEquiMap, plot_dims = c('lat', 'lon'),
          filled.continents = FALSE, 
          toptitle = 'Near-surface temperature Nov.',
          titles = paste(2000:2005))
+ }
 }
diff --git a/man/VizMatrix.Rd b/man/VizMatrix.Rd
index 0705246f3531e0671a9a7933219abf5831d431de..016603c5f6b90041ba3d91a572143b0524bc5790 100644
--- a/man/VizMatrix.Rd
+++ b/man/VizMatrix.Rd
@@ -101,7 +101,8 @@ colors instead of numbers.
 }
 \examples{
 #Example with random data
-VizMatrix(var = matrix(rnorm(n = 120, mean = 0.3), 10, 12),
+ \dontrun{
+esviz:::VizMatrix(var = matrix(rnorm(n = 120, mean = 0.3), 10, 12),
            cols = c('white','#fef0d9','#fdd49e','#fdbb84','#fc8d59',
                      '#e34a33','#b30000', '#7f0000'),
            brks = c(-1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1),
@@ -109,5 +110,6 @@ VizMatrix(var = matrix(rnorm(n = 120, mean = 0.3), 10, 12),
            xtitle = "Forecast time (month)", ytitle = "Start date",
            xlabels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", 
                        "Aug", "Sep", "Oct", "Nov", "Dec"))
+ }
 
 }
diff --git a/man/VizMostLikelyQuantileMap.Rd b/man/VizMostLikelyQuantileMap.Rd
index 7b8ee9a253a6f9d84c1e943b3bd8a380e7d3a4fe..c15ab05578e73672c930c75f7a2f84a663ef44e5 100644
--- a/man/VizMostLikelyQuantileMap.Rd
+++ b/man/VizMostLikelyQuantileMap.Rd
@@ -79,47 +79,17 @@ VizMostLikelyQuantileMap(list(a, b, c), lons, lats,
 }
 
 # More complex example
-n_lons <- 40
-n_lats <- 20
-n_timesteps <- 100
-n_bins <- 4
 
 # 1. Generation of sample data
-lons <- seq(0, 359.5, length = n_lons)
-lats <- seq(-89.5, 89.5, length = n_lats)
-
-# This function builds a 3-D gaussian at a specified point in the map.
-make_gaussian <- function(lon, sd_lon, lat, sd_lat) {
- w <- outer(lons, lats, function(x, y) dnorm(x, lon, sd_lon) * dnorm(y, lat, sd_lat))
- min_w <- min(w)
- w <- w - min_w
- w <- w / max(w)
- w <- t(w)
- names(dim(w)) <- c('lat', 'lon')
- w
-}
-
-# This function generates random time series (with values ranging 1 to 5) 
-# according to 2 input weights.
-gen_data <- function(w1, w2, n) {
- r <- sample(1:5, n, 
-             prob = c(.05, .9 * w1, .05, .05, .9 * w2), 
-             replace = TRUE)
- r <- r + runif(n, -0.5, 0.5)
- dim(r) <- c(time = n)
- r
-}
-
-# We build two 3-D gaussians.
-w1 <- make_gaussian(120, 80, 20, 30)
-w2 <- make_gaussian(260, 60, -10, 40)
+lons <- seq(0, 359.5, length = 40)
+lats <- seq(-89.5, 89.5, length = 20)
 
-# We generate sample data (with dimensions time, lat, lon) according 
-# to the generated gaussians
-sample_data <- multiApply::Apply(list(w1, w2), NULL, 
-                                gen_data, n = n_timesteps)$output1
+# Generate sample data (with dimensions time, lat, lon)
+sample_data <- sample(1:5, 50 * 20 * 40, replace = TRUE)
+dim(sample_data) <- c(time = 50, lat = 20, lon = 40)
 
 # 2. Binning sample data
+n_bins <- 4
 prob_thresholds <- 1:n_bins / n_bins
 prob_thresholds <- prob_thresholds[1:(n_bins - 1)]
 thresholds <- quantile(sample_data, prob_thresholds)
@@ -127,7 +97,6 @@ thresholds <- quantile(sample_data, prob_thresholds)
 binning <- function(x, thresholds) {
  n_samples <- length(x)
  n_bins <- length(thresholds) + 1
-
  thresholds <- c(thresholds, max(x))
  result <- 1:n_bins
  lower_threshold <- min(x) - 1
@@ -135,20 +104,17 @@ binning <- function(x, thresholds) {
    result[i] <- sum(x > lower_threshold & x <= thresholds[i]) / n_samples
    lower_threshold <- thresholds[i]
  }
-
  dim(result) <- c(bin = n_bins)
  result
 }
+bins <- apply(sample_data, 2:3, binning, thresholds)
+names(dim(bins))[1] <- "bin"
 
-bins <- multiApply::Apply(sample_data, 'time', binning, thresholds)$output1
-
-# 3. Plotting most likely quantile/bin
 \dontrun{
 VizMostLikelyQuantileMap(bins, lons, lats, 
                         toptitle = 'Most likely quantile map',
                         bar_titles = paste('\% of belonging to', letters[1:n_bins]),
-                        mask = 1 - (w1 + w2 / max(c(w1, w2))), 
-                        brks = 20, width = 10, height = 8)
+                        brks = 20, width = 10, height = 10)
 }
 }
 \seealso{
diff --git a/man/VizRobinson.Rd b/man/VizRobinson.Rd
index 2da516256ccf7512833253ac705af1c343834dd0..632074413a548c004d9c8e35a6c504bc5a0d18c4 100644
--- a/man/VizRobinson.Rd
+++ b/man/VizRobinson.Rd
@@ -10,7 +10,7 @@ VizRobinson(
   lat,
   lon_dim = NULL,
   lat_dim = NULL,
-  target_proj = 54030,
+  target_proj = "ESRI:54030",
   legend = "bar",
   style = "point",
   dots = NULL,
@@ -172,6 +172,7 @@ data <- array(rep(seq(-10, 10, length.out = 181), 360) + rnorm(360),
 dots <- data
 dots[which(dots < 4 & dots > -4)] <- 0
 dots[which(dots != 0)] <- 1
+ \dontrun{
 VizRobinson(data, lon = 0:359, lat = -90:90, dots = dots,
            brks = seq(-10, 10, length.out = 11),
            toptitle = 'synthetic example', vertical = FALSE,
@@ -182,5 +183,5 @@ VizRobinson(data, lon = 0:359, lat = -90:90, mask = dots, legend = 'ggplot2',
           color_fun = ClimPalette("purpleorange"), colNA = 'green',
           toptitle = 'synthetic example', caption = 'Mollweide Projection',
           width = 8, height = 6)
-
+ }
 }
diff --git a/man/VizSection.Rd b/man/VizSection.Rd
index c11f2c910976505c534c966cbce49ed7fdfb4bad..bc2947b13342b88ec916f870f86116d16109bf06 100644
--- a/man/VizSection.Rd
+++ b/man/VizSection.Rd
@@ -84,6 +84,6 @@ Plot a (longitude,depth) or (latitude,depth) section.
 # Synthetic data
 data <- array(rep(seq(25, 10, length.out = 7), each = 21) - rnorm(147),
              dim = c(lat = 21, depth = 7))
-VizSection(data, horiz = 0:20, depth = seq(0, 300, length.out = 7),
+esviz:::VizSection(data, horiz = 0:20, depth = seq(0, 300, length.out = 7),
           toptitle = 'Temperature cross-section', units = "degC")
 }
diff --git a/man/VizVsLTime.Rd b/man/VizVsLTime.Rd
index f8905175adeae244d4b37f5b0ac4023264cee25c..7ff5ad47239f4ce486a4f4196726fb95b2b28532 100644
--- a/man/VizVsLTime.Rd
+++ b/man/VizVsLTime.Rd
@@ -127,7 +127,7 @@ input_cor[, 1, ] <- corr_ano$conf.lower[, 1, 1, ]
 input_cor[, 2, ] <- corr_ano$corr[, 1, 1, ]
 input_cor[, 3, ] <- corr_ano$conf.upper[, 1, 1, ]
 input_cor[, 4, ] <- corr_ano$p.val[, 1, 1, ]
-VizVsLTime(input_cor, toptitle = "Correlation",
+esviz:::VizVsLTime(input_cor, toptitle = "Correlation",
           monini = 11, limits = c(-1, 2), listexp = 'SEAS5', 
           listobs = 'ERA5', biglab = FALSE, hlines = c(-1, 0, 1))
 
diff --git a/vignettes/map_proj_equi.md b/vignettes/map_proj_equi.md
index ed5519725a3ff785df77c09e7488db07d57b8705..9bef2227a4338dba758a9207d4f68a8556c20d12 100644
--- a/vignettes/map_proj_equi.md
+++ b/vignettes/map_proj_equi.md
@@ -145,7 +145,14 @@ VizEquiMap(exp1[1, 1, , ], lon = lons, lat = lats, filled.continents = FALSE,
 ## VizAnimateMap()
 
 `VizAnimateMap()` generates animation of maps from `VizEquiMap()` or 
-`VizStereoMap()` and saves them as .gif file . The input array should have dimensions as follow: 
+`VizStereoMap()` and saves them as .gif file . 
+Note that VizAnimateMap() is not exported, so we need to specify the fu
+nction name first.
+```r
+VizAnimateMap <- esviz:::VizAnimateMap
+```
+
+The input array should have dimensions as follow: 
 ```r
   c(n. of experiments/observations, n. of observations, 3, n. of forecast times, n. of latitudes, n. of longitudes)
 ```
diff --git a/vignettes/map_proj_stereo.md b/vignettes/map_proj_stereo.md
index b163b6007cc9bfc7341a2b5eab177cead9947f6b..66d5a3a7bfbfc42a336480d5762410a563bdd48e 100644
--- a/vignettes/map_proj_stereo.md
+++ b/vignettes/map_proj_stereo.md
@@ -90,7 +90,14 @@ VizLayout(VizStereoMap, c('lat', 'lon'), exp1_clim[, lats_np_ind, ],
 
 ## VizAnimateMap()
 
-`VizAnimateMap()` can generate animation with the plots generated by `VizStereoMap()`. The input array should have dimensions as follow: 
+`VizAnimateMap()` can generate animation with the plots generated by `VizStereoMap()`.
+Note that VizAnimateMap() is not exported, so we need to specify the fu
+nction name first.
+```r
+VizAnimateMap <- esviz:::VizAnimateMap
+```
+
+The input array should have dimensions as follow: 
 ```r
   c(n. of experiments/observations, n. of observations, 3, n. of forecast times, n. of latitudes, n. of longitudes)
 ```
diff --git a/vignettes/time_series.md b/vignettes/time_series.md
index 1e62b8e4c2aeb82766082ea9388a95f8744c7120..e6d6a0c9b6a00bf1c418e4058c61835d4ff49868 100644
--- a/vignettes/time_series.md
+++ b/vignettes/time_series.md
@@ -13,6 +13,16 @@ vignette: >
 
 In this vignette, we will go through the functions that plot the time series type of figures: 
 `VizAno()`, `VizClim()`, `VizVsLtime()`, `Viz2VarsVsLTime()`, `VizBoxWhisker()`, `VizACC()`. 
+Note that all the functions here are not exported, so we need to specify the function names first.
+
+```r
+VizAno <- esviz:::VizAno
+VizClim <- esviz:::VizClim
+VizVsLtime <- esviz:::VizVsLtime
+Viz2VarsVsLTime <- esviz:::Viz2VarsVsLTime
+VizBoxWhisker <- esviz:::VizBoxWhisker
+VizACC <- esviz:::VizACC
+```
 
 We will use the data loaded in [Data Loading](data_loading.md). 
 Please follow that vignette first to get the data if you want to replicate the results in this vignette. Load European region data with `region = "euro"` in the loading script.