diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index ccaf0809ed7097283f5eaa8c02955282f902a5e7..4dd045cee587be221dc7a3641bbf326bb377c167 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -40,3 +40,16 @@ unit-test-subseasonal:  # This job runs in the test stage.
     - echo "Running subseasonal unit tests..."
     - Rscript ./tests/test_subseasonal.R
 
+unit-test-utilities:  # This job runs in the test stage.
+  stage: test
+  script:
+    - echo "Loading modules..."
+    - module load R/4.1.2-foss-2015a-bare
+    - module load CDO/1.9.8-foss-2015a
+    - module load GEOS/3.7.2-foss-2015a-Python-3.7.3
+    - module load GDAL/2.2.1-foss-2015a
+    - module load PROJ/4.8.0-foss-2015a
+    - module list
+    - echo "Running subseasonal unit tests..."
+    - Rscript ./tests/test_utils.R
+
diff --git a/modules/Aggregation/Aggregation.R b/modules/Aggregation/Aggregation.R
index d5c09faca9f7a44d8a650fee2bc286f88fabb69f..e4cc6b978de82951127a63164b672156c27c0d92 100644
--- a/modules/Aggregation/Aggregation.R
+++ b/modules/Aggregation/Aggregation.R
@@ -36,7 +36,7 @@ Aggregation <- function(recipe, data) {
                       } else {
                         result <- data[x]
                       }}, simplify = TRUE)
