diff --git a/NAMESPACE b/NAMESPACE index 1531631ff8198bd6ece591b6e0cbc9516336852a..a6a423293b3fcd6149386ec9e68b85f21501a845 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,13 +4,64 @@ export(ClimColors) export(ClimPalette) export(ColorBarContinuous) export(ColorBarDiscrete) +export(PlotCombinedMap) +export(PlotEquiMap) +export(PlotForecastPDF) +export(PlotMostLikelyQuantileMap) +export(PlotPDFsOLE) export(PlotRobinson) +export(PlotTriangles4Categories) +export(PlotWeeklyClim) export(ShapeToMask) +import(RColorBrewer) import(cowplot) import(easyNCDF) import(ggplot2) +import(graphics) +import(maps) +import(multiApply) import(rnaturalearth) +import(scales) import(sf) +importFrom(CSTools,SplitDim) +importFrom(ClimProjDiags,Subset) +importFrom(RColorBrewer,brewer.pal) +importFrom(data.table,CJ) +importFrom(data.table,data.table) +importFrom(data.table,setkey) +importFrom(grDevices,adjustcolor) +importFrom(grDevices,bmp) importFrom(grDevices,col2rgb) importFrom(grDevices,colorRampPalette) +importFrom(grDevices,dev.cur) +importFrom(grDevices,dev.new) +importFrom(grDevices,dev.off) +importFrom(grDevices,gray) +importFrom(grDevices,hcl) +importFrom(grDevices,jpeg) +importFrom(grDevices,pdf) +importFrom(grDevices,png) +importFrom(grDevices,postscript) importFrom(grDevices,rgb) +importFrom(grDevices,svg) +importFrom(grDevices,tiff) +importFrom(graphics,axis) +importFrom(graphics,box) +importFrom(graphics,image) +importFrom(graphics,layout) +importFrom(graphics,mtext) +importFrom(graphics,par) +importFrom(graphics,plot) +importFrom(graphics,plot.new) +importFrom(graphics,points) +importFrom(graphics,polygon) +importFrom(graphics,text) +importFrom(graphics,title) +importFrom(lubridate,ymd) +importFrom(maps,map) +importFrom(plyr,.) +importFrom(plyr,dlply) +importFrom(reshape2,melt) +importFrom(s2dv,InsertDim) +importFrom(s2dv,MeanDims) +importFrom(stats,cor) diff --git a/R/ColorBarContinuous.R b/R/ColorBarContinuous.R index d19b6f567de1a0abdf482955bb2d9c67cb0dca5b..b3c74df949e7a930b9b74afe211b3b47054b8990 100644 --- a/R/ColorBarContinuous.R +++ b/R/ColorBarContinuous.R @@ -288,7 +288,7 @@ ColorBarContinuous <- function(brks = NULL, cols = NULL, vertical = TRUE, triangle_ends <- triangle_ends } } - if (plot) { + if (plot && !is.null(var_limits)) { if ((bar_limits[1] >= var_limits[1]) && !triangle_ends[1]) { warning("There are variable values smaller or equal to the lower limit ", "of the colour bar and the lower triangle end has been ", diff --git a/R/PlotCombinedMap.R b/R/PlotCombinedMap.R new file mode 100644 index 0000000000000000000000000000000000000000..871a9c838a7fccd68f513f7a6b6a9022c7a60ef3 --- /dev/null +++ b/R/PlotCombinedMap.R @@ -0,0 +1,529 @@ +#'Plot Multiple Lon-Lat Variables In a Single Map According to a Decision Function +#' +#'@description Plot a number a two dimensional matrices with (longitude, +#'latitude) dimensions on a single map with the cylindrical equidistant +#'latitude and longitude projection. +#' +#'@author Nicolau Manubens, \email{nicolau.manubens@bsc.es} +#'@author Veronica Torralba, \email{veronica.torralba@bsc.es} +#' +#'@param maps List of matrices to plot, each with (longitude, latitude) +#' dimensions, or 3-dimensional array with the dimensions (longitude, latitude, +#' map). Dimension names are required. +#'@param lon Vector of longitudes. Must match the length of the corresponding +#' dimension in 'maps'. +#'@param lat Vector of latitudes. Must match the length of the corresponding +#' dimension in 'maps'. +#'@param map_select_fun Function that selects, for each grid point, which value +#' to take among all the provided maps. This function receives as input a +#' vector of values for a same grid point for all the provided maps, and must +#' return a single selected value (not its index!) or NA. For example, the +#' \code{min} and \code{max} functions are accepted. +#'@param display_range Range of values to be displayed for all the maps. This +#' must be a numeric vector c(range min, range max). The values in the +#' parameter 'maps' can go beyond the limits specified in this range. If the +#' selected value for a given grid point (according to 'map_select_fun') falls +#' outside the range, it will be coloured with 'col_unknown_map'. +#'@param map_dim Optional name for the dimension of 'maps' along which the +#' multiple maps are arranged. Only applies when 'maps' is provided as a +#' 3-dimensional array. Takes the value 'map' by default. +#'@param brks Colour levels to be sent to PlotEquiMap. This parameter is +#' optional and adjusted automatically by the function. +#'@param cols List of vectors of colours to be sent to PlotEquiMap for the +#' colour bar of each map. This parameter is optional and adjusted +#' automatically by the function (up to 5 maps). The colours provided for each +#' colour bar will be automatically interpolated to match the number of breaks. +#' Each item in this list can be named, and the name will be used as title for +#' the corresponding colour bar (equivalent to the parameter 'bar_titles'). +#'@param col_unknown_map Colour to use to paint the grid cells for which a map +#' is not possible to be chosen according to 'map_select_fun' or for those +#' values that go beyond 'display_range'. Takes the value 'white' by default. +#'@param mask Optional numeric array with dimensions (latitude, longitude), with +#' values in the range [0, 1], indicating the opacity of the mask over each +#' grid point. Cells with a 0 will result in no mask, whereas cells with a 1 +#' will result in a totally opaque superimposed pixel coloured in 'col_mask'. +#'@param col_mask Colour to be used for the superimposed mask (if specified in +#' 'mask'). Takes the value 'grey' by default. +#'@param dots Array of same dimensions as 'var' or with dimensions +#' c(n, dim(var)), where n is the number of dot/symbol layers to add to the +#' plot. A value of TRUE at a grid cell will draw a dot/symbol on the +#' corresponding square of the plot. By default all layers provided in 'dots' +#' are plotted with dots, but a symbol can be specified for each of the +#' layers via the parameter 'dot_symbol'. +#'@param bar_titles Optional vector of character strings providing the titles to +#' be shown on top of each of the colour bars. +#'@param legend_scale Scale factor for the size of the colour bar labels. Takes +#' 1 by default. +#'@param cex_bar_titles Scale factor for the sizes of the bar titles. Takes 1.5 +#' by default. +#'@param plot_margin Numeric vector of length 4 for the margin sizes in the +#' following order: bottom, left, top, and right. If not specified, use the +#' default of par("mar"), c(5.1, 4.1, 4.1, 2.1). Used as 'margin_scale' in +#' PlotEquiMap. +#'@param fileout File where to save the plot. If not specified (default) a +#' graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp +#' and tiff +#'@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 drawleg Where to draw the common colour bar. Can take values TRUE, +#' FALSE or:\cr +#' 'up', 'u', 'U', 'top', 't', 'T', 'north', 'n', 'N'\cr +#' 'down', 'd', 'D', 'bottom', 'b', 'B', 'south', 's', 'S' (default)\cr +#' 'right', 'r', 'R', 'east', 'e', 'E'\cr +#' 'left', 'l', 'L', 'west', 'w', 'W' +#'@param return_leg A logical value indicating if the color bars information +#' should be returned by the function. If TRUE, the function doesn't plot the +#' color bars but still creates the layout with color bar areas, and the +#' arguments for GradientCatsColorBar() or ColorBarContinuous() will be +#' returned. It is convenient for users to adjust the color bars manually. The +#' default is FALSE, the color bars will be plotted directly. +#'@param ... Additional parameters to be passed on to \code{PlotEquiMap}. +#' +#'@examples +#'# Simple example +#'x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200 +#'a <- x * 0.6 +#'b <- (1 - x) * 0.6 +#'c <- 1 - (a + b) +#'lons <- seq(0, 359.5, length = 20) +#'lats <- seq(-89.5, 89.5, length = 10) +#'\dontrun{ +#'PlotCombinedMap(list(a, b, c), lons, lats, +#' toptitle = 'Maximum map', +#' map_select_fun = max, +#' display_range = c(0, 1), +#' bar_titles = paste('% of belonging to', c('a', 'b', 'c')), +#' brks = 20, width = 12, height = 10) +#'} +#' +#'Lon <- c(0:40, 350:359) +#'Lat <- 51:26 +#'data <- rnorm(51 * 26 * 3) +#'dim(data) <- c(map = 3, lon = 51, lat = 26) +#'mask <- sample(c(0,1), replace = TRUE, size = 51 * 26) +#'dim(mask) <- c(lat = 26, lon = 51) +#'\dontrun{ +#'PlotCombinedMap(data, lon = Lon, lat = Lat, map_select_fun = max, +#' display_range = range(data), mask = mask, +#' width = 14, height = 10) +#'} +#' +#'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap} +#' +#'@importFrom maps map +#'@importFrom graphics box image layout mtext par plot.new +#'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off +#' hcl jpeg pdf png postscript svg tiff +#'@export +PlotCombinedMap <- function(maps, lon, lat, + map_select_fun, display_range, + map_dim = 'map', + brks = NULL, cols = NULL, + bar_limits = NULL, triangle_ends = c(F, F), col_inf = NULL, col_sup = NULL, + col_unknown_map = 'white', + mask = NULL, col_mask = 'grey', + dots = NULL, + bar_titles = NULL, legend_scale = 1, + cex_bar_titles = 1.5, + plot_margin = NULL, bar_extra_margin = c(2, 0, 2, 0), + fileout = NULL, width = 8, height = 5, + size_units = 'in', res = 100, drawleg = T, return_leg = FALSE, + ...) { + args <- list(...) + + # 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 + } + + # Check probs + error <- FALSE + # Change list into an array + if (is.list(maps)) { + if (length(maps) < 1) { + stop("Parameter 'maps' must be of length >= 1 if provided as a list.") + } + check_fun <- function(x) { + is.numeric(x) && (length(dim(x)) == 2) + } + if (!all(sapply(maps, check_fun))) { + error <- TRUE + } + ref_dims <- dim(maps[[1]]) + equal_dims <- all(sapply(maps, function(x) identical(dim(x), ref_dims))) + if (!equal_dims) { + stop("All arrays in parameter 'maps' must have the same dimension ", + "sizes and names when 'maps' is provided as a list of arrays.") + } + num_maps <- length(maps) + maps <- unlist(maps) + dim(maps) <- c(ref_dims, map = num_maps) + map_dim <- 'map' + } + if (!is.numeric(maps)) { + error <- TRUE + } + if (is.null(dim(maps))) { + error <- TRUE + } + if (length(dim(maps)) != 3) { + error <- TRUE + } + if (error) { + stop("Parameter 'maps' must be either a numeric array with 3 dimensions ", + " or a list of numeric arrays of the same size with the 'lon' and ", + "'lat' dimensions.") + } + dimnames <- names(dim(maps)) + + # Check map_dim + if (is.character(map_dim)) { + if (is.null(dimnames)) { + stop("Specified a dimension name in 'map_dim' but no dimension names provided ", + "in 'maps'.") + } + map_dim <- which(dimnames == map_dim) + if (length(map_dim) < 1) { + stop("Dimension 'map_dim' not found in 'maps'.") + } else { + map_dim <- map_dim[1] + } + } else if (!is.numeric(map_dim)) { + stop("Parameter 'map_dim' must be either a numeric value or a ", + "dimension name.") + } + if (length(map_dim) != 1) { + stop("Parameter 'map_dim' must be of length 1.") + } + map_dim <- round(map_dim) + + # Work out lon_dim and lat_dim + lon_dim <- NULL + if (!is.null(dimnames)) { + lon_dim <- which(dimnames %in% c('lon', 'longitude'))[1] + } + if (length(lon_dim) < 1) { + lon_dim <- (1:3)[-map_dim][1] + } + lon_dim <- round(lon_dim) + + lat_dim <- NULL + if (!is.null(dimnames)) { + lat_dim <- which(dimnames %in% c('lat', 'latitude'))[1] + } + if (length(lat_dim) < 1) { + lat_dim <- (1:3)[-map_dim][2] + } + lat_dim <- round(lat_dim) + + # Check lon + if (!is.numeric(lon)) { + stop("Parameter 'lon' must be a numeric vector.") + } + if (length(lon) != dim(maps)[lon_dim]) { + stop("Parameter 'lon' does not match the longitude dimension in 'maps'.") + } + + # Check lat + if (!is.numeric(lat)) { + stop("Parameter 'lat' must be a numeric vector.") + } + if (length(lat) != dim(maps)[lat_dim]) { + stop("Parameter 'lat' does not match the longitude dimension in 'maps'.") + } + + # Check map_select_fun + if (is.numeric(map_select_fun)) { + if (length(dim(map_select_fun)) != 2) { + stop("Parameter 'map_select_fun' must be an array with dimensions ", + "'lon' and 'lat' if provided as an array.") + } + if (!identical(dim(map_select_fun), dim(maps)[-map_dim])) { + stop("The dimensions 'lon' and 'lat' in the 'map_select_fun' array must ", + "have the same size, name and order as in the 'maps' parameter.") + } + } + if (!is.function(map_select_fun)) { + stop("The parameter 'map_select_fun' must be a function or a numeric array.") + } + + # Generate the desired brks and cols. Only nmap, brks, cols, bar_limits, and + # bar_titles matter here because plot = F. + var_limits_maps <- range(maps, na.rm = TRUE) + if (is.null(bar_limits)) bar_limits <- display_range + nmap <- dim(maps)[map_dim] + colorbar <- GradientCatsColorBar(nmap = nmap, + brks = brks, cols = cols, vertical = FALSE, + subsampleg = NULL, bar_limits = bar_limits, + var_limits = var_limits_maps, + triangle_ends = triangle_ends, col_inf = col_inf, col_sup = col_sup, plot = FALSE, draw_separators = TRUE, + bar_titles = bar_titles, title_scale = cex_bar_titles, label_scale = legend_scale * 1.5, + extra_margin = bar_extra_margin) + + # Check legend_scale + if (!is.numeric(legend_scale)) { + stop("Parameter 'legend_scale' must be numeric.") + } + + # Check col_unknown_map + if (!is.character(col_unknown_map)) { + stop("Parameter 'col_unknown_map' must be a character string.") + } + + # Check col_mask + if (!is.character(col_mask)) { + stop("Parameter 'col_mask' must be a character string.") + } + + # Check mask + if (!is.null(mask)) { + if (!is.numeric(mask)) { + stop("Parameter 'mask' must be numeric.") + } + if (length(dim(mask)) != 2) { + stop("Parameter 'mask' must have two dimensions.") + } + if ((dim(mask)[1] != dim(maps)[lat_dim]) || + (dim(mask)[2] != dim(maps)[lon_dim])) { + stop("Parameter 'mask' must have dimensions c(lat, lon).") + } + } + # Check dots + if (!is.null(dots)) { + if (length(dim(dots)) != 2) { + stop("Parameter 'mask' must have two dimensions.") + } + if ((dim(dots)[1] != dim(maps)[lat_dim]) || + (dim(dots)[2] != dim(maps)[lon_dim])) { + stop("Parameter 'mask' must have dimensions c(lat, lon).") + } + } + + #---------------------- + # Identify the most likely map + #---------------------- + #TODO: Consider col_inf + if (!is.null(colorbar$col_inf[[1]])) { + .warning("Lower triangle is not supported now. Please contact maintainer if you have this need.") + } + if (!is.null(colorbar$col_sup[[1]])) { + + brks_norm <- vector('list', length = nmap) + range_width <- vector('list', length = nmap) + slightly_tune_val <- vector('list', length = nmap) + for (ii in 1:nmap) { + brks_norm[[ii]] <- seq(0, 1, length.out = length(colorbar$brks[[ii]]) + 1) # add one break for col_sup + slightly_tune_val[[ii]] <- brks_norm[[ii]][2] / (length(brks_norm[[ii]]) * 2) + range_width[[ii]] <- diff(range(colorbar$brks[[ii]])) + } + ml_map <- apply(maps, c(lat_dim, lon_dim), function(x) { + if (any(is.na(x))) { + res <- NA + } else { + res <- which(x == map_select_fun(x)) + if (length(res) > 0) { + res <- res_ind <- res[1] + if (map_select_fun(x) < display_range[1] || map_select_fun(x) > display_range[2]) { + res <- -0.5 + } else { + if (map_select_fun(x) > tail(colorbar$brks[[res_ind]], 1)) { # col_sup + res <- res + 1 - slightly_tune_val[[res_ind]] + } else { + res <- res + ((map_select_fun(x) - colorbar$brks[[res_ind]][1]) / + range_width[[res_ind]] * ((length(brks_norm[[res_ind]]) - 2) / (length(brks_norm[[res_ind]]) - 1))) + if (map_select_fun(x) == colorbar$brks[[res_ind]][1]) { + res <- res + slightly_tune_val[[res_ind]] + } + } + } + } else { + res <- -0.5 + } + } + res + }) + + } else { + + brks_norm <- vector('list', length = nmap) + range_width <- vector('list', length = nmap) + slightly_tune_val <- vector('list', length = nmap) + for (ii in 1:nmap) { + brks_norm[[ii]] <- seq(0, 1, length.out = length(colorbar$brks[[ii]])) + slightly_tune_val[[ii]] <- brks_norm[[ii]][2] / (length(brks_norm[[ii]]) * 2) + range_width[[ii]] <- diff(range(colorbar$brks[[ii]])) + } + ml_map <- apply(maps, c(lat_dim, lon_dim), function(x) { + if (any(is.na(x))) { + res <- NA + } else { + res <- which(x == map_select_fun(x)) + if (length(res) > 0) { + res <- res_ind <- res[1] + if (map_select_fun(x) < display_range[1] || + map_select_fun(x) > display_range[2]) { + res <- -0.5 + } else { + res <- res + ((map_select_fun(x) - colorbar$brks[[res_ind]][1]) / + range_width[[res_ind]]) + if (map_select_fun(x) == colorbar$brks[[res_ind]][1]) { + res <- res + slightly_tune_val[[res_ind]] + } + } + } else { + res <- -0.5 + } + } + res + }) + } + + nlat <- length(lat) + nlon <- length(lon) + + #---------------------- + # Set latitudes from minimum to maximum + #---------------------- + if (lat[1] > lat[nlat]) { + lat <- lat[nlat:1] + indices <- list(nlat:1, TRUE) + ml_map <- do.call("[", c(list(x = ml_map), indices)) + if (!is.null(mask)){ + mask <- mask[nlat:1, ] + } + if (!is.null(dots)){ + dots <- dots[nlat:1,] + } + } + + #---------------------- + # Set layout and parameters + #---------------------- + # 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) + } + #NOTE: I think plot.new() is not necessary in any case. +# plot.new() + #TODO: Don't hardcoded. Let users decide. + par(font.main = 1) + # If colorbars need to be plotted, re-define layout. + if (drawleg) { + layout(matrix(c(rep(1, nmap),2:(nmap + 1)), 2, nmap, byrow = TRUE), heights = c(6, 1.5)) + } + + #---------------------- + # Set colors and breaks and then PlotEquiMap + #---------------------- + if (!is.null(colorbar$col_sup[[1]])) { + tcols <- c(col_unknown_map, colorbar$cols[[1]], colorbar$col_sup[[1]]) + tbrks <- c(-1, brks_norm[[1]] + rep(1, each = length(brks_norm[[1]]))) + for (k in 2:nmap) { + tcols <- append(tcols, c(col_unknown_map, colorbar$cols[[k]], colorbar$col_sup[[k]])) + tbrks <- c(tbrks, brks_norm[[k]] + rep(k, each = length(brks_norm[[k]]))) + } + } else { # original code + tcols <- c(col_unknown_map, colorbar$cols[[1]]) + tbrks <- c(-1, brks_norm[[1]] + rep(1, each = length(brks_norm[[1]]))) + for (k in 2:nmap) { + tcols <- append(tcols, c(col_unknown_map, colorbar$cols[[k]])) + tbrks <- c(tbrks, brks_norm[[k]] + rep(k, each = length(brks_norm[[k]]))) + } + } + + if (is.null(plot_margin)) { + plot_margin <- c(5, 4, 4, 2) + 0.1 # default of par()$mar + } + + PlotEquiMap(var = ml_map, lon = lon, lat = lat, + brks = tbrks, cols = tcols, drawleg = FALSE, + filled.continents = FALSE, dots = dots, margin_scale = plot_margin, ...) + + #---------------------- + # Add overplot on top + #---------------------- + if (!is.null(mask)) { + dims_mask <- dim(mask) + latb <- sort(lat, index.return = TRUE) + dlon <- lon[2:dims_mask[2]] - lon[1:(dims_mask[2] - 1)] + wher <- which(dlon > (mean(dlon) + 1)) + if (length(wher) > 0) { + lon[(wher + 1):dims_mask[2]] <- lon[(wher + 1):dims_mask[2]] - 360 + } + lonb <- sort(lon, index.return = TRUE) + + cols_mask <- sapply(seq(from = 0, to = 1, length.out = 10), + function(x) adjustcolor(col_mask, alpha.f = x)) + image(lonb$x, latb$x, t(mask)[lonb$ix, latb$ix], + axes = FALSE, col = cols_mask, + breaks = seq(from = 0, to = 1, by = 0.1), + xlab='', ylab='', add = TRUE, xpd = TRUE) + if (!exists('coast_color')) { + coast_color <- 'black' + } + if (min(lon) < 0) { + map('world', interior = FALSE, add = TRUE, lwd = 1, col = coast_color) # Low resolution world map (lon -180 to 180). + } else { + map('world2', interior = FALSE, add = TRUE, lwd = 1, col = coast_color) # Low resolution world map (lon 0 to 360). + } + box() + } + + #---------------------- + # Add colorbars + #---------------------- + if ('toptitle' %in% names(args)) { + size_title <- 1 + if ('title_scale' %in% names(args)) { + size_title <- args[['title_scale']] + } + old_mar <- par('mar') + old_mar[3] <- old_mar[3] - (2 * size_title + 1) + par(mar = old_mar) + } + + if (drawleg & !return_leg) { + GradientCatsColorBar(nmap = dim(maps)[map_dim], + brks = colorbar$brks, cols = colorbar$cols, vertical = FALSE, + subsampleg = NULL, bar_limits = bar_limits, + var_limits = var_limits_maps, + triangle_ends = triangle_ends, col_inf = colorbar$col_inf, col_sup = colorbar$col_sup, + plot = TRUE, draw_separators = TRUE, + bar_titles = bar_titles, title_scale = cex_bar_titles, label_scale = legend_scale * 1.5, + extra_margin = bar_extra_margin) + } + + if (!return_leg) { + # If the graphic was saved to file, close the connection with the device + if (!is.null(fileout)) dev.off() + } + + if (return_leg) { + tmp <- list(nmap = dim(maps)[map_dim], + brks = colorbar$brks, cols = colorbar$cols, vertical = FALSE, + subsampleg = NULL, bar_limits = bar_limits, + var_limits = var_limits_maps, + triangle_ends = triangle_ends, col_inf = colorbar$col_inf, col_sup = colorbar$col_sup, + plot = TRUE, draw_separators = TRUE, + bar_titles = bar_titles, title_scale = cex_bar_titles, label_scale = legend_scale * 1.5, + extra_margin = bar_extra_margin) + .warning("The device is not off yet. Use dev.off() after plotting the color bars.") + return(tmp) + #NOTE: The device is not off! Can keep plotting the color bars. + } + +} + diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R new file mode 100644 index 0000000000000000000000000000000000000000..086880954d81e11f1e5679a9b53dd516971a2978 --- /dev/null +++ b/R/PlotForecastPDF.R @@ -0,0 +1,579 @@ +#'Plot one or multiple ensemble forecast pdfs for the same event +#' +#'@author Llorenç Lledó \email{llledo@bsc.es} +#'@description This function plots the probability distribution function of +#'several ensemble forecasts. Separate panels are used to plot forecasts valid +#'or initialized at different times or by different models or even at different +#'locations. Probabilities for tercile categories are computed, plotted in +#'colors and annotated. An asterisk marks the tercile with higher probabilities. +#'Probabilities for extreme categories (above P90 and below P10) can also be +#'included as hatched areas. Individual ensemble members can be plotted as +#'jittered points. The observed value is optionally shown as a diamond. +#' +#'@param fcst A dataframe or array containing all the ensember members for each +#' forecast. If \code{'fcst'} is an array, it should have two labelled +#' dimensions, and one of them should be \code{'members'}. If \code{'fcsts'} is +#' a data.frame, each column shoul be a separate forecast, with the rows beeing +#' the different ensemble members. +#'@param tercile.limits An array or vector with P33 and P66 values that define +#' the tercile categories for each panel. Use an array of dimensions +#' (nforecasts,2) to define different terciles for each forecast panel, or a +#' vector with two elements to reuse the same tercile limits for all forecast +#' panels. +#'@param extreme.limits (optional) An array or vector with P10 and P90 values +#' that define the extreme categories for each panel. Use an array of +#' (nforecasts,2) to define different extreme limits for each forecast panel, +#' or a vector with two elements to reuse the same tercile limits for all +#' forecast panels. (Default: extreme categories are not shown). +#'@param obs (optional) A vector providing the observed values for each forecast +#' panel or a single value that will be reused for all forecast panels. +#' (Default: observation is not shown). +#'@param plotfile (optional) A filename (pdf, png...) where the plot will be +#' saved. (Default: the plot is not saved). +#'@param title A string with the plot title. +#'@param var.name A string with the variable name and units. +#'@param fcst.names (optional) An array of strings with the titles of each +#' individual forecast. +#'@param add.ensmemb Either to add the ensemble members \code{'above'} (default) +#' or \code{'below'} the pdf, or not (\code{'no'}). +#'@param color.set A selection of predefined color sets: use \code{'ggplot'} +#' (default) for blue/green/red, \code{'s2s4e'} for blue/grey/orange, +#' \code{'hydro'} for yellow/gray/blue (suitable for precipitation and +#' inflows) or the \code{"vitigeoss"} color set. +#'@param memb_dim A character string indicating the name of the member +#' dimension. +#' +#'@return A ggplot object containing the plot. +#'@examples +#'fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2)) +#'PlotForecastPDF(fcsts,c(-1,1)) +#'@importFrom data.table data.table +#'@importFrom data.table CJ +#'@importFrom data.table setkey +#'@import ggplot2 +#'@importFrom reshape2 melt +#'@importFrom plyr . +#'@importFrom plyr dlply +#'@importFrom s2dv InsertDim +#'@export +PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, + plotfile = NULL, title = "Set a title", var.name = "Varname (units)", + fcst.names = NULL, add.ensmemb = c("above", "below", "no"), + color.set = c("ggplot", "s2s4e", "hydro", "vitigeoss"), + memb_dim = 'member') { + value <- init <- extremes <- x <- ymin <- ymax <- tercile <- NULL + y <- xend <- yend <- yjitter <- MLT <- lab.pos <- NULL + ggColorHue <- function(n) { + hues <- seq(15, 375, length = n + 1) + hcl(h = hues, l = 65, c = 100)[1:n] + } + #------------------------ + # Define color sets + #------------------------ + color.set <- match.arg(color.set) + if (color.set == "s2s4e") { + colorFill <- c("#FF764D", "#b5b5b5", "#33BFD1") # AN, N, BN fill colors + colorHatch <- c("indianred3", "deepskyblue3") # AP90, BP10 line colors + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("red", "blue") # AP90, BP10 text colors + } else if (color.set == "hydro") { + colorFill <- c("#41CBC9", "#b5b5b5", "#FFAB38") + colorHatch <- c("deepskyblue3", "darkorange1") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("blue", "darkorange3") + } else if (color.set == "ggplot") { + colorFill <- ggColorHue(3) + colorHatch <- c("indianred3", "deepskyblue3") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("red", "blue") + } else if (color.set == "vitigeoss") { + colorFill <- rev(c("#007be2", "#acb2b5", "#f40000")) + colorHatch <- rev(c("#211b79", "#ae0003")) + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- colorHatch + } else { + stop("Parameter 'color.set' should be one of ggplot/s2s4e/hydro") + } + #------------------------ + # Check input arguments + #------------------------ + add.ensmemb <- match.arg(add.ensmemb) + #------------------------ + # Check fcst type and convert to data.frame if needed + #------------------------ + if (is.array(fcst)) { + if (!memb_dim %in% names(dim(fcst)) | length(dim(fcst)) != 2) { + stop("Parameter 'fcst' should be a two-dimensional array with labelled dimensions and one of them should be for member. The name of this dimension can be adjusted with 'memb_dim'.") + } + dim.members <- which(names(dim(fcst)) == memb_dim) + if (dim.members == 1) { + fcst.df <- data.frame(fcst) + } else { + fcst.df <- data.frame(t(fcst)) + } + } else if (is.data.frame(fcst)) { + fcst.df <- fcst + } else { + stop("Parameter 'fcst' should be an array or a data.frame") + } + npanels <- ncol(fcst.df) + #------------------------ + # Check observations + #------------------------ + if (!is.null(obs)) { + if (!is.vector(obs)){ + stop("Observations should be a vector") + } else if (!length(obs) %in% c(1, npanels)) { + stop("The number of observations should equal one or the number of forecasts") + } + } + #------------------------ + # Check tercile limits + #------------------------ + if (is.vector(tercile.limits)) { + if (length(tercile.limits) != 2) { + stop("Provide two tercile limits") + } + tercile.limits <- InsertDim(tercile.limits, 1, npanels, name = "new") + } else if (is.array(tercile.limits)) { + if (length(dim(tercile.limits)) == 2) { + if (dim(tercile.limits)[2] != 2) { + stop("Provide two tercile limits for each panel") + } + if (dim(tercile.limits)[1] != npanels) { + stop("The number of tercile limits does not match the number of forecasts provided") + } + } else { + stop("Tercile limits should have two dimensions") + } + } else { + stop("Tercile limits should be a vector of length two or an array of dimension (nfcsts,2)") + } + # check consistency of tercile limits + if (any(tercile.limits[, 1] >= tercile.limits[, 2])) { + stop("Inconsistent tercile limits") + } + #------------------------ + # Check extreme limits + #------------------------ + if (!is.null(extreme.limits)) { + if (is.vector(extreme.limits)) { + if (length(extreme.limits) != 2) { + stop("Provide two extreme limits") + } + extreme.limits <- InsertDim(extreme.limits, 1, npanels, name = "new") + } else if (is.array(extreme.limits)) { + if (length(dim(extreme.limits)) == 2) { + if (dim(extreme.limits)[2] != 2) { + stop("Provide two extreme limits for each panel") + } + if (dim(extreme.limits)[1] != npanels) { + stop("The number of extreme limits does not match the number of forecasts provided") + } + } else { + stop("extreme limits should have two dimensions") + } + } else { + stop("Extreme limits should be a vector of length two or an array of dimensions (nfcsts,2)") + } + # Check that extreme limits are consistent with tercile limits + if (any(extreme.limits[, 1] >= tercile.limits[, 1])) { + stop("Inconsistent lower extreme limits") + } + if (any(extreme.limits[, 2] <= tercile.limits[, 2])) { + stop("Inconsistent higher extreme limits") + } + } + #------------------------ + # Set proper fcst names + #------------------------ + if (!is.null(fcst.names)) { + if (length(fcst.names) != dim(fcst.df)[2]) { + stop("Parameter 'fcst.names' should be an array with as many names as distinct forecasts provided") + } + colnames(fcst.df) <- factor(fcst.names, levels = fcst.names) + } + #------------------------ + # Produce a first plot with the pdf for each init in a panel + #------------------------ + melt.df <- reshape2::melt(fcst.df, variable.name = "init", id.vars = NULL) + plot <- ggplot(melt.df, aes(x = value)) + + geom_density(alpha = 1, na.rm = T) + + coord_flip() + facet_wrap(~init, strip.position = "top", nrow = 1) + + xlim(range(c(obs, density(melt.df$value, na.rm = T)$x))) + ggp <- ggplot_build(plot) + #------------------------ + # Gather the coordinates of the plots together with init and corresponding + # terciles + #------------------------ + tmp.df <- ggp$data[[1]][, c("x", "ymin", "ymax", "PANEL")] + if (!is.null(ggp$layout$layout)) { + tmp.df$init <- ggp$layout$layout$init[as.numeric(tmp.df$PANEL)] + } else if (!is.null(ggp$layout$panel_layout)) { + tmp.df$init <- ggp$layout$panel_layout$init[as.numeric(tmp.df$PANEL)] + } else { + stop("Cannot find PANELS in ggp object") + } + tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL, 1], + "Below normal", ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL, 2], + "Normal", "Above normal")), + levels = c("Above normal", "Normal", "Below normal")) + #------------------------ + # Get the height and width of a panel + #------------------------ + pan.width <- diff(range(tmp.df$x)) + pan.height <- max(tmp.df$ymax) + magic.ratio <- 9 * pan.height/pan.width + #------------------------ + # Compute hatch coordinates for extremes + #------------------------ + if (!is.null(extreme.limits)) { + tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL, 1], + "Below P10", ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL, 2], + "Normal", "Above P90")), + levels = c("Above P90", "Normal", "Below P10")) + hatch.ls <- dlply(tmp.df, .(init, extremes), function(x) { + # close the polygon + tmp.df2 <- data.frame(x = c(x$x, max(x$x), min(x$x)), y = c(x$ymax, 0, + 0)) + # compute the hatches for this polygon + hatches <- .polygon.fullhatch(tmp.df2$x, tmp.df2$y, angle = 60, density = 10, + width.units = pan.width, height.units = pan.height) + # add bottom segment + end1 <- data.frame(x = x$x[1], y = x$ymax[1], xend = x$x[1], yend = 0) + # add top segment + end2 <- data.frame(x = x$x[length(x$x)], y = x$ymax[length(x$x)], xend = x$x[length(x$x)], + yend = 0) + return(rbind(hatches, end1, end2)) + }) + attr <- attr(hatch.ls, "split_labels") + for (i in 1:length(hatch.ls)) { + hatch.ls[[i]] <- cbind(hatch.ls[[i]], attr[i, ], row.names = NULL) + } + hatch.df <- do.call("rbind", hatch.ls) + # Compute max y for each extreme category + max.ls <- dlply(tmp.df, .(init, extremes), function(x) { + data.frame(y = min(0.85 * pan.height, max(x$ymax))) + }) + attr <- attr(max.ls, "split_labels") + for (i in 1:length(max.ls)) { + max.ls[[i]] <- cbind(max.ls[[i]], attr[i, ], row.names = NULL) + } + max.df <- do.call("rbind", max.ls) + } + #------------------------ + # Compute jitter space for ensemble members + #------------------------ + if (add.ensmemb != "no") { + jitter.df <- reshape2::melt(data.frame(dlply(melt.df, .(init), function(x) { + .jitter.ensmemb(sort(x$value, na.last = T), pan.width / 100) + }), check.names = F), value.name = "yjitter", variable.name = "init", id.vars = NULL) + jitter.df$x <- reshape2::melt(data.frame(dlply(melt.df, .(init), function(x) { + sort(x$value, na.last = T) + })), value.name = "x", id.vars = NULL)$x + } + #------------------------ + # Get y coordinates for observed x values, using a cool data.table feature: merge + # to nearest value + #------------------------ + if (!is.null(obs)) { + tmp.dt <- data.table(tmp.df, key = c("init", "x")) + obs.dt <- data.table(init = factor(colnames(fcst.df), levels = colnames(fcst.df)), + value = rep(obs, dim(fcst.df)[2])) + setkey(obs.dt, init, value) + obs.xy <- tmp.dt[obs.dt, roll = "nearest"] + } + #------------------------ + # Fill each pdf with different colors for the terciles + #------------------------ + plot <- plot + + geom_ribbon(data = tmp.df, + aes(x = x, ymin = ymin, ymax = ymax, fill = tercile), + alpha = 0.7) + #------------------------ + # Add hatches for extremes + #------------------------ + if (!is.null(extreme.limits)) { + if (nrow(hatch.df[hatch.df$extremes != "Normal", ]) == 0) { + warning("The provided extreme categories are outside the plot bounds. The extremes will not be drawn.") + extreme.limits <- NULL + } else { + plot <- plot + + geom_segment(data = hatch.df[hatch.df$extremes != "Normal", ], + aes(x = x, y = y, + xend = xend, yend = yend, color = extremes)) + } + } + #------------------------ + # Add obs line + #------------------------ + if (!is.null(obs)) { + plot <- plot + + geom_vline(data = obs.dt, + aes(xintercept = value), + linetype = "dashed", color = colorObs) + } + #------------------------ + # Add ensemble members + #------------------------ + if (add.ensmemb == "below") { + plot <- plot + + # this adds a grey box for ensmembers + geom_rect(aes(xmin = -Inf, xmax = Inf, + ymin = -Inf, ymax = -pan.height / 10), + fill = "gray95", color = "black", width = 0.2) + + # this adds the ensemble members + geom_point(data = jitter.df, + aes(x = x, + y = -pan.height / 10 - magic.ratio * yjitter, + shape = "Ensemble members"), + color = "black", fill = colorMember, alpha = 1) + + } else if (add.ensmemb == "above") { + plot <- plot + + geom_point(data = jitter.df, + aes(x = x, + y = 0.7 * magic.ratio * yjitter, + shape = "Ensemble members"), + color = "black", fill = colorMember, alpha = 1) + + } + #------------------------ + # Add obs diamond + #------------------------ + if (!is.null(obs)) { + plot <- plot + + # this adds the obs diamond + geom_point(data = obs.xy, + aes(x = x, y = ymax, size = "Observation"), + shape = 23, color = "black", fill = colorObs) + } + #------------------------ + # Compute probability for each tercile and identify MLT + #------------------------ + tmp.dt <- data.table(tmp.df) + pct <- tmp.dt[, .(pct = integrate(approxfun(x, ymax), lower = min(x), upper = max(x))$value), + by = .(init, tercile)] + # include potentially missing groups + pct <- merge(pct, CJ(init = factor(levels(pct$init), levels = levels(pct$init)), + tercile = factor(c("Below normal", "Normal", "Above normal"), + levels = c("Above normal", "Normal", "Below normal"))), + by = c("init", "tercile"), all.y = T) + pct[is.na(pct),"pct"] <- 0 + tot <- pct[, .(tot = sum(pct)), by = init] + pct <- merge(pct, tot, by = "init") + pct$pct <- round(100 * pct$pct / pct$tot, 0) + pct$MLT <- pct[, .(MLT = pct == max(pct)), by = init]$MLT + pct <- pct[order(init, tercile)] + pct$lab.pos <- as.vector(apply(tercile.limits, 1, function(x) {c(max(x), mean(x), min(x))})) + #------------------------ + # Compute probability for extremes + #------------------------ + if (!is.null(extreme.limits)) { + pct2 <- tmp.dt[, .(pct = integrate(approxfun(x, ymax), lower = min(x), upper = max(x))$value), + by = .(init, extremes)] + # include potentially missing groups + pct2 <- merge(pct2, CJ(init = factor(levels(pct2$init), levels = levels(pct2$init)), + extremes = factor(c("Below P10", "Normal", "Above P90"), + levels = c("Above P90", "Normal", "Below P10"))), + by = c("init", "extremes"), all.y=T) + pct2[is.na(pct),"pct"] <- 0 + tot2 <- pct2[, .(tot = sum(pct)), by = init] + pct2 <- merge(pct2, tot2, by = "init") + pct2$pct <- round(100 * pct2$pct / pct2$tot, 0) + pct2 <- pct2[order(init, extremes)] + pct2$lab.pos <- as.vector(apply(extreme.limits, 1, function(x) {c(x[2], NA, x[1])})) + pct2 <- merge(pct2, max.df, by = c("init", "extremes"), all.x = T) + } + #------------------------ + # Add probability labels for terciles + #------------------------ + if (add.ensmemb == "above") { + labpos <- -0.2 * pan.height + vjust <- 0 + } else { + labpos <- 0 + vjust <- -0.5 + } + plot <- plot + + geom_text(data = pct, + aes(x = lab.pos, y = labpos, label = paste0(pct, "%"), + hjust = as.integer(tercile) * -1.5 + 3.5), + vjust = vjust, angle = -90, size = 3.2) + + geom_text(data = pct[MLT == T, ], + aes(x = lab.pos, y = labpos, label = "*", + hjust = as.integer(tercile) * -3.5 + 9), + vjust = 0.1, angle = -90, size = 7, color = "black") + #------------------------ + # Add probability labels for extremes + #------------------------ + if (!is.null(extreme.limits)) { + plot <- plot + + geom_text(data = pct2[extremes != "Normal", ], + aes(x = lab.pos, y = 0.9 * y, label = paste0(pct, "%"), + hjust = as.integer(extremes) * -1.5 + 3.5), + vjust = -0.5, angle = -90, size = 3.2, + color = rep(colorLab, dim(fcst.df)[2])) + } + #------------------------ + # Finish all theme and legend details + #------------------------ + plot <- plot + + theme_minimal() + + scale_fill_manual(name = "Probability of\nterciles", + values = colorFill, drop = F) + + scale_color_manual(name = "Probability of\nextremes", + values = colorHatch) + + scale_shape_manual(name = "Ensemble\nmembers", + values = c(21)) + + scale_size_manual(name = "Observation", + values = c(3)) + + labs(x = var.name, + y = "Probability density\n(total area=1)", + title = title) + + theme(axis.text.x = element_blank(), + panel.grid.minor.x = element_blank(), + legend.key.size = unit(0.3, "in"), + panel.border = element_rect(fill = NA, color = "black"), + strip.background = element_rect(colour = "black", fill = "gray80"), + panel.spacing = unit(0.2, "in"), + panel.grid.major.x = element_line(color = "grey93"), + panel.background = element_rect(fill = "white"), + plot.background = element_rect(fill = "white", color = NA)) + + guides(fill = guide_legend(order = 1), + color = guide_legend(order = 2), + shape = guide_legend(order = 3, label = F), + size = guide_legend(order = 4, label = F)) + #------------------------ + # Save to plotfile if needed, and return plot + #------------------------ + if (!is.null(plotfile)) { + ggsave(plotfile, plot) + } + return(plot) +} +.jitter.ensmemb <- function(x, thr = 0.1) { + # Idea: start with first level. Loop all points, and if distance to last point in + # the level is more than a threshold, include the point to the level. Otherwise + # keep the point for another round. Do one round in each direction to avoid + # uggly patterns. + if (is.unsorted(x, na.rm = T)) { + stop("Provide a sorted array!") + } + lev <- x * 0 + lev[is.na(lev)] <- -1 + level <- 1 + while (any(lev == 0)) { + last <- -1/0 + for (i in 1:length(x)) { + if (lev[i] != 0) { + next + } + if (x[i] - last > thr) { + lev[i] <- level + last <- x[i] + } + } + level <- level + 1 + last <- 1/0 + for (i in seq(length(x), 1, -1)) { + if (lev[i] != 0) { + next + } + if (last - x[i] > thr) { + lev[i] <- level + last <- x[i] + } + } + level <- level + 1 + } + lev[lev == -1] <- NA + return(lev * thr * sqrt(3)/2) +} +.polygon.onehatch <- function(x, y, x0, y0, xd, yd, fillOddEven = F) { + halfplane <- as.integer(xd * (y - y0) - yd * (x - x0) <= 0) + cross <- halfplane[-1L] - halfplane[-length(halfplane)] + does.cross <- cross != 0 + if (!any(does.cross)) { + return() + } + x1 <- x[-length(x)][does.cross] + y1 <- y[-length(y)][does.cross] + x2 <- x[-1L][does.cross] + y2 <- y[-1L][does.cross] + t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1))/(xd * (y2 - y1) - yd * + (x2 - x1))) + o <- order(t) + tsort <- t[o] + crossings <- cumsum(cross[does.cross][o]) + if (fillOddEven) { + crossings <- crossings%%2 + } + drawline <- crossings != 0 + lx <- x0 + xd * tsort + ly <- y0 + yd * tsort + lx1 <- lx[-length(lx)][drawline] + ly1 <- ly[-length(ly)][drawline] + lx2 <- lx[-1L][drawline] + ly2 <- ly[-1L][drawline] + return(data.frame(x = lx1, y = ly1, xend = lx2, yend = ly2)) +} +.polygon.fullhatch <- function(x, y, density, angle, width.units, height.units, inches = c(5, + 1)) { + x <- c(x, x[1L]) + y <- c(y, y[1L]) + angle <- angle%%180 + upi <- c(width.units, height.units)/inches + if (upi[1L] < 0) { + angle <- 180 - angle + } + if (upi[2L] < 0) { + angle <- 180 - angle + } + upi <- abs(upi) + xd <- cos(angle/180 * pi) * upi[1L] + yd <- sin(angle/180 * pi) * upi[2L] + hatch.ls <- list() + i <- 1 + if (angle < 45 || angle > 135) { + if (angle < 45) { + first.x <- max(x) + last.x <- min(x) + } else { + first.x <- min(x) + last.x <- max(x) + } + y.shift <- upi[2L]/density/abs(cos(angle/180 * pi)) + x0 <- 0 + y0 <- floor((min(y) - first.x * yd/xd)/y.shift) * y.shift + y.end <- max(y) - last.x * yd/xd + while (y0 < y.end) { + hatch.ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) + i <- i + 1 + y0 <- y0 + y.shift + } + } else { + if (angle < 90) { + first.y <- max(y) + last.y <- min(y) + } else { + first.y <- min(y) + last.y <- max(y) + } + x.shift <- upi[1L]/density/abs(sin(angle/180 * pi)) + x0 <- floor((min(x) - first.y * xd/yd)/x.shift) * x.shift + y0 <- 0 + x.end <- max(x) - last.y * xd/yd + while (x0 < x.end) { + hatch.ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) + i <- i + 1 + x0 <- x0 + x.shift + } + } + return(do.call("rbind", hatch.ls)) +} + diff --git a/R/PlotMostLikelyQuantileMap.R b/R/PlotMostLikelyQuantileMap.R new file mode 100644 index 0000000000000000000000000000000000000000..a43280591ff64749a0176b17652364d0e2124211 --- /dev/null +++ b/R/PlotMostLikelyQuantileMap.R @@ -0,0 +1,222 @@ +#'Plot Maps of Most Likely Quantiles +#' +#'@author Veronica Torralba, \email{veronica.torralba@bsc.es}, Nicolau Manubens, +#'\email{nicolau.manubens@bsc.es} +#'@description This function receives as main input (via the parameter +#'\code{probs}) a collection of longitude-latitude maps, each containing the +#'probabilities (from 0 to 1) of the different grid cells of belonging to a +#'category. As many categories as maps provided as inputs are understood to +#'exist. The maps of probabilities must be provided on a common rectangular +#'regular grid, and a vector with the longitudes and a vector with the latitudes +#'of the grid must be provided. The input maps can be provided in two forms, +#'either as a list of multiple two-dimensional arrays (one for each category) or +#'as a three-dimensional array, where one of the dimensions corresponds to the +#'different categories. +#' +#'@param probs A list of bi-dimensional arrays with the named dimensions +#' 'latitude' (or 'lat') and 'longitude' (or 'lon'), with equal size and in the +#' same order, or a single tri-dimensional array with an additional dimension +#' (e.g. 'bin') for the different categories. The arrays must contain +#' probability values between 0 and 1, and the probabilities for all categories +#' of a grid cell should not exceed 1 when added. +#'@param lon A numeric vector with the longitudes of the map grid, in the same +#' order as the values along the corresponding dimension in \code{probs}. +#'@param lat A numeric vector with the latitudes of the map grid, in the same +#' order as the values along the corresponding dimension in \code{probs}. +#'@param cat_dim The name of the dimension along which the different categories +#' are stored in \code{probs}. This only applies if \code{probs} is provided in +#' the form of 3-dimensional array. The default expected name is 'bin'. +#'@param bar_titles Vector of character strings with the names to be drawn on +#' top of the color bar for each of the categories. As many titles as +#' categories provided in \code{probs} must be provided. +#'@param col_unknown_cat Character string with a colour representation of the +#' colour to be used to paint the cells for which no category can be clearly +#' assigned. Takes the value 'white' by default. +#'@param drawleg Where to draw the common colour bar. Can take values TRUE, +#' FALSE or:\cr +#' 'up', 'u', 'U', 'top', 't', 'T', 'north', 'n', 'N'\cr +#' 'down', 'd', 'D', 'bottom', 'b', 'B', 'south', 's', 'S' (default)\cr +#' 'right', 'r', 'R', 'east', 'e', 'E'\cr +#' 'left', 'l', 'L', 'west', 'w', 'W' +#'@param ... Additional parameters to be sent to \code{PlotCombinedMap} and +#' \code{PlotEquiMap}. +#'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap} +#'@examples +#'# Simple example +#'x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200 +#'a <- x * 0.6 +#'b <- (1 - x) * 0.6 +#'c <- 1 - (a + b) +#'lons <- seq(0, 359.5, length = 20) +#'lats <- seq(-89.5, 89.5, length = 10) +#'\dontrun{ +#'PlotMostLikelyQuantileMap(list(a, b, c), lons, lats, +#' toptitle = 'Most likely tercile map', +#' bar_titles = paste('% of belonging to', c('a', 'b', 'c')), +#' brks = 20, width = 10, height = 8) +#'} +#' +#'# 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) +#' +#'# 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 +#' +#'# 2. Binning sample data +#'prob_thresholds <- 1:n_bins / n_bins +#'prob_thresholds <- prob_thresholds[1:(n_bins - 1)] +#'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 +#' for (i in 1:n_bins) { +#' result[i] <- sum(x > lower_threshold & x <= thresholds[i]) / n_samples +#' lower_threshold <- thresholds[i] +#' } +#' +#' dim(result) <- c(bin = n_bins) +#' result +#'} +#' +#'bins <- multiApply::Apply(sample_data, 'time', binning, thresholds)$output1 +#' +#'# 3. Plotting most likely quantile/bin +#'\dontrun{ +#'PlotMostLikelyQuantileMap(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) +#'} +#'@importFrom maps map +#'@importFrom graphics box image layout mtext par plot.new +#'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off hcl jpeg pdf png postscript svg tiff +#'@export +PlotMostLikelyQuantileMap <- function(probs, lon, lat, cat_dim = 'bin', + bar_titles = NULL, + col_unknown_cat = 'white', drawleg = T, + ...) { + # Check probs + error <- FALSE + if (is.list(probs)) { + if (length(probs) < 1) { + stop("Parameter 'probs' must be of length >= 1 if provided as a list.") + } + check_fun <- function(x) { + is.numeric(x) && (length(dim(x)) == 2) + } + if (!all(sapply(probs, check_fun))) { + error <- TRUE + } + ref_dims <- dim(probs[[1]]) + equal_dims <- all(sapply(probs, function(x) identical(dim(x), ref_dims))) + if (!equal_dims) { + stop("All arrays in parameter 'probs' must have the same dimension ", + "sizes and names when 'probs' is provided as a list of arrays.") + } + num_probs <- length(probs) + probs <- unlist(probs) + dim(probs) <- c(ref_dims, map = num_probs) + cat_dim <- 'map' + } + if (!is.numeric(probs)) { + error <- TRUE + } + if (is.null(dim(probs))) { + error <- TRUE + } + if (length(dim(probs)) != 3) { + error <- TRUE + } + if (error) { + stop("Parameter 'probs' must be either a numeric array with 3 dimensions ", + " or a list of numeric arrays of the same size with the 'lon' and ", + "'lat' dimensions.") + } + dimnames <- names(dim(probs)) + + # Check cat_dim + if (is.character(cat_dim)) { + if (is.null(dimnames)) { + stop("Specified a dimension name in 'cat_dim' but no dimension names provided ", + "in 'probs'.") + } + cat_dim <- which(dimnames == cat_dim) + if (length(cat_dim) < 1) { + stop("Dimension 'cat_dim' not found in 'probs'.") + } + cat_dim <- cat_dim[1] + } else if (!is.numeric(cat_dim)) { + stop("Parameter 'cat_dim' must be either a numeric value or a ", + "dimension name.") + } + if (length(cat_dim) != 1) { + stop("Parameter 'cat_dim' must be of length 1.") + } + cat_dim <- round(cat_dim) + nprobs <- dim(probs)[cat_dim] + + # Check bar_titles + if (is.null(bar_titles)) { + if (nprobs == 3) { + bar_titles <- c("Below normal (%)", "Normal (%)", "Above normal (%)") + } else if (nprobs == 5) { + bar_titles <- c("Low (%)", "Below normal (%)", + "Normal (%)", "Above normal (%)", "High (%)") + } else { + bar_titles <- paste0("Cat. ", 1:nprobs, " (%)") + } + } + + minimum_value <- ceiling(1 / nprobs * 10 * 1.1) * 10 + + # By now, the PlotCombinedMap function is included below in this file. + # In the future, PlotCombinedMap will be part of s2dverification and will + # be properly imported. + PlotCombinedMap(probs * 100, lon, lat, map_select_fun = max, + display_range = c(minimum_value, 100), + map_dim = cat_dim, + bar_titles = bar_titles, + col_unknown_map = col_unknown_cat, + drawleg = drawleg, ...) +} \ No newline at end of file diff --git a/R/PlotPDFsOLE.R b/R/PlotPDFsOLE.R new file mode 100644 index 0000000000000000000000000000000000000000..bf95abb76c745410447d0cdf59c22f34b2509231 --- /dev/null +++ b/R/PlotPDFsOLE.R @@ -0,0 +1,259 @@ +#'Plotting two probability density gaussian functions and the optimal linear +#'estimation (OLE) as result of combining them. +#' +#'@author Eroteida Sanchez-Garcia - AEMET, //email{esanchezg@aemet.es} +#' +#'@description This function plots two probability density gaussian functions +#'and the optimal linear estimation (OLE) as result of combining them. +#' +#'@param pdf_1 A numeric array with a dimension named 'statistic', containg +#' two parameters: mean' and 'standard deviation' of the first gaussian pdf +#' to combining. +#'@param pdf_2 A numeric array with a dimension named 'statistic', containg +#' two parameters: mean' and 'standard deviation' of the second gaussian pdf +#' to combining. +#'@param nsigma (optional) A numeric value for setting the limits of X axis. +#' (Default nsigma = 3). +#'@param legendPos (optional) A character value for setting the position of the +#' legend ("bottom", "top", "right" or "left")(Default 'bottom'). +#'@param legendSize (optional) A numeric value for setting the size of the +#' legend text. (Default 1.0). +#'@param plotfile (optional) A filename where the plot will be saved. +#' (Default: the plot is not saved). +#'@param width (optional) A numeric value indicating the plot width in +#' units ("in", "cm", or "mm"). (Default width = 30). +#'@param height (optional) A numeric value indicating the plot height. +#' (Default height = 15). +#'@param units (optional) A character value indicating the plot size +#' unit. (Default units = 'cm'). +#'@param dpi (optional) A numeric value indicating the plot resolution. +#' (Default dpi = 300). +#' +#'@return PlotPDFsOLE() returns a ggplot object containing the plot. +#' +#'@examples +#'# Example 1 +#'pdf_1 <- c(1.1,0.6) +#'attr(pdf_1, "name") <- "NAO1" +#'dim(pdf_1) <- c(statistic = 2) +#'pdf_2 <- c(1,0.5) +#'attr(pdf_2, "name") <- "NAO2" +#'dim(pdf_2) <- c(statistic = 2) +#' +#'PlotPDFsOLE(pdf_1, pdf_2) +#'@import ggplot2 +#'@export +PlotPDFsOLE <- function(pdf_1, pdf_2, nsigma = 3, legendPos = 'bottom', + legendSize = 1.0, plotfile = NULL, width = 30, + height = 15, units = "cm", dpi = 300) { + y <- type <- NULL + + if(!is.null(plotfile)){ + if (!is.numeric(dpi)) { + stop("Parameter 'dpi' must be numeric.") + } + if (length(dpi) > 1) { + warning("Parameter 'dpi' has length greater than 1 and ", + "only the first element will be used.") + dpi <- dpi[1] + } + if (!is.character(units)) { + stop("Parameter 'units' must be character") + } + if (length(units) > 1) { + warning("Parameter 'units' has length greater than 1 and ", + "only the first element will be used.") + units <- units[1] + } + if(!(units %in% c("in", "cm", "mm"))) { + stop("Parameter 'units' must be equal to 'in', 'cm' or 'mm'.") + } + if (!is.numeric(height)) { + stop("Parameter 'height' must be numeric.") + } + if (length(height) > 1) { + warning("Parameter 'height' has length greater than 1 and ", + "only the first element will be used.") + height <- height[1] + } + if (!is.numeric(width)) { + stop("Parameter 'width' must be numeric.") + } + if (length(width) > 1) { + warning("Parameter 'width' has length greater than 1 and ", + "only the first element will be used.") + width <- width[1] + } + if (!is.character(plotfile)) { + stop("Parameter 'plotfile' must be a character string ", + "indicating the path and name of output png file.") + } + } + if (!is.character(legendPos)) { + stop("Parameter 'legendPos' must be character") + } + if(!(legendPos %in% c("bottom", "top", "right", "left"))) { + stop("Parameter 'legendPos' must be equal to 'bottom', 'top', 'right' or 'left'.") + } + if (!is.numeric(legendSize)) { + stop("Parameter 'legendSize' must be numeric.") + } + if (!is.numeric(nsigma)) { + stop("Parameter 'nsigma' must be numeric.") + } + if (length(nsigma) > 1) { + warning("Parameter 'nsigma' has length greater than 1 and ", + "only the first element will be used.") + nsigma <- nsigma[1] + } + if (!is.array(pdf_1)) { + stop("Parameter 'pdf_1' must be an array.") + } + if (!is.array(pdf_2)) { + stop("Parameter 'pdf_2' must be an array.") + } + if (!is.numeric(pdf_1)) { + stop("Parameter 'pdf_1' must be a numeric array.") + } + if (!is.numeric(pdf_2)) { + stop("Parameter 'pdf_2' must be a numeric array.") + } + if (is.null(names(dim(pdf_1))) || + is.null(names(dim(pdf_2)))) { + stop("Parameters 'pdf_1' and 'pdf_2' ", + "should have dimmension names.") + } + if(!('statistic' %in% names(dim(pdf_1)))) { + stop("Parameter 'pdf_1' must have dimension 'statistic'.") + } + if(!('statistic' %in% names(dim(pdf_2)))) { + stop("Parameter 'pdf_2' must have dimension 'statistic'.") + } + if (length(dim(pdf_1)) != 1) { + stop("Parameter 'pdf_1' must have only dimension 'statistic'.") + } + if (length(dim(pdf_2)) != 1) { + stop("Parameter 'pdf_2' must have only dimension 'statistic'.") + } + if ((dim(pdf_1)['statistic'] != 2) || (dim(pdf_2)['statistic'] != 2)) { + stop("Length of dimension 'statistic'", + "of parameter 'pdf_1' and 'pdf_2' must be equal to 2.") + } + if(!is.null(attr(pdf_1, "name"))){ + if(!is.character(attr(pdf_1, "name"))){ + stop("The 'name' attribute of parameter 'pdf_1' must be a character ", + "indicating the name of the variable of parameter 'pdf_1'.") + } + } + if(!is.null(attr(pdf_2, "name"))){ + if(!is.character(attr(pdf_2, "name"))){ + stop("The 'name' attribute of parameter 'pdf_2' must be a character ", + "indicating the name of the variable of parameter 'pdf_2'.") + } + } + if(is.null(attr(pdf_1, "name"))){ + name1 <- "variable 1" + } else { + name1 <- attr(pdf_1, "name") + } + if(is.null(attr(pdf_2, "name"))){ + name2 <- "Variable 2" + } else { + name2 <- attr(pdf_2, "name") + } + + #----------------------------------------------------------------------------- + # Set parameters of gaussian distributions (mean and sd) + #----------------------------------------------------------------------------- + mean1 <- pdf_1[1] + sigma1 <- pdf_1[2] + mean2 <- pdf_2[1] + sigma2 <- pdf_2[2] + pdfBest <- CombinedPDFs(pdf_1, pdf_2) + meanBest <- pdfBest[1] + sigmaBest <- pdfBest[2] + + + #----------------------------------------------------------------------------- + # Plot the gaussian distributions + #----------------------------------------------------------------------------- + nameBest <- paste0(name1, " + ", name2) + graphicTitle <- "OPTIMAL LINEAR ESTIMATION" + xlimSup <- max(nsigma * sigmaBest + meanBest, nsigma * sigma1 + mean1, + nsigma * sigma2 + mean2) + xlimInf <- min(-nsigma * sigmaBest+meanBest, - nsigma * sigma1 + mean1, + -nsigma * sigma2 + mean2) + # deltax <- 0.02 + deltax <- (xlimSup - xlimInf) / 10000 + + x <- seq(xlimInf, xlimSup, deltax) + df1 <- data.frame(x = x, y = dnorm(x, mean = mean1, sd = sigma1), + type = name1) + df2 <- data.frame(x = x, y = dnorm(x, mean = mean2, sd = sigma2), + type = name2) + df3 <- data.frame(x = x, y = dnorm(x, mean = meanBest, sd = sigmaBest), + type = nameBest) + df123 <- rbind(df1, df2, df3) + label1 <- paste0(name1, ": N(mean=",round(mean1, 2), ", sd=", round(sigma1, 2), + ")") + label2 <- paste0(name2, ": N(mean=",round(mean2, 2), ", sd=", round(sigma2, 2), + ")") + labelBest <- paste0(nameBest, ": N(mean=",round(meanBest,2), ", sd=", + round(sigmaBest, 2), ")") + cols <- c("#DC3912", "#13721A", "#1F5094") + names(cols) <- c(name1, name2, nameBest) + g <- ggplot(df123) + geom_line(aes(x, y, colour = type), size = rel(1.2)) + + g <- g + scale_colour_manual(values = cols, + limits = c(name1, name2, nameBest), + labels = c(label1, label2, labelBest)) + + g <- g + theme(plot.title=element_text(size=rel(1.1), colour="black", + face= "bold"), + axis.text.x = element_text(size=rel(1.2)), + axis.text.y = element_text(size=rel(1.2)), + axis.title.x = element_blank(), + legend.title = element_blank(), + legend.position = legendPos, + legend.text = element_text(face = "bold", size=rel(legendSize))) + + g <- g + ggtitle(graphicTitle) + g <- g + labs(y="probability", size=rel(1.9)) + g <- g + stat_function(fun = dnorm_limit, args = list(mean=mean1, sd=sigma1), + fill = cols[name1], alpha=0.2, geom="area") + g <- g + stat_function(fun = dnorm_limit, args = list(mean=mean2, sd=sigma2), + fill = cols[name2], alpha=0.2, geom="area") + g <- g + stat_function(fun = dnorm_limit, args = list(mean=meanBest, + sd=sigmaBest), + fill = cols[nameBest], alpha=0.2, geom="area") + + + #----------------------------------------------------------------------------- + # Save to plotfile if needed, and return plot + #----------------------------------------------------------------------------- + if (!is.null(plotfile)) { + ggsave(plotfile, g, width = width, height = height, units = units, dpi = dpi) + } + return(g) +} + +# Auxiliar function to plot +CombinedPDFs <- function(pdf_1, pdf_2) { + mean_1 <- pdf_1[1] + sigma_1 <- pdf_1[2] + mean_2 <- pdf_2[1] + sigma_2 <- pdf_2[2] + a_1 <- (sigma_2^2)/((sigma_1^2)+(sigma_2^2)) + a_2 <- (sigma_1^2)/((sigma_1^2)+(sigma_2^2)) + pdf_mean <- a_1*mean_1 + a_2*mean_2 + pdf_sigma <- sqrt((sigma_1^2)*(sigma_2^2)/((sigma_1^2)+(sigma_2^2))) + data <- c(pdf_mean, pdf_sigma) + dim(data) <- c(statistic = 2) + return(data) +} + +dnorm_limit <- function(x,mean,sd){ + y <- dnorm(x,mean,sd) + y[x mean+sd] <- NA + return(y) +} diff --git a/R/PlotTriangles4Categories.R b/R/PlotTriangles4Categories.R new file mode 100644 index 0000000000000000000000000000000000000000..501536c72df82b181efe858c15fe0a1dcdc9d5d7 --- /dev/null +++ b/R/PlotTriangles4Categories.R @@ -0,0 +1,296 @@ +#'Function to convert any 3-d numerical array to a grid of coloured triangles. +#' +#'This function converts a 3-d numerical data array into a coloured +#'grid with triangles. It is useful for a slide or article to present tabular +#'results as colors instead of numbers. This can be used to compare the outputs +#'of two or four categories (e.g. modes of variability, clusters, or forecast +#'systems). +#' +#'@param data Array with three named dimensions: 'dimx', 'dimy', 'dimcat', +#' containing the values to be displayed in a coloured image with triangles. +#'@param brks A vector of the color bar intervals. The length must be one more +#' than the parameter 'cols'. Use ColorBar() to generate default values. +#'@param cols A vector of valid colour identifiers for color bar. The length +#' must be one less than the parameter 'brks'. Use ColorBar() to generate +#' default values. +#'@param toptitle A string of the title of the grid. Set NULL as default. +#'@param sig_data Logical array with the same dimensions as 'data' to add layers +#' to the plot. A value of TRUE at a grid cell will draw a dot/symbol on the +#' corresponding triangle of the plot. Set NULL as default. +#'@param pch_sig Symbol to be used to represent sig_data. Takes 18 +#' (diamond) by default. See 'pch' in par() for additional accepted options. +#'@param col_sig Colour of the symbol to represent sig_data. +#'@param cex_sig Parameter to increase/reduce the size of the symbols used +#' to represent sig_data. +#'@param xlab A logical value (TRUE) indicating if xlabels should be plotted +#'@param ylab A logical value (TRUE) indicating if ylabels should be plotted +#'@param xlabels A vector of labels of the x-axis The length must be +#' length of the col of parameter 'data'. Set the sequence from 1 to the +#' length of the row of parameter 'data' as default. +#'@param xtitle A string of title of the x-axis. Set NULL as default. +#'@param ylabels A vector of labels of the y-axis The length must be +#' length of the row of parameter 'data'. Set the sequence from 1 to the +#' length of the row of parameter 'data' as default. +#'@param ytitle A string of title of the y-axis. Set NULL as default. +#'@param legend A logical value to decide to draw the color bar legend or not. +#' Set TRUE as default. +#'@param lab_legend A vector of labels indicating what is represented in each +#'category (i.e. triangle). Set the sequence from 1 to the length of +#' the categories (2 or 4). +#'@param cex_leg A number to indicate the increase/reductuion of the lab_legend +#' used to represent sig_data. +#'@param col_leg Color of the legend (triangles). +#'@param cex_axis A number to indicate the increase/reduction of the axis labels. +#'@param fileout A string of full directory path and file name indicating where +#' to save the plot. If not specified (default), a graphics device will pop up. +#'@param mar A numerical vector of the form c(bottom, left, top, right) which +#' gives the number of lines of margin to be specified on the four sides of the +#' plot. +#'@param size_units A string indicating the units of the size of the device +#' (file or window) to plot in. Set 'px' as default. See ?Devices and the +#' creator function of the corresponding device. +#'@param res A positive number indicating resolution of the device (file or +#' window) to plot in. See ?Devices and the creator function of the +#' corresponding device. +#'@param figure.width a numeric value to control the width of the plot. +#'@param ... The additional parameters to be passed to function +#' ColorBarContinuous() in for color legend creation. +#'@return A figure in popup window by default, or saved to the specified path. +#' +#'@author History:\cr +#'1.0 - 2020-10 (V.Torralba, \email{veronica.torralba@bsc.es}) - Original code +#' +#'@examples +#'# Example with random data +#'arr1 <- array(runif(n = 4 * 5 * 4, min = -1, max = 1), dim = c(4,5,4)) +#'names(dim(arr1)) <- c('dimx', 'dimy', 'dimcat') +#'arr2 <- array(TRUE, dim = dim(arr1)) +#'arr2[which(arr1 < 0.3)] <- FALSE +#'PlotTriangles4Categories(data = arr1, +#' cols = c('white','#fef0d9','#fdd49e','#fdbb84','#fc8d59'), +#' brks = c(-1, 0, 0.1, 0.2, 0.3, 0.4), +#' lab_legend = c('NAO+', 'BL','AR','NAO-'), +#' xtitle = "Target month", ytitle = "Lead time", +#' xlabels = c("Jan", "Feb", "Mar", "Apr")) +#'@importFrom grDevices dev.new dev.off dev.cur +#'@importFrom graphics plot points polygon text title axis +#'@importFrom RColorBrewer brewer.pal +#'@importFrom ClimProjDiags Subset +#'@export +PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, + toptitle = NULL, sig_data = NULL, + pch_sig = 18, col_sig = 'black', + cex_sig = 1, xlab = TRUE, ylab = TRUE, + xlabels = NULL, xtitle = NULL, + ylabels = NULL, ytitle = NULL, + legend = TRUE, lab_legend = NULL, + cex_leg = 1, col_leg = 'black', + cex_axis = 1.5, mar = c(5, 4, 0, 0), + fileout = NULL, size_units = 'px', + res = 100, figure.width = 1, ...) { + # Checking the dimensions + if (length(dim(data)) != 3) { + stop("Parameter 'data' must be an array with three dimensions.") + } + + if (any(is.na(data))){ + stop("Parameter 'data' cannot contain NAs.") + } + + if (is.null(names(dim(data)))) { + stop("Parameter 'data' must be an array with named dimensions.") + }else{ + if (!any(names(dim(data)) == 'dimx') | !any(names(dim(data)) == 'dimy') | + !any(names(dim(data)) == 'dimcat')) { + stop("Parameter 'data' should contain 'dimx', 'dimy' and 'dimcat' dimension names. ") + } + } + if (!is.vector(mar) & length(mar) != 4) { + stop("Parameter 'mar' must be a vector of length 4.") + } + if (!is.null(sig_data)) { + if (!is.logical(sig_data)) { + stop("Parameter 'sig_data' array must be logical.")} + else if (length(dim(sig_data)) != 3) { + stop("Parameter 'sig_data' must be an array with three dimensions.") + }else if (any(dim(sig_data) != dim(data))){ + stop("Parameter 'sig_data' must be an array with the same dimensions as 'data'.") + }else if(!is.null(names(dim(sig_data)))) { + if (any(names(dim(sig_data)) != names(dim(data)))) { + stop("Parameter 'sig_data' must be an array with the same named dimensions as 'data'.")} + } + } + + if (dim(data)['dimcat'] != 4 && dim(data)['dimcat'] != 2) { + stop( + "Parameter 'data' should contain a dimcat dimension with length equals + to two or four as only two or four categories can be plotted.") + } + + # Checking what is available and generating missing information + if (!is.null(lab_legend) && + length(lab_legend) != 4 && length(lab_legend) != 2) { + stop("Parameter 'lab_legend' should contain two or four names.") + } + + datadim <- dim(data) + nrow <- dim(data)['dimy'] + ncol <- dim(data)['dimx'] + ncat <- dim(data)['dimcat'] + + # 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 = 80 * ncol * figure.width, + height = 80 * nrow, + units = size_units, res = res) + saveToFile <- deviceInfo$fun + fileout <- deviceInfo$files + } + + # 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 = 8 * figure.width, height = 5) + } + + if (is.null(xlabels)){ + xlabels = 1:ncol + } + if (is.null(ylabels)){ + ylabels = 1:nrow + } + + if (!is.null(brks) && !is.null(cols)) { + if (length(brks) != length(cols) + 1) { + stop("The length of the parameter 'brks' must be one more than 'cols'.") + } + } + if (is.null(brks)) { + brks <- seq(min(data, na.rm = T), max(data, na.rm = T), length.out = 9) + } + if (is.null(cols)) { + cols <- rev(brewer.pal(length(brks) - 1, 'RdBu')) + } + + # The colours for each triangle/category are defined + data_cat <- array(cols[length(cols)], dim = datadim) + names(dim(data_cat)) <- names(dim(data)) + for (i in (length(cols) - 1):1) { + data_cat[data < brks[i + 1]] <- cols[i] + } + + if(legend){ + layout(matrix(c(1, 2, 1, 3), 2, 2, byrow = T), + widths = c(10, 3.4), heights = c(10, 3.5)) + par(oma = c(1, 1, 1, 1), mar = mar) + if(is.null(lab_legend)) { + lab_legend = 1:ncat + } + } + + plot(ncol, nrow, xlim = c(0, ncol), ylim=c(0, nrow), xaxs = "i", yaxs = 'i', type = "n", + xaxt = "n", yaxt = "n", ann = F, axes = F) + + box(col = 'black',lwd = 1) + + if (! is.null(toptitle)){ + title(toptitle, cex = 1.5) + } + + if (!is.null(xtitle)){ + mtext(side = 1, text = xtitle, line = 4, cex = 1.5) + } + if (!is.null(ytitle)){ + mtext(side = 2, text = ytitle, line = 2.5, cex = 1.5) + } + + if (xlab){ + axis(1, at =(1:ncol) - 0.5, las = 2, labels = xlabels, cex.axis = cex_axis) + } + if (ylab){ + axis(2, at = (1:nrow) - 0.5, las = 2, labels = ylabels, cex.axis = cex_axis) + } + + + #The triangles are plotted + for(p in 1:ncol){ + for(l in 1:nrow){ + if (ncat == 4){ + coord_triangl <- list(xs=list(c(p-1, p-0.5, p-1),c(p-1, p-0.5, p),c(p, p-0.5, p),c(p-1, p-0.5, p)), + ys=list( c(l-1, -0.5+l, l), c(l-1, -0.5+l, l-1),c(l-1, -0.5+l, l),c(l, -0.5+l, l))) + + coord_sig <- list(x=c(p-0.75,p-0.5,p-0.25,p-0.5),y=c(l-0.5,l-0.75,l-0.5,l-0.25)) + } + + if (ncat==2){ + coord_triangl <- list(xs=list(c(p-1, p, p-1),c(p-1, p, p)), + ys=list(c(l-1, l, l),c(l-1,l-1, l))) + coord_sig <- list(x=c(p-(2/3),p-(1/3)),y=c(l-(1/3),l-(2/3))) + } + for (n in 1:ncat) { + polygon(coord_triangl$xs[[n]], + coord_triangl$ys[[n]], + col = Subset( + data_cat, + along = c('dimcat', 'dimx', 'dimy'), + indices = list(n, p, l))) + if (!is.null(sig_data) && + Subset(sig_data,along = c('dimcat', 'dimx', 'dimy'), + indices = list(n, p, l))) { + points( + x = coord_sig$x[n], + y = coord_sig$y[n], + pch = pch_sig, + cex = cex_sig, + col = col_sig + ) + } + } + } + } + + # legend + + if(legend){ + # Colorbar + par(mar=c(0,0,0,0)) + ColorBarContinuous(brks = brks, cols = cols, vertical = T, draw_ticks = T, draw_separators = T, + # extra_margin = c(0,0,2.5,0),label_scale = 1.5,...) + extra_margin = c( 0, 0, 0, 0), label_scale = 1.5, ...) + + par(mar = c(0.5, 2.5, 0.5, 2.5)) + plot(1, 1, xlim = c(0, 1), ylim =c(0, 1), xaxs = "i", yaxs = 'i', type = "n", + xaxt = "n", yaxt = "n", ann = F, axes = F) + + box(col = col_leg) + p = l = 1 + if (ncat == 4){ + coord_triangl <- list(xs = list(c(p -1, p - 0.5, p - 1), c(p - 1, p - 0.5, p), + c(p, p - 0.5, p), c(p - 1, p - 0.5, p)), + ys = list(c(l - 1, - 0.5 + l, l), c(l - 1, - 0.5 + l, l - 1), + c(l - 1, - 0.5 + l, l), c(l, - 0.5 + l, l))) + + coord_sig <- list(x = c(p - 0.75, p - 0.5, p - 0.25, p - 0.5), + y = c(l - 0.5, l - 0.75, l - 0.5, l - 0.25)) + } + + if (ncat==2){ + coord_triangl<- list(xs=list(c(p-1, p, p),c(p-1, p, p-1)), + ys=list( c(l-1,l-1, l), c(l-1, l, l))) + coord_sig<- list(x=c(p-(2/3),p-(1/3)),y=c(l-(1/3),l-(2/3))) + } + for (n in 1:ncat) { + polygon(coord_triangl$xs[[n]], + coord_triangl$ys[[n]],border=col_leg) + text(x=coord_sig$x[[n]],y=coord_sig$y[[n]],labels = lab_legend[n],cex=cex_leg,col=col_leg) + + } + } + + # If the graphic was saved to file, close the connection with the device + if (!is.null(fileout)) dev.off() +} diff --git a/R/PlotWeeklyClim.R b/R/PlotWeeklyClim.R new file mode 100644 index 0000000000000000000000000000000000000000..4fa829bce06477006e25134eba2224b32267afa9 --- /dev/null +++ b/R/PlotWeeklyClim.R @@ -0,0 +1,315 @@ +#'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. +#' +#'@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 ylim A numeric vector of length two providing limits of the scale. +#' Use NA to refer to the existing minimum or maximum. For more information, +#' see 'ggplot2' documentation of 'scale_y_continuous' parameter. +#'@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 ggplot2 +#'@import RColorBrewer +#'@import scales +#'@importFrom ClimProjDiags Subset +#'@importFrom s2dv MeanDims +#'@importFrom CSTools SplitDim +#'@importFrom lubridate ymd +#'@export +PlotWeeklyClim <- function(data, first_date, ref_period, last_date = NULL, + data_years = NULL, time_dim = 'time', + sdate_dim = 'sdate', ylim = NULL, + 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") + if (!all(col_names %in% names(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 <- s2dv::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 <- lubridate::ymd(first_date) + target_year <- lubridate::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 <- lubridate::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") + } + # ylim + if (is.character(ylim)) { + warning("Parameter 'ylim' can't be a character string, it will not be used.") + ylim <- NULL + } + + 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 <- ClimProjDiags::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 <- ClimProjDiags::Subset(data_subset, dims_subset, + as.list(rep(1, length(dims_subset))), drop = TRUE) + } + # observed daily data creation + daily <- ClimProjDiags::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 <- ClimProjDiags::Subset(data_subset, along = sdate_dim, + indices = indexes_reference_period) + } + + ## Weekly aggregations for reference period + weekly_aggre <- CSTools::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 <- s2dv::MeanDims(weekly_aggre, time_dim) + weekly_clim <- s2dv::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)) + + ## Prepare observations from 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])) + + 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)) + + scale_y_continuous(limits = ylim) + + # 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) + } +} \ No newline at end of file diff --git a/R/ShapeToMask.R b/R/ShapeToMask.R index 8c91675e6e436bc77ff7c5eac3959ea037edbe3b..7bf8de45efdcf6bea8db58f698b77a3de36531b3 100644 --- a/R/ShapeToMask.R +++ b/R/ShapeToMask.R @@ -185,8 +185,8 @@ ShapeToMask <- function(shp_file, ref_grid = NULL, ## Method 1: ref_grid is a netCDF file if (is.null(lat_dim) | is.null(lon_dim)) { var_names <- easyNCDF::NcReadVarNames(ref_grid) - lat_dim <- var_names[which(var_names %in% s2dv:::.KnownLatNames())] - lon_dim <- var_names[which(var_names %in% s2dv:::.KnownLonNames())] + lat_dim <- var_names[which(var_names %in% .KnownLatNames())] + lon_dim <- var_names[which(var_names %in% .KnownLonNames())] } latlon <- NcToArray(ref_grid, vars_to_read = c(lat_dim, lon_dim)) lat <- NcToArray(ref_grid, vars_to_read = lat_dim)[1, ] @@ -198,9 +198,8 @@ ShapeToMask <- function(shp_file, ref_grid = NULL, stop("If 'ref_grid' is a list, it must have two items for longitude and latitude.") } if (is.null(lat_dim) | is.null(lon_dim)) { - # NOTE: the names come from s2dv:::.KnownLonNames and .KnownLatNames - lon_known_names <- c(s2dv:::.KnownLonNames(), 'lons') - lat_known_names <- c(s2dv:::.KnownLatNames(), 'lats') + lon_known_names <- c(.KnownLonNames(), 'lons') + lat_known_names <- c(.KnownLatNames(), 'lats') lon_dim <- lon_known_names[which(lon_known_names %in% names(ref_grid))] lat_dim <- lat_known_names[which(lat_known_names %in% names(ref_grid))] diff --git a/R/Utils.R b/R/Utils.R index 383a46f855ca3d8790cce574809909b9e879b935..b4874d06a4f56382f2bc020ff455bdc524d74a46 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -99,3 +99,146 @@ list(fun = saveToFile, files = fileout) } +#Draws Color Bars for Categories +#A wrapper of ColorBar to generate multiple color bars for different +#categories, and each category has different color set. +#Draws Color Bars for Categories +#A wrapper of ColorBarContinuous to generate multiple color bars for different +#categories, and each category has different color set. +GradientCatsColorBar <- function(nmap, brks = NULL, cols = NULL, vertical = TRUE, subsampleg = NULL, + bar_limits, var_limits = NULL, + triangle_ends = NULL, col_inf = NULL, col_sup = NULL, plot = TRUE, + draw_separators = FALSE, + bar_titles = NULL, title_scale = 1, label_scale = 1, extra_margin = rep(0, 4), + ...) { + + # bar_limits: a vector of 2 or a list + if (!is.list(bar_limits)) { + if (!is.numeric(bar_limits) || length(bar_limits) != 2) { + stop("Parameter 'bar_limits' must be a numeric vector of length 2 or a list containing that.") + } + # turn into list + bar_limits <- rep(list(bar_limits), nmap) + } else { + if (any(!sapply(bar_limits, is.numeric)) || any(sapply(bar_limits, length) != 2)) { + stop("Parameter 'bar_limits' must be a numeric vector of length 2 or a list containing that.") + } + if (length(bar_limits) != nmap) { + stop("Parameter 'bar_limits' must have the length of 'nmap'.") + } + } + # Check brks + if (!is.list(brks)) { + if (is.null(brks)) { + brks <- 5 + } else if (!is.numeric(brks)) { + stop("Parameter 'brks' must be a numeric vector.") + } + # Turn it into list + brks <- rep(list(brks), nmap) + } else { + if (length(brks) != nmap) { + stop("Parameter 'brks' must have the length of 'nmap'.") + } + } + for (i_map in 1:nmap) { + if (length(brks[[i_map]]) == 1) { + brks[[i_map]] <- seq(from = bar_limits[[i_map]][1], to = bar_limits[[i_map]][2], length.out = brks[[i_map]]) + } + } + + # Check cols + col_sets <- list(c("#A1D99B", "#74C476", "#41AB5D", "#238B45"), + c("#6BAED6FF", "#4292C6FF", "#2171B5FF", "#08519CFF"), + c("#FFEDA0FF", "#FED976FF", "#FEB24CFF", "#FD8D3CFF"), + c("#FC4E2AFF", "#E31A1CFF", "#BD0026FF", "#800026FF"), + c("#FCC5C0", "#FA9FB5", "#F768A1", "#DD3497")) + if (is.null(cols)) { + if (length(col_sets) >= nmap) { + chosen_sets <- 1:nmap + chosen_sets <- chosen_sets + floor((length(col_sets) - length(chosen_sets)) / 2) + } else { + chosen_sets <- array(1:length(col_sets), nmap) + } + cols <- col_sets[chosen_sets] + + # Set triangle_ends, col_sup, col_inf + #NOTE: The "col" input of ColorBar() later is not NULL (since we determine it here) + # so ColorBar() cannot decide these parameters for us. + #NOTE: Here, col_inf and col_sup are prior to triangle_ends, which is consistent with ColorBar(). + #TODO: Make triangle_ends a list + if (is.null(triangle_ends)) { + if (!is.null(var_limits)) { + triangle_ends <- c(FALSE, FALSE) + #TODO: bar_limits is a list + if (bar_limits[1] >= var_limits[1] | !is.null(col_inf)) { + triangle_ends[1] <- TRUE + if (is.null(col_inf)) { + col_inf <- lapply(cols, head, 1) + cols <- lapply(cols, '[', -1) + } + } + if (bar_limits[2] < var_limits[2] | !is.null(col_sup)) { + triangle_ends[2] <- TRUE + if (is.null(col_sup)) { + col_sup <- lapply(cols, tail, 1) + cols <- lapply(cols, '[', -length(cols[[1]])) + } + } + } else { + triangle_ends <- c(!is.null(col_inf), !is.null(col_sup)) + } + } else { # triangle_ends has values + if (triangle_ends[1] & is.null(col_inf)) { + col_inf <- lapply(cols, head, 1) + cols <- lapply(cols, '[', -1) + } + if (triangle_ends[2] & is.null(col_sup)) { + col_sup <- lapply(cols, tail, 1) + cols <- lapply(cols, '[', -length(cols[[1]])) + } + } + + } else { + if (!is.list(cols)) { + stop("Parameter 'cols' must be a list of character vectors.") + } + if (!all(sapply(cols, is.character))) { + stop("Parameter 'cols' must be a list of character vectors.") + } + if (length(cols) != nmap) { + stop("Parameter 'cols' must be a list of the same length as 'nmap'.") + } + } + for (i_map in 1:length(cols)) { + if (length(cols[[i_map]]) != (length(brks[[i_map]]) - 1)) { + cols[[i_map]] <- grDevices::colorRampPalette(cols[[i_map]])(length(brks[[i_map]]) - 1) + } + } + + # Check bar_titles + if (is.null(bar_titles)) { + if (nmap == 3) { + bar_titles <- c("Below normal (%)", "Normal (%)", "Above normal (%)") + } else if (nmap == 5) { + bar_titles <- c("Low (%)", "Below normal (%)", + "Normal (%)", "Above normal (%)", "High (%)") + } else { + bar_titles <- paste0("Cat. ", 1:nmap, " (%)") + } + } + + if (plot) { + for (k in 1:nmap) { + ColorBarContinuous(brks = brks[[k]], cols = cols[[k]], vertical = FALSE, subsampleg = subsampleg, + bar_limits = bar_limits[[k]], #var_limits = var_limits, + triangle_ends = triangle_ends, col_inf = col_inf[[k]], col_sup = col_sup[[k]], plot = TRUE, + draw_separators = draw_separators, + title = bar_titles[[k]], title_scale = title_scale, + label_scale = label_scale, extra_margin = extra_margin) + } + } else { + return(list(brks = brks, cols = cols, col_inf = col_inf, col_sup = col_sup)) + } + +} diff --git a/man/PlotCombinedMap.Rd b/man/PlotCombinedMap.Rd new file mode 100644 index 0000000000000000000000000000000000000000..d9ec5badcb218b249348b3a51a261f16e12006a4 --- /dev/null +++ b/man/PlotCombinedMap.Rd @@ -0,0 +1,184 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PlotCombinedMap.R +\name{PlotCombinedMap} +\alias{PlotCombinedMap} +\title{Plot Multiple Lon-Lat Variables In a Single Map According to a Decision Function} +\usage{ +PlotCombinedMap( + maps, + lon, + lat, + map_select_fun, + display_range, + map_dim = "map", + brks = NULL, + cols = NULL, + bar_limits = NULL, + triangle_ends = c(F, F), + col_inf = NULL, + col_sup = NULL, + col_unknown_map = "white", + mask = NULL, + col_mask = "grey", + dots = NULL, + bar_titles = NULL, + legend_scale = 1, + cex_bar_titles = 1.5, + plot_margin = NULL, + bar_extra_margin = c(2, 0, 2, 0), + fileout = NULL, + width = 8, + height = 5, + size_units = "in", + res = 100, + drawleg = T, + return_leg = FALSE, + ... +) +} +\arguments{ +\item{maps}{List of matrices to plot, each with (longitude, latitude) +dimensions, or 3-dimensional array with the dimensions (longitude, latitude, +map). Dimension names are required.} + +\item{lon}{Vector of longitudes. Must match the length of the corresponding +dimension in 'maps'.} + +\item{lat}{Vector of latitudes. Must match the length of the corresponding +dimension in 'maps'.} + +\item{map_select_fun}{Function that selects, for each grid point, which value +to take among all the provided maps. This function receives as input a +vector of values for a same grid point for all the provided maps, and must +return a single selected value (not its index!) or NA. For example, the +\code{min} and \code{max} functions are accepted.} + +\item{display_range}{Range of values to be displayed for all the maps. This +must be a numeric vector c(range min, range max). The values in the +parameter 'maps' can go beyond the limits specified in this range. If the +selected value for a given grid point (according to 'map_select_fun') falls +outside the range, it will be coloured with 'col_unknown_map'.} + +\item{map_dim}{Optional name for the dimension of 'maps' along which the +multiple maps are arranged. Only applies when 'maps' is provided as a +3-dimensional array. Takes the value 'map' by default.} + +\item{brks}{Colour levels to be sent to PlotEquiMap. This parameter is +optional and adjusted automatically by the function.} + +\item{cols}{List of vectors of colours to be sent to PlotEquiMap for the +colour bar of each map. This parameter is optional and adjusted +automatically by the function (up to 5 maps). The colours provided for each +colour bar will be automatically interpolated to match the number of breaks. +Each item in this list can be named, and the name will be used as title for +the corresponding colour bar (equivalent to the parameter 'bar_titles').} + +\item{col_unknown_map}{Colour to use to paint the grid cells for which a map +is not possible to be chosen according to 'map_select_fun' or for those +values that go beyond 'display_range'. Takes the value 'white' by default.} + +\item{mask}{Optional numeric array with dimensions (latitude, longitude), with +values in the range [0, 1], indicating the opacity of the mask over each +grid point. Cells with a 0 will result in no mask, whereas cells with a 1 +will result in a totally opaque superimposed pixel coloured in 'col_mask'.} + +\item{col_mask}{Colour to be used for the superimposed mask (if specified in +'mask'). Takes the value 'grey' by default.} + +\item{dots}{Array of same dimensions as 'var' or with dimensions +c(n, dim(var)), where n is the number of dot/symbol layers to add to the +plot. A value of TRUE at a grid cell will draw a dot/symbol on the +corresponding square of the plot. By default all layers provided in 'dots' +are plotted with dots, but a symbol can be specified for each of the +layers via the parameter 'dot_symbol'.} + +\item{bar_titles}{Optional vector of character strings providing the titles to +be shown on top of each of the colour bars.} + +\item{legend_scale}{Scale factor for the size of the colour bar labels. Takes +1 by default.} + +\item{cex_bar_titles}{Scale factor for the sizes of the bar titles. Takes 1.5 +by default.} + +\item{plot_margin}{Numeric vector of length 4 for the margin sizes in the +following order: bottom, left, top, and right. If not specified, use the +default of par("mar"), c(5.1, 4.1, 4.1, 2.1). Used as 'margin_scale' in +PlotEquiMap.} + +\item{fileout}{File where to save the plot. If not specified (default) a +graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp +and tiff} + +\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{drawleg}{Where to draw the common colour bar. Can take values TRUE, +FALSE or:\cr +'up', 'u', 'U', 'top', 't', 'T', 'north', 'n', 'N'\cr +'down', 'd', 'D', 'bottom', 'b', 'B', 'south', 's', 'S' (default)\cr +'right', 'r', 'R', 'east', 'e', 'E'\cr +'left', 'l', 'L', 'west', 'w', 'W'} + +\item{return_leg}{A logical value indicating if the color bars information +should be returned by the function. If TRUE, the function doesn't plot the +color bars but still creates the layout with color bar areas, and the +arguments for GradientCatsColorBar() or ColorBarContinuous() will be +returned. It is convenient for users to adjust the color bars manually. The +default is FALSE, the color bars will be plotted directly.} + +\item{...}{Additional parameters to be passed on to \code{PlotEquiMap}.} +} +\description{ +Plot a number a two dimensional matrices with (longitude, +latitude) dimensions on a single map with the cylindrical equidistant +latitude and longitude projection. +} +\examples{ +# Simple example +x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200 +a <- x * 0.6 +b <- (1 - x) * 0.6 +c <- 1 - (a + b) +lons <- seq(0, 359.5, length = 20) +lats <- seq(-89.5, 89.5, length = 10) +\dontrun{ +PlotCombinedMap(list(a, b, c), lons, lats, + toptitle = 'Maximum map', + map_select_fun = max, + display_range = c(0, 1), + bar_titles = paste('\% of belonging to', c('a', 'b', 'c')), + brks = 20, width = 12, height = 10) +} + +Lon <- c(0:40, 350:359) +Lat <- 51:26 +data <- rnorm(51 * 26 * 3) +dim(data) <- c(map = 3, lon = 51, lat = 26) +mask <- sample(c(0,1), replace = TRUE, size = 51 * 26) +dim(mask) <- c(lat = 26, lon = 51) +\dontrun{ +PlotCombinedMap(data, lon = Lon, lat = Lat, map_select_fun = max, + display_range = range(data), mask = mask, + width = 14, height = 10) +} + +} +\seealso{ +\code{PlotCombinedMap} and \code{PlotEquiMap} +} +\author{ +Nicolau Manubens, \email{nicolau.manubens@bsc.es} + +Veronica Torralba, \email{veronica.torralba@bsc.es} +} diff --git a/man/PlotEquiMap.Rd b/man/PlotEquiMap.Rd new file mode 100644 index 0000000000000000000000000000000000000000..bc470e24f1e2a4eefe4931cd58d0e4e2c9d47771 --- /dev/null +++ b/man/PlotEquiMap.Rd @@ -0,0 +1,408 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PlotEquiMap.R +\name{PlotEquiMap} +\alias{PlotEquiMap} +\title{Maps A Two-Dimensional Variable On A Cylindrical Equidistant Projection} +\usage{ +PlotEquiMap( + var, + lon, + lat, + varu = NULL, + varv = NULL, + toptitle = NULL, + sizetit = NULL, + units = NULL, + brks = NULL, + cols = NULL, + bar_limits = NULL, + triangle_ends = NULL, + col_inf = NULL, + col_sup = NULL, + colNA = NULL, + color_fun = ClimPalette(), + square = TRUE, + filled.continents = NULL, + filled.oceans = FALSE, + country.borders = FALSE, + coast_color = NULL, + coast_width = 1, + lake_color = NULL, + shapefile = NULL, + shapefile_color = NULL, + shapefile_lwd = 1, + contours = NULL, + brks2 = NULL, + contour_lwd = 0.5, + contour_color = "black", + contour_lty = 1, + contour_draw_label = TRUE, + contour_label_scale = 1, + dots = NULL, + dot_symbol = 4, + dot_size = 1, + arr_subsamp = floor(length(lon)/30), + arr_scale = 1, + arr_ref_len = 15, + arr_units = "m/s", + arr_scale_shaft = 1, + arr_scale_shaft_angle = 1, + axelab = TRUE, + labW = FALSE, + lab_dist_x = NULL, + lab_dist_y = NULL, + degree_sym = FALSE, + intylat = 20, + intxlon = 20, + xlonshft = 0, + ylatshft = 0, + xlabels = NULL, + ylabels = NULL, + axes_tick_scale = 1, + axes_label_scale = 1, + drawleg = TRUE, + subsampleg = NULL, + bar_extra_labels = NULL, + draw_bar_ticks = TRUE, + draw_separators = FALSE, + triangle_ends_scale = 1, + bar_label_digits = 4, + bar_label_scale = 1, + units_scale = 1, + bar_tick_scale = 1, + bar_extra_margin = rep(0, 4), + boxlim = NULL, + boxcol = "purple2", + boxlwd = 5, + margin_scale = rep(1, 4), + title_scale = 1, + numbfig = NULL, + fileout = NULL, + width = 8, + height = 5, + size_units = "in", + res = 100, + ... +) +} +\arguments{ +\item{var}{Array with the values at each cell of a grid on a regular +rectangular or gaussian grid. The array is expected to have two +dimensions: c(latitude, longitude). Longitudes can be in ascending or +descending order and latitudes in any order. It can contain NA values +(coloured with 'colNA'). Arrays with dimensions c(longitude, latitude) +will also be accepted but 'lon' and 'lat' will be used to disambiguate so +this alternative is not appropriate for square arrays. It is allowed that +the positions of the longitudinal and latitudinal coordinate dimensions +are interchanged.} + +\item{lon}{Numeric vector of longitude locations of the cell centers of the +grid of 'var', in ascending or descending order (same as 'var'). Expected +to be regularly spaced, within either of the ranges [-180, 180] or +[0, 360]. Data for two adjacent regions split by the limits of the +longitude range can also be provided, e.g. \code{lon = c(0:50, 300:360)} +('var' must be provided consitently).} + +\item{lat}{Numeric vector of latitude locations of the cell centers of the +grid of 'var', in any order (same as 'var'). Expected to be from a regular +rectangular or gaussian grid, within the range [-90, 90].} + +\item{varu}{Array of the zonal component of wind/current/other field with +the same dimensions as 'var'. It is allowed that the positions of the +longitudinal and latitudinal coordinate dimensions are interchanged.} + +\item{varv}{Array of the meridional component of wind/current/other field +with the same dimensions as 'var'. It is allowed that the positions of the +longitudinal and latitudinal coordinate dimensions are interchanged.} + +\item{toptitle}{Top title of the figure, scalable with parameter +'title_scale'.} + +\item{sizetit}{Scale factor for the figure top title provided in parameter +'toptitle'. Deprecated. Use 'title_scale' instead.} + +\item{units}{Title at the top of the colour bar, most commonly the units of +the variable provided in parameter 'var'.} + +\item{brks, cols, bar_limits, triangle_ends}{Usually only providing 'brks' is +enough to generate the desired colour bar. These parameters allow to +define n breaks that define n - 1 intervals to classify each of the values +in 'var'. The corresponding grid cell of a given value in 'var' will be +coloured in function of the interval it belongs to. These parameters are +sent to \code{ColorBar()} to generate the breaks and colours. Additional +colours for values beyond the limits of the colour bar are also generated +and applied to the plot if 'bar_limits' or 'brks' and 'triangle_ends' are +properly provided to do so. See ?ColorBar for a full explanation.} + +\item{col_inf, col_sup, colNA}{Colour identifiers to colour the values in +'var' that go beyond the extremes of the colour bar and to colour NA +values, respectively. 'colNA' takes attr(cols, 'na_color') if available by +default, where cols is the parameter 'cols' if provided or the vector of +colors returned by 'color_fun'. If not available, it takes 'pink' by +default. 'col_inf' and 'col_sup' will take the value of 'colNA' if not +specified. See ?ColorBar for a full explanation on 'col_inf' and 'col_sup'.} + +\item{color_fun, subsampleg, bar_extra_labels, draw_bar_ticks}{Set of +parameters to control the visual aspect of the drawn colour bar +(1/3). See ?ColorBar for a full explanation.} + +\item{square}{Logical value to choose either to draw a coloured square for +each grid cell in 'var' (TRUE; default) or to draw contour lines and fill +the spaces in between with colours (FALSE). In the latter case, +'filled.continents' will take the value FALSE if not specified.} + +\item{filled.continents}{Colour to fill in drawn projected continents. +Takes the value gray(0.5) by default or, if 'square = FALSE', takes the +value FALSE. If set to FALSE, continents are not filled in.} + +\item{filled.oceans}{A logical value or the color name to fill in drawn +projected oceans. The default value is FALSE. If it is TRUE, the default +colour is "light blue".} + +\item{country.borders}{A logical value indicating if the country borders +should be plotted (TRUE) or not (FALSE). It only works when +'filled.continents' is FALSE. The default value is FALSE.} + +\item{coast_color}{Colour of the coast line of the drawn projected continents. +Takes the value gray(0.5) by default.} + +\item{coast_width}{Line width of the coast line of the drawn projected +continents. Takes the value 1 by default.} + +\item{lake_color}{Colour of the lake or other water body inside continents. +The default value is NULL.} + +\item{shapefile}{A character string of the path to a .rds file or a list +object containinig shape file data. If it is a .rds file, it should contain +a list. The list should contains 'x' and 'y' at least, which indicate the +location of the shape. The default value is NULL.} + +\item{shapefile_color}{Line color of the shapefile.} + +\item{shapefile_lwd}{Line width of the shapefile. The default value is 1.} + +\item{contours}{Array of same dimensions as 'var' to be added to the plot +and displayed with contours. Parameter 'brks2' is required to define the +magnitude breaks for each contour curve. Disregarded if 'square = FALSE'. +It is allowed that the positions of the longitudinal and latitudinal +coordinate dimensions are interchanged.} + +\item{brks2}{Vector of magnitude breaks where to draw contour curves for the +array provided in 'contours' or if 'square = FALSE'.} + +\item{contour_lwd}{Line width of the contour curves provided via 'contours' +and 'brks2', or if 'square = FALSE'.} + +\item{contour_color}{Line color of the contour curves provided via 'contours' +and 'brks2', or if 'square = FALSE'.} + +\item{contour_lty}{Line type of the contour curves. Takes 1 (solid) by +default. See help on 'lty' in par() for other accepted values.} + +\item{contour_draw_label}{A logical value indicating whether to draw the +contour labels or not. The default value is TRUE.} + +\item{contour_label_scale}{Scale factor for the superimposed labels when +drawing contour levels.} + +\item{dots}{Array of same dimensions as 'var' or with dimensions +c(n, dim(var)), where n is the number of dot/symbol layers to add to the +plot. A value of TRUE at a grid cell will draw a dot/symbol on the +corresponding square of the plot. By default all layers provided in 'dots' +are plotted with dots, but a symbol can be specified for each of the +layers via the parameter 'dot_symbol'. It is allowed that the positions of +the longitudinal and latitudinal coordinate dimensions are interchanged.} + +\item{dot_symbol}{Single character/number or vector of characters/numbers +that correspond to each of the symbol layers specified in parameter 'dots'. +If a single value is specified, it will be applied to all the layers in +'dots'. Takes 15 (centered square) by default. See 'pch' in par() for +additional accepted options.} + +\item{dot_size}{Scale factor for the dots/symbols to be plotted, specified +in 'dots'. If a single value is specified, it will be applied to all +layers in 'dots'. Takes 1 by default.} + +\item{arr_subsamp}{Subsampling factor to select a subset of arrows in +'varu' and 'varv' to be drawn. Only one out of arr_subsamp arrows will +be drawn. Takes 1 by default.} + +\item{arr_scale}{Scale factor for drawn arrows from 'varu' and 'varv'. +Takes 1 by default.} + +\item{arr_ref_len}{Length of the refence arrow to be drawn as legend at the +bottom of the figure (in same units as 'varu' and 'varv', only affects the +legend for the wind or variable in these arrays). Defaults to 15.} + +\item{arr_units}{Units of 'varu' and 'varv', to be drawn in the legend. +Takes 'm/s' by default.} + +\item{arr_scale_shaft}{Parameter for the scale of the shaft of the arrows +(which also depend on the number of figures and the arr_scale parameter). +Defaults to 1.} + +\item{arr_scale_shaft_angle}{Parameter for the scale of the angle of the +shaft of the arrows (which also depend on the number of figure and the +arr_scale parameter). Defaults to 1.} + +\item{axelab}{Whether to draw longitude and latitude axes or not. +TRUE by default.} + +\item{labW}{Whether to label the longitude axis with a 'W' instead of minus +for negative values. Defaults to FALSE.} + +\item{lab_dist_x}{A numeric of the distance of the longitude labels to the +box borders. The default value is NULL and is automatically adjusted by +the function.} + +\item{lab_dist_y}{A numeric of the distance of the latitude labels to the +box borders. The default value is NULL and is automatically adjusted by +the function.} + +\item{degree_sym}{A logical indicating whether to include degree symbol +(30° N) or not (30N; default).} + +\item{intylat}{Interval between latitude ticks on y-axis, in degrees. +Defaults to 20.} + +\item{intxlon}{Interval between latitude ticks on x-axis, in degrees. +Defaults to 20.} + +\item{xlonshft}{A numeric of the degrees to shift the latitude ticks. The +default value is 0.} + +\item{ylatshft}{A numeric of the degrees to shift the longitude ticks. The +default value is 0.} + +\item{xlabels}{A vector of character string of the custumized x-axis labels. +The values should correspond to each tick, which is decided by the longitude +and parameter 'intxlon'. The default value is NULL and the labels will be +automatically generated.} + +\item{ylabels}{A vector of character string of the custumized y-axis labels. +The values should correspond to each tick, which is decided by the latitude +and parameter 'intylat'. The default value is NULL and the labels will be +automatically generated.} + +\item{axes_tick_scale}{Scale factor for the tick lines along the longitude +and latitude axes.} + +\item{axes_label_scale}{Scale factor for the labels along the longitude +and latitude axes.} + +\item{drawleg}{Whether to plot a color bar (legend, key) or not. Defaults to +TRUE. It is not possible to plot the colour bar if 'add = TRUE'. Use +ColorBar() and the return values of PlotEquiMap() instead.} + +\item{draw_separators, triangle_ends_scale, bar_label_digits}{Set of +parameters to control the visual aspect of the drawn colour bar +(2/3). See ?ColorBar for a full explanation.} + +\item{bar_label_scale, units_scale, bar_tick_scale, bar_extra_margin}{Set of +parameters to control the visual aspect of the drawn colour bar (3/3). +See ?ColorBar for a full explanation.} + +\item{boxlim}{Limits of a box to be added to the plot, in degrees: +c(x1, y1, x2, y2). A list with multiple box specifications can also be +provided.} + +\item{boxcol}{Colour of the box lines. A vector with a colour for each of +the boxes is also accepted. Defaults to 'purple2'.} + +\item{boxlwd}{Line width of the box lines. A vector with a line width for +each of the boxes is also accepted. Defaults to 5.} + +\item{margin_scale}{Scale factor for the margins around the map plot, with +the format c(y1, x1, y2, x2). Defaults to rep(1, 4). If drawleg = TRUE, +then margin_scale[1] is subtracted 1 unit.} + +\item{title_scale}{Scale factor for the figure top title. Defaults to 1.} + +\item{numbfig}{Number of figures in the layout the plot will be put into. +A higher numbfig will result in narrower margins and smaller labels, +axe labels, ticks, thinner lines, ... Defaults to 1.} + +\item{fileout}{File where to save the plot. If not specified (default) a +graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, +bmp and tiff.} + +\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{\dots}{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 font font.axis font.lab font.main font.sub lend +lheight ljoin lmitre mex mfcol mfrow mfg mkh omd omi page pch pin plt +pty smo srt tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog \cr +For more information about the parameters see `par`.} +} +\value{ +\item{brks}{ + Breaks used for colouring the map (and legend if drawleg = TRUE). +} +\item{cols}{ + Colours used for colouring the map (and legend if drawleg = TRUE). Always + of length length(brks) - 1. +} +\item{col_inf}{ + Colour used to draw the lower triangle end in the colour bar (NULL if not + drawn at all). +} +\item{col_sup}{ + Colour used to draw the upper triangle end in the colour bar (NULL if not + drawn at all). +} +} +\description{ +Map longitude-latitude array (on a regular rectangular or gaussian grid) +on a cylindrical equidistant latitude and longitude projection with coloured +grid cells. Only the region for which data has been provided is displayed. +A colour bar (legend) can be plotted and adjusted. It is possible to draw +superimposed arrows, dots, symbols, contour lines and boxes. A number of +options is provided to adjust the position, size and colour of the +components. Some parameters are provided to add and adjust the masks that +include continents, oceans, and lakes. This plot function is compatible with +figure layouts if colour bar is disabled. +} +\examples{ +# See examples on Load() to understand the first lines in this example + \dontrun{ +data_path <- system.file('sample_data', package = 's2dv') +expA <- list(name = 'experiment', path = file.path(data_path, + 'model/$EXP_NAME$/$STORE_FREQ$_mean/$VAR_NAME$_3hourly', + '$VAR_NAME$_$START_DATE$.nc')) +obsX <- list(name = 'observation', path = file.path(data_path, + '$OBS_NAME$/$STORE_FREQ$_mean/$VAR_NAME$', + '$VAR_NAME$_$YEAR$$MONTH$.nc')) + +# Now we are ready to use Load(). +startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101') +sampleData <- Load('tos', list(expA), list(obsX), startDates, + leadtimemin = 1, leadtimemax = 4, output = 'lonlat', + latmin = 27, latmax = 48, lonmin = -12, lonmax = 40) + } + \dontshow{ +startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101') +sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'), + c('observation'), startDates, + leadtimemin = 1, + leadtimemax = 4, + output = 'lonlat', + latmin = 27, latmax = 48, + lonmin = -12, lonmax = 40) + } +PlotEquiMap(sampleData$mod[1, 1, 1, 1, , ], sampleData$lon, sampleData$lat, + toptitle = 'Predicted sea surface temperature for Nov 1960 from 1st Nov', + title_scale = 0.5) +} diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd new file mode 100644 index 0000000000000000000000000000000000000000..31e3001fce60c33ce81730189eec426c55f5c5ce --- /dev/null +++ b/man/PlotForecastPDF.Rd @@ -0,0 +1,84 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PlotForecastPDF.R +\name{PlotForecastPDF} +\alias{PlotForecastPDF} +\title{Plot one or multiple ensemble forecast pdfs for the same event} +\usage{ +PlotForecastPDF( + fcst, + tercile.limits, + extreme.limits = NULL, + obs = NULL, + plotfile = NULL, + title = "Set a title", + var.name = "Varname (units)", + fcst.names = NULL, + add.ensmemb = c("above", "below", "no"), + color.set = c("ggplot", "s2s4e", "hydro", "vitigeoss"), + memb_dim = "member" +) +} +\arguments{ +\item{fcst}{A dataframe or array containing all the ensember members for each +forecast. If \code{'fcst'} is an array, it should have two labelled +dimensions, and one of them should be \code{'members'}. If \code{'fcsts'} is +a data.frame, each column shoul be a separate forecast, with the rows beeing +the different ensemble members.} + +\item{tercile.limits}{An array or vector with P33 and P66 values that define +the tercile categories for each panel. Use an array of dimensions +(nforecasts,2) to define different terciles for each forecast panel, or a +vector with two elements to reuse the same tercile limits for all forecast +panels.} + +\item{extreme.limits}{(optional) An array or vector with P10 and P90 values +that define the extreme categories for each panel. Use an array of +(nforecasts,2) to define different extreme limits for each forecast panel, +or a vector with two elements to reuse the same tercile limits for all +forecast panels. (Default: extreme categories are not shown).} + +\item{obs}{(optional) A vector providing the observed values for each forecast +panel or a single value that will be reused for all forecast panels. +(Default: observation is not shown).} + +\item{plotfile}{(optional) A filename (pdf, png...) where the plot will be +saved. (Default: the plot is not saved).} + +\item{title}{A string with the plot title.} + +\item{var.name}{A string with the variable name and units.} + +\item{fcst.names}{(optional) An array of strings with the titles of each +individual forecast.} + +\item{add.ensmemb}{Either to add the ensemble members \code{'above'} (default) +or \code{'below'} the pdf, or not (\code{'no'}).} + +\item{color.set}{A selection of predefined color sets: use \code{'ggplot'} +(default) for blue/green/red, \code{'s2s4e'} for blue/grey/orange, +\code{'hydro'} for yellow/gray/blue (suitable for precipitation and +inflows) or the \code{"vitigeoss"} color set.} + +\item{memb_dim}{A character string indicating the name of the member +dimension.} +} +\value{ +A ggplot object containing the plot. +} +\description{ +This function plots the probability distribution function of +several ensemble forecasts. Separate panels are used to plot forecasts valid +or initialized at different times or by different models or even at different +locations. Probabilities for tercile categories are computed, plotted in +colors and annotated. An asterisk marks the tercile with higher probabilities. +Probabilities for extreme categories (above P90 and below P10) can also be +included as hatched areas. Individual ensemble members can be plotted as +jittered points. The observed value is optionally shown as a diamond. +} +\examples{ +fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2)) +PlotForecastPDF(fcsts,c(-1,1)) +} +\author{ +Llorenç Lledó \email{llledo@bsc.es} +} diff --git a/man/PlotMostLikelyQuantileMap.Rd b/man/PlotMostLikelyQuantileMap.Rd new file mode 100644 index 0000000000000000000000000000000000000000..0dde63ff0480fd076769cc77b48a14ac9ec563d2 --- /dev/null +++ b/man/PlotMostLikelyQuantileMap.Rd @@ -0,0 +1,160 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PlotMostLikelyQuantileMap.R +\name{PlotMostLikelyQuantileMap} +\alias{PlotMostLikelyQuantileMap} +\title{Plot Maps of Most Likely Quantiles} +\usage{ +PlotMostLikelyQuantileMap( + probs, + lon, + lat, + cat_dim = "bin", + bar_titles = NULL, + col_unknown_cat = "white", + drawleg = T, + ... +) +} +\arguments{ +\item{probs}{A list of bi-dimensional arrays with the named dimensions +'latitude' (or 'lat') and 'longitude' (or 'lon'), with equal size and in the +same order, or a single tri-dimensional array with an additional dimension +(e.g. 'bin') for the different categories. The arrays must contain +probability values between 0 and 1, and the probabilities for all categories +of a grid cell should not exceed 1 when added.} + +\item{lon}{A numeric vector with the longitudes of the map grid, in the same +order as the values along the corresponding dimension in \code{probs}.} + +\item{lat}{A numeric vector with the latitudes of the map grid, in the same +order as the values along the corresponding dimension in \code{probs}.} + +\item{cat_dim}{The name of the dimension along which the different categories +are stored in \code{probs}. This only applies if \code{probs} is provided in +the form of 3-dimensional array. The default expected name is 'bin'.} + +\item{bar_titles}{Vector of character strings with the names to be drawn on +top of the color bar for each of the categories. As many titles as +categories provided in \code{probs} must be provided.} + +\item{col_unknown_cat}{Character string with a colour representation of the +colour to be used to paint the cells for which no category can be clearly +assigned. Takes the value 'white' by default.} + +\item{drawleg}{Where to draw the common colour bar. Can take values TRUE, +FALSE or:\cr +'up', 'u', 'U', 'top', 't', 'T', 'north', 'n', 'N'\cr +'down', 'd', 'D', 'bottom', 'b', 'B', 'south', 's', 'S' (default)\cr +'right', 'r', 'R', 'east', 'e', 'E'\cr +'left', 'l', 'L', 'west', 'w', 'W'} + +\item{...}{Additional parameters to be sent to \code{PlotCombinedMap} and +\code{PlotEquiMap}.} +} +\description{ +This function receives as main input (via the parameter +\code{probs}) a collection of longitude-latitude maps, each containing the +probabilities (from 0 to 1) of the different grid cells of belonging to a +category. As many categories as maps provided as inputs are understood to +exist. The maps of probabilities must be provided on a common rectangular +regular grid, and a vector with the longitudes and a vector with the latitudes +of the grid must be provided. The input maps can be provided in two forms, +either as a list of multiple two-dimensional arrays (one for each category) or +as a three-dimensional array, where one of the dimensions corresponds to the +different categories. +} +\examples{ +# Simple example +x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200 +a <- x * 0.6 +b <- (1 - x) * 0.6 +c <- 1 - (a + b) +lons <- seq(0, 359.5, length = 20) +lats <- seq(-89.5, 89.5, length = 10) +\dontrun{ +PlotMostLikelyQuantileMap(list(a, b, c), lons, lats, + toptitle = 'Most likely tercile map', + bar_titles = paste('\% of belonging to', c('a', 'b', 'c')), + brks = 20, width = 10, height = 8) +} + +# 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) + +# 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 + +# 2. Binning sample data +prob_thresholds <- 1:n_bins / n_bins +prob_thresholds <- prob_thresholds[1:(n_bins - 1)] +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 + for (i in 1:n_bins) { + result[i] <- sum(x > lower_threshold & x <= thresholds[i]) / n_samples + lower_threshold <- thresholds[i] + } + + dim(result) <- c(bin = n_bins) + result +} + +bins <- multiApply::Apply(sample_data, 'time', binning, thresholds)$output1 + +# 3. Plotting most likely quantile/bin +\dontrun{ +PlotMostLikelyQuantileMap(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) +} +} +\seealso{ +\code{PlotCombinedMap} and \code{PlotEquiMap} +} +\author{ +Veronica Torralba, \email{veronica.torralba@bsc.es}, Nicolau Manubens, +\email{nicolau.manubens@bsc.es} +} diff --git a/man/PlotPDFsOLE.Rd b/man/PlotPDFsOLE.Rd new file mode 100644 index 0000000000000000000000000000000000000000..e2c6606eb20822b0b29d40c99411dbb0de451710 --- /dev/null +++ b/man/PlotPDFsOLE.Rd @@ -0,0 +1,74 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PlotPDFsOLE.R +\name{PlotPDFsOLE} +\alias{PlotPDFsOLE} +\title{Plotting two probability density gaussian functions and the optimal linear +estimation (OLE) as result of combining them.} +\usage{ +PlotPDFsOLE( + pdf_1, + pdf_2, + nsigma = 3, + legendPos = "bottom", + legendSize = 1, + plotfile = NULL, + width = 30, + height = 15, + units = "cm", + dpi = 300 +) +} +\arguments{ +\item{pdf_1}{A numeric array with a dimension named 'statistic', containg +two parameters: mean' and 'standard deviation' of the first gaussian pdf +to combining.} + +\item{pdf_2}{A numeric array with a dimension named 'statistic', containg +two parameters: mean' and 'standard deviation' of the second gaussian pdf +to combining.} + +\item{nsigma}{(optional) A numeric value for setting the limits of X axis. +(Default nsigma = 3).} + +\item{legendPos}{(optional) A character value for setting the position of the +legend ("bottom", "top", "right" or "left")(Default 'bottom').} + +\item{legendSize}{(optional) A numeric value for setting the size of the +legend text. (Default 1.0).} + +\item{plotfile}{(optional) A filename where the plot will be saved. +(Default: the plot is not saved).} + +\item{width}{(optional) A numeric value indicating the plot width in +units ("in", "cm", or "mm"). (Default width = 30).} + +\item{height}{(optional) A numeric value indicating the plot height. +(Default height = 15).} + +\item{units}{(optional) A character value indicating the plot size +unit. (Default units = 'cm').} + +\item{dpi}{(optional) A numeric value indicating the plot resolution. +(Default dpi = 300).} +} +\value{ +PlotPDFsOLE() returns a ggplot object containing the plot. +} +\description{ +This function plots two probability density gaussian functions +and the optimal linear estimation (OLE) as result of combining them. +} +\examples{ +# Example 1 +pdf_1 <- c(1.1,0.6) +attr(pdf_1, "name") <- "NAO1" +dim(pdf_1) <- c(statistic = 2) +pdf_2 <- c(1,0.5) +attr(pdf_2, "name") <- "NAO2" +dim(pdf_2) <- c(statistic = 2) + +PlotPDFsOLE(pdf_1, pdf_2) +} +\author{ +Eroteida Sanchez-Garcia - AEMET, //email{esanchezg@aemet.es} +} diff --git a/man/PlotTriangles4Categories.Rd b/man/PlotTriangles4Categories.Rd new file mode 100644 index 0000000000000000000000000000000000000000..62a076ef4277a81af0540d4296fa451d18f4dfd8 --- /dev/null +++ b/man/PlotTriangles4Categories.Rd @@ -0,0 +1,136 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PlotTriangles4Categories.R +\name{PlotTriangles4Categories} +\alias{PlotTriangles4Categories} +\title{Function to convert any 3-d numerical array to a grid of coloured triangles.} +\usage{ +PlotTriangles4Categories( + data, + brks = NULL, + cols = NULL, + toptitle = NULL, + sig_data = NULL, + pch_sig = 18, + col_sig = "black", + cex_sig = 1, + xlab = TRUE, + ylab = TRUE, + xlabels = NULL, + xtitle = NULL, + ylabels = NULL, + ytitle = NULL, + legend = TRUE, + lab_legend = NULL, + cex_leg = 1, + col_leg = "black", + cex_axis = 1.5, + mar = c(5, 4, 0, 0), + fileout = NULL, + size_units = "px", + res = 100, + figure.width = 1, + ... +) +} +\arguments{ +\item{data}{Array with three named dimensions: 'dimx', 'dimy', 'dimcat', +containing the values to be displayed in a coloured image with triangles.} + +\item{brks}{A vector of the color bar intervals. The length must be one more +than the parameter 'cols'. Use ColorBar() to generate default values.} + +\item{cols}{A vector of valid colour identifiers for color bar. The length +must be one less than the parameter 'brks'. Use ColorBar() to generate +default values.} + +\item{toptitle}{A string of the title of the grid. Set NULL as default.} + +\item{sig_data}{Logical array with the same dimensions as 'data' to add layers +to the plot. A value of TRUE at a grid cell will draw a dot/symbol on the +corresponding triangle of the plot. Set NULL as default.} + +\item{pch_sig}{Symbol to be used to represent sig_data. Takes 18 +(diamond) by default. See 'pch' in par() for additional accepted options.} + +\item{col_sig}{Colour of the symbol to represent sig_data.} + +\item{cex_sig}{Parameter to increase/reduce the size of the symbols used +to represent sig_data.} + +\item{xlab}{A logical value (TRUE) indicating if xlabels should be plotted} + +\item{ylab}{A logical value (TRUE) indicating if ylabels should be plotted} + +\item{xlabels}{A vector of labels of the x-axis The length must be +length of the col of parameter 'data'. Set the sequence from 1 to the +length of the row of parameter 'data' as default.} + +\item{xtitle}{A string of title of the x-axis. Set NULL as default.} + +\item{ylabels}{A vector of labels of the y-axis The length must be +length of the row of parameter 'data'. Set the sequence from 1 to the +length of the row of parameter 'data' as default.} + +\item{ytitle}{A string of title of the y-axis. Set NULL as default.} + +\item{legend}{A logical value to decide to draw the color bar legend or not. +Set TRUE as default.} + +\item{lab_legend}{A vector of labels indicating what is represented in each +category (i.e. triangle). Set the sequence from 1 to the length of +the categories (2 or 4).} + +\item{cex_leg}{A number to indicate the increase/reductuion of the lab_legend +used to represent sig_data.} + +\item{col_leg}{Color of the legend (triangles).} + +\item{cex_axis}{A number to indicate the increase/reduction of the axis labels.} + +\item{mar}{A numerical vector of the form c(bottom, left, top, right) which +gives the number of lines of margin to be specified on the four sides of the +plot.} + +\item{fileout}{A string of full directory path and file name indicating where +to save the plot. If not specified (default), a graphics device will pop up.} + +\item{size_units}{A string indicating the units of the size of the device +(file or window) to plot in. Set 'px' as default. See ?Devices and the +creator function of the corresponding device.} + +\item{res}{A positive number indicating resolution of the device (file or +window) to plot in. See ?Devices and the creator function of the +corresponding device.} + +\item{figure.width}{a numeric value to control the width of the plot.} + +\item{...}{The additional parameters to be passed to function +ColorBarContinuous() in for color legend creation.} +} +\value{ +A figure in popup window by default, or saved to the specified path. +} +\description{ +This function converts a 3-d numerical data array into a coloured +grid with triangles. It is useful for a slide or article to present tabular +results as colors instead of numbers. This can be used to compare the outputs +of two or four categories (e.g. modes of variability, clusters, or forecast +systems). +} +\examples{ +# Example with random data +arr1 <- array(runif(n = 4 * 5 * 4, min = -1, max = 1), dim = c(4,5,4)) +names(dim(arr1)) <- c('dimx', 'dimy', 'dimcat') +arr2 <- array(TRUE, dim = dim(arr1)) +arr2[which(arr1 < 0.3)] <- FALSE +PlotTriangles4Categories(data = arr1, + cols = c('white','#fef0d9','#fdd49e','#fdbb84','#fc8d59'), + brks = c(-1, 0, 0.1, 0.2, 0.3, 0.4), + lab_legend = c('NAO+', 'BL','AR','NAO-'), + xtitle = "Target month", ytitle = "Lead time", + xlabels = c("Jan", "Feb", "Mar", "Apr")) +} +\author{ +History:\cr +1.0 - 2020-10 (V.Torralba, \email{veronica.torralba@bsc.es}) - Original code +} diff --git a/man/PlotWeeklyClim.Rd b/man/PlotWeeklyClim.Rd new file mode 100644 index 0000000000000000000000000000000000000000..3e064e8db54ad73d1cca4ccebb00ec0a1591d067 --- /dev/null +++ b/man/PlotWeeklyClim.Rd @@ -0,0 +1,119 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PlotWeeklyClim.R +\name{PlotWeeklyClim} +\alias{PlotWeeklyClim} +\title{Plots the observed weekly means and climatology of a timeseries data} +\usage{ +PlotWeeklyClim( + data, + first_date, + ref_period, + last_date = NULL, + data_years = NULL, + time_dim = "time", + sdate_dim = "sdate", + ylim = NULL, + title = NULL, + subtitle = NULL, + ytitle = NULL, + legend = TRUE, + palette = "Blues", + fileout = NULL, + device = NULL, + width = 8, + height = 6, + units = "in", + dpi = 300 +) +} +\arguments{ +\item{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'.} + +\item{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.} + +\item{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'.} + +\item{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.} + +\item{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.} + +\item{time_dim}{A character string indicating the daily time dimension name. +The default value is 'time'.} + +\item{sdate_dim}{A character string indicating the start year dimension name. +The default value is 'sdate'.} + +\item{ylim}{A numeric vector of length two providing limits of the scale. +Use NA to refer to the existing minimum or maximum. For more information, +see 'ggplot2' documentation of 'scale_y_continuous' parameter.} + +\item{title}{The text for the top title of the plot. It is NULL by default.} + +\item{subtitle}{The text for the subtitle of the plot. It is NULL bu default.} + +\item{ytitle}{Character string to be drawn as y-axis title. It is NULL by +default.} + +\item{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.} + +\item{palette}{A palette name from the R Color Brewer’s package. The default +value is 'Blues'.} + +\item{fileout}{A character string indicating the file name where to save the +plot. If not specified (default) a graphics device will pop up.} + +\item{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).} + +\item{width}{A numeric value of the plot width in units ("in", "cm", "mm", or +"px"). It is set to 8 by default.} + +\item{height}{A numeric value of the plot height in units ("in", "cm", "mm", +or "px"). It is set to 6 by default.} + +\item{units}{Units of the size of the device (file or window) to plot in. +Inches (’in’) by default.} + +\item{dpi}{A numeric value of the plot resolution. It is set to 300 by +default.} +} +\value{ +A ggplot object containing the plot. +} +\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. +} +\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", ")")) + +}