Newer
Older
#'Plots the observed weekly means and climatology of a timeseries data
#'
#'@description This function plots the observed weekly means and climatology of
#'a timeseries data using ggplot package. It compares the weekly climatology in
#'a specified period (reference period) to the observed conditions during the
#'target period analyzed in the case study (included in the reference period).
#'
#'@param data A multidimensional array with named dimensions with at least sdate
#' and time dimensions containing observed daily data. It can also be a
#' dataframe with computed percentiles as input for ggplot. If it's a
#' dataframe, it must contain the following column names: 'week', 'clim',
#' 'p10', 'p90', 'p33', 'p66', 'week_mean', 'day' and 'data'.
#'@param first_date The first date of the observed values of timeseries. It can
#' be of class 'Date', 'POSIXct' or a character string in the format
#' 'yyyy-mm-dd'. If parameter 'data_years' is not provided, it must be a date
#' included in the reference period.
#'@param last_date Optional parameter indicating the last date of the target
#' period of the daily timeseries. It can be of class 'Date', 'POSIXct' or a
#' character string in the format 'yyyy-mm-dd'. If it is NULL, the last date of
#' the daily timeseries will be set as the last date of 'data'. As the data is
#' plotted by weeks, only full groups of 7 days will be plotted. If the last
#' date of the timeseries is not a multiple of 7 days, the last week will
#' not be plotted.
#'@param ref_period A vector of numeric values indicating the years of the
#' reference period. If parameter 'data_years' is not specified, it must
#' be of the same length of dimension 'sdate_dim' of parameter 'data'.
#'@param data_years A vector of numeric values indicating the years of the
#' data. It must be of the same length of dimension 'sdate_dim' of parameter
#' 'data'. It is optional, if not specified, all the years will be used as the
#' target period.
#'@param time_dim A character string indicating the daily time dimension name.
#' The default value is 'time'.
#'@param sdate_dim A character string indicating the start year dimension name.
#' The default value is 'sdate'.
#'@param title The text for the top title of the plot. It is NULL by default.
#'@param subtitle The text for the subtitle of the plot. It is NULL bu default.
#'@param ytitle Character string to be drawn as y-axis title. It is NULL by
#' default.
#'@param legend A logical value indicating whether a legend should be included
#' in the plot. If it is TRUE or NA, the legend will be included. If it is
#' FALSE, the legend will not be included. It is TRUE by default.
#'@param palette A palette name from the R Color Brewer’s package. The default
#' value is 'Blues'.
#'@param fileout A character string indicating the file name where to save the
#' plot. If not specified (default) a graphics device will pop up.
#'@param device A character string indicating the device to use. Can either be
#' a device function (e.g. png), or one of "eps", "ps", "tex" (pictex),
#' "pdf", "jpeg", "tiff", "png", "bmp", "svg" or "wmf" (windows only).
#'@param width A numeric value of the plot width in units ("in", "cm", "mm", or
#' "px"). It is set to 8 by default.
#'@param height A numeric value of the plot height in units ("in", "cm", "mm",
#' or "px"). It is set to 6 by default.
#'@param units Units of the size of the device (file or window) to plot in.
#' Inches (’in’) by default.
#'@param dpi A numeric value of the plot resolution. It is set to 300 by
#' default.
#'
#'@return A ggplot object containing the plot.
#'
#'@examples
#'data <- array(rnorm(49*20*3, 274), dim = c(time = 49, sdate = 20, member = 3))
#'PlotWeeklyClim(data = data, first_date = '2002-08-09',
#' last_date = '2002-09-15', ref_period = 2010:2019,
#' data_years = 2000:2019, time_dim = 'time', sdate_dim = 'sdate',
#' title = "Observed weekly means and climatology",
#' subtitle = "Target years: 2010 to 2019",
#' ytitle = paste0('tas', " (", "deg.C", ")"))
#'
#'@import multiApply
#'@import lubridate
#'@import ggplot2
#'@import RColorBrewer
#'@import scales
#'@importFrom ClimProjDiags Subset
#'@importFrom s2dv MeanDims
#'@export
PlotWeeklyClim <- function(data, first_date, ref_period, last_date = NULL,
data_years = NULL, time_dim = 'time',
sdate_dim = 'sdate', title = NULL, subtitle = NULL,
ytitle = NULL, legend = TRUE,
palette = "Blues", fileout = NULL, device = NULL,
width = 8, height = 6, units = 'in', dpi = 300) {
# Check input arguments
# data
if (is.array(data)) {
if (is.null(names(dim(data)))) {
stop("Parameter 'data' must have named dimensions.")
}
is_array <- TRUE
} else if (is.data.frame(data)) {
col_names <- c("week", "clim", "p10", "p90", "p33", "p66",
"week_mean", "day", "data")
stop(paste0("If parameter 'data' is a data frame, it must contain the ",
"following column names: 'week', 'clim', 'p10', 'p90', 'p33', ",
"'p66', 'week_mean', 'day' and 'data'."))
}
is_array <- FALSE
} else {
stop("Parameter 'data' must be an array or a data frame.")
}
if (is_array) {
# time_dim
if (!is.character(time_dim)) {
stop("Parameter 'time_dim' must be a character string.")
}
if (!all(time_dim %in% names(dim(data)))) {
stop("Parameter 'time_dim' is not found in 'data' dimension.")
}
if (dim(data)[time_dim] < 7) {
stop(paste0("Parameter 'data' must have the dimension 'time_dim' of ",
"length equal or grater than 7 to compute the weekly means."))
}
# sdate_dim
if (!is.character(sdate_dim)) {
stop("Parameter 'sdate_dim' must be a character string.")
}
if (!sdate_dim %in% names(dim(data))) {
warning(paste0("Parameter 'sdate_dim' is not found in 'data' dimension. ",
"A dimension of length 1 has been added."))
data <- InsertDim(data, 1, lendim = 1, name = sdate_dim)
}
# legend
if (!is.logical(legend)) {
stop("Parameter 'legend' must be a logical value.")
}
if (is.na(legend)) {
legend <- TRUE
} else if (legend) {
legend <- NA
}
# ref_period (1)
if (!is.numeric(ref_period)) {
stop("Parameter 'ref_period' must be numeric.")
}
# first_date
if ((!inherits(first_date, "POSIXct") & !inherits(first_date, "Date")) &&
(!is.character(first_date) | nchar(first_date) != 10)) {
stop(paste0("Parameter 'first_date' must be a character string ",
"indicating the date in the format 'yyyy-mm-dd', 'POSIXct' ",
"or 'Dates' class."))
first_date <- ymd(first_date)
target_year <- year(first_date)
taget_year_outside_reference <- FALSE
# data_years
if (!is.null(data_years)) {
if (!is.numeric(data_years)) {
stop("Parameter 'data_years' must be numeric.")
}
if (length(data_years) != dim(data)[sdate_dim]) {
stop(paste0("Parameter 'data_years' must have the same length as ",
"the dimension '", sdate_dim, "' of 'data'."))
}
if (!all(ref_period %in% data_years)) {
stop(paste0("The 'ref_period' must be included in the 'data_years' ",
"period."))
}
if (!any(target_year %in% data_years)) {
stop(paste0("Parameter 'first_date' must be a date included ",
"in the 'data_years' period."))
}
taget_year_outside_reference <- TRUE
} else {
# ref_period (2)
if (length(ref_period) != dim(data)[sdate_dim]) {
stop(paste0("Parameter 'ref_period' must have the same length as the ",
"dimension '", sdate_dim ,"' of 'data' if ",
"'data_years' is not provided."))
}
if (!any(target_year %in% ref_period)) {
stop(paste0("If parameter 'data_years' is NULL, parameter 'first_date' ",
"must be a date included in the 'ref_period' period."))
}
data_years <- ref_period
}
# last_date
if (!is.null(last_date)) {
if ((!inherits(last_date, "POSIXct") & !inherits(last_date, "Date")) &&
(!is.character(last_date) | nchar(last_date) != 10)) {
stop(paste0("Parameter 'last_date' must be a character string ",
"indicating the date in the format 'yyyy-mm-dd', 'POSIXct' ",
"or 'Dates' class."))
}
last_date <- ymd(last_date)
dates <- seq(first_date, last_date, by = "1 day")
if (length(dates) > dim(data)[time_dim]) {
warning(paste0("Parameter 'last_date' is greater than the last date ",
"of 'data'. The last date of 'data' will be used."))
dates <- seq(first_date, first_date + days(dim(data)[time_dim]-1), by = "1 day")
}
} else {
dates <- seq(first_date, first_date + days(dim(data)[time_dim]-1), by = "1 day")
}
index_first_date <- which(dates == first_date)
index_last_date <- length(dates) - (length(dates) %% 7)
last_date <- dates[index_last_date]
## Data preparation
# subset time_dim for weeks
data_subset <- Subset(data, along = time_dim,
indices = index_first_date:index_last_date)
# remove other dimensions
dims_subset <- names(dim(data_subset))[which(!names(dim(data_subset)) %in% c(time_dim, sdate_dim))]
if (!identical(dims_subset, character(0))) {
data_subset <- Subset(data_subset, dims_subset, as.list(rep(1, length(dims_subset))), drop = TRUE)
}
# observed daily data creation
daily <- Subset(data_subset, along = sdate_dim,
indices = which(data_years == target_year),
drop = TRUE)
if (taget_year_outside_reference) {
indexes_reference_period <- which(data_years %in% ref_period)
# remove values outside reference period for computing the means
data_subset <- Subset(data_subset, along = sdate_dim,
indices = indexes_reference_period)
}
## Weekly aggregations for reference period
weekly_aggre <- SplitDim(data_subset, split_dim = time_dim,
indices = sort(rep(1:(length(index_first_date:index_last_date)/7), 7)),
new_dim_name = 'week')
weekly_means <- MeanDims(weekly_aggre, time_dim)
weekly_clim <- MeanDims(weekly_means, sdate_dim)
weekly_p10 <- Apply(weekly_means, target_dims = sdate_dim,
fun = function(x) {quantile(x, 0.10)})$output1
weekly_p90 <- Apply(weekly_means, target_dims = sdate_dim,
fun = function(x) {quantile(x, 0.90)})$output1
weekly_p33 <- Apply(weekly_means, target_dims = sdate_dim,
fun = function(x) {quantile(x, 0.33)})$output1
weekly_p66 <- Apply(weekly_means, target_dims = sdate_dim,
fun = function(x) {quantile(x, 0.66)})$output1
clim <- p10 <- p90 <- p33 <- p66 <- NULL
weekly_data <- data.frame(clim = as.vector(weekly_clim),
p10 = as.vector(weekly_p10),
p90 = as.vector(weekly_p90),
p33 = as.vector(weekly_p33),
p66 = as.vector(weekly_p66),
week = 1:(length(index_first_date:index_last_date)/7))
## observations for target year
daily_data <- data.frame(day = seq(first_date, last_date, by = "1 day"),
data = daily,
week = sort(rep(1:(length(index_first_date:index_last_date)/7), 7)))
week_mean <- aggregate(data ~ week, data = daily_data, mean)
weekly_data <- cbind(weekly_data, week_mean$data)
colnames(weekly_data)[7] <- 'week_mean'
all <- merge(weekly_data, daily_data, by = 'week')
} else {
all <- data
}
# 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"),
alpha = 0.7, show.legend = legend) + # extremes clim
geom_ribbon(aes(ymin = p33, ymax = p66, group = week, fill = "p33-p66"),
alpha = 0.7, show.legend = legend) + # terciles clim
geom_line(aes(y = clim, group = week, color = "climatological mean",
linetype = "climatological mean"),
alpha = 1.0, size = 0.7, show.legend = legend) + # mean clim
geom_line(aes(y = data, color = "observed daily mean",
linetype = "observed daily mean"),
alpha = 1, size = 0.2, show.legend = legend) + # daily evolution
geom_line(aes(y = week_mean, group = week, color = "observed weekly mean",
linetype = "observed weekly mean"),
alpha = 1, size = 0.7, show.legend = legend) + # weekly evolution
theme_bw() + ylab(ytitle) + xlab(NULL) +
ggtitle(title, subtitle = subtitle) +
scale_x_date(breaks = seq(min(all$day), max(all$day), by = "7 days"),
minor_breaks = NULL, expand = c(0.03, 0.03),
labels = date_format("%d %b %Y")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid.major = element_line(size = 0.5, linetype = 'solid',
colour = "gray92"),
panel.grid.minor = element_line(size = 0.25, linetype = 'solid',
colour = "gray92"),
legend.spacing = unit(-0.2, "cm")) +
scale_fill_manual(name = NULL,
values = c("p10-p90" = cols[3], "p33-p66" = cols[4])) +
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
scale_color_manual(name = NULL, values = c("climatological mean" = cols[5],
"observed daily mean" = "grey20",
"observed weekly mean" = "black")) +
scale_linetype_manual(name = NULL, values = c("climatological mean" = "solid",
"observed daily mean" = "dashed",
"observed weekly mean" = "solid"),
guide = guide_legend(override.aes = list(lwd = c(0.7, 0.2, 0.7)))) +
guides(fill = guide_legend(order = 1))
# Return the ggplot object
if (is.null(fileout)) {
return(p)
} else {
ggsave(filename = fileout, plot = p, device = device, height = height,
width = width, units = units, dpi = dpi)
}
}