-      }                   
+      }
     } else if (!is.null(custom)) {
       # calCUlate aggregations from fmonth1 to ftmonth2
       # for decadal, it makes sense to get all seasons and annual means too
diff --git a/modules/Aggregation/R/agg_ini_end.R b/modules/Aggregation/R/agg_ini_end.R
index 563745bd89687fb6d795116a0b1091e79921f69e..bb3a410ac2221b38ab5cb50dfd0557cf723da939 100644
--- a/modules/Aggregation/R/agg_ini_end.R
+++ b/modules/Aggregation/R/agg_ini_end.R
@@ -3,7 +3,7 @@ agg_ini_end <- function(x, ini, end, indices = NULL, method, na.rm ,ncores) {
   if (!is.null(ini) && is.null(indices)) {
     # list of vectors for the indices from ini to end pairs
     indices <- lapply(1:length(ini), function(x) { 
-                    ini[x]:end[x]})
+                      ini[x]:end[x]})
     plotting_attr <- list(ini_ftime = ini, end_ftime = end)
   } else {
     # take the first and last element of each indices list for time_bounds saving
diff --git a/modules/Visualization/R/get_plot_time_labels.R b/modules/Visualization/R/get_plot_time_labels.R
new file mode 100644
index 0000000000000000000000000000000000000000..2dfc420244cd58e4d6493698be4d1fa2d757120a
--- /dev/null
+++ b/modules/Visualization/R/get_plot_time_labels.R
@@ -0,0 +1,76 @@
+# This function returns the day, month and year labels for forecast plots.
+# If the data is time-aggregated (time_bounds is not NULL); then the time
+# labels reflect the aggregated period, and will include the years to avoid
+# displaying incorrect or repetitive information.
+# dates: the array of POSIXt Dates in the s2dv_cube
+# horizon: The forecast horizon: 'subseasonal', 'seasonal' or 'decadal'.
+# init_date: Initialization date as POSIXt object
+# start_date: start date 
+# i_syear: The forecast year
+# time_bounds: the array of POSIXt time_bounds from the s2dv_cube
+
+.get_plot_time_labels <- function(dates, horizon, init_date, start_date,
+                                  i_syear, time_bounds = NULL) {
+  years <- lubridate::year(dates[1, 1, which(start_date == i_syear), ])
+  if (is.null(attributes(time_bounds))) {
+    if (horizon == "subseasonal") {
+      time_labels <- dates[1, 1, which(start_date == i_syear), ]
+      monday <- ymd_hms(time_labels) - days(wday(ymd_hms(time_labels),
+                                                 week_start = 1) - 1)
+      sunday <- monday + days(6)
+      time_labels <- paste0("Valid from ", format(monday,"%d-%m"), " to ",
+                            format(sunday, "%d-%m"), " of ")
+    } else {
+      time_labels <- lubridate::month(dates[1, 1, which(start_date == i_syear), ],
+                                      label = T, abb = F)
+    }
+  } else {
+    if (length(attributes(time_bounds)$plotting_attr) > 1) {
+      years <- NULL
+      start <- time_bounds$start[1, 1, which(start_date == i_syear), ]
+      end <- time_bounds$end[1, 1, which(start_date == i_syear), ]
+      time_labels <-
+        unlist(lapply(1:length(start), function(i) {
+          ftime_ini <- attributes(time_bounds)$plotting_attr$ini[i]
+          ftime_end <- attributes(time_bounds)$plotting_attr$end[i]
+          # labels for file name:
+          ftime <- paste0(ftime_ini, "-", ftime_end)
+          # title names:
+          if (horizon %in% c("seasonal", "decadal")) {
+            # start <- Subset(time_bounds$start,
+            #                 along = c("syear", "time"),
+            #                 indices = list(syear = i_syear, time = i),
+            #                 drop = TRUE)
+            # end <- Subset(time_bounds$end,
+            #               along = c("syear", "time"),
+            #               indices = list(syear = i_syear, time = i),
+            #               drop = TRUE)
+            year_ini <- year(start[i])
+            month_ini <- month.name[month(start[i])]
+            year_end <- year(end[i])
+            month_end <- month.name[month(end[i])]
+            if (year_ini == year_end) {
+              toptitle <- paste(month_ini, "to", month_end, year_end)
+            } else {
+              toptitle <- paste(month_ini, year_ini, "to",
+                                month_end, year_end)
+            }
+          } else if (horizon == "subseasonal") {
+            ## set to make ftime_ini be a Monday
+            ## if init_week always Thursday, the 2nd addend can be simplified to +4
+            ftime_ini <- ymd(init_date) +
+                         ((8 - wday(ymd(init_date), week_start = 1)) %% 7) +
+                         weeks(ftime_ini - 1)
+            ## set to make ftime_end be a Sunday
+            ftime_end <- ymd(init_date) +
+                         ((8 - wday(ymd(init_date), week_start = 1)) %% 7) +
+                         weeks(ftime_end - 1) + 6
+            toptitle <- paste("Valid from", ftime_ini, "to", ftime_end)
+          }
+        }))
+    } else {
+      time_labels <- attributes(time_bounds)$plotting_attr[[1]]
+    }
+  }
+  return(list(time_labels = time_labels, years = years))
+}
diff --git a/modules/Visualization/R/plot_ensemble_mean.R b/modules/Visualization/R/plot_ensemble_mean.R
index 19521f7dc99ab4e568bb83ed36085ce331e83684..953aafd1615d37345962b679e1c8223e4b3cc523 100644
--- a/modules/Visualization/R/plot_ensemble_mean.R
+++ b/modules/Visualization/R/plot_ensemble_mean.R
@@ -114,59 +114,16 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
         var_dots <- as.numeric(var_dots <= 0)
         dim(var_dots) <- dim_dots
       }
