From f5bf433b755e19b1f1a00edac03c303a259d5e3c Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Thu, 21 Sep 2023 09:56:25 +0200 Subject: [PATCH 1/6] Add CSTools Plotting functions except PlotCombinedMap and PlotMostLikelyQuantileMap --- NAMESPACE | 29 ++ R/PlotForecastPDF.R | 577 ++++++++++++++++++++++++++++++++ R/PlotPDFsOLE.R | 259 ++++++++++++++ R/PlotTriangles4Categories.R | 437 ++++++++++++++++++++++++ R/PlotWeeklyClim.R | 315 +++++++++++++++++ R/Utils.R | 141 ++++++++ man/PlotForecastPDF.Rd | 84 +++++ man/PlotPDFsOLE.Rd | 74 ++++ man/PlotTriangles4Categories.Rd | 136 ++++++++ man/PlotWeeklyClim.Rd | 119 +++++++ 10 files changed, 2171 insertions(+) create mode 100644 R/PlotForecastPDF.R create mode 100644 R/PlotPDFsOLE.R create mode 100644 R/PlotTriangles4Categories.R create mode 100644 R/PlotWeeklyClim.R create mode 100644 man/PlotForecastPDF.Rd create mode 100644 man/PlotPDFsOLE.Rd create mode 100644 man/PlotTriangles4Categories.Rd create mode 100644 man/PlotWeeklyClim.Rd diff --git a/NAMESPACE b/NAMESPACE index a92f986..c767e23 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,11 +4,40 @@ export(ClimColors) export(ClimPalette) export(ColorBarContinuous) export(ColorBarDiscrete) +export(PlotForecastPDF) +export(PlotPDFsOLE) export(PlotRobinson) +export(PlotTriangles4Categories) +export(PlotWeeklyClim) +import(RColorBrewer) import(cowplot) import(ggplot2) +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,col2rgb) importFrom(grDevices,colorRampPalette) +importFrom(grDevices,dev.cur) +importFrom(grDevices,dev.new) +importFrom(grDevices,dev.off) importFrom(grDevices,rgb) +importFrom(graphics,axis) +importFrom(graphics,plot) +importFrom(graphics,points) +importFrom(graphics,polygon) +importFrom(graphics,text) +importFrom(graphics,title) +importFrom(lubridate,ymd) +importFrom(plyr,.) +importFrom(plyr,dlply) +importFrom(reshape2,melt) +importFrom(s2dv,ColorBar) +importFrom(s2dv,InsertDim) +importFrom(s2dv,MeanDims) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R new file mode 100644 index 0000000..7143294 --- /dev/null +++ b/R/PlotForecastPDF.R @@ -0,0 +1,577 @@ +#'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")) + + 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/PlotPDFsOLE.R b/R/PlotPDFsOLE.R new file mode 100644 index 0000000..bf95abb --- /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 0000000..8cf775a --- /dev/null +++ b/R/PlotTriangles4Categories.R @@ -0,0 +1,437 @@ +#'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 ColorBar() in +#' s2dv 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 s2dv ColorBar +#'@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)) + ColorBar(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() +} +.SelectDevice <- function(fileout, width, height, units, res) { + # This function is used in the plot functions to check the extension of the + # files where the graphics will be stored and select the right R device to + # save them. + # If the vector of filenames ('fileout') has files with different + # extensions, then it will only accept the first one, changing all the rest + # of the filenames to use that extension. + + # We extract the extension of the filenames: '.png', '.pdf', ... + ext <- regmatches(fileout, regexpr("\\.[a-zA-Z0-9]*$", fileout)) + + if (length(ext) != 0) { + # If there is an extension specified, select the correct device + ## units of width and height set to accept inches + if (ext[1] == ".png") { + saveToFile <- function(fileout) { + png(filename = fileout, width = width, height = height, res = res, units = units) + } + } else if (ext[1] == ".jpeg") { + saveToFile <- function(fileout) { + jpeg(filename = fileout, width = width, height = height, res = res, units = units) + } + } else if (ext[1] %in% c(".eps", ".ps")) { + saveToFile <- function(fileout) { + postscript(file = fileout, width = width, height = height) + } + } else if (ext[1] == ".pdf") { + saveToFile <- function(fileout) { + pdf(file = fileout, width = width, height = height) + } + } else if (ext[1] == ".svg") { + saveToFile <- function(fileout) { + svg(filename = fileout, width = width, height = height) + } + } else if (ext[1] == ".bmp") { + saveToFile <- function(fileout) { + bmp(filename = fileout, width = width, height = height, res = res, units = units) + } + } else if (ext[1] == ".tiff") { + saveToFile <- function(fileout) { + tiff(filename = fileout, width = width, height = height, res = res, units = units) + } + } else { + .warning("file extension not supported, it will be used '.eps' by default.") + ## In case there is only one filename + fileout[1] <- sub("\\.[a-zA-Z0-9]*$", ".eps", fileout[1]) + ext[1] <- ".eps" + saveToFile <- function(fileout) { + postscript(file = fileout, width = width, height = height) + } + } + # Change filenames when necessary + if (any(ext != ext[1])) { + .warning(paste0("some extensions of the filenames provided in 'fileout' are not ", ext[1],". The extensions are being converted to ", ext[1], ".")) + fileout <- sub("\\.[a-zA-Z0-9]*$", ext[1], fileout) + } + } else { + # Default filenames when there is no specification + .warning("there are no extensions specified in the filenames, default to '.eps'") + fileout <- paste0(fileout, ".eps") + saveToFile <- postscript + } + + # return the correct function with the graphical device, and the correct + # filenames + list(fun = saveToFile, files = fileout) +} + +.warning <- function(...) { + # Function to use the 'warning' R function with our custom settings + # Default: no call information, indent to 0, exdent to 3, + # collapse to \n + args <- list(...) + + ## In case we need to specify warning arguments + if (!is.null(args[["call."]])) { + call <- args[["call."]] + } else { + ## Default: don't show info about the call where the warning came up + call <- FALSE + } + if (!is.null(args[["immediate."]])) { + immediate <- args[["immediate."]] + } else { + ## Default value in warning function + immediate <- FALSE + } + if (!is.null(args[["noBreaks."]])) { + noBreaks <- args[["noBreaks."]] + } else { + ## Default value warning function + noBreaks <- FALSE + } + if (!is.null(args[["domain"]])) { + domain <- args[["domain"]] + } else { + ## Default value warning function + domain <- NULL + } + args[["call."]] <- NULL + args[["immediate."]] <- NULL + args[["noBreaks."]] <- NULL + args[["domain"]] <- NULL + + ## To modify strwrap indent and exdent arguments + if (!is.null(args[["indent"]])) { + indent <- args[["indent"]] + } else { + indent <- 0 + } + if (!is.null(args[["exdent"]])) { + exdent <- args[["exdent"]] + } else { + exdent <- 3 + } + args[["indent"]] <- NULL + args[["exdent"]] <- NULL + + ## To modify paste collapse argument + if (!is.null(args[["collapse"]])) { + collapse <- args[["collapse"]] + } else { + collapse <- "\n!" + } + args[["collapse"]] <- NULL + + ## Warning tag + if (!is.null(args[["tag"]])) { + tag <- args[["tag"]] + } else { + tag <- "! Warning: " + } + args[["tag"]] <- NULL + + warning(paste0(tag, paste(strwrap( + args, indent = indent, exdent = exdent + ), collapse = collapse)), call. = call, immediate. = immediate, + noBreaks. = noBreaks, domain = domain) +} + diff --git a/R/PlotWeeklyClim.R b/R/PlotWeeklyClim.R new file mode 100644 index 0000000..4fa829b --- /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/Utils.R b/R/Utils.R index 94797e7..598ea8d 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -99,3 +99,144 @@ list(fun = saveToFile, files = fileout) } +#Draws Color Bars for Categories +#A wrapper of s2dv::ColorBar 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) { +#TODO: Add s2dv:: + ColorBar(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)) + } + +} \ No newline at end of file diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd new file mode 100644 index 0000000..31e3001 --- /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/PlotPDFsOLE.Rd b/man/PlotPDFsOLE.Rd new file mode 100644 index 0000000..e2c6606 --- /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 0000000..8356da9 --- /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 ColorBar() in +s2dv 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 0000000..3e064e8 --- /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", ")")) + +} -- GitLab From 602ad28c431de4a7399b98d1c4f73ba15150b077 Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Thu, 21 Sep 2023 10:50:51 +0200 Subject: [PATCH 2/6] Remove GradientCatsColorBar from Utils --- R/Utils.R | 141 ------------------------------------------------------ 1 file changed, 141 deletions(-) diff --git a/R/Utils.R b/R/Utils.R index 598ea8d..94797e7 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -99,144 +99,3 @@ list(fun = saveToFile, files = fileout) } -#Draws Color Bars for Categories -#A wrapper of s2dv::ColorBar 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) { -#TODO: Add s2dv:: - ColorBar(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)) - } - -} \ No newline at end of file -- GitLab From 047bdd1bc2e45bba784eadfabc538b9a86fe2d59 Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Thu, 28 Sep 2023 17:42:06 +0200 Subject: [PATCH 3/6] Correct background color of PlotForecastPDF --- R/PlotForecastPDF.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 7143294..0868809 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -441,7 +441,9 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N 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.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), -- GitLab From 5aa7da0a731f47d42b61bf98aadaafcbcf05cc06 Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Thu, 28 Sep 2023 17:43:28 +0200 Subject: [PATCH 4/6] Add PlotCombinedMap and PlotMostLikelyQuantileMap --- R/PlotCombinedMap.R | 529 ++++++++++++++++++++++++++++++++++ R/PlotMostLikelyQuantileMap.R | 222 ++++++++++++++ 2 files changed, 751 insertions(+) create mode 100644 R/PlotCombinedMap.R create mode 100644 R/PlotMostLikelyQuantileMap.R diff --git a/R/PlotCombinedMap.R b/R/PlotCombinedMap.R new file mode 100644 index 0000000..d7d5f49 --- /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 +#' s2dv::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 ColorBar() 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 s2dv PlotEquiMap ColorBar +#'@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)) { + .SelectDevice <- utils::getFromNamespace(".SelectDevice", "s2dv") + 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 <- display_range[2] - display_range[1] #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) + } + 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[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) - display_range[1]) / range_width + if (map_select_fun(x) == display_range[1]) { + res <- res + slightly_tune_val + } + } + } 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/PlotMostLikelyQuantileMap.R b/R/PlotMostLikelyQuantileMap.R new file mode 100644 index 0000000..18e4a8b --- /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, ...) +} -- GitLab From dfa719c641bc85be705d306d5f421ad1ab52fda8 Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Tue, 3 Oct 2023 18:07:14 +0200 Subject: [PATCH 5/6] Add PlotCombinedMap, PlotMostLilelyQuantileMap, PlotTriangles --- NAMESPACE | 19 +++- R/PlotCombinedMap.R | 2 - R/PlotTriangles4Categories.R | 3 +- R/Utils.R | 140 +++++++++++++++++++++++ man/PlotCombinedMap.Rd | 184 +++++++++++++++++++++++++++++++ man/PlotMostLikelyQuantileMap.Rd | 160 +++++++++++++++++++++++++++ 6 files changed, 503 insertions(+), 5 deletions(-) create mode 100644 man/PlotCombinedMap.Rd create mode 100644 man/PlotMostLikelyQuantileMap.Rd diff --git a/NAMESPACE b/NAMESPACE index c767e23..a1504d4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,7 +4,9 @@ export(ClimColors) export(ClimPalette) export(ColorBarContinuous) export(ColorBarDiscrete) +export(PlotCombinedMap) export(PlotForecastPDF) +export(PlotMostLikelyQuantileMap) export(PlotPDFsOLE) export(PlotRobinson) export(PlotTriangles4Categories) @@ -22,22 +24,37 @@ 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,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,ColorBar) importFrom(s2dv,InsertDim) importFrom(s2dv,MeanDims) diff --git a/R/PlotCombinedMap.R b/R/PlotCombinedMap.R index d7d5f49..642f4fd 100644 --- a/R/PlotCombinedMap.R +++ b/R/PlotCombinedMap.R @@ -117,7 +117,6 @@ #' #'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap} #' -#'@importFrom s2dv PlotEquiMap ColorBar #'@importFrom maps map #'@importFrom graphics box image layout mtext par plot.new #'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off @@ -142,7 +141,6 @@ PlotCombinedMap <- function(maps, lon, lat, # If there is any filenames to store the graphics, process them # to select the right device if (!is.null(fileout)) { - .SelectDevice <- utils::getFromNamespace(".SelectDevice", "s2dv") deviceInfo <- .SelectDevice(fileout = fileout, width = width, height = height, units = size_units, res = res) saveToFile <- deviceInfo$fun diff --git a/R/PlotTriangles4Categories.R b/R/PlotTriangles4Categories.R index 8cf775a..463830e 100644 --- a/R/PlotTriangles4Categories.R +++ b/R/PlotTriangles4Categories.R @@ -75,7 +75,6 @@ #'@importFrom grDevices dev.new dev.off dev.cur #'@importFrom graphics plot points polygon text title axis #'@importFrom RColorBrewer brewer.pal -#'@importFrom s2dv ColorBar #'@importFrom ClimProjDiags Subset #'@export PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, @@ -259,7 +258,7 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, if(legend){ # Colorbar par(mar=c(0,0,0,0)) - ColorBar(brks = brks, cols = cols, vertical = T, draw_ticks = T, draw_separators = T, + 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, ...) diff --git a/R/Utils.R b/R/Utils.R index 94797e7..665faf6 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -99,3 +99,143 @@ 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. +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)) + } + +} \ No newline at end of file diff --git a/man/PlotCombinedMap.Rd b/man/PlotCombinedMap.Rd new file mode 100644 index 0000000..c734f71 --- /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 +s2dv::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 ColorBar() 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/PlotMostLikelyQuantileMap.Rd b/man/PlotMostLikelyQuantileMap.Rd new file mode 100644 index 0000000..0dde63f --- /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} +} -- GitLab From 73e59a855076b1f94dff439e8e6eb37faa1b2e8a Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Wed, 4 Oct 2023 16:39:12 +0200 Subject: [PATCH 6/6] Add PlotCombinedMap and PlotMostLikelyQuantileMap from CSTools and remove s2dv unneeded dependencies --- NAMESPACE | 7 +- R/ColorBarContinuous.R | 2 +- R/PlotCombinedMap.R | 22 +- R/PlotMostLikelyQuantileMap.R | 2 +- R/PlotTriangles4Categories.R | 148 +----------- R/ShapeToMask.R | 9 +- R/Utils.R | 7 +- man/PlotCombinedMap.Rd | 8 +- man/PlotEquiMap.Rd | 408 ++++++++++++++++++++++++++++++++ man/PlotTriangles4Categories.Rd | 4 +- 10 files changed, 447 insertions(+), 170 deletions(-) create mode 100644 man/PlotEquiMap.Rd diff --git a/NAMESPACE b/NAMESPACE index 23ad0a3..a6a4232 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,17 +5,20 @@ export(ClimPalette) export(ColorBarContinuous) export(ColorBarDiscrete) export(PlotCombinedMap) +export(PlotEquiMap) export(PlotForecastPDF) export(PlotMostLikelyQuantileMap) export(PlotPDFsOLE) export(PlotRobinson) export(PlotTriangles4Categories) export(PlotWeeklyClim) -import(RColorBrewer) export(ShapeToMask) +import(RColorBrewer) import(cowplot) import(easyNCDF) import(ggplot2) +import(graphics) +import(maps) import(multiApply) import(rnaturalearth) import(scales) @@ -33,6 +36,7 @@ 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) @@ -60,3 +64,4 @@ 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 d19b6f5..b3c74df 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 index 642f4fd..871a9c8 100644 --- a/R/PlotCombinedMap.R +++ b/R/PlotCombinedMap.R @@ -59,7 +59,7 @@ #'@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 -#' s2dv::PlotEquiMap. +#' 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 @@ -81,9 +81,9 @@ #'@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 ColorBar() 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. +#' 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 @@ -357,11 +357,12 @@ PlotCombinedMap <- function(maps, lon, lat, } else { brks_norm <- vector('list', length = nmap) - range_width <- display_range[2] - display_range[1] #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))) { @@ -369,15 +370,16 @@ PlotCombinedMap <- function(maps, lon, lat, } else { res <- which(x == map_select_fun(x)) if (length(res) > 0) { - res <- res[1] + 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) - display_range[1]) / range_width - if (map_select_fun(x) == display_range[1]) { - res <- res + slightly_tune_val - } + 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 diff --git a/R/PlotMostLikelyQuantileMap.R b/R/PlotMostLikelyQuantileMap.R index 18e4a8b..a432805 100644 --- a/R/PlotMostLikelyQuantileMap.R +++ b/R/PlotMostLikelyQuantileMap.R @@ -219,4 +219,4 @@ PlotMostLikelyQuantileMap <- function(probs, lon, lat, cat_dim = 'bin', bar_titles = bar_titles, col_unknown_map = col_unknown_cat, drawleg = drawleg, ...) -} +} \ No newline at end of file diff --git a/R/PlotTriangles4Categories.R b/R/PlotTriangles4Categories.R index 463830e..501536c 100644 --- a/R/PlotTriangles4Categories.R +++ b/R/PlotTriangles4Categories.R @@ -53,8 +53,8 @@ #' 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 ColorBar() in -#' s2dv for color legend creation. +#'@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 @@ -259,8 +259,8 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, # 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, ...) + # 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", @@ -294,143 +294,3 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, # If the graphic was saved to file, close the connection with the device if (!is.null(fileout)) dev.off() } -.SelectDevice <- function(fileout, width, height, units, res) { - # This function is used in the plot functions to check the extension of the - # files where the graphics will be stored and select the right R device to - # save them. - # If the vector of filenames ('fileout') has files with different - # extensions, then it will only accept the first one, changing all the rest - # of the filenames to use that extension. - - # We extract the extension of the filenames: '.png', '.pdf', ... - ext <- regmatches(fileout, regexpr("\\.[a-zA-Z0-9]*$", fileout)) - - if (length(ext) != 0) { - # If there is an extension specified, select the correct device - ## units of width and height set to accept inches - if (ext[1] == ".png") { - saveToFile <- function(fileout) { - png(filename = fileout, width = width, height = height, res = res, units = units) - } - } else if (ext[1] == ".jpeg") { - saveToFile <- function(fileout) { - jpeg(filename = fileout, width = width, height = height, res = res, units = units) - } - } else if (ext[1] %in% c(".eps", ".ps")) { - saveToFile <- function(fileout) { - postscript(file = fileout, width = width, height = height) - } - } else if (ext[1] == ".pdf") { - saveToFile <- function(fileout) { - pdf(file = fileout, width = width, height = height) - } - } else if (ext[1] == ".svg") { - saveToFile <- function(fileout) { - svg(filename = fileout, width = width, height = height) - } - } else if (ext[1] == ".bmp") { - saveToFile <- function(fileout) { - bmp(filename = fileout, width = width, height = height, res = res, units = units) - } - } else if (ext[1] == ".tiff") { - saveToFile <- function(fileout) { - tiff(filename = fileout, width = width, height = height, res = res, units = units) - } - } else { - .warning("file extension not supported, it will be used '.eps' by default.") - ## In case there is only one filename - fileout[1] <- sub("\\.[a-zA-Z0-9]*$", ".eps", fileout[1]) - ext[1] <- ".eps" - saveToFile <- function(fileout) { - postscript(file = fileout, width = width, height = height) - } - } - # Change filenames when necessary - if (any(ext != ext[1])) { - .warning(paste0("some extensions of the filenames provided in 'fileout' are not ", ext[1],". The extensions are being converted to ", ext[1], ".")) - fileout <- sub("\\.[a-zA-Z0-9]*$", ext[1], fileout) - } - } else { - # Default filenames when there is no specification - .warning("there are no extensions specified in the filenames, default to '.eps'") - fileout <- paste0(fileout, ".eps") - saveToFile <- postscript - } - - # return the correct function with the graphical device, and the correct - # filenames - list(fun = saveToFile, files = fileout) -} - -.warning <- function(...) { - # Function to use the 'warning' R function with our custom settings - # Default: no call information, indent to 0, exdent to 3, - # collapse to \n - args <- list(...) - - ## In case we need to specify warning arguments - if (!is.null(args[["call."]])) { - call <- args[["call."]] - } else { - ## Default: don't show info about the call where the warning came up - call <- FALSE - } - if (!is.null(args[["immediate."]])) { - immediate <- args[["immediate."]] - } else { - ## Default value in warning function - immediate <- FALSE - } - if (!is.null(args[["noBreaks."]])) { - noBreaks <- args[["noBreaks."]] - } else { - ## Default value warning function - noBreaks <- FALSE - } - if (!is.null(args[["domain"]])) { - domain <- args[["domain"]] - } else { - ## Default value warning function - domain <- NULL - } - args[["call."]] <- NULL - args[["immediate."]] <- NULL - args[["noBreaks."]] <- NULL - args[["domain"]] <- NULL - - ## To modify strwrap indent and exdent arguments - if (!is.null(args[["indent"]])) { - indent <- args[["indent"]] - } else { - indent <- 0 - } - if (!is.null(args[["exdent"]])) { - exdent <- args[["exdent"]] - } else { - exdent <- 3 - } - args[["indent"]] <- NULL - args[["exdent"]] <- NULL - - ## To modify paste collapse argument - if (!is.null(args[["collapse"]])) { - collapse <- args[["collapse"]] - } else { - collapse <- "\n!" - } - args[["collapse"]] <- NULL - - ## Warning tag - if (!is.null(args[["tag"]])) { - tag <- args[["tag"]] - } else { - tag <- "! Warning: " - } - args[["tag"]] <- NULL - - warning(paste0(tag, paste(strwrap( - args, indent = indent, exdent = exdent - ), collapse = collapse)), call. = call, immediate. = immediate, - noBreaks. = noBreaks, domain = domain) -} - diff --git a/R/ShapeToMask.R b/R/ShapeToMask.R index 8c91675..7bf8de4 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 74a741c..b4874d0 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -102,6 +102,9 @@ #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, @@ -229,7 +232,7 @@ GradientCatsColorBar <- function(nmap, brks = NULL, cols = NULL, vertical = TRUE 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, + 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) @@ -238,4 +241,4 @@ GradientCatsColorBar <- function(nmap, brks = NULL, cols = NULL, vertical = TRUE return(list(brks = brks, cols = cols, col_inf = col_inf, col_sup = col_sup)) } -} \ No newline at end of file +} diff --git a/man/PlotCombinedMap.Rd b/man/PlotCombinedMap.Rd index c734f71..d9ec5ba 100644 --- a/man/PlotCombinedMap.Rd +++ b/man/PlotCombinedMap.Rd @@ -104,7 +104,7 @@ 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 -s2dv::PlotEquiMap.} +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 @@ -133,9 +133,9 @@ FALSE or:\cr \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 ColorBar() 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.} +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}.} } diff --git a/man/PlotEquiMap.Rd b/man/PlotEquiMap.Rd new file mode 100644 index 0000000..bc470e2 --- /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/PlotTriangles4Categories.Rd b/man/PlotTriangles4Categories.Rd index 8356da9..62a076e 100644 --- a/man/PlotTriangles4Categories.Rd +++ b/man/PlotTriangles4Categories.Rd @@ -104,8 +104,8 @@ corresponding device.} \item{figure.width}{a numeric value to control the width of the plot.} -\item{...}{The additional parameters to be passed to function ColorBar() in -s2dv for color legend creation.} +\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. -- GitLab