diff --git a/DESCRIPTION b/DESCRIPTION index f4636f1c8d43fe5fb6f430d1c263ac35eda868cb..ddef70fefe6b088eaa5d0b6865155245e618ce5f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -17,7 +17,11 @@ Imports: maps, mapproj, methods, - ClimProjDiags + ClimProjDiags, + sf, + ggplot2, + rnaturalearth, + cowplot Suggests: testthat License: GPL-3 diff --git a/NAMESPACE b/NAMESPACE index e7fb6d08f9645a15cd177a0932ed53fe00aff161..a92f986203957e98011e8230c1f2ef46054085be 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,7 +2,13 @@ export(ClimColors) export(ClimPalette) +export(ColorBarContinuous) export(ColorBarDiscrete) +export(PlotRobinson) +import(cowplot) +import(ggplot2) +import(rnaturalearth) +import(sf) importFrom(grDevices,col2rgb) importFrom(grDevices,colorRampPalette) importFrom(grDevices,rgb) diff --git a/R/ColorBarContinuous.R b/R/ColorBarContinuous.R new file mode 100644 index 0000000000000000000000000000000000000000..5d2e79b5f4dc15179e87f9154c99a712732fd75c --- /dev/null +++ b/R/ColorBarContinuous.R @@ -0,0 +1,592 @@ +#'Draws a Continuous Color Bar +#' +#'Generates a color bar to use as colouring function for map plots and +#'optionally draws it (horizontally or vertically) to be added to map +#'multipanels or plots. It is possible to draw triangles at the ends of the +#'colour bar to represent values that go beyond the range of interest. A +#'number of options is provided to adjust the colours and the position and +#'size of the components. The drawn colour bar spans a whole figure region +#'and is compatible with figure layouts.\cr\cr +#'The generated colour bar consists of a set of breaks that define the +#'length(brks) - 1 intervals to classify each of the values in each of the +#'grid cells of a two-dimensional field. The corresponding grid cell of a +#'given value of the field will be coloured in function of the interval it +#'belongs to.\cr\cr +#'The only mandatory parameters are 'var_limits' or 'brks' (in its second +#'format, see below). +#' +#'@param brks Can be provided in two formats: +#'\itemize{ +#' \item{A single value with the number of breaks to be generated +#' automatically, between the minimum and maximum specified in 'var_limits' +#' (both inclusive). Hence the parameter 'var_limits' is mandatory if 'brks' +#' is provided with this format. If 'bar_limits' is additionally provided, +#' values only between 'bar_limits' will be generated. The higher the value +#' of 'brks', the smoother the plot will look.} +#' \item{A vector with the actual values of the desired breaks. Values will +#' be reordered by force to ascending order. If provided in this format, no +#' other parameters are required to generate/plot the colour bar.} +#'} +#' This parameter is optional if 'var_limits' is specified. If 'brks' not +#' specified but 'cols' is specified, it will take as value length(cols) + 1. +#' If 'cols' is not specified either, 'brks' will take 21 as value. +#'@param cols Vector of length(brks) - 1 valid colour identifiers, for each +#' interval defined by the breaks. This parameter is optional and will be +#' filled in with a vector of length(brks) - 1 colours generated with the +#' function provided in 'color_fun' (\code{clim.colors} by default).\cr 'cols' +#' can have one additional colour at the beginning and/or at the end with the +#' aim to colour field values beyond the range of interest represented in the +#' colour bar. If any of these extra colours is provided, parameter +#' 'triangle_ends' becomes mandatory in order to disambiguate which of the +#' ends the colours have been provided for. +#'@param vertical TRUE/FALSE for vertical/horizontal colour bar +#' (disregarded if plot = FALSE). +#'@param subsampleg The first of each subsampleg breaks will be ticked on the +#' colorbar. Takes by default an approximation of a value that yields a +#' readable tick arrangement (extreme breaks always ticked). If set to 0 or +#' lower, no labels are drawn. See the code of the function for details or +#' use 'extra_labels' for customized tick arrangements. +#'@param bar_limits Vector of two numeric values with the extremes of the +#' range of values represented in the colour bar. If 'var_limits' go beyond +#' this interval, the drawing of triangle extremes is triggered at the +#' corresponding sides, painted in 'col_inf' and 'col_sup'. Either of them +#' can be set as NA and will then take as value the corresponding extreme in +#' 'var_limits' (hence a triangle end won't be triggered for these sides). +#' Takes as default the extremes of 'brks' if available, else the same values +#' as 'var_limits'. +#'@param var_limits Vector of two numeric values with the minimum and maximum +#' values of the field to represent. These are used to know whether to draw +#' triangle ends at the extremes of the colour bar and what colour to fill +#' them in with. If not specified, take the same value as the extremes of +#' 'brks'. Hence the parameter 'brks' is mandatory if 'var_limits' is not +#' specified. +#'@param triangle_ends Vector of two logical elements, indicating whether to +#' force the drawing of triangle ends at each of the extremes of the colour +#' bar. This choice is automatically made from the provided 'brks', +#' 'bar_limits', 'var_limits', 'col_inf' and 'col_sup', but the behaviour +#' can be manually forced to draw or not to draw the triangle ends with this +#' parameter. If 'cols' is provided, 'col_inf' and 'col_sup' will take +#' priority over 'triangle_ends' when deciding whether to draw the triangle +#' ends or not. +#'@param col_inf Colour to fill the inferior triangle end with. Useful if +#' specifying colours manually with parameter 'cols', to specify the colour +#' and to trigger the drawing of the lower extreme triangle, or if 'cols' is +#' not specified, to replace the colour automatically generated by ColorBar(). +#'@param col_sup Colour to fill the superior triangle end with. Useful if +#' specifying colours manually with parameter 'cols', to specify the colour +#' and to trigger the drawing of the upper extreme triangle, or if 'cols' is +#' not specified, to replace the colour automatically generated by ColorBar(). +#'@param color_fun Function to generate the colours of the color bar. Must +#' take an integer and must return as many colours. The returned colour vector +#' can have the attribute 'na_color', with a colour to draw NA values. This +#' parameter is set by default to ClimPalette(). +#'@param plot Logical value indicating whether to only compute its breaks and +#' colours (FALSE) or to also draw it on the current device (TRUE). +#'@param draw_ticks Whether to draw ticks for the labels along the colour bar +#' (TRUE) or not (FALSE). TRUE by default. Disregarded if 'plot = FALSE'. +#'@param draw_separators Whether to draw black lines in the borders of each of +#' the colour rectancles of the colour bar (TRUE) or not (FALSE). FALSE by +#' default. Disregarded if 'plot = FALSE'. +#'@param triangle_ends_scale Scale factor for the drawn triangle ends of the +#' colour bar, if drawn at all. Takes 1 by default (rectangle triangle +#' proportional to the thickness of the colour bar). Disregarded if +#' 'plot = FALSE'. +#'@param extra_labels Numeric vector of extra labels to draw along axis of +#' the colour bar. The number of provided decimals will be conserved. +#' Disregarded if 'plot = FALSE'. +#'@param title Title to draw on top of the colour bar, most commonly with the +#' units of the represented field in the neighbour figures. Empty by default. +#'@param title_scale Scale factor for the 'title' of the colour bar. +#' Takes 1 by default. +#'@param label_scale Scale factor for the labels of the colour bar. +#' Takes 1 by default. +#'@param tick_scale Scale factor for the length of the ticks of the labels +#' along the colour bar. Takes 1 by default. +#'@param extra_margin Extra margins to be added around the colour bar, +#' in the format c(y1, x1, y2, x2). The units are margin lines. Takes +#' rep(0, 4) by default. +#'@param label_digits Number of significant digits to be displayed in the +#' labels of the colour bar, usually to avoid too many decimal digits +#' overflowing the figure region. This does not have effect over the labels +#' provided in 'extra_labels'. Takes 4 by default. +#'@param ... Arguments to be passed to the method. Only accepts the following +#' graphical parameters:\cr adj ann ask bg bty cex.lab cex.main cex.sub cin +#' col.axis col.lab col.main col.sub cra crt csi cxy err family fg fig fin +#' font font.axis font.lab font.main font.sub lend lheight ljoin lmitre lty +#' lwd mai mex mfcol mfrow mfg mkh oma omd omi page pch pin plt pty smo srt +#' tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog.\cr For more +#' information about the parameters see `par`. +#' +#'@return +#'\item{brks}{ +#' Breaks used for splitting the range in intervals. +#'} +#'\item{cols}{ +#' Colours generated for each of the length(brks) - 1 intervals. +#' 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). +#'} +#' +#'@examples +#'cols <- c("dodgerblue4", "dodgerblue1", "forestgreen", "yellowgreen", "white", +#' "white", "yellow", "orange", "red", "saddlebrown") +#'lims <- seq(-1, 1, 0.2) +#'ColorBarContinuous(lims, cols) +#'@importFrom grDevices col2rgb rgb +#'@export +ColorBarContinuous <- function(brks = NULL, cols = NULL, vertical = TRUE, + subsampleg = NULL, bar_limits = NULL, var_limits = NULL, + triangle_ends = NULL, col_inf = NULL, col_sup = NULL, + color_fun = ClimPalette(), plot = TRUE, + draw_ticks = TRUE, draw_separators = FALSE, + triangle_ends_scale = 1, extra_labels = NULL, + title = NULL, title_scale = 1, + label_scale = 1, tick_scale = 1, + extra_margin = rep(0, 4), label_digits = 4, ...) { + # Required checks + if ((is.null(brks) || length(brks) < 2) && is.null(bar_limits) && is.null(var_limits)) { + stop("At least one of 'brks' with the desired breaks, 'bar_limits' or ", + "'var_limits' must be provided to generate the colour bar.") + } + + # Check brks + if (!is.null(brks)) { + if (!is.numeric(brks)) { + stop("Parameter 'brks' must be numeric if specified.") + } else if (length(brks) > 1) { + reorder <- sort(brks, index.return = TRUE) + if (!is.null(cols)) { + cols <- cols[reorder$ix[which(reorder$ix <= length(cols))]] + } + brks <- reorder$x + } + } + + # Check bar_limits + if (!is.null(bar_limits)) { + if (!(all(is.na(bar_limits) | is.numeric(bar_limits)) && (length(bar_limits) == 2))) { + stop("Parameter 'bar_limits' must be a vector of two numeric elements or NAs.") + } + } + + # Check var_limits + if (!is.null(var_limits)) { + if (!(is.numeric(var_limits) && (length(var_limits) == 2))) { + stop("Parameter 'var_limits' must be a numeric vector of length 2.") + } else if (anyNA(var_limits)) { + stop("Parameter 'var_limits' must not contain NA values.") + } else if (any(is.infinite(var_limits))) { + stop("Parameter 'var_limits' must not contain infinite values.") + } + } + + # Check cols + if (!is.null(cols)) { + if (!is.character(cols)) { + stop("Parameter 'cols' must be a vector of character strings.") + } else if (any(!sapply(cols, .IsColor))) { + stop("Parameter 'cols' must contain valid colour identifiers.") + } + } + + # Check color_fun + if (!is.function(color_fun)) { + stop("Parameter 'color_fun' must be a colour-generator function.") + } + + # Check integrity among brks, bar_limits and var_limits + if (is.null(brks) || (length(brks) < 2)) { + if (is.null(brks)) { + if (is.null(cols)) { + brks <- 21 + } else { + brks <- length(cols) + 1 + } + } + if (is.null(bar_limits) || anyNA(bar_limits)) { + # var_limits is defined + if (is.null(bar_limits)) { + bar_limits <- c(NA, NA) + } + half_width <- 0.5 * (var_limits[2] - var_limits[1]) / (brks - 1) + bar_limits[which(is.na(bar_limits))] <- c(var_limits[1] - half_width, var_limits[2] + half_width)[which(is.na(bar_limits))] + brks <- seq(bar_limits[1], bar_limits[2], length.out = brks) + } else if (is.null(var_limits)) { + # bar_limits is defined + var_limits <- bar_limits + half_width <- 0.5 * (var_limits[2] - var_limits[1]) / (brks - 1) + brks <- seq(bar_limits[1], bar_limits[2], length.out = brks) + var_limits[1] <- var_limits[1] + half_width / 50 + } else { + # both bar_limits and var_limits are defined + brks <- seq(bar_limits[1], bar_limits[2], length.out = brks) + } + } else if (is.null(bar_limits)) { + if (is.null(var_limits)) { + # brks is defined + bar_limits <- c(head(brks, 1), tail(brks, 1)) + var_limits <- bar_limits + half_width <- 0.5 * (var_limits[2] - var_limits[1]) / (length(brks) - 1) + var_limits[1] <- var_limits[1] + half_width / 50 + } else { + # brks and var_limits are defined + bar_limits <- c(head(brks, 1), tail(brks, 1)) + } + } else { + # brks and bar_limits are defined + # or + # brks, bar_limits and var_limits are defined + if (head(brks, 1) != bar_limits[1] || tail(brks, 1) != bar_limits[2]) { + stop("Parameters 'brks' and 'bar_limits' are inconsistent.") + } + } + + # Check col_inf + if (!is.null(col_inf)) { + if (!.IsColor(col_inf)) { + stop("Parameter 'col_inf' must be a valid colour identifier.") + } + } + + # Check col_sup + if (!is.null(col_sup)) { + if (!.IsColor(col_sup)) { + stop("Parameter 'col_sup' must be a valid colour identifier.") + } + } + + # Check triangle_ends + if (!is.null(triangle_ends) && (!is.logical(triangle_ends) || length(triangle_ends) != 2)) { + stop("Parameter 'triangle_ends' must be a logical vector with two elements.") + } + teflc <- triangle_ends_from_limit_cols <- c(!is.null(col_inf), !is.null(col_sup)) + if (is.null(triangle_ends) && is.null(col_inf) && is.null(col_sup)) { + triangle_ends <- c(FALSE, FALSE) + if (bar_limits[1] >= var_limits[1]) { + triangle_ends[1] <- TRUE + } + if (bar_limits[2] < var_limits[2]) { + triangle_ends[2] <- TRUE + } + } else if (!is.null(triangle_ends) && is.null(col_inf) && is.null(col_sup)) { + triangle_ends <- triangle_ends + } else if (is.null(triangle_ends) && (!is.null(col_inf) || !is.null(col_sup))) { + triangle_ends <- teflc + } else if (any(teflc != triangle_ends)) { + if (!is.null(brks) && length(brks) > 1 && !is.null(cols) && length(cols) >= length(brks)) { + triangle_ends <- teflc + } else if (!is.null(cols)) { + triangle_ends <- teflc + } else { + triangle_ends <- triangle_ends + } + } + if (plot) { + 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 ", + "disabled. These will be painted in the colour for NA values.") + } + if ((bar_limits[2] < var_limits[2]) && !triangle_ends[2]) { + .warning("There are variable values greater than the higher limit ", + "of the colour bar and the higher triangle end has been ", + "disabled. These will be painted in the colour for NA values.") + } + } + + # Generate colours if needed + if (is.null(cols)) { + cols <- color_fun(length(brks) - 1 + sum(triangle_ends)) + attr_bk <- attributes(cols) + if (triangle_ends[1]) { + if (is.null(col_inf)) col_inf <- head(cols, 1) + cols <- cols[-1] + } + if (triangle_ends[2]) { + if (is.null(col_sup)) col_sup <- tail(cols, 1) + cols <- cols[-length(cols)] + } + attributes(cols) <- attr_bk + } else if ((length(cols) != (length(brks) - 1))) { + stop("Incorrect number of 'brks' and 'cols'. There must be one more break than the number of colours.") + } + + # Check vertical + if (!is.logical(vertical)) { + stop("Parameter 'vertical' must be TRUE or FALSE.") + } + + # Check extra_labels + if (is.null(extra_labels)) { + extra_labels <- numeric(0) + } + if (!is.numeric(extra_labels)) { + stop("Parameter 'extra_labels' must be numeric.") + } else { + if (any(extra_labels > bar_limits[2]) || any(extra_labels < bar_limits[1])) { + stop("Parameter 'extra_labels' must not contain ticks beyond the color bar limits.") + } + } + extra_labels <- sort(extra_labels) + + # Check subsampleg + primes <- function(x) { + # Courtesy of Chase. See http://stackoverflow.com/questions/6424856/r-function-for-returning-all-factors + x <- as.integer(x) + div <- seq_len(abs(x)) + factors <- div[x %% div == 0L] + factors <- list(neg = -factors, pos = factors) + return(factors) + } + remove_final_tick <- FALSE + added_final_tick <- TRUE + if (is.null(subsampleg)) { + subsampleg <- 1 + while (length(brks) / subsampleg > 15 - 1) { + next_factor <- primes((length(brks) - 1) / subsampleg)$pos + next_factor <- next_factor[length(next_factor) - ifelse(length(next_factor) > 2, 1, 0)] + subsampleg <- subsampleg * next_factor + } + if (subsampleg > (length(brks) - 1) / 4) { + subsampleg <- max(1, round(length(brks) / 4)) + extra_labels <- c(extra_labels, bar_limits[2]) + added_final_tick <- TRUE + if ((length(brks) - 1) %% subsampleg < (length(brks) - 1) / 4 / 2) { + remove_final_tick <- TRUE + } + } + } else if (!is.numeric(subsampleg)) { + stop("Parameter 'subsampleg' must be numeric.") + } + subsampleg <- round(subsampleg) + draw_labels <- TRUE + if ((subsampleg) < 1) { + draw_labels <- FALSE + } + + # Check plot + if (!is.logical(plot)) { + stop("Parameter 'plot' must be logical.") + } + + # Check draw_separators + if (!is.logical(draw_separators)) { + stop("Parameter 'draw_separators' must be logical.") + } + + # Check triangle_ends_scale + if (!is.numeric(triangle_ends_scale)) { + stop("Parameter 'triangle_ends_scale' must be numeric.") + } + + # Check draw_ticks + if (!is.logical(draw_ticks)) { + stop("Parameter 'draw_ticks' must be logical.") + } + + # Check title + if (is.null(title)) { + title <- '' + } + if (!is.character(title)) { + stop("Parameter 'title' must be a character string.") + } + + # Check title_scale + if (!is.numeric(title_scale)) { + stop("Parameter 'title_scale' must be numeric.") + } + + # Check label_scale + if (!is.numeric(label_scale)) { + stop("Parameter 'label_scale' must be numeric.") + } + + # Check tick_scale + if (!is.numeric(tick_scale)) { + stop("Parameter 'tick_scale' must be numeric.") + } + + # Check extra_margin + if (!is.numeric(extra_margin) || length(extra_margin) != 4) { + stop("Parameter 'extra_margin' must be a numeric vector of length 4.") + } + + # Check label_digits + if (!is.numeric(label_digits)) { + stop("Parameter 'label_digits' must be numeric.") + } + label_digits <- round(label_digits) + + # Process the user graphical parameters that may be passed in the call + ## Graphical parameters to exclude + excludedArgs <- c("cex", "cex.axis", "col", "lab", "las", "mar", "mgp", "new", "ps") + userArgs <- .FilterUserGraphicArgs(excludedArgs, ...) + + # + # Plotting colorbar + # ~~~~~~~~~~~~~~~~~~~ + # + if (plot) { + pars_to_save <- c('mar', 'cex', names(userArgs), 'mai', 'mgp', 'las', 'xpd') + saved_pars <- par(pars_to_save) + par(mar = c(0, 0, 0, 0), cex = 1) + image(1, 1, t(t(1)), col = rgb(0, 0, 0, 0), axes = FALSE, xlab = '', ylab = '') + # Get the availale space + figure_size <- par('fin') + cs <- par('csi') + # This allows us to assume we always want to plot horizontally + if (vertical) { + figure_size <- rev(figure_size) + } +# pannel_to_redraw <- par('mfg') +# .SwitchToFigure(pannel_to_redraw[1], pannel_to_redraw[2]) + # Load the user parameters + par(new = TRUE) + par(userArgs) + # Set up color bar plot region + margins <- c(0.0, 0, 0.0, 0) + cex_title <- 1 * title_scale + cex_labels <- 0.9 * label_scale + cex_ticks <- -0.3 * tick_scale + spaceticklab <- max(-cex_ticks, 0) + if (vertical) { + margins[1] <- margins[1] + (1.2 * cex_labels * 3 + spaceticklab) * cs + margins <- margins + extra_margin[c(4, 1:3)] * cs + } else { + margins[1] <- margins[1] + (1.2 * cex_labels * 1 + spaceticklab) * cs + margins <- margins + extra_margin * cs + } + if (title != '') { + margins[3] <- margins[3] + (1.0 * cex_title) * cs + } + margins[3] <- margins[3] + sqrt(figure_size[2] / (margins[1] + margins[3])) * + figure_size[2] / 6 * ifelse(title != '', 0.5, 0.8) + # Set side margins + margins[2] <- margins[2] + figure_size[1] / 16 + margins[4] <- margins[4] + figure_size[1] / 16 + triangle_ends_prop <- 1 / 32 * triangle_ends_scale + triangle_ends_cex <- triangle_ends_prop * figure_size[2] + if (triangle_ends[1]) { + margins[2] <- margins[2] + triangle_ends_cex + } + if (triangle_ends[2]) { + margins[4] <- margins[4] + triangle_ends_cex + } + ncols <- length(cols) + # Set up the points of triangles + # Compute the proportion of horiz. space occupied by one plot unit + prop_unit <- (1 - (margins[2] + margins[4]) / figure_size[1]) / ncols + # Convert triangle height to plot inits + triangle_height <- triangle_ends_prop / prop_unit + left_triangle <- list(x = c(1, 1 - triangle_height, 1) - 0.5, + y = c(1.4, 1, 0.6)) + right_triangle <- list(x = c(ncols, ncols + triangle_height, ncols) + 0.5, + y = c(1.4, 1, 0.6)) + # Draw the color squares and title + if (vertical) { + par(mai = c(margins[2:4], margins[1]), + mgp = c(0, spaceticklab + 0.2, 0), las = 1) + d <- 4 + image(1, 1:ncols, t(1:ncols), axes = FALSE, col = cols, + xlab = '', ylab = '') + title(ylab = title, line = cex_title * (0.2 + 0.1), cex.lab = cex_title) + # Draw top and bottom border lines + lines(c(0.6, 0.6), c(1 - 0.5, ncols + 0.5)) + lines(c(1.4, 1.4), c(1 - 0.5, ncols + 0.5)) + # Rotate triangles + names(left_triangle) <- rev(names(left_triangle)) + names(right_triangle) <- rev(names(right_triangle)) + } else { + # The term - cex_labels / 4 * (3 / cex_labels - 1) was found by + # try and error + par(mai = margins, + mgp = c(0, cex_labels / 2 + spaceticklab + - cex_labels / 4 * (3 / cex_labels - 1), 0), + las = 1) + d <- 1 + image(1:ncols, 1, t(t(1:ncols)), axes = FALSE, col = cols, + xlab = '', ylab = '') + title(title, line = cex_title * (0.2 + 0.1), cex.main = cex_title) + # Draw top and bottom border lines + lines(c(1 - 0.5, ncols + 0.5), c(0.6, 0.6)) + lines(c(1 - 0.5, ncols + 0.5), c(1.4, 1.4)) + tick_length <- -0.4 + } + # Draw the triangles + par(xpd = TRUE) + if (triangle_ends[1]) { + # Draw left triangle + polygon(left_triangle$x, left_triangle$y, col = col_inf, border = NA) + lines(left_triangle$x, left_triangle$y) + } + if (triangle_ends[2]) { + # Draw right triangle + polygon(right_triangle$x, right_triangle$y, col = col_sup, border = NA) + lines(right_triangle$x, right_triangle$y) + } + par(xpd = FALSE) + + # Put the separators + if (vertical) { + if (draw_separators) { + for (i in 1:(ncols - 1)) { + lines(c(0.6, 1.4), c(i, i) + 0.5) + } + } + if (draw_separators || is.null(col_inf)) { + lines(c(0.6, 1.4), c(0.5, 0.5)) + } + if (draw_separators || is.null(col_sup)) { + lines(c(0.6, 1.4), c(ncols + 0.5, ncols + 0.5)) + } + } else { + if (draw_separators) { + for (i in 1:(ncols - 1)) { + lines(c(i, i) + 0.5, c(0.6, 1.4)) + } + } + if (draw_separators || is.null(col_inf)) { + lines(c(0.5, 0.5), c(0.6, 1.4)) + } + if (draw_separators || is.null(col_sup)) { + lines(c(ncols + 0.5, ncols + 0.5), c(0.6, 1.4)) + } + } + # Put the ticks + plot_range <- length(brks) - 1 + var_range <- tail(brks, 1) - head(brks, 1) + extra_labels_at <- ((extra_labels - head(brks, 1)) / var_range) * plot_range + 0.5 + at <- seq(1, length(brks), subsampleg) + labels <- brks[at] + # Getting rid of next-to-last tick if too close to last one + if (remove_final_tick) { + at <- at[-length(at)] + labels <- labels[-length(labels)] + } + labels <- signif(labels, label_digits) + if (added_final_tick) { + extra_labels[length(extra_labels)] <- signif(tail(extra_labels, 1), label_digits) + } + at <- at - 0.5 + at <- c(at, extra_labels_at) + labels <- c(labels, extra_labels) + tick_reorder <- sort(at, index.return = TRUE) + at <- tick_reorder$x + if (draw_labels) { + labels <- labels[tick_reorder$ix] + } else { + labels <- FALSE + } + axis(d, at = at, tick = draw_ticks, labels = labels, cex.axis = cex_labels, tcl = cex_ticks) + par(saved_pars) + } + invisible(list(brks = brks, cols = cols, col_inf = col_inf, col_sup = col_sup)) +} diff --git a/R/PlotRobinson.R b/R/PlotRobinson.R new file mode 100644 index 0000000000000000000000000000000000000000..dc9b36d4619fff23241fae6e8df1b20b430bae6c --- /dev/null +++ b/R/PlotRobinson.R @@ -0,0 +1,552 @@ +#'Plot map in Robinson or other projections +#' +#'Transform a regular grid longitude-latitude data to a different projection and +#'plot the map. The target projection must be a valid CRS string, preferrably be +#'EPSG or ESRI code; check \link[sf]{st_crs} for more explanation. This function +#'is mainly tested for Robinson projection (ESRI:54030), but it can work with +#'other projection types in theory.\n +#'The map can be plotted by points or polygon. A legend can be plotted as either +#'a color bar or a discrete ggplot legend. Dots can be drawn on top of the data, +#'which can be used for significance test. A mask can be added to not plot the +#'data specified. A number of options is provided to adjust aesthetics, like +#'position, size, colors, etc. +#' +#'@param data A numeric array with longitude and latitude dimensions. The grid +#' should be regular grid. It can contain NA values. +#'@param lon A numeric vector of longitude locations of the cell centers of the +#' grid of 'data'. Expected to be regularly spaced, within the range of either +#' [-180, 180] or [0, 360]. +#'@param lat A numeric vector of latitude locations of the cell centers of the +#' grid of 'data'. Expected to be regularly spaced, within the range [-90, 90] +#' of ascending or descending order. +#'@param lon_dim A character string indicating the longitude dimension name in +#' 'data'. If it is NULL, the function tries to find the name in +#' \code{s2dv:::.KnownLonNames}. The default value is NULL. +#'@param lat_dim A character string indicating the latitude dimension name in +#' 'data'. If it is NULL, the function tries to find the name in +#' \code{s2dv:::.KnownLatNames}. The default value is NULL. +#'@param target_proj A character string indicating the target projection. It +#' should be a valid crs string. The default projection is Robinson +#' (ESRI:54030). Note that the character string may work differently depending +#' on PROJ and GDAL module version. +#'@param legend A character string indicating the legend style. It can be 's2dv' +#' (color bar by \code{ColorBarContinuous()}), 'ggplot2' (discrete legend by +#' ggplot2), or NULL (no legend), +#'@param style A character string indicating the plotting style. It can be +#' 'point' or 'polygon'. The default value is 'point'. Note that 'polygon' may +#' be time- and memory-consuming for global or high-resolution data. +#'@param dots An array with the same dimensions as 'data' of [0, 1] or logical +#' indicating the grids to plot dots. The value 0 or FALSE is the point to be +#' dotted. +#'@param mask An array with the same dimensions as 'data' of [0, 1] or logical +#' indicating the grids to not plot data. The value 0 or FALSE is the point not +#' to be plotted. +#'@param brks,cols,bar_limits,triangle_ends Usually only providing 'brks' is +#' enough to generate the desired color bar. These parameters allow to +#' define n breaks that define n - 1 intervals to classify each of the values +#' in 'data'. The corresponding grid cell of a given value in 'data' will be +#' colored in function of the interval it belongs to. These parameters are +#' sent to \code{ColorBarContinuous()} to generate the breaks and colours. +#' Additional colors 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 ?ColorBarContinuous for +#' a full explanation. +#'@param col_inf,col_sup,colNA Colour identifiers to color the values that +#' excess the extremes of the color bar and to color NAs, 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'. 'col_inf' and 'col_sup' will take the value of 'colNA' if not +#' specified. See ?ColorBarContinuous for a full explanation. +#'@param color_fun,bar_extra_margin Set of +#' parameters to control the visual aspect of the drawn colour bar +#' (1/3). See ?ColorBarContinuous for a full explanation. +#'@param vertical A logical value indicating the direction of colorbar if +#' parameter 'legend' is 's2dv'. The default value is TRUE. +#'@param toptitle A character string of the top title of the figure, scalable +#' with parameter 'title_size'. +#'@param caption A character string of the caption located at left-bottom of the +#' plot. +#'@param units A character string of the data units, which is the title of the +#' legend. +#'@param crop_coastlines A named numeric vector [lonmin, lonmax, latmin, latmax] +#' indicating the region to plot coastlines. Note that the longitude range +#' cannot exceed 180 degrees. +#'@param point_size A number of the size of the data points if "style = 'point'". +#' The default is 'auto' and the function tries to find the appropriate size. +#'@param title_size A number of the size of the top title. The default is 16. +#'@param dots_size A number of the size of the dots. The default is 0.5. +#'@param dots_shape A number indicating the dot shape recognized by parameter +#' 'shape' in \code{geom_point()}. +#'@param coastlines_width A number indicating the width of the coastlines. +#'@param fileout A character string of the path to save the plot. If not +#' specified (default), a graphic device will pop up. The extension should be +#' accepted by \code{ggsave()}. +#'@param width A number of the plot width, in the units specified in parameter +#' 'size_units'. The default is 8. +#'@param height A number of the plot height, in the units specified in parameter +#' 'size_units'. The default is 4. +#'@param size_units A character string of the units of the size of the device +#' (file or window) to plot in. The default is 'in' (inches). See ?ggsave and +#' ?Devices for details of the corresponding device. +#'@param res Resolution of the device (file or window) to plot in. The default +#' value is 300. See ?ggsave 'dpi' and ?Devices for details of the +#' corresponding device. +#' +#'@return A map plot with speficied projection, either in pop-up window or a +#' saved file. +#' +#'@examples +#'data <- array(rep(seq(-10, 10, length.out = 181), 360) + rnorm(360), +#' dim = c(lat = 181, lon = 360)) +#'dots <- data +#'dots[which(dots < 4 & dots > -4)] <- 0 +#'dots[which(dots != 0)] <- 1 +#'PlotRobinson(data, lon = 0:359, lat = -90:90, dots = dots, +#' brks = seq(-10, 10, length.out = 11), +#' toptitle = 'synthetic example', vertical = F, +#' caption = 'Robinson Global\ns2dv::PlotRobinson example', +#' bar_extra_margin = c(0, 1, 0, 1), width = 8, height = 6) +#'PlotRobinson(data, lon = 0:359, lat = -90:90, mask = dots, legend = 'ggplot2', +#' target_proj = "+proj=moll", brks = seq(-10, 10, length.out = 11), +#' color_fun = ClimPalette("purpleorange"), colNA = 'green', +#' toptitle = 'synthetic example', +#' caption = 'Mollweide Global\ns2dv::PlotRobinson example', +#' width = 8, height = 6) +#' +#'@import sf ggplot2 rnaturalearth cowplot +#'@export +PlotRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL, + target_proj = 54030, legend = 's2dv', style = 'point', + dots = NULL, mask = NULL, brks = NULL, cols = NULL, bar_limits = NULL, + triangle_ends = NULL, col_inf = NULL, col_sup = NULL, colNA = NULL, + color_fun = ClimPalette(), bar_extra_margin = rep(0, 4), vertical = TRUE, + toptitle = NULL, caption = NULL, units = NULL, crop_coastlines = NULL, + point_size = "auto", title_size = 16, dots_size = 0.5, + dots_shape = 47, coastlines_width = 0.3, + fileout = NULL, width = 8, height = 4, size_units = "in", + res = 300) { + + # Sanity check + # data + data <- drop(data) + if (length(dim(data)) != 2) { + stop("Parameter 'data' must have two dimensions.") + } + dims <- dim(data) + # lon, lon_dim + if (is.null(lon_dim)) { + lon_dim <- names(dims)[names(dims) %in% .KnownLonNames()] + if (identical(lon_dim, character(0))) { + stop("Cannot find known longitude name in data dimension. Please define parameter 'lon_dim'.") + } + } + if (is.unsorted(lon)) { + warning("Parameter 'lon' should be sorted to guarantee the correct result.") + } + # lat, lat_dim + if (is.null(lat_dim)) { + lat_dim <- names(dims)[names(dims) %in% .KnownLatNames()] + if (identical(lat_dim, character(0))) { + stop("Cannot find known latitude name in data dimension. Please define parameter 'lat_dim'.") + } + } + if (!all(names(dims) %in% c(lat_dim, lon_dim))) { + stop("Dimensions names in paramter 'data' should match 'lat_dim' and 'lon_dim.") + } + if (length(lon) != dims[lon_dim]) { + stop("Length of parameter 'lon' should match longitude dimension in 'data'.") + } + if (length(lat) != dims[lat_dim]) { + stop("Length of parameter 'lat' should match latitude dimension in 'data'.") + } + + # Reorder data + data <- aperm(data, match(names(dim(data)), c(lon_dim, lat_dim))) + + # Make lat always from 90 to -90 + sort_lat <- FALSE + if (!is.unsorted(lat)) { + lat <- rev(lat) + data <- ClimProjDiags::Subset(data, along = lat_dim, indices = seq(length(lat), 1, -1)) + sort_lat <- TRUE + } + + # original_proj: it can only be regular grid now + original_proj <- st_crs(4326) + # tartget_proj + if (is.null(target_proj)) { + stop("Parameter 'target_proj' cannot be NULL.") + } else { + target_proj_tmp <- st_crs(target_proj) + if (is.na(target_proj_tmp)) { + warning(paste0("Try ESRI code: ESRI:", target_proj)) + target_proj <- st_crs(paste0("ESRI:", target_proj)) + } else { + target_proj <- target_proj_tmp + } + } + + # legend + if (!is.null(legend) && (!legend %in% c('s2dv', 'ggplot2'))) { + stop("Parameter 'legend' must be NULL, ggplot2 or s2dv.") + } + # style + if (!style %in% c('point', 'polygon') || length(style) != 1) { + stop("Parameter 'style' must be 'point' or 'polygon'.") + } + if (style == 'polygon') { + # polygon is slow for global map (and may be wrong) Confirm if users want to proceed + if ((abs(diff(range(lon))) > 350 & abs(diff(range(lat))) > 175) | + (prod(dim(data)) >= (180 * 360))) { + if (!isTRUE(utils::askYesNo("The region seems to be global and style 'polygon' is chosen. It may be time- and memory-consuming to plot the map. Are you sure that you want to continue?"))) { + return(invisible()) + } + } + } + # dots + if (!is.null(dots)) { + dots <- drop(dots) + if (!is.array(dots) || any(!names(dim(dots)) %in% c(lon_dim, lat_dim))) { + stop("Parameter 'dots' must have two dimensions named as longitude and latitude dimensions in 'data'.") + } else { + dots <- aperm(dots, match(names(dim(dots)), c(lon_dim, lat_dim))) + } + if (!identical(dim(dots), dim(data))) { + stop("Parameter 'dots' must have the same dimensions as 'data'.") + } else if (is.numeric(dots)) { + if (all(dots %in% c(0, 1))) { + dots <- array(as.logical(dots), dim = dim(dots)) + } else { + stop("Parameter 'dots' must have only TRUE/FALSE or 0/1.") + } + } else if (is.logical(dots)) { + if (!all(dots %in% c(T, F))) { + stop("Parameter 'dots' must have only TRUE/FALSE or 0/1.") + } + } else { + stop("Parameter 'dots' must be a logical or numerical array.") + } + } + # mask + if (!is.null(mask)) { + mask <- drop(mask) + if (!is.array(mask) || any(!names(dim(mask)) %in% c(lon_dim, lat_dim))) { + stop("Parameter 'mask' must have two dimensions named as longitude and latitude dimensions in 'data'.") + } else { + mask <- aperm(mask, match(names(dim(mask)), c(lon_dim, lat_dim))) + } + if (!identical(dim(mask), dim(data))) { + stop("Parameter 'mask' must have the same dimensions as 'data'.") + } else if (is.numeric(mask)) { + if (all(mask %in% c(0, 1))) { + mask <- array(as.logical(mask), dim = dim(mask)) + } else { + stop("Parameter 'mask' must have only TRUE/FALSE or 0/1.") + } + } else if (is.logical(mask)) { + if (!all(mask %in% c(T, F))) { + stop("Parameter 'mask' must have only TRUE/FALSE or 0/1.") + } + } else { + stop("Parameter 'mask' must be a logical or numerical array.") + } + } + + # Color bar + ## Check: brks, cols, bar_limits, color_fun, bar_extra_margin, units + ## Build: brks, cols, bar_limits, col_inf, col_sup + var_limits <- c(min(data, na.rm = TRUE), max(data, na.rm = TRUE)) + colorbar <- ColorBarContinuous(brks = brks, cols = cols, vertical = vertical, subsampleg = NULL, + bar_limits = bar_limits, var_limits = var_limits, triangle_ends = triangle_ends, + col_inf = col_inf, col_sup = col_sup, color_fun = color_fun, + plot = FALSE, draw_ticks = TRUE, draw_separators = FALSE, + triangle_ends_scale = 1, extra_labels = NULL, + title = units, title_scale = 1, # units_scale + label_scale = 1, tick_scale = 1, #bar_tick_scale + extra_margin = bar_extra_margin, label_digits = 4) + brks <- colorbar$brks + cols <- colorbar$cols + col_inf <- colorbar$col_inf + col_sup <- colorbar$col_sup + bar_limits <- c(head(brks, 1), tail(brks, 1)) + # colNA + if (is.null(colNA)) { + if ('na_color' %in% names(attributes(cols))) { + colNA <- attr(cols, 'na_color') + if (!.IsColor(colNA)) { + stop("The 'na_color' provided as attribute of the colour vector must be a valid colour identifier.") + } + } else { + colNA <- 'pink' + } + } else if (!.IsColor(colNA)) { + stop("Parameter 'colNA' must be a valid colour identifier.") + } + # toptitle + if (!is.null(toptitle) && !is.character(toptitle)) { + stop("Parameter 'toptitle' must be a character string.") + } + # caption + if (!is.null(caption) && !is.character(caption)) { + stop("Parameter 'caption' must be a character string.") + } + # crop_coastlines + if (!is.null(crop_coastlines)) { + # if crop_coastlines doesn't have name, [lonmin, lonmax, latmin, latmax] + if (is.null(names(crop_coastlines))) { + names(crop_coastlines) <- c("lonmin", "lonmax", "latmin", "latmax") + } else if (!identical(sort(names(crop_coastlines)), sort(c("latmax", "latmin", "lonmax", "lonmin")))) { + stop("Parameter 'crop_coastlines' needs to have names 'latmax', 'latmin', 'lonmax', 'lonmin'.") + } + } + + # point_size + if (point_size == 'auto') { + # 360x181 with default setting, 0.05 + point_size <- round(0.05 * (360 * 181) / (length(lon) * length(lat)), 2) + } else if (!is.numeric(point_size)) { + stop("Parameter 'point_size' must be a number.") + } + # + +#================================================================= + + # Adapt s2dv ColorBar parameters to ggplot plot + # If legend is NULL, still tune with s2dv legend way + if (is.null(legend) || legend == 's2dv') { + # the colorbar triangle color. If it is NULL (no triangle plotted), use colNA + col_inf_image <- ifelse(is.null(col_inf), colNA, col_inf) + col_sup_image <- ifelse(is.null(col_sup), colNA, col_sup) + cols_ggplot <- c(col_inf_image, cols, col_sup_image) + + # Add triangles to brks + brks_ggplot <- brks + if (max(data, na.rm = T) > tail(brks, 1)) { + brks_ggplot <- c(brks_ggplot, max(data, na.rm = T)) + } else { + brks_ggplot <- c(brks_ggplot, tail(brks, 1) + diff(tail(brks, 2))) + } + if (min(data, na.rm = T) < brks[1]) { + brks_ggplot <- c(min(data, na.rm = T), brks_ggplot) + } else { + brks_ggplot <- c(brks[1] - diff(brks[1:2]), brks_ggplot) + } + + } else { # ggplot2 legend + brks_ggplot <- brks + cols_ggplot <- cols + } + + # Build data dataframe + lonlat_df <- data.frame(lon = rep(as.vector(lon), length(lat)), + lat = sort(rep(as.vector(lat), length(lon)), decreasing = TRUE)) + data_df <- lonlat_df %>% + dplyr::mutate(dat = as.vector(data)) + + lonlat_df_ori <- NULL + # Remove the points where mask = FALSE + if (!is.null(mask)) { + # Save original lonlat_df to plot with expected region + lonlat_df_ori <- st_as_sf(lonlat_df, coords = c("lon", "lat"), crs = original_proj) + lonlat_df_ori <- st_transform(lonlat_df_ori, crs = target_proj) + lonlat_df_ori <- as.data.frame(st_coordinates(lonlat_df_ori)) + names(lonlat_df_ori) <- c('long', 'lat') + + if (sort_lat) { + mask <- ClimProjDiags::Subset(mask, along = lat_dim, indices = seq(length(lat), 1, -1)) + } + mask_df <- data.frame(lon = rep(as.vector(lon), length(lat)), + lat = sort(rep(as.vector(lat), length(lon)), decreasing = TRUE), + mask = as.vector(mask)) + data_df <- data_df[mask_df$mask == TRUE, ] + lonlat_df <- data_df[, 1:2] + } + + #NOTE: if target_proj = "ESRI:54030", Nord3v2 has different behavior from hub and ws!! + data_df <- st_as_sf(data_df, coords = c("lon", "lat"), crs = original_proj) + data_df <- st_transform(data_df, crs = target_proj) + data_df <- data_df %>% + dplyr::mutate(long = st_coordinates(data_df)[, 1], + lat = st_coordinates(data_df)[, 2]) + + # Re-project dots + if (!is.null(dots)) { + if (sort_lat) { + dots <- ClimProjDiags::Subset(dots, along = lat_dim, indices = seq(length(lat), 1, -1)) + } + dots_df <- data.frame(lon = rep(as.vector(lon), length(lat)), + lat = sort(rep(as.vector(lat), length(lon)), decreasing = TRUE), + dot = as.vector(dots)) + + dots_df <- st_as_sf(dots_df, coords = c("lon", "lat"), crs = original_proj) + dots_df <- st_transform(dots_df, crs = target_proj) + dots_df <- dots_df %>% + dplyr::mutate(long = st_coordinates(dots_df)[, 1], + lat = st_coordinates(dots_df)[, 2]) + dots_df <- subset(dots_df, dot == FALSE) + } + + # coastlines + coastlines <- rnaturalearth::ne_coastline(scale = "medium", returnclass = "sf") + ## crop the coastlines to the desired range + if (!is.null(crop_coastlines)) { + suppressWarnings({ + coastlines <- st_crop(coastlines, + xmin = as.numeric(crop_coastlines['lonmin']), + xmax = as.numeric(crop_coastlines['lonmax']), + ymin = as.numeric(crop_coastlines['latmin']), + ymax = as.numeric(crop_coastlines['latmax'])) + }) + } + coastlines <- st_transform(coastlines, crs = target_proj) + + if (style == 'polygon') { + # Calculate polygon points from regular lat/lon + #NOTE: The original grid must be regular grid with same space + d_lon <- abs(lon[2] - lon[1]) / 2 + d_lat <- abs(lat[2] - lat[1]) / 2 + lon_poly <- lat_poly <- rep(NA, 4 * dim(lonlat_df)[1]) + for (ii in 1:dim(lonlat_df)[1]) { + lon_poly[(ii*4-3):(ii*4)] <- c(lonlat_df$lon[ii] - d_lon, lonlat_df$lon[ii] + d_lon, + lonlat_df$lon[ii] + d_lon, lonlat_df$lon[ii] - d_lon) + lat_poly[(ii*4-3):(ii*4)] <- c(lonlat_df$lat[ii] - d_lat, lonlat_df$lat[ii] - d_lat, + lonlat_df$lat[ii] + d_lat, lonlat_df$lat[ii] + d_lat) + } + # # To prevent out-of-global lon + # lon_poly[which(lon_poly >= 180)] <- 179.9 + # lon_poly[which(lon_poly < -180)] <- -180 + # To prevent out-of-global lat + lat_poly[which(lat_poly > 90)] <- 90 + lat_poly[which(lat_poly < -90)] <- -90 + + lonlat_df <- data.frame(lon = lon_poly, lat = lat_poly) + # Transfer lon/lat to projection + proj_lonlat <- st_as_sf(lonlat_df, coords = c("lon", "lat"), crs = original_proj) + #NOTE: if target_proj = "ESRI:54030", on Nord3v2, st_transform has lon and lat swapped! + proj_lonlat <- st_transform(proj_lonlat, crs = target_proj) + lonlat_df_proj <- st_coordinates(proj_lonlat) + + # Use id to create groups for each polygon + ids <- factor(paste0("id_", 1:dim(data_df)[1])) + values <- data.frame(id = ids, + value = data_df$dat) + positions <- data.frame(id = rep(ids, each = 4), + x = lonlat_df_proj[, 1], + y = lonlat_df_proj[, 2]) + datapoly <- merge(values, positions, by = "id") + datapoly <- st_as_sf(datapoly, coords = c("x", "y"), crs = target_proj) + datapoly <- datapoly %>% + dplyr::group_by(id) %>% + dplyr::summarise() %>% #NOTE: VERY SLOW if plot global + dplyr::mutate(value = values[order(values$id), ]$value) %>% + st_cast("POLYGON") %>% + st_convex_hull() # maintain outer polygen (no bowtie shape) + } + + # Plots + if (style == 'polygon') { + res_p <- ggplot(data = data_df) + #NOTE: must be data_df? + geom_sf(data = datapoly, + aes(col = cut(value, breaks = brks_ggplot, include.lowest = T), + fill = cut(value, breaks = brks_ggplot, include.lowest = T))) + } else if (style == 'point') { + res_p <- ggplot(data = data_df) + + geom_point(aes(x = long, y = lat, + col = cut(dat, breaks = brks_ggplot, include.lowest = T)), + #NOTE: These two lines make point size vary with lat + #size = point_size / (data_df$lat / min(data_df$lat))) + + #size = (sort(rep(as.vector(lat), length(lon))) / max(lat)) * point_size) + + size = point_size) + } + + if (is.null(lonlat_df_ori)) { + coord_sf_lim <- c(range(data_df$long), range(data_df$lat)) + } else { + coord_sf_lim <- c(range(lonlat_df_ori$long), range(lonlat_df_ori$lat)) + } + res_p <- res_p + + geom_sf(data = coastlines, colour ='black', size = coastlines_width) + + # Remove background grid and lat/lon label; add white background + theme_void() + theme(plot.background = element_rect(fill = 'white', colour = 'white')) + + # crop the projection + coord_sf(xlim = coord_sf_lim[1:2], ylim = coord_sf_lim[3:4], + expand = TRUE, datum = target_proj) + + if (!is.null(dots)) { + res_p <- res_p + geom_point(data = dots_df, aes(x = long, y = lat), + shape = dots_shape, size = dots_size) + #NOTE: This line makes point size vary with lat + #size = dots_size / (dots_df$lat / min(dots_df$lat))) + } + + if (identical(legend, 'ggplot2')) { + if (style == 'polygon') { + res_p <- res_p + scale_colour_manual(values = cols_ggplot, + aesthetics = c("colour", "fill"), + drop = FALSE, na.value = colNA) + + guides(fill = guide_legend(title = units, override.aes = list(size = 1)), + color = "none") + } else if (style == 'point') { + res_p <- res_p + scale_colour_manual(values = cols_ggplot, + drop = FALSE, na.value = colNA) + + guides(colour = guide_legend(title = units, override.aes = list(size = 1))) + } + + } else { # s2dv or NULL + if (style == 'polygon') { + res_p <- res_p + scale_colour_manual(values = cols_ggplot, + aesthetics = c("colour", "fill"), + drop = FALSE, na.value = colNA) + } else if (style == 'point') { + res_p <- res_p + scale_colour_manual(values = cols_ggplot, + drop = FALSE, na.value = colNA) + } + # Remove ggplot legend + res_p <- res_p + theme(legend.position = "none", plot.margin = margin(0.5, 0, 0, 0, 'cm')) + } + + if (!is.null(toptitle)) { + res_p <- res_p + ggtitle(toptitle) + + theme(plot.title = element_text(size = title_size, hjust = 0.5, vjust = 3)) + } + if (!is.null(caption)) { + res_p <- res_p + labs(caption = caption) + + theme(plot.caption = element_text(hjust = 0, vjust = 1, margin = margin(0, 0, 0, 0, 'cm'))) + } + + # s2dv legend fun to put in cowplot::plot_grid + if (identical(legend, 's2dv')) { + fun_legend <- function() { + if (vertical) { + par(mar = c(7.1, 2.2, 7.1, 3.1), mgp = c(3, 1, 0)) + } else { + par(mar = c(1.1, 1.2, 0.1, 1.1), mgp = c(3, 1, 0)) + } + ColorBarContinuous(brks = brks, cols = cols, vertical = vertical, subsampleg = NULL, + bar_limits = bar_limits, var_limits = var_limits, triangle_ends = triangle_ends, + col_inf = col_inf, col_sup = col_sup, color_fun = color_fun, + plot = TRUE, draw_ticks = TRUE, draw_separators = FALSE, + triangle_ends_scale = 1, extra_labels = NULL, + title = units, title_scale = 1, # units_scale + label_scale = 1, tick_scale = 1, #bar_tick_scale + extra_margin = bar_extra_margin, label_digits = 4) + } + if (vertical) { + res_p <- cowplot::plot_grid(res_p, fun_legend, rel_widths = c(6, 1)) + } else { + res_p <- cowplot::plot_grid(res_p, fun_legend, rel_heights = c(5, 1), ncol = 1) + } + res_p <- res_p + theme(plot.background = element_rect(fill = "white", colour = "white")) + } + + if (!is.null(fileout)) { + ext <- regmatches(fileout, regexpr("[a-zA-Z0-9]*$", fileout)) + ggsave(fileout, res_p, width = width, height = height, dpi = res, units = size_units, + device = ext) + } else { # pop-up window + dev.new(units = size_units, res = res, width = width, height = height) + res_p + } + +} + diff --git a/R/Utils.R b/R/Utils.R index f262e0a9e50261ca18aca16fa45d5f5ce1397735..94797e7da89ff9fd0be0e61fe3aa7d362b14981a 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -1,3 +1,11 @@ +.KnownLonNames <- function() { + known_lon_names <- c('lon', 'longitude', 'x', 'i', 'nav_lon') +} + +.KnownLatNames <- function() { + known_lat_names <- c('lat', 'latitude', 'y', 'j', 'nav_lat') +} + .IsColor <- function(x) { res <- try(col2rgb(x), silent = TRUE) return(!"try-error" %in% class(res)) diff --git a/man/ColorBarContinuous.Rd b/man/ColorBarContinuous.Rd new file mode 100644 index 0000000000000000000000000000000000000000..78d6aa526f7b7394c20ee288b773a9081ddfb79d --- /dev/null +++ b/man/ColorBarContinuous.Rd @@ -0,0 +1,194 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ColorBarContinuous.R +\name{ColorBarContinuous} +\alias{ColorBarContinuous} +\title{Draws a Continuous Color Bar} +\usage{ +ColorBarContinuous( + brks = NULL, + cols = NULL, + vertical = TRUE, + subsampleg = NULL, + bar_limits = NULL, + var_limits = NULL, + triangle_ends = NULL, + col_inf = NULL, + col_sup = NULL, + color_fun = ClimPalette(), + plot = TRUE, + draw_ticks = TRUE, + draw_separators = FALSE, + triangle_ends_scale = 1, + extra_labels = NULL, + title = NULL, + title_scale = 1, + label_scale = 1, + tick_scale = 1, + extra_margin = rep(0, 4), + label_digits = 4, + ... +) +} +\arguments{ +\item{brks}{Can be provided in two formats: +\itemize{ + \item{A single value with the number of breaks to be generated + automatically, between the minimum and maximum specified in 'var_limits' + (both inclusive). Hence the parameter 'var_limits' is mandatory if 'brks' + is provided with this format. If 'bar_limits' is additionally provided, + values only between 'bar_limits' will be generated. The higher the value + of 'brks', the smoother the plot will look.} + \item{A vector with the actual values of the desired breaks. Values will + be reordered by force to ascending order. If provided in this format, no + other parameters are required to generate/plot the colour bar.} +} + This parameter is optional if 'var_limits' is specified. If 'brks' not + specified but 'cols' is specified, it will take as value length(cols) + 1. + If 'cols' is not specified either, 'brks' will take 21 as value.} + +\item{cols}{Vector of length(brks) - 1 valid colour identifiers, for each +interval defined by the breaks. This parameter is optional and will be +filled in with a vector of length(brks) - 1 colours generated with the +function provided in 'color_fun' (\code{clim.colors} by default).\cr 'cols' +can have one additional colour at the beginning and/or at the end with the +aim to colour field values beyond the range of interest represented in the +colour bar. If any of these extra colours is provided, parameter +'triangle_ends' becomes mandatory in order to disambiguate which of the +ends the colours have been provided for.} + +\item{vertical}{TRUE/FALSE for vertical/horizontal colour bar +(disregarded if plot = FALSE).} + +\item{subsampleg}{The first of each subsampleg breaks will be ticked on the +colorbar. Takes by default an approximation of a value that yields a +readable tick arrangement (extreme breaks always ticked). If set to 0 or +lower, no labels are drawn. See the code of the function for details or +use 'extra_labels' for customized tick arrangements.} + +\item{bar_limits}{Vector of two numeric values with the extremes of the +range of values represented in the colour bar. If 'var_limits' go beyond +this interval, the drawing of triangle extremes is triggered at the +corresponding sides, painted in 'col_inf' and 'col_sup'. Either of them +can be set as NA and will then take as value the corresponding extreme in +'var_limits' (hence a triangle end won't be triggered for these sides). +Takes as default the extremes of 'brks' if available, else the same values +as 'var_limits'.} + +\item{var_limits}{Vector of two numeric values with the minimum and maximum +values of the field to represent. These are used to know whether to draw +triangle ends at the extremes of the colour bar and what colour to fill +them in with. If not specified, take the same value as the extremes of +'brks'. Hence the parameter 'brks' is mandatory if 'var_limits' is not +specified.} + +\item{triangle_ends}{Vector of two logical elements, indicating whether to +force the drawing of triangle ends at each of the extremes of the colour +bar. This choice is automatically made from the provided 'brks', +'bar_limits', 'var_limits', 'col_inf' and 'col_sup', but the behaviour +can be manually forced to draw or not to draw the triangle ends with this +parameter. If 'cols' is provided, 'col_inf' and 'col_sup' will take +priority over 'triangle_ends' when deciding whether to draw the triangle +ends or not.} + +\item{col_inf}{Colour to fill the inferior triangle end with. Useful if +specifying colours manually with parameter 'cols', to specify the colour +and to trigger the drawing of the lower extreme triangle, or if 'cols' is +not specified, to replace the colour automatically generated by ColorBar().} + +\item{col_sup}{Colour to fill the superior triangle end with. Useful if +specifying colours manually with parameter 'cols', to specify the colour +and to trigger the drawing of the upper extreme triangle, or if 'cols' is +not specified, to replace the colour automatically generated by ColorBar().} + +\item{color_fun}{Function to generate the colours of the color bar. Must +take an integer and must return as many colours. The returned colour vector +can have the attribute 'na_color', with a colour to draw NA values. This +parameter is set by default to ClimPalette().} + +\item{plot}{Logical value indicating whether to only compute its breaks and +colours (FALSE) or to also draw it on the current device (TRUE).} + +\item{draw_ticks}{Whether to draw ticks for the labels along the colour bar +(TRUE) or not (FALSE). TRUE by default. Disregarded if 'plot = FALSE'.} + +\item{draw_separators}{Whether to draw black lines in the borders of each of +the colour rectancles of the colour bar (TRUE) or not (FALSE). FALSE by +default. Disregarded if 'plot = FALSE'.} + +\item{triangle_ends_scale}{Scale factor for the drawn triangle ends of the +colour bar, if drawn at all. Takes 1 by default (rectangle triangle +proportional to the thickness of the colour bar). Disregarded if +'plot = FALSE'.} + +\item{extra_labels}{Numeric vector of extra labels to draw along axis of +the colour bar. The number of provided decimals will be conserved. +Disregarded if 'plot = FALSE'.} + +\item{title}{Title to draw on top of the colour bar, most commonly with the +units of the represented field in the neighbour figures. Empty by default.} + +\item{title_scale}{Scale factor for the 'title' of the colour bar. +Takes 1 by default.} + +\item{label_scale}{Scale factor for the labels of the colour bar. +Takes 1 by default.} + +\item{tick_scale}{Scale factor for the length of the ticks of the labels +along the colour bar. Takes 1 by default.} + +\item{extra_margin}{Extra margins to be added around the colour bar, +in the format c(y1, x1, y2, x2). The units are margin lines. Takes +rep(0, 4) by default.} + +\item{label_digits}{Number of significant digits to be displayed in the +labels of the colour bar, usually to avoid too many decimal digits +overflowing the figure region. This does not have effect over the labels +provided in 'extra_labels'. Takes 4 by default.} + +\item{...}{Arguments to be passed to the method. Only accepts the following + graphical parameters:\cr adj ann ask bg bty cex.lab cex.main cex.sub cin + col.axis col.lab col.main col.sub cra crt csi cxy err family fg fig fin + font font.axis font.lab font.main font.sub lend lheight ljoin lmitre lty + lwd mai mex mfcol mfrow mfg mkh oma omd omi page pch pin plt pty smo srt + tck 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 splitting the range in intervals. +} +\item{cols}{ + Colours generated for each of the length(brks) - 1 intervals. + 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{ +Generates a color bar to use as colouring function for map plots and +optionally draws it (horizontally or vertically) to be added to map +multipanels or plots. It is possible to draw triangles at the ends of the +colour bar to represent values that go beyond the range of interest. A +number of options is provided to adjust the colours and the position and +size of the components. The drawn colour bar spans a whole figure region +and is compatible with figure layouts.\cr\cr +The generated colour bar consists of a set of breaks that define the +length(brks) - 1 intervals to classify each of the values in each of the +grid cells of a two-dimensional field. The corresponding grid cell of a +given value of the field will be coloured in function of the interval it +belongs to.\cr\cr +The only mandatory parameters are 'var_limits' or 'brks' (in its second +format, see below). +} +\examples{ +cols <- c("dodgerblue4", "dodgerblue1", "forestgreen", "yellowgreen", "white", + "white", "yellow", "orange", "red", "saddlebrown") +lims <- seq(-1, 1, 0.2) +ColorBarContinuous(lims, cols) +} diff --git a/man/ColorBarDiscrete.Rd b/man/ColorBarDiscrete.Rd index 82dcc7ca8414e6396c3115a8100965cb1271790a..9ba59d9721fa9c5cdbd236b489fb49db690a9b55 100644 --- a/man/ColorBarDiscrete.Rd +++ b/man/ColorBarDiscrete.Rd @@ -11,7 +11,7 @@ ColorBarDiscrete( subsampleg = NULL, bar_limits = NULL, var_limits = NULL, - color_fun = clim.palette(), + color_fun = ClimPalette(), plot = TRUE, draw_ticks = FALSE, draw_separators = TRUE, @@ -81,7 +81,7 @@ specified.} \item{color_fun}{Function to generate the colours of the color bar. Must take an integer and must return as many colours. The returned colour vector can have the attribute 'na_color', with a colour to draw NA values. This -parameter is set by default to clim.palette().} +parameter is set by default to ClimPalette().} \item{plot}{Logical value indicating whether to only compute its breaks and colours (FALSE) or to also draw it on the current device (TRUE).} diff --git a/man/PlotRobinson.Rd b/man/PlotRobinson.Rd new file mode 100644 index 0000000000000000000000000000000000000000..cccf0f58b5bdb5a250fe3865304190659934d2ce --- /dev/null +++ b/man/PlotRobinson.Rd @@ -0,0 +1,187 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PlotRobinson.R +\name{PlotRobinson} +\alias{PlotRobinson} +\title{Plot map in Robinson or other projections} +\usage{ +PlotRobinson( + data, + lon, + lat, + lon_dim = NULL, + lat_dim = NULL, + target_proj = 54030, + legend = "s2dv", + style = "point", + dots = NULL, + mask = NULL, + brks = NULL, + cols = NULL, + bar_limits = NULL, + triangle_ends = NULL, + col_inf = NULL, + col_sup = NULL, + colNA = NULL, + color_fun = ClimPalette(), + bar_extra_margin = rep(0, 4), + vertical = TRUE, + toptitle = NULL, + caption = NULL, + units = NULL, + crop_coastlines = NULL, + point_size = "auto", + title_size = 16, + dots_size = 0.5, + dots_shape = 47, + coastlines_width = 0.3, + fileout = NULL, + width = 8, + height = 4, + size_units = "in", + res = 300 +) +} +\arguments{ +\item{data}{A numeric array with longitude and latitude dimensions. The grid +should be regular grid. It can contain NA values.} + +\item{lon}{A numeric vector of longitude locations of the cell centers of the +grid of 'data'. Expected to be regularly spaced, within the range of either +[-180, 180] or [0, 360].} + +\item{lat}{A numeric vector of latitude locations of the cell centers of the +grid of 'data'. Expected to be regularly spaced, within the range [-90, 90] +of ascending or descending order.} + +\item{lon_dim}{A character string indicating the longitude dimension name in +'data'. If it is NULL, the function tries to find the name in +\code{s2dv:::.KnownLonNames}. The default value is NULL.} + +\item{lat_dim}{A character string indicating the latitude dimension name in +'data'. If it is NULL, the function tries to find the name in +\code{s2dv:::.KnownLatNames}. The default value is NULL.} + +\item{target_proj}{A character string indicating the target projection. It +should be a valid crs string. The default projection is Robinson +(ESRI:54030). Note that the character string may work differently depending +on PROJ and GDAL module version.} + +\item{legend}{A character string indicating the legend style. It can be 's2dv' +(color bar by \code{ColorBarContinuous()}), 'ggplot2' (discrete legend by +ggplot2), or NULL (no legend),} + +\item{style}{A character string indicating the plotting style. It can be +'point' or 'polygon'. The default value is 'point'. Note that 'polygon' may +be time- and memory-consuming for global or high-resolution data.} + +\item{dots}{An array with the same dimensions as 'data' of [0, 1] or logical +indicating the grids to plot dots. The value 0 or FALSE is the point to be +dotted.} + +\item{mask}{An array with the same dimensions as 'data' of [0, 1] or logical +indicating the grids to not plot data. The value 0 or FALSE is the point not +to be plotted.} + +\item{brks, cols, bar_limits, triangle_ends}{Usually only providing 'brks' is +enough to generate the desired color bar. These parameters allow to +define n breaks that define n - 1 intervals to classify each of the values +in 'data'. The corresponding grid cell of a given value in 'data' will be +colored in function of the interval it belongs to. These parameters are +sent to \code{ColorBarContinuous()} to generate the breaks and colours. +Additional colors 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 ?ColorBarContinuous for +a full explanation.} + +\item{col_inf, col_sup, colNA}{Colour identifiers to color the values that +excess the extremes of the color bar and to color NAs, 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'. 'col_inf' and 'col_sup' will take the value of 'colNA' if not +specified. See ?ColorBarContinuous for a full explanation.} + +\item{color_fun, bar_extra_margin}{Set of +parameters to control the visual aspect of the drawn colour bar +(1/3). See ?ColorBarContinuous for a full explanation.} + +\item{vertical}{A logical value indicating the direction of colorbar if +parameter 'legend' is 's2dv'. The default value is TRUE.} + +\item{toptitle}{A character string of the top title of the figure, scalable +with parameter 'title_size'.} + +\item{caption}{A character string of the caption located at left-bottom of the +plot.} + +\item{units}{A character string of the data units, which is the title of the +legend.} + +\item{crop_coastlines}{A named numeric vector [lonmin, lonmax, latmin, latmax] +indicating the region to plot coastlines. Note that the longitude range +cannot exceed 180 degrees.} + +\item{point_size}{A number of the size of the data points if "style = 'point'". +The default is 'auto' and the function tries to find the appropriate size.} + +\item{title_size}{A number of the size of the top title. The default is 16.} + +\item{dots_size}{A number of the size of the dots. The default is 0.5.} + +\item{dots_shape}{A number indicating the dot shape recognized by parameter +'shape' in \code{geom_point()}.} + +\item{coastlines_width}{A number indicating the width of the coastlines.} + +\item{fileout}{A character string of the path to save the plot. If not +specified (default), a graphic device will pop up. The extension should be +accepted by \code{ggsave()}.} + +\item{width}{A number of the plot width, in the units specified in parameter +'size_units'. The default is 8.} + +\item{height}{A number of the plot height, in the units specified in parameter +'size_units'. The default is 4.} + +\item{size_units}{A character string of the units of the size of the device +(file or window) to plot in. The default is 'in' (inches). See ?ggsave and +?Devices for details of the corresponding device.} + +\item{res}{Resolution of the device (file or window) to plot in. The default +value is 300. See ?ggsave 'dpi' and ?Devices for details of the +corresponding device.} +} +\value{ +A map plot with speficied projection, either in pop-up window or a + saved file. +} +\description{ +Transform a regular grid longitude-latitude data to a different projection and +plot the map. The target projection must be a valid CRS string, preferrably be +EPSG or ESRI code; check \link[sf]{st_crs} for more explanation. This function +is mainly tested for Robinson projection (ESRI:54030), but it can work with +other projection types in theory.\n +The map can be plotted by points or polygon. A legend can be plotted as either +a color bar or a discrete ggplot legend. Dots can be drawn on top of the data, +which can be used for significance test. A mask can be added to not plot the +data specified. A number of options is provided to adjust aesthetics, like +position, size, colors, etc. +} +\examples{ +data <- array(rep(seq(-10, 10, length.out = 181), 360) + rnorm(360), + dim = c(lat = 181, lon = 360)) +dots <- data +dots[which(dots < 4 & dots > -4)] <- 0 +dots[which(dots != 0)] <- 1 +PlotRobinson(data, lon = 0:359, lat = -90:90, dots = dots, + brks = seq(-10, 10, length.out = 11), + toptitle = 'synthetic example', vertical = F, + caption = 'Robinson Global\ns2dv::PlotRobinson example', + bar_extra_margin = c(0, 1, 0, 1), width = 8, height = 6) +PlotRobinson(data, lon = 0:359, lat = -90:90, mask = dots, legend = 'ggplot2', + target_proj = "+proj=moll", brks = seq(-10, 10, length.out = 11), + color_fun = ClimPalette("purpleorange"), colNA = 'green', + toptitle = 'synthetic example', + caption = 'Mollweide Global\ns2dv::PlotRobinson example', + width = 8, height = 6) + +} diff --git a/tests/testthat/_snaps/PlotRobinson/PlotRobinson_1a.png b/tests/testthat/_snaps/PlotRobinson/PlotRobinson_1a.png new file mode 100644 index 0000000000000000000000000000000000000000..bb12bad6df3b712142afcbd07398bbc4c8ec86ee Binary files /dev/null and b/tests/testthat/_snaps/PlotRobinson/PlotRobinson_1a.png differ diff --git a/tests/testthat/_snaps/PlotRobinson/PlotRobinson_1b.png b/tests/testthat/_snaps/PlotRobinson/PlotRobinson_1b.png new file mode 100644 index 0000000000000000000000000000000000000000..31a42335dc96c816b19ba26f072b2ff44f4b0bd6 Binary files /dev/null and b/tests/testthat/_snaps/PlotRobinson/PlotRobinson_1b.png differ diff --git a/tests/testthat/_snaps/PlotRobinson/PlotRobinson_1c.png b/tests/testthat/_snaps/PlotRobinson/PlotRobinson_1c.png new file mode 100644 index 0000000000000000000000000000000000000000..6becc0fb0f9f168c7150cd899cb7a7e44e24cfe6 Binary files /dev/null and b/tests/testthat/_snaps/PlotRobinson/PlotRobinson_1c.png differ diff --git a/tests/testthat/_snaps/PlotRobinson/PlotRobinson_2a.png b/tests/testthat/_snaps/PlotRobinson/PlotRobinson_2a.png new file mode 100644 index 0000000000000000000000000000000000000000..811ceb7dfe2b862c9e2772cdff963c02d2fd4084 Binary files /dev/null and b/tests/testthat/_snaps/PlotRobinson/PlotRobinson_2a.png differ diff --git a/tests/testthat/test-PlotRobinson.R b/tests/testthat/test-PlotRobinson.R new file mode 100644 index 0000000000000000000000000000000000000000..c16d1bc599bad572b9cc591f415c8f369d646e19 --- /dev/null +++ b/tests/testthat/test-PlotRobinson.R @@ -0,0 +1,157 @@ +#=============================================================== +# NOTE: The figures are generated with the following environment: +#- On workstation +#- R/4.1.2-foss-2015a-bare +#- GDAL/2.2.1-foss-2015a +#- PROJ/4.8.0-foss-2015a +#- GEOS/3.7.2-foss-2015a-Python-2.7.9 +#=============================================================== + +# data1: global +set.seed(19) +data1 <- array(rep(seq(-10, 10, length.out = 181), 360) + rnorm(360), + dim = c(sdate = 1, lat = 181, lon = 360)) +mask1 <- drop(data1) +mask1[which(mask1 > 9 & mask1 < -9)] <- 0 +mask1[which(mask1 != 0)] <- 1 +dots1 <- mask1 +lon1_1 <- -180:179 +lon1_2 <- 0:359 +lat1 <- -90:90 + +data1_NA <- data1; data1_NA[10000:15000] <- NA + +# data2: Europe +data2 <- data1[,21:61, 161:221] +lon2 <- lon1_1[161:221] +lat2 <- 70:30 +dots2 <- array(t(dots1[21:61, 161:221]), dim = c(lon = length(lon2), lat = length(lat2))) + +#-------------------------------------------------------------------- + +test_that("1. Input checks", { + +# data +expect_error( + PlotRobinson(array(data1*2, dim = c(sdate = 2, lat = 181, lon = 360))), + "Parameter 'data' must have two dimensions." +) +# lon, lon_dim +tmp <- data1; dim(tmp) <- c(sdate = 1, latt = 181, lonn = 360) +expect_error( + PlotRobinson(tmp), + "Cannot find known longitude name in data dimension. Please define parameter 'lon_dim'." +) +expect_error( + PlotRobinson(data1, lon = 1:10), + "Length of parameter 'lon' should match longitude dimension in 'data'." +) +# target_proj +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, target_proj = NULL), + "Parameter 'target_proj' cannot be NULL." +) +# legend +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, legend = 1), + "Parameter 'legend' must be NULL, ggplot2 or s2dv." +) +# style +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, style = 'line'), + "Parameter 'style' must be 'point' or 'polygon'." +) +# dots +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, dots = c(1, 0)), + "Parameter 'dots' must have two dimensions named as longitude and latitude dimensions in 'data'." +) +tmp <- dots1; dim(tmp) <- c(sdate = 1, lat = 360, lon = 181) +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, dots = tmp), + "Parameter 'dots' must have the same dimensions as 'data'." +) +tmp <- dots1; tmp[1] <- NA +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, dots = tmp), + "Parameter 'dots' must have only TRUE/FALSE or 0/1." +) +# colNA +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, colNA = 'fair'), + "Parameter 'colNA' must be a valid colour identifier." +) +# toptitle +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, toptitle = 1:10), + "Parameter 'toptitle' must be a character string." +) +# caption +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, caption = 1:10), + "Parameter 'caption' must be a character string." +) +# crop_coastlines +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, crop_coastlines = c(x1 = 1, x2 = 10, y1 = -10, y2 = 10)), + "Parameter 'crop_coastlines' needs to have names 'latmax', 'latmin', 'lonmax', 'lonmin'." +) +# point_size +expect_error( + PlotRobinson(data1, lon = lon1_1, lat = lat1, point_size = 'small'), + "Parameter 'point_size' must be a number." +) + +}) + +#-------------------------------------------------------- +# PlotRobinson +save_PlotRobinson <- function(...) { + path <- tempfile(fileext = ".png") + do.call(PlotRobinson, list(..., fileout = path)) + path +} + +test_that("2. data1: Global", { + +# up triangle end only +expect_snapshot_file( + save_PlotRobinson(data1, lon = lon1_1, lat = lat1, dots = dots1, brks = seq(-12.6, 12.6, 1), vertical = F, toptitle = 'PlotRobinson', caption = 'unit test: 1a'), + name = 'PlotRobinson_1a.png' +) +# NAs +expect_snapshot_file( + save_PlotRobinson(data1_NA, lon = lon1_1, lat = lat1, dots = dots1, brks = seq(-13, 13, 2), toptitle = 'PlotRobinson', caption = 'unit test: 1b', bar_extra_margin = c(5, 0, 5, 0)), + name = 'PlotRobinson_1b.png' +) + +expect_snapshot_file( + save_PlotRobinson(data1, lon = lon1_1, lat = lat1, mask = mask1, brks = seq(-11, 11, 2), color_fun = ClimPalette('purpleorange'), colNA = 'red', legend = 'ggplot2', toptitle = 'PlotRobinson', caption = 'unit test: 1c'), + name = 'PlotRobinson_1c.png' +) + + +}) + + +test_that("3. data2: Europe", { + +expect_snapshot_file( + save_PlotRobinson(data2, lon = lon2, lat = lat2, dots = dots2, target_proj = 102014, + brks = seq(-10, 0, 2), cols = c("#FDD49E", "#FDBB84", "#FC8D59", "#E34A33", "#B30000"), col_inf = "#FEF0D9", + toptitle = 'PlotRobinson', caption = 'unit test: 1a\n projection: Lambert Europe (ESRI:102014)', point_size = 1.1, title_size = 12, + width = 6, vertical = F, bar_extra_margin = c(0.5, 4, 0.5, 4), crop_coastlines = c(latmin = 27, latmax = 70, lonmin = -20, lonmax = 40)), + name = 'PlotRobinson_2a.png' +) + +#polygon +#expect_snapshot_file( +# save_PlotRobinson(data2, lon = lon2, lat = lat2, dots = dots2, target_proj = 102014, +# brks = seq(-10, 0, 2), cols = c("#FDD49E", "#FDBB84", "#FC8D59", "#E34A33", "#B30000"), col_inf = "#FEF0D9", +# toptitle = 'PlotRobinson', caption = 'unit test: 1a\n projection: Lambert Europe (ESRI:102014)', point_size = 1.1, title_size = 12, +# width = 6, vertical = F, bar_extra_margin = c(0.5, 4, 0.5, 4), crop_coastlines = c(latmin = 27, latmax = 70, lonmin = -20, lonmax = 40), style = 'polygon'), +# name = 'PlotRobinson_2a.png' +#) + + +}) diff --git a/vignettes/Figures/vis_proj_eu_1.png b/vignettes/Figures/vis_proj_eu_1.png new file mode 100644 index 0000000000000000000000000000000000000000..7275ea88d10278c12702bb3419e4aa7802839af4 Binary files /dev/null and b/vignettes/Figures/vis_proj_eu_1.png differ diff --git a/vignettes/Figures/vis_proj_eu_2.png b/vignettes/Figures/vis_proj_eu_2.png new file mode 100644 index 0000000000000000000000000000000000000000..dae2c4d5b6d8380365e623bc39c18867e92c129b Binary files /dev/null and b/vignettes/Figures/vis_proj_eu_2.png differ diff --git a/vignettes/Figures/vis_proj_global_1.png b/vignettes/Figures/vis_proj_global_1.png new file mode 100644 index 0000000000000000000000000000000000000000..dbd18c15d35392f1993e20fe9cee8c053f1cc6cf Binary files /dev/null and b/vignettes/Figures/vis_proj_global_1.png differ diff --git a/vignettes/Figures/vis_proj_global_2.png b/vignettes/Figures/vis_proj_global_2.png new file mode 100644 index 0000000000000000000000000000000000000000..e82903bc2c2b136b0c9c49618d8c35e1162f3294 Binary files /dev/null and b/vignettes/Figures/vis_proj_global_2.png differ diff --git a/vignettes/Figures/vis_proj_np_1.png b/vignettes/Figures/vis_proj_np_1.png new file mode 100644 index 0000000000000000000000000000000000000000..a3f4a1995290c4813729d4a9f152aff49a5d5286 Binary files /dev/null and b/vignettes/Figures/vis_proj_np_1.png differ diff --git a/vignettes/Figures/vis_proj_np_2.png b/vignettes/Figures/vis_proj_np_2.png new file mode 100644 index 0000000000000000000000000000000000000000..b2a7e49ed4726fd2278e8e6345289aca39b32d9a Binary files /dev/null and b/vignettes/Figures/vis_proj_np_2.png differ diff --git a/vignettes/visualization_projection.md b/vignettes/visualization_projection.md new file mode 100644 index 0000000000000000000000000000000000000000..432b2fda53e393aecfa3eb203590099c83f2be05 --- /dev/null +++ b/vignettes/visualization_projection.md @@ -0,0 +1,237 @@ +--- +author: "An-Chi Ho" +date: "2023-05-03" +output: html_document +vignette: > + %\VignetteEngine{R.rsp::md} + %\VignetteIndexEntry{Visualisation} + %\usepackage[utf8]{inputenc} +--- + +# Visualization of different projections + +This vignette shows how to use function `PlotRobinson()`, which can plot maps with the projection of your choice. Since packages `sf` and `rnaturalearth` are required by the function, make sure you have modules `PROJ`, `GEOS`, `GDAL` loaded. +The projection types and IDs can be found on https://epsg.io/. + +```r +## Load packages +library(sf) +library(ggplot2) +library(rnaturalearth) +library(RColorBrewer) +library(cowplot) +library(startR) + +#TODO: Change to load the package +## Source function +source("/esarchive/scratch/aho/tmp/PlotRobinson.R") +## Source unexported s2dv functions used in PlotRobinson +.KnownLonNames <- s2dv:::.KnownLonNames +.KnownLatNames <- s2dv:::.KnownLatNames +.IsColor <- s2dv:::.IsColor +.warning <- s2dv:::.warning + +``` +Specify where you want to save the figures later. +```r +fileout_path <- "/esarchive/scratch/" +``` + +## Global map + +### 1. Load data + +We use `startR::Start` to load system5c3s monthly temperature data from esarchive. To focus on visualization, we simply load one file with one ensemble and one time step. + +```r +repos <- "/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc" +data <- Start(dat = list(list(name = 'system5c3s', path = repos)), + var = 'tas', + sdate = '20170101', + ensemble = indices(1), + time = indices(1), + lon = 'all', + lon_reorder = CircularSort(-180, 180), + lat = 'all', + lat_reorder = Sort(decreasing = T), + synonims = list(lat = c('lat', 'latitude'), + lon = c('lon', 'longitude')), + return_vars = list(time = 'sdate', longitude = 'dat', latitude = 'dat'), + retrieve = T) +lon <- attributes(data)$Variable$system5c3s$lon +lat <- attributes(data)$Variable$system5c3s$lat +data <- data - 273.15 # K to degC +``` + +### 2. Create mask + +We can create a logical or [0, 1] array to serve as mask in the plot. +```r +mask <- drop(data) +mask[which(mask > -20 & mask < -10)] <- 0 +mask[which(mask != 0)] <- 1 +``` + +### 3. Create colorbar breaks and colors + +Create the adequete colorbar elements for the plot +```r +## have exceeded values +brks1 <- seq(-40, 40, 10) +cols1 <- brewer.pal(9, 'PuOr')[-9] +## cover all values +#brks1 <- seq(-ceiling(max(abs(range(data)))), ceiling(max(abs(range(data)))), length.out = 11) +``` +### 4. Plot + +You can find more projection types through the URL above. `PlotRobinson()` is mainly designed for Robinson projection, but it is possible to use other projections as well. +E.g., Mollweide projection is "ESRI:54009" or "ESRI:53009". +Note that different PROJ and GDAL module versions have different CRS codes. If on Nord3v2, please use "ESRI:53030" or "+proj=robin" to get the correct Robinson projection. + +There are two types of legend, 'ggplot2' or 's2dv'. You can also set `legend` as `NULL` to not plot legend. + +```r +## legend = 'ggplot2' +PlotRobinson(data, lon, lat, brks = brks1, cols = NULL, legend = 'ggplot2', + target_proj = 54030, style = 'point', + toptitle = "system5c3s monthly tas Jan 2017", caption = "Projection: Robinson", + dots = mask, units = 'degC', colNA = 'red', + width = 10, height = 6) + +## legend = 's2dv' +PlotRobinson(data, lon, lat, brks = brks1, cols = cols1, legend = 's2dv', + target_proj = 54030, style = 'point', + toptitle = "system5c3s monthly tas Jan 2017", caption = "Projection: Robinson", + dots = mask, units = "degC", col_inf = 'brown', + width = 10, height = 5, + bar_extra_margin = c(0, 5, 0, 5), + vertical = F) +``` + + + + + + +## European region + +### 1. Subset data + +Crop the european area from the global data we loaded above. + +```r +data_eu <- ClimProjDiags::Subset(data, c('lat', 'lon'), list(21:61, 161:221)) +lon_eu <- lon[161:221] +lat_eu <- lat[21:61] +``` +### 2. Crop mask + +```r +mask_eu <- ClimProjDiags::Subset(mask, c('lat', 'lon'), list(21:61, 161:221)) +``` + +### 3. Create colorbar breaks + +```r +## have exceeded values +brks_eu <- seq(-18, 18, length.out = 11) +## cover all values +#brks_eu <- seq(-ceiling(max(abs(range(data_eu)))), ceiling(max(abs(range(data_eu)))), length.out = 11) +``` + +### 4. Plot +In the previous plots, we use `dots` to mark the masked area. This time, we use parameter `mask` to remove the unwanted data. The removed points are left as blank. + +The parameter `crop_coastlines` helps to adjust the map coastline plotting. You can test and find the best range. Note that this parameter is not suitable for the case that longitude range exceeds half of the earth (i.e., 180 degrees.) + +We can also choose to plot the map by polygon instead of point. Note that polygon is relatively **slow** so it takes a long time for global map, and it **does not work on our workstation** for some technical reason. + +Note that `target_proj = 102014` does not work well on Nord3v2. Use PROJ.4 string `target_proj = "+proj=lcc +lat_0=30 +lon_0=10 +lat_1=43 +lat_2=62 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs +type=crs"` instead. + +```r +## style = point +PlotRobinson(data_eu, lon_eu, lat_eu, brks = brks_eu, cols = NULL, legend = 's2dv', + target_proj = 102014, style = 'point', + toptitle = "system5c3s monthly tas Jan 2017", + caption = paste0("Projection: Lambert Europe\nESRI:102014\n", "Europe"), + mask = mask_eu, units = "degC", + width = 8, height = 6, point_size = 1.3, + crop_coastlines = c(lonmin = -20, lonmax = 40, latmin = 27, latmax = 70), + triangle_ends = c(T, T), bar_extra_margin = c(0.4, 4, 0.4, 4), + vertical = F) + +## style = polygon +PlotRobinson(data_eu, lon_eu, lat_eu, brks = brks_eu, cols = NULL, legend = 's2dv', + target_proj = 102014, style = 'polygon', + toptitle = "system5c3s monthly tas Jan 2017", + caption = paste0("Projection: Lambert Europe\nESRI:102014\n", "Europe"), + mask = mask_eu, units = "degC", + width = 8, height = 6, + crop_coastlines = c(lonmin = -20, lonmax = 40, latmin = 27, latmax = 70), + triangle_ends = c(T, T), bar_extra_margin = c(0.4, 4, 0.4, 4), + vertical = F) +``` + + + + + + +## North Pole region + +### 1. Subset data + +Crop the north pole area from the global data we loaded above. +```r +data_np <- ClimProjDiags::Subset(data, c('lat'), list(1:31)) +lon_np <- lon +lat_np <- lat[1:31] +``` +### 2. Crop mask + +```r +mask_np <- ClimProjDiags::Subset(mask, c('lat'), list(1:31)) +``` + +### 3. Create colorbar breaks + +```r +## have exceeded values +brks_np <- seq(-40, 0, length.out = 11) +## cover all values +#brks_np <- seq(-ceiling(max(abs(range(data_np)))), ceiling(max(abs(range(data_np)))), length.out = 11) +``` + +### 4. Plot + +#TODO: Update PlotStereoMap to the new function +We compare `PlotRobinson()` with `PlotStereoMap()`, which is specifically for the stereographic projection and it uses polygon to plot. +Note that we cannot use `crop_coastlines` here in `PlotRobinson()` because the longitude range is global. +It is recommended saving the figures for these two plots (especially `PlotStereoMap()`) because polygon plotting is time-consuming when it is in the pop-up window. + +The projection used here is Arctic Polar Stereographic (EPSG:3995). +```r +#NOTE: A bit slow, be patience. +PlotRobinson(data_np, lon_np, lat_np, brks = brks_np, cols = NULL, legend = 's2dv', + target_proj = 3995, style = 'polygon', + toptitle = "system5c3s monthly tas Jan 2017", + caption = paste0("Projection: WGS 84 / Arctic Polar Stereographic\nEPSG:3995\n", "North Pole"), + dots = mask_np, dots_size = 1.5, + units = "degC", color_fun = clim.palette("orangepurple"), + width = 8, height = 8, + triangle_ends = c(T, T), bar_extra_margin = c(0.4, 4, 0.4, 4), + vertical = F, fileout = paste0(fileout_path, '/PlotRobinson_np.png')) + +# Compare with PlotStereoMap +#NOTE: Save the figure instead of plotting in pop-up window because it is slow. +PlotStereoMap(data_np, lon = lon_np, lat = lat_np, brks = brks_np, cols = NULL, + color_fun = clim.palette("orangepurple"), + toptitle = "system5c3s monthly tas Jan 2017 (PlotStereoMap)", + width = 8, height = 8, filled.continent = F, + dots = !mask_np, units = "degC", title_scale = 0.8, + fileout = paste0(fileout_path, '/PlotStereoMap_np.png')) +``` + + + +