-      if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
-        toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
-                           "\n", "Forecast Ensemble Mean / ", "Issued on ", 
-                           format(ymd(start_date), "%d-%m-%Y"))
-      } else {
-        toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
-                           "\n", "Forecast Ensemble Mean / ", "Init.: ", i_syear)
-      }
-      if (is.null(attributes(fcst$attrs$time_bounds))) {
-        if (tolower(recipe$Analysis$Horizon) %in% c("seasonal", "decadal")) {
-          time_labels <- lubridate::month(fcst$attrs$Dates[1, 1, which(start_date == i_syear), ],
-                                          label = T, abb = F)
-        } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
-          time_labels <- fcst$attrs$Dates[1, 1, which(start_date == i_syear), ] 
-          monday <- ymd_hms(time_labels) - days(wday(ymd_hms(time_labels), week_start = 1) - 1)
-          sunday <- monday + days(6)
-          time_labels <- paste0("Valid from ", format(monday,"%d-%m"), " to ", format(sunday, "%d-%m"))
-        }
-      } else {
-        if (length(attributes(fcst$attrs$time_bounds)$plotting_attr) > 1) {
-          time_labels <-
-            unlist(lapply(1:length(fcst$attrs$time_bounds$start), function(i) {
-              ftime_ini <- attributes(fcst$attrs$time_bounds)$plotting_attr$ini[i]
-              ftime_end <- attributes(fcst$attrs$time_bounds)$plotting_attr$end[i]
-              # labels for file name:
-              ftime <- paste0(ftime_ini, "-", ftime_end)
-              # title names:
-              if (tolower(recipe$Analysis$Horizon) == "seasonal") {
-                ftime_ini <- init_date + ftime_ini - 1
-                ftime_ini <- ifelse(ftime_ini > 12, ftime_ini - 12, ftime_ini)
-                ftime_ini <- month.name[ftime_ini]
-                ftime_end <- init_date + ftime_end - 1
-                ftime_end <- ifelse(ftime_end > 12, ftime_end - 12, ftime_end)
-                ftime_end <- month.name[ftime_end]
-                toptitle <- paste(ftime_ini, "to", ftime_end)
-              } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
-                ## set to make ftime_ini be a Monday
-                ## if init_week always Thursday, the 2nd addend can be simplified to +4
-                ftime_ini <- ymd(init_date) + 
-                             ((8 - wday(ymd(init_date), week_start = 1)) %% 7) + 
-                             weeks(ftime_ini)
-                ## set to make ftime_end be a Sunday
-                ftime_end <- ymd(init_week) + 
-                             ((8 - wday(ymd(init_date), week_start = 1)) %% 7) + 
-                             weeks(ftime_end) + 6
-                toptitle <- paste("Valid from", ftime_ini, "to", ftime_end)
-                }            
-              })) 
-        } else {
-          time_labels <- attributes(fcst$attrs$time_bounds)$plotting_attr[[1]]
-        }
-      }
-      years <- lubridate::year(fcst$attrs$Dates[1, 1, which(start_date == i_syear), ])
+
+      # Plot title labels look different depending on horizon and aggregation
+      labels <- .get_plot_time_labels(dates = fcst$attrs$Dates,
+                                      horizon = recipe$Analysis$Horizon,
+                                      init_date = init_date,
+                                      start_date = start_date,
+                                      i_syear = i_syear,
+                                      time_bounds = fcst$attrs$time_bounds)
+      time_labels <- labels$time_labels
+      years <- labels$years
       
       if (recipe$Analysis$Workflow$Visualization$multi_panel) { 
         # Define name of output file and titles
@@ -175,6 +132,14 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
         } else {
           titles <- NULL
         }
