From 5d451445f3c4ee6f356b5aa4b9b8b97326b923e7 Mon Sep 17 00:00:00 2001 From: vagudets Date: Wed, 28 Aug 2024 16:23:43 +0200 Subject: [PATCH 1/3] Set first forecast time as January for decadal models initialized in November --- modules/Loading/R/helper_loading_decadal.R | 4 +--- modules/Loading/R/load_decadal.R | 20 ++++++++++++++++---- 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/modules/Loading/R/helper_loading_decadal.R b/modules/Loading/R/helper_loading_decadal.R index 65058a6b..95f29c8c 100644 --- a/modules/Loading/R/helper_loading_decadal.R +++ b/modules/Loading/R/helper_loading_decadal.R @@ -20,11 +20,9 @@ get_daily_time_ind <- function(ftimemin, ftimemax, initial_month, sdates, calend total_months <- ifelse(total_months > 12, total_months %% 12, total_months) total_months <- ifelse(total_months == 0, 12, total_months) - if (calendar == '360-day') { + if (calendar == '360_day') { daily_time_ind <- ((ftimemin - 1) * 30 + 1): (((ftimemin - 1) * 30) + (ftimemax - ftimemin + 1) * 30) - } else { - total_days_sum <- sum(days_in_month(total_months)) if (ftimemin == 1) { daily_time_ind <- 1:total_days_sum diff --git a/modules/Loading/R/load_decadal.R b/modules/Loading/R/load_decadal.R index 21df4223..80d6e474 100644 --- a/modules/Loading/R/load_decadal.R +++ b/modules/Loading/R/load_decadal.R @@ -38,12 +38,24 @@ load_decadal <- function(recipe) { sdates_hcst <- as.numeric(recipe$Analysis$Time$hcst_start):as.numeric(recipe$Analysis$Time$hcst_end) #1960:2015 sdates_fcst <- recipe$Analysis$Time$fcst + # Some models are initialized in November, but all models should load + # data starting in January, so the first two months are ignored for + # November initializations. + ## TODO: Consider May initialized models + ftime_min <- as.numeric(recipe$Analysis$Time$ftime_min) + ftime_max <- as.numeric(recipe$Analysis$Time$ftime_max) + initial_month <- as.numeric(archive$System[[exp.name]]$initial_month) + + if (initial_month == 11) { + ftime_min <- ftime_min + 2 + ftime_max <- ftime_max + 2 + } if (store.freq == "monthly_mean") { - time_ind <- (as.numeric(recipe$Analysis$Time$ftime_min):as.numeric(recipe$Analysis$Time$ftime_max)) + time_ind <- ftime_min:ftime_max } else if (store.freq %in% c("daily", "daily_mean")) { - time_ind <- get_daily_time_ind(ftimemin = as.numeric(recipe$Analysis$Time$ftime_min), - ftimemax = as.numeric(recipe$Analysis$Time$ftime_max), - initial_month = archive$System[[exp.name]]$initial_month, + time_ind <- get_daily_time_ind(ftimemin = ftime_min, + ftimemax = ftime_max, + initial_month = initial_month, sdates = sdates_hcst, calendar = archive$System[[exp.name]]$calendar) store.freq <- "daily_mean" -- GitLab From 2441f85fef2789993945e56f2bfba298614937ab Mon Sep 17 00:00:00 2001 From: vagudets Date: Wed, 28 Aug 2024 16:24:23 +0200 Subject: [PATCH 2/3] Bugfix: Add conditions for decadal plots and fix forecast time definition --- modules/Visualization/R/plot_ensemble_mean.R | 4 +- modules/Visualization/R/plot_metrics.R | 44 +++++++++++++++++++- 2 files changed, 44 insertions(+), 4 deletions(-) diff --git a/modules/Visualization/R/plot_ensemble_mean.R b/modules/Visualization/R/plot_ensemble_mean.R index 08a894ea..19521f7d 100644 --- a/modules/Visualization/R/plot_ensemble_mean.R +++ b/modules/Visualization/R/plot_ensemble_mean.R @@ -123,7 +123,7 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o "\n", "Forecast Ensemble Mean / ", "Init.: ", i_syear) } if (is.null(attributes(fcst$attrs$time_bounds))) { - if (tolower(recipe$Analysis$Horizon) == "seasonal") { + 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") { @@ -251,7 +251,7 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o toptitle <- paste0(system_name, " / ", str_to_title(var_long_name), "\n", "Ensemble Mean / ", - time_labels[i], + time_labels[i], " ", years[i], " / Start date: ", i_syear) } diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R index 15ce05e6..8904ed57 100644 --- a/modules/Visualization/R/plot_metrics.R +++ b/modules/Visualization/R/plot_metrics.R @@ -52,9 +52,12 @@ plot_metrics <- function(recipe, data_cube, metrics, weeks <- paste0(0, ftime_min:ftime_max) # This week_label appears on the name of the file. It's just the start date. week_label <- recipe$Analysis$Time$sdate - } else { ## ? + } else { # Decadal init_month <- 1 init_week <- 1 + months <- lubridate::month(Subset(data_cube$attrs$Dates, + "syear", indices = 1), + label = T, abb = F,locale = "en_GB") } if (!is.null(recipe$Analysis$Workflow$Visualization$projection)) { @@ -358,8 +361,45 @@ plot_metrics <- function(recipe, data_cube, metrics, year(ymd(start_date)))) # "/ valid week", format(weeks[i], # "%Y-%m-%d"), "/", hcst_period) + } else if (recipe$Analysis$Horizon == "decadal") { + # Case without time aggregation: + if (is.null(attributes(data_cube$attrs$time_bounds))) { + forecast_time <- match(months[i], month.name) + floor((i - 1)/12)*12 + forecast_time <- sprintf("%02d", forecast_time) + # Define plot title + toptitle <- paste(system_name, "/", + str_to_title(var_long_name), + "\n", display_name, "/", months[i], "/", + hcst_period) + } else { + if (length(attributes(data_cube$attrs$time_bounds)$plotting_attr) > 1) { + forecast_time_ini <- attributes(data_cube$attrs$time_bounds)$plotting_attr$ini[i] + forecast_time_end <- attributes(data_cube$attrs$time_bounds)$plotting_attr$end[i] + # labels for file name: + forecast_time <- paste0(forecast_time_ini, "-", forecast_time_end) + # title names: + 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 <- month.name[forecast_time_end] + toptitle <- paste(system_name, "/", + str_to_title(var_long_name), + "\n", display_name, "/", forecast_time_ini, "to", + forecast_time_end, "/", + hcst_period) + } else { + forecast_time <- attributes(data_cube$attrs$time_bounds)$plotting_attr[[1]][i] + toptitle <- paste(system_name, "/", + str_to_title(var_long_name), + "\n", display_name, "/", + forecast_time, "/", + hcst_period) + } + } } else { - warning("Plotting decadal?") + warning("Unknown time horizon?") } # Modify base arguments base_args[[1]] <- metric[i, , ] -- GitLab From 6530fa651c75655e146ce444f5e2ecd56f5de4fe Mon Sep 17 00:00:00 2001 From: vagudets Date: Wed, 28 Aug 2024 16:25:01 +0200 Subject: [PATCH 3/3] Update decadal daily unit test --- tests/recipes/recipe-decadal_daily_1.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/recipes/recipe-decadal_daily_1.yml b/tests/recipes/recipe-decadal_daily_1.yml index 88b87622..f9a343a4 100644 --- a/tests/recipes/recipe-decadal_daily_1.yml +++ b/tests/recipes/recipe-decadal_daily_1.yml @@ -18,8 +18,8 @@ Analysis: hcst_start: 1990 hcst_end: 1992 season: 'Annual' - ftime_min: 3 - ftime_max: 5 + ftime_min: 1 + ftime_max: 3 Region: latmin: 10 #-90 latmax: 20 #90 -- GitLab