+        if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+          toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
+                             "\n", "Forecast Ensemble Mean / ", "Issued on ", 
+                             format(ymd(start_date), "%d-%m-%Y"))
+        } else {
+          toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
+                             "\n", "Forecast Ensemble Mean / ", "Init.: ", i_syear)
+        }
         # Plots
         output_configuration <- output_conf$Multipanel$forecast_ensemble_mean
         base_args <- list(fun = "PlotEquiMap",
@@ -246,7 +211,7 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
                                "\n", "Ensemble Mean / ",
                                "Issued on ", 
                                format(ymd(start_date), "%d-%m-%Y"),
-                               "\n", time_labels[i], " of ",  years[i])
+                               "\n", time_labels[i], years[i])
           } else {
             toptitle <- paste0(system_name, " / ",
                                str_to_title(var_long_name),
diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index 8904ed575fd2a0f4189614ac9aca181f05d29fd6..6ffddafc0d8a21ab144c86446f460cf62f8a192e 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -328,11 +328,13 @@ plot_metrics <- function(recipe, data_cube, metrics,
                   # labels for file name:
                   forecast_time <- paste0(forecast_time_ini, "-", forecast_time_end)
                   # title names:
+                  ## TODO: There's probably a better way to do this, like using
+                  ## time_bounds$start and time_bounds$end directly.
                   forecast_time_ini <- init_month + forecast_time_ini - 1
                   forecat_time_ini <- ifelse(forecast_time_ini > 12, forecast_time_ini - 12, forecast_time_ini) 
                   forecast_time_ini <- month.name[forecast_time_ini]
                   forecast_time_end <- init_month + forecast_time_end - 1
-                  forecat_time_end <- ifelse(forecast_time_end > 12, forecast_time_end - 12, forecast_time_end)
+                  forecast_time_end <- ifelse(forecast_time_end > 12, forecast_time_end - 12, forecast_time_end)
                   forecast_time_end <- month.name[forecast_time_end]
                   toptitle <- paste(system_name, "/",
                                     str_to_title(var_long_name),
diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index 3787a1a9c67e6a04c9d36a03e6ade36346ee5bb0..9d03347f23f7225c7a73f1f957a2341caed0c881 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -122,68 +122,30 @@ plot_most_likely_terciles <- function(recipe,
       } else {
         var_dots <- NULL
       }
-      # Define top title: system name, variable, forecast date
-      if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
-        toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
-                           "\n", "Most Likely Tercile / ", "Issued on ",
-                           format(ymd(start_date), "%d-%m-%Y"))
-      } else {
-        toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
-                           "\n", "Most Likely Tercile / ", "Initialization: ",
-                           i_syear)
-      }
-      # Define labels for forecast times
-      if (tolower(recipe$Analysis$Horizon == "subseasonal")) {
-        ## TODO: time_bounds case for subseasonal data
-        time_labels <- fcst$attrs$Dates[1, 1, which(start_date == i_syear), ]
-        monday <- ymd_hms(time_labels) - days(wday(ymd_hms(time_labels), week_start = 1) - 1)
-        sunday <- monday + days(6)
-        time_labels <- paste0("Valid from ", format(monday,"%d-%m"), " to ", format(sunday, "%d-%m"))
-      } else {
-        if (is.null(attributes(fcst$attrs$time_bounds))) {
-          ## TODO: Does this need to be redefined for subseasonal?
-          time_labels <- lubridate::month(fcst$attrs$Dates[1, 1, which(start_date == i_syear), ],
-                                          label = T, abb = F)
-        } else {
-          if (length(attributes(fcst$attrs$time_bounds)$plotting_attr) > 1) {
-            time_labels <- 
-              unlist(lapply(1:length(fcst$attrs$time_bounds$start), function(i) {
-                ftime_ini <- attributes(fcst$attrs$time_bounds)$plotting_attr$ini[i]
-                ftime_end <- attributes(fcst$attrs$time_bounds)$plotting_attr$end[i]
-                # labels for file name:
-                ftime <- paste0(ftime_ini, "-", ftime_end)
-                # title names:
-                if (tolower(recipe$Analysis$Horizon) == "seasonal") {
-                  ftime_ini <- init_date + ftime_ini - 1
-                  ftime_ini <- ifelse(ftime_ini > 12, ftime_ini - 12, ftime_ini)
-                  ftime_ini <- month.name[ftime_ini]
-                  ftime_end <- init_date + ftime_end - 1
-                  ftime_end <- ifelse(ftime_end > 12, ftime_end - 12, ftime_end)
-                  ftime_end <- month.name[ftime_end]
-                  toptitle <- paste(ftime_ini, "to", ftime_end)
-                } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
-                  ## set to make ftime_ini be a Monday
-                  ## if init_date always Thursday, the 2nd addend can be simplified to +4
-                  ftime_ini <- ymd(init_date) +
-                               ((8 - wday(ymd(init_date), week_start = 1)) %% 7) +
-                               weeks(ftime_ini)
-                  ## set to make ftime_end be a Sunday
-                  ftime_end <- ymd(init_date) +
-                               ((8 - wday(ymd(init_date), week_start = 1)) %% 7) +
-                               weeks(ftime_end) + 6
-                  toptitle <- paste("Valid from", ftime_ini, "to", ftime_end)
-                  }
-                }))
-          } else {
-            time_labels <- attributes(fcst$attrs$time_bounds)$plotting_attr[[1]]
-          }
-        }
-      }
-      years <- lubridate::year(fcst$attrs$Dates[1, 1, which(start_date == i_syear), ])
+
+      # Plot title labels look different depending on horizon and aggregation
+      labels <- .get_plot_time_labels(dates = fcst$attrs$Dates,
+                                      horizon = recipe$Analysis$Horizon,
+                                      init_date = init_date,
+                                      start_date = start_date,
+                                      i_syear = i_syear,
+                                      time_bounds = fcst$attrs$time_bounds)
+      time_labels <- labels$time_labels
+      years <- labels$years
       
       if (recipe$Analysis$Workflow$Visualization$multi_panel) {
         ## TODO: Ensure this works for daily and sub-daily cases
         titles <- as.vector(time_labels)
+        # Define top title: system name, variable, forecast date
+        if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+          toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
+                             "\n", "Most Likely Tercile / ", "Issued on ",
+                             format(ymd(start_date), "%d-%m-%Y"))
+        } else {
+          toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
+                             "\n", "Most Likely Tercile / ", "Initialization: ",
+                             i_syear)
+        }
         # Plots
         ## NOTE: PlotLayout() and PlotMostLikelyQuantileMap() are still being worked
         ## on. This option does not work with mask or dots for now. 
@@ -247,7 +209,8 @@ plot_most_likely_terciles <- function(recipe,
               forecast_time <- sprintf("%02d", forecast_time)
             } else {
               ## TODO: Check this
-              forecast_time <- sprintf("%02d", time_labels[i])
+              # forecast_time <- sprintf("%02d", time_labels[i])
+              forecast_time <- time_labels[i]
             }
           } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
             forecast_time <- sprintf("%02d", i)
@@ -282,7 +245,7 @@ plot_most_likely_terciles <- function(recipe,
                                "\n", "Most Likely Tercile / ",
                                "Issued on ",
                                format(ymd(start_date), "%d-%m-%Y"),
-                               "\n", time_labels[i], " of ",  years[i])
+                               "\n", time_labels[i], years[i])
           } else {
             toptitle <- paste0(system_name, " / ",
                                str_to_title(var_long_name),
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index 52ce506875362ab4e7f04bdd5e96aa3e493a8810..be46f68e53b670de5ff5a3094075be473d4707f9 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -6,6 +6,7 @@
 
 source("modules/Visualization/R/plot_metrics.R")
 source("modules/Visualization/R/get_proj_code.R")
+source("modules/Visualization/R/get_plot_time_labels.R")
 ## TODO: Remove after the next s2dv release
 source("modules/Visualization/R/tmp/PlotRobinson.R")
 source("modules/Visualization/R/plot_most_likely_terciles_map.R")
diff --git a/tests/test_utils.R b/tests/test_utils.R
new file mode 100644
index 0000000000000000000000000000000000000000..87769761e559475068c3c2070f92d0e7f16e8cce
--- /dev/null
+++ b/tests/test_utils.R
@@ -0,0 +1,9 @@
+library(testthat)
+
+path_testthat <- file.path('./tests/testthat/')
+files_testthat <- list.files('./tests/testthat/', pattern = 'test-utils')
+
+for (i_file in 1:length(files_testthat)) {
+  source(paste0('./tests/testthat/', files_testthat[i_file]))
+}
+
diff --git a/tests/testthat/test-utils-get_plot_time_labels.R b/tests/testthat/test-utils-get_plot_time_labels.R
new file mode 100644
index 0000000000000000000000000000000000000000..686bdfb1e30f486f3ab57c8071e60381e6fe6f3b
--- /dev/null
+++ b/tests/testthat/test-utils-get_plot_time_labels.R
@@ -0,0 +1,269 @@
+context("get_plot_time_labels.R")
+
+source("tools/libs.R")
+source("modules/Visualization/R/get_plot_time_labels.R")
+###########################################
+
+test_that("1. Subseasonal, time-aggregated data", {
+  horizon <- "subseasonal"
+  init_date <- "20231221"
+  start_date <- init_date
+  i_syear <- init_date
+
+  time_bounds <- list(start = as.POSIXct("2023-12-21 09:00:00 UTC"),
+                      end = as.POSIXct("2024-01-11 09:00:00 UTC"))
+  attr(time_bounds, "plotting_attr") <- list(ini_ftime = 1,
+                                             end_ftime = 4)
+
+  dates <- as.POSIXct("2023-12-21 09:00:00 UTC")
+  dim(dates) <- c(sday = 1, sweek = 1, syear = 1, time = 1)
+  time_bounds[] <- lapply(time_bounds, 
+                          function(x) {
+                            dim(x) <- dim(dates)
+                            return(x)
+                          })
+
+  labels <- .get_plot_time_labels(
+    dates = dates,
+    horizon = horizon,
+    init_date = init_date,
+    start_date = start_date,
+    i_syear = i_syear,
+    time_bounds = time_bounds
+  )
+
+  # Results
+  expect_equal(
+    labels$years,
+    NULL
+  )
+  expect_equal(
+    labels$time_labels[1],
+    "Valid from 2023-12-25 to 2024-01-21"
+  )
+})
+
+test_that("2. Subseasonal, week-by-week data", {
+  horizon <- "subseasonal"
+  init_date <- "20231214"
+  start_date <- init_date
+  i_syear <- init_date
+
+  time_bounds <- NULL
+
+  dates <- as.POSIXct(c("2023-12-21 09:00:00 UTC", "2023-12-28 09:00:00 UTC",
+                        "2024-01-04 09:00:00 UTC", "2024-01-11 09:00:00 UTC"))
+  dim(dates) <- c(sday = 1, sweek = 1, syear = 1, time = 4)
+
+  labels <- .get_plot_time_labels(
+    dates = dates,
+    horizon = horizon,
+    init_date = init_date,
+    start_date = start_date,
+    i_syear = i_syear,
+    time_bounds = time_bounds
+  )
+
+  # Results
+  expect_equal(
+    labels$years,
+    c(2023, 2023, 2024, 2024)
+  )
+  expect_equal(
+    labels$time_labels[1],
+    "Valid from 18-12 to 24-12 of "
+  )
+  expect_equal(
+    labels$time_labels[2],
+    "Valid from 25-12 to 31-12 of "
+  )
+  expect_equal(
+    labels$time_labels[3],
+    "Valid from 01-01 to 07-01 of "
+  )
+  expect_equal(
+    labels$time_labels[4],
+    "Valid from 08-01 to 14-01 of "
+  )
+})
+
+test_that("3. Seasonal, time-aggregated data", {
+  horizon <- "seasonal"
+  init_date <- 10
+  start_date <- "20231001"
+  i_syear <- "20231001"
+
+  time_bounds <- list(start = c(as.POSIXct("2023-10-16 12:00:00 UTC"),
+                                as.POSIXct("2023-12-16 09:00:00 UTC")),
+                      end = c(as.POSIXct("2023-12-16 09:00:00 UTC"),
+                              as.POSIXct("2024-03-16 09:00:00 UTC")))
+  attr(time_bounds, "plotting_attr") <- list(ini_ftime = c(1, 3),
+                                             end_ftime = c(3, 6)
+  )
+
+  dates <- as.POSIXct(c("2023-10-16 12:00:00 UTC", "2023-12-16 09:00:00 UTC"))
+  dim(dates) <- c(sday = 1, sweek = 1, syear = 1, time = 2)
+  time_bounds[] <- lapply(time_bounds, 
+                          function(x) {
+                            dim(x) <- dim(dates)
+                            return(x)
+                          })
+
+  labels <- .get_plot_time_labels(
+    dates = dates,
+    horizon = horizon,
+    init_date = init_date,
+    start_date = start_date,
+    i_syear = i_syear,
+    time_bounds = time_bounds
+  )
+
+  # Results
+  expect_equal(
+    labels$years,
+    NULL
+  )
+  expect_equal(
+    labels$time_labels[1],
+    "October to December 2023"
+  )
+  expect_equal(
+    labels$time_labels[2],
+    "December 2023 to March 2024"
+  )
+})
+
+test_that("4. Seasonal, month-by-month data", {
+  horizon <- "seasonal"
+  init_date <- 10
+  start_date <- "20231001"
+  i_syear <- "20231001"
+
+  time_bounds <- NULL
+
+  dates <- as.POSIXct(c("2023-10-16 12:00:00 UTC", "2023-11-15 12:00:00 UTC",
+                        "2023-12-16 09:00:00 UTC", "2024-01-15 12:00:00 UTC"))
+  dim(dates) <- c(sday = 1, sweek = 1, syear = 1, time = 4)
+
+  labels <- .get_plot_time_labels(
+    dates = dates,
+    horizon = horizon,
+    init_date = init_date,
+    start_date = start_date,
+    i_syear = i_syear,
+    time_bounds = time_bounds
+  )
+
+  # Results
+  expect_equal(
+    labels$years,
+    c(2023, 2023, 2023, 2024)
+  )
+  expect_equal(
+    as.vector(labels$time_labels),
+    c("October", "November", "December", "January")
+  )
+})
+
+test_that("5. Decadal, time-aggregated data", {
+
+  horizon <- "decadal"
+  init_date <- 1
+  start_date <- c("2020", "2021")
+
+  time_bounds <- list(start = as.POSIXct(c("2020-11-16 UTC", "2021-11-16 UTC",
+                                           "2020-11-16 UTC", "2021-11-16 UTC")),
+                      end = as.POSIXct(c("2021-02-15 00:00:00 UTC",
+                                         "2022-02-15 00:00:00 UTC",
+                                         "2020-12-16 12:00:00 UTC",
+                                         "2021-12-16 12:00:00 UTC")))
+  attr(time_bounds, "plotting_attr") <- list(ini_ftime = c(1, 1),
+                                             end_ftime = c(4, 2)
+  )
+
+  dates <- as.POSIXct(c("2020-11-16 UTC", "2021-11-16 UTC", "2020-11-16 UTC", "2021-11-16 UTC"))
+  dim(dates) <- c(sday = 1, sweek = 1, syear = 2, time = 2)
+  time_bounds[] <- lapply(time_bounds, 
+                          function(x) {
+                            dim(x) <- dim(dates)
+                            return(x)
+                          })
+
+  labels <- vector(mode = "list", length = 2)
+  for (i in 1:length(start_date)) {
+    labels[[i]] <- .get_plot_time_labels(
+      dates = dates,
+      horizon = horizon,
+      init_date = init_date,
+      start_date = start_date,
+      i_syear = start_date[i],
+      time_bounds = time_bounds
+    )
+  }
+
+  # Results
+  expect_equal(
+    labels[[1]]$years,
+    labels[[2]]$years,
+    NULL
+  )
+
+  expect_equal(
+    labels[[1]]$time_labels[1],
+    "November 2020 to February 2021"
+  )
+  expect_equal(
+    labels[[1]]$time_labels[2],
+    "November to December 2020"
+  )
+  expect_equal(
+    labels[[2]]$time_labels[1],
+    "November 2021 to February 2022"
+  )
+  expect_equal(
+    labels[[2]]$time_labels[2],
+    "November to December 2021"
+  )
+})
+
+
+test_that("6. Decadal, month-by-month data", {
+
+  horizon <- "decadal"
+  init_date <- 1
+  start_date <- c("2020", "2021")
+
+  dates <- as.POSIXct(c("2020-11-16 00:00:00 UTC", "2021-11-16 00:00:00 UTC",
+                        "2020-12-16 12:00:00 UTC", "2021-12-16 12:00:00 UTC",
+                        "2021-01-16 12:00:00 UTC", "2022-01-16 12:00:00 UTC",
+                        "2021-02-15 00:00:00 UTC", "2022-02-15 00:00:00 UTC"))
+  dim(dates) <- c(sday = 1, sweek = 1, syear = 2, time = 4)
+
+  labels <- vector(mode = "list", length = 2)
+  for (i in 1:length(start_date)) {
+    labels[[i]] <- .get_plot_time_labels(
+      dates = dates,
+      horizon = horizon,
+      init_date = init_date,
+      start_date = start_date,
+      i_syear = start_date[i],
+      time_bounds = NULL
+    )
+  }
+
+  # Results
+  expect_equal(
+    labels[[1]]$years,
+    c(2020, 2020, 2021, 2021)
+  )
+  expect_equal(
+    labels[[2]]$years,
+    c(2021, 2021, 2022, 2022)
+  )
+
+  expect_equal(
+    as.vector(labels[[1]]$time_labels),
+    as.vector(labels[[2]]$time_labels),
+    c("November", "December", "January", "February")
+  )
+})