diff --git a/DESCRIPTION b/DESCRIPTION index 32c52d1..1ff5461 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,9 +1,9 @@ Package: visualizeR Depends: R(>= 3.1.0), - sm + sm, + transformeR(>= 1.2.0) Imports: - transformeR, methods, fields, RColorBrewer, @@ -15,15 +15,21 @@ Imports: SpecsVerification, verification, sp, - lattice + lattice, + latticeExtra, + magrittr, + padr, + grDevices, + data.table Suggests: loadeR.ECOMS, loadeR, + easyVerification, knitr Type: Package Title: Visualizing and Communicating Uncertainty in Seasonal Climate Prediction -Version: 1.0.0 -Date: 2017-10-20 +Version: 1.1.0 +Date: 2017-12-20 Authors@R: as.person(c( "Santander Meteorology Group [ctb]", "Maria Dolores Frias [aut, cre]", diff --git a/NAMESPACE b/NAMESPACE index d41bd02..ed9b31e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,15 +3,23 @@ export(MrQuantile) export(QuantileProbs) export(bubblePlot) +export(clim2sgdf) export(detrend.data) +export(map.lines) +export(map.stippling) export(reliabilityCategories) export(rocss.fun) export(seasMean) +export(skillMap) export(spatialMean) +export(spatialPlot) export(spreadPlot) +export(temporalPlot) export(tercileBarplot) +export(tercileMap) export(tercilePlot) import(lattice) +import(latticeExtra) import(methods) import(sm) import(verification) @@ -20,8 +28,12 @@ importFrom(RColorBrewer,brewer.pal) importFrom(RCurl,getURL) importFrom(SpecsVerification,Auc) importFrom(abind,abind) +importFrom(abind,asub) +importFrom(data.table,rleid) importFrom(fields,image.plot) importFrom(grDevices,col2rgb) +importFrom(grDevices,colorRampPalette) +importFrom(grDevices,colors) importFrom(grDevices,gray) importFrom(grDevices,grey.colors) importFrom(grDevices,rgb) @@ -40,29 +52,50 @@ importFrom(graphics,points) importFrom(graphics,polygon) importFrom(graphics,text) importFrom(graphics,title) +importFrom(lattice,panel.grid) +importFrom(lattice,panel.polygon) +importFrom(lattice,panel.xyplot) +importFrom(lattice,xyplot) +importFrom(magrittr,"%>%") importFrom(mapplots,add.pie) importFrom(mapplots,draw.pie) +importFrom(padr,pad) importFrom(scales,alpha) +importFrom(sp,GridTopology) +importFrom(sp,Line) +importFrom(sp,Lines) importFrom(sp,Polygon) importFrom(sp,Polygons) +importFrom(sp,SpatialGridDataFrame) +importFrom(sp,SpatialLines) importFrom(sp,SpatialPoints) +importFrom(sp,SpatialPointsDataFrame) importFrom(sp,SpatialPolygons) importFrom(sp,over) +importFrom(sp,spplot) importFrom(stats,coef) importFrom(stats,complete.cases) importFrom(stats,filter) importFrom(stats,lm) +importFrom(stats,median) importFrom(stats,na.omit) +importFrom(stats,qnorm) importFrom(stats,quantile) +importFrom(transformeR,aggregateGrid) importFrom(transformeR,array3Dto2Dmat) +importFrom(transformeR,climatology) importFrom(transformeR,draw.world.lines) importFrom(transformeR,getCoordinates) importFrom(transformeR,getDim) importFrom(transformeR,getGrid) +importFrom(transformeR,getShape) importFrom(transformeR,getYearsAsINDEX) importFrom(transformeR,interpGrid) +importFrom(transformeR,isRegular) importFrom(transformeR,makeMultiGrid) importFrom(transformeR,mat2Dto3Darray) importFrom(transformeR,redim) importFrom(transformeR,subsetGrid) +importFrom(utils,head) importFrom(utils,packageDescription) +importFrom(utils,tail) diff --git a/NEWS b/NEWS index ee407de..ad8f478 100644 --- a/NEWS +++ b/NEWS @@ -1,7 +1,8 @@ -visualizeR 1.0.0 +visualizeR 1.1.0 ================= -* Internal adaptation of reliabilityCategories to be consistent with the other functions. -* Implemented the ROCSS significance in tercilePlot and tercileBarplot. -* Added CITATION file to Environmental Modelling & Software. +* New function `spatialPlot` (essentially a renaming of `transformeR::plotClimatology`, the latter planned to be deprecated in future `transformeR` releases) +* New plotting functions for seasonal forecast products: + * `tercileMap` for forecast tercile maps (based on C3S seasonal forecast graphical products, http://climate.copernicus.eu/s/charts/c3s_seasonal/) + * `skillMap`: A convenient wrapper of `SpatialPlot` and `map.stippling` for plotting skill maps * Other minor changes and documentation updates. diff --git a/R/auxiliaryFun.R b/R/auxiliaryFun.R index b55f210..5f20ca0 100644 --- a/R/auxiliaryFun.R +++ b/R/auxiliaryFun.R @@ -210,6 +210,7 @@ getSeason.S4 <- function(obj) { #' @return The ROC area skill score and the significance (TRUE or FALSE) #' @author M. D. Frias \email{mariadolores.frias@@unican.es} and J. Fernandez #' @importFrom SpecsVerification Auc +#' @importFrom stats qnorm #' @export rocss.fun <- function (obs, pred, conf.level = 0.95){ no.nan <- complete.cases(obs, pred) diff --git a/R/reliabilityCategories.R b/R/reliabilityCategories.R index 17d0930..0a9912b 100644 --- a/R/reliabilityCategories.R +++ b/R/reliabilityCategories.R @@ -31,6 +31,8 @@ #' @param n.bins (optional): number of probability bins considered. By default n.bins = 10 #' @param n.boot number of samples considered for bootstrapping. By default n.boot = 100 #' @param conf.level Confidence interval for the reliability line. By default \code{conf.level} = 0.75 (two sided), as in Weisheimer et al. 2014 +#' @param na.rate Allowed proportion of NA values in each region. Regions with proportions higher than na.rate +#' are excuded from the analysis. Default is 0.75. #' @param diagrams Logical (default = TRUE). Plotting results. #' @param cex0 numeric (default is 0.5). Minimum size of the points shown in the reliability diagrams, i.e. size of the point #' for the minimum n frequency (n = 1) (see parameter \code{n.bins}. The sizes for points that correspond to n > 1 @@ -108,6 +110,7 @@ reliabilityCategories <- function(hindcast, n.bins = 10, n.boot = 100, conf.level = 0.75, + na.rate = 0.75, diagrams = TRUE, cex0 = 0.5, cex.scale = 20, @@ -166,7 +169,8 @@ reliabilityCategories <- function(hindcast, for(l in 1:length(regs)){ o <- over(spoints, regs[l,]) w <- which(!is.na(o)) - if(length(w)>0){ + if(1-length(w)/length(o) > na.rate) w <- NULL + if(length(w) > 0) { ob <- ob.full[, w, drop = F] se <- array(dim = c(nmem, ntime, length(w))) for(i in 1:nmem){ @@ -265,8 +269,6 @@ reliabilityCategories <- function(hindcast, at = c(1, 2, 2.75, 3.25, 4, 5), labels = c("dangerously useless", "not useful","marginally useful", "marginally useful +","still useful","perfect")))) - - print(pc) }else{ x1 <- 1/n.events y1 <- 1/n.events @@ -289,7 +291,7 @@ reliabilityCategories <- function(hindcast, # Customized Lattice Example - xyp <- xyplot(y~x|z, par.strip = list(lines = 1), ylim = c(0,1), strip = strip.custom(fg = rgb(0,0,0,0), strip.names = c(T,F), strip.levels = c(F,T), factor.levels = labels), + pc <- xyplot(y~x|z, par.strip = list(lines = 1), ylim = c(0,1), strip = strip.custom(fg = rgb(0,0,0,0), strip.names = c(T,F), strip.levels = c(F,T), factor.levels = labels), scales=list(x = list(at = seq(0,1,round(1/n.bins, digits = 2)), labels = seq(0,1,round(1/n.bins, digits = 2))), y = list(at = seq(0,1,round(1/n.bins, digits = 2)), @@ -336,7 +338,7 @@ reliabilityCategories <- function(hindcast, layout = layout, xlab = "Predicted probability", ylab= "Observed frequency", main= list(cex = 1, font = 1, label = sprintf("n = %d years x %d points", nyear, npoint))) - print(xyp) + # update(xyp, par.settings = list(fontsize = list(text = 8, points = 10))) ######################################################### @@ -394,6 +396,7 @@ reliabilityCategories <- function(hindcast, ########################################################### } + print(pc) } result.grid <- mg attr(result.grid$Data, "dimensions") <- c("cat", "var", "member", "time", "lat", "lon") diff --git a/R/skillMap.R b/R/skillMap.R new file mode 100644 index 0000000..88bddbd --- /dev/null +++ b/R/skillMap.R @@ -0,0 +1,144 @@ +# skillMap.R Skill maps for seasonal forecasts +# +# Copyright (C) 2017 Santander Meteorology Group (http://www.meteo.unican.es) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +#' @title A wrapper of \code{spatialPlot} for the creation of verification maps for seasonal forecast systems +#' +#' @description A wrapper of \code{spatialPlot} for the creation of verification maps for seasonal forecast systems. +#' It provides a convenient interface for \code{\link{map.stippling}} +#' +#' +#' @param easyVeriGrid A climatological grid with a verification output. See details +#' @param stippling A key-value list of arguments passed to \code{map.stippling}: +#' \code{threshold} and \code{condition}. Ignored by default, returning a map without stippling. +#' @param stippling.point.options Default to \code{NULL}. Further graphical arguments passed to points +#' (e.g. \code{cex}, \code{pch} etc.) +#' @param backdrop.theme See \code{\link{spatialPlot}} +#' @param title Title of the plot +#' @param ... Further graphical options passed to \code{\link{spatialPlot}}. +#' +#' @details The function applies the \code{\link{spatialPlot}} function, that in turn uses \code{lattice-methods}. +#' +#' \strong{Graphical options} +#' +#' Some examples of specific map graphical options are available in the help of function \code{\link[sp]{spplot}}. +#' In addtion, fine-tuning of the resulting plots can be obtained using the arguments of \pkg{lattice} plots. For +#' an overview, see the help of function \code{\link[lattice]{xyplot}}. +#' +#' @seealso The bridging function \code{\link[transformeR]{easyVeri2grid}} from package \pkg{transformeR} allows for +#' the conversion of verification outputs from package \pkg{easyVerification} to the \code{climate4R} data structure +#' used by the function. +#' +#' Many different aspects of the plot can be controlled passing the relevant arguments to +#' \code{\link[sp]{spplot}}. +#' +#' @return A lattice plot of class \dQuote{trellis}. +#' +#' @author J. Bedia +#' @export +#' @examples \dontrun{ +#' +#' # The package 'easyVerification' will be used to calculate the RPSS: +#' +#' require(easyVerification) +#' require(transformeR) +#' # First of all, a data subset is done, considering a spatial domain centered on the North Atlantic: +#' tas.cfs2 <- subsetGrid(tas.cfs, lonLim = c(-100, 40), latLim = c(-5, 75)) +#' # The same is done with the reanalysis dataset: +#' tas.ncep2 <- subsetGrid(tas.ncep, lonLim = c(-100, 40), latLim = c(-5, 75)) +#' # In the next step, the reanalysis data are interpolated to the hindcast grid: +#' tas.ncep2.int <- interpGrid(tas.ncep2, new.coordinates = getGrid(tas.cfs2), method = "nearest") +#' # We compute the Ranked Probabiloity SKill Score using veriApply, and a cross-validation strategy: +#' ev <- easyVerification::veriApply(verifun = "EnsRpss", +#' fcst = tas.cfs2$Data, +#' obs = tas.ncep2.int$Data, +#' prob = 1:2/3, +#' tdim = 2, +#' ensdim = 1, +#' parallel = TRUE, +#' ncpus = 3, +#' strategy = "crossval") +#' # The bridging function 'easyVeri2grid' converts the object returned by 'veriApply' +#' # to a 'climate4R' climatological grid: +#' +#' easyVeriGrid <- easyVeri2grid(ev$skillscore, +#' obs.grid = tas.ncep2.int, +#' verifun = "EnsRpss") +#' +#' # A basic RPSS map, using the spatialPlot defaults: +#' skillMap(easyVeriGrid = easyVeriGrid, backdrop.theme = "coastline") +#' +#' # Stippling. Mark significant RPSS values at the 95% ci +#' thresh <- ev$skillscore.sd*qnorm(0.95) +#' +#' skillMap(easyVeriGrid = easyVeriGrid, +#' stippling = list(threshold = th, condition = "GT"), # GT = greater than +#' stippling.point.options = list(pch = 19, cex = .2, col = "black"), +#' backdrop.theme = "coastline") +#' +#' # Further customization: For instance a more elaborated title, colorblind-friendly palette etc.: +#' +#' cb.colors <- colorRampPalette(rev(RColorBrewer::brewer.pal(11, "Spectral"))) +#' skillMap(easyVeriGrid = easyVeriGrid, +#' stippling = list(threshold = th, condition = "GT"), +#' stippling.point.options = list(pch = 19, cex = .2, col = "black"), +#' backdrop.theme = "coastline", +#' at = seq(-0.7, 0.7, .05), +#' set.max = 0.7, set.min = -0.7, +#' scales = list(draw = TRUE, alternating = 3), +#' colorkey = list(space = "bottom"), +#' col.regions = cb.colors(51), +#' title = paste("Ranked Probability Skill Score", +#' "Forecasting System: CFSv2 - 24 members", +#' "November Initializations", +#' "Observing System: NCEP/NCAR Reanalysis 1", +#' "t2m DJF 1983-2010", sep = "\n") +#' ) +#' } + +skillMap <- function(easyVeriGrid, + stippling = list(threshold = NULL, condition = NULL), + stippling.point.options = NULL, + backdrop.theme = NULL, + title = NULL, + ...) { + arg.list <- list(...) + arg.list[["grid"]] <- easyVeriGrid + if (!is.null(stippling$threshold)) { + stippling[["clim"]] <- easyVeriGrid + if (!is.null(stippling.point.options)) { + stippling <- c(stippling, stippling.point.options) + } + stip.list <- list(do.call("map.stippling", stippling)) + if ("sp.layout" %in% names(arg.list)) { + arg.list[["sp.layout"]] <- c(arg.list[["sp.layout"]], stip.list) + } else { + arg.list[["sp.layout"]] <- stip.list + } + } + if (!is.null(backdrop.theme)) { + arg.list[["backdrop.theme"]] <- backdrop.theme + } + if (!is.null(title)) { + arg.list[["main"]] <- list(label = title, + cex = .9, + col = "blue", + font = 1, + just = "left", x = 0.05) + } + map <- do.call("spatialPlot", arg.list) + print(map) +} diff --git a/R/spatialPlot.R b/R/spatialPlot.R new file mode 100644 index 0000000..00f6f9a --- /dev/null +++ b/R/spatialPlot.R @@ -0,0 +1,442 @@ +# spatialPlot.R Lattice plot methods for climatological grids +# +# Copyright (C) 2017 Santander Meteorology Group (http://www.meteo.unican.es) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + + + +#' @title Lattice plot methods for climatological grids +#' @description A wrapper for the lattice (trellis) plot methods for spatial data in \code{sp::spplot} +#' @param grid Input grid +#' @param backdrop.theme Reference geographical lines to be added to the plot. See Details. +#' @param set.min Numeric value indicating an absolute minimum value (default to \code{NULL}). All grid values below this are mapped to \code{set.min}. See details. +#' @param set.max Same as \code{set.min} argument, but to force a ceiling. +#' @param lonCenter Value of the longitude to be centered in the plot. +#' @param ... Further arguments passed to \code{spplot} +#' @details The function applies the \code{\link[sp]{spplot}} method after conversion of the climatological map(s) to a +#' \code{SpatialGridDataFrame}. +#' +#' The \code{set.min} and \code{set.max} options are useful in order to preserve adequate ranges for map representation, +#' avoiding the influence of extreme values. Note that this is different than setting a range of values with an +#' interval using the \code{at} argument. The latter choice, that overrides \code{set.min} and \code{set.max}, +#' leaves blank grid points for outlying values. +#' +#' \strong{Multigrids} +#' +#' Multigrids of climatologies can be created using \code{makeMultiGrid} +#' for trellis visualization of different variables, or for instance, for the comparison of +#' raw and corrected/downscaled scenarios side to side. In case of multimember multigrids, +#' the function will internally compute the ensemble mean of each variable in the multigrid +#' for representation (with a message). +#' +#' \strong{Backdrop theme} +#' +#' Current implemented options are \code{"none"} and \code{"coastline"}, which contains +#' a simplied vector theme delineating the world coastlines. Any other themes can be introduced +#' by the user using the \code{sp.layout} options in \code{spplot}. +#' +#' \strong{Controlling graphical parameters} +#' +#' Many different aspects of the map can be controlled passing the relevant arguments to +#' \code{spplot}. Fine control of graphical parameters for the trellis display can +#' be also controlled using \code{\link[lattice]{trellis.par.set}}. +#' +#' Some examples of specific map graphical options are available in the help of function \code{\link[sp]{spplot}}. +#' In addtion, fine-tuning of the resulting plots can be obtained using the arguments of \pkg{lattice} plots. For +#' an overview, see the help of function \code{\link[lattice]{xyplot}}. +# +#' @return As spplot, \code{spatialPlot} returns a lattice plot of class \dQuote{trellis}. +#' If you fail to \dQuote{see} it, explicitly call \code{print(spatialPlot(...))}. +#' +#' @references \itemize{ +#' \item Bivand, R.S., Pebesma, E.J., Gomez-Rubio, V., 2013. Applied Spatial Data Analysis with R, 2nd ed, useR! Springer, NY. +#' \item For some graticulate customization examples, visit the \emph{sp Gallery}: \url{https://edzer.github.io/sp/} +#' } +#' +#' @importFrom abind abind asub +#' @importFrom sp spplot SpatialGridDataFrame SpatialPointsDataFrame GridTopology SpatialPoints +#' @importFrom grDevices colorRampPalette +#' @importFrom utils tail head +#' +#' +#' @export +#' +#' @author J. Bedia +#' @seealso \code{\link{climatology}} for details on climatology calculation. +#' \code{\link{map.stippling}}, for adding a custom point layer on top of the map. +#' \code{\link{map.lines}}, to add lines and polygons to climatological maps +#' Also see \code{\link[sp]{spplot}} in package \pkg{sp} for further information on plotting capabilities and options +#' @examples +#' require(transformeR) +#' data("CFS_Iberia_tas") +#' # Climatology is computed: +#' clim <- climatology(CFS_Iberia_tas, by.member = TRUE) +#' spatialPlot(clim) +#' # Geographical lines can be added using the argument 'backdrop.theme': +#' spatialPlot(clim, backdrop.theme = "coastline") +#' spatialPlot(clim, backdrop.theme = "coastline", center = 180) +#' +#' # Further arguments can be passed to 'spplot'... +#' +#' # ... a subset of members to be displayed, using 'zcol': +#' spatialPlot(clim, +#' backdrop.theme = "coastline", +#' zcol = 1:4) +#' +#' # ... regional focuses (e.g. the Iberian Peninsula). +#' spatialPlot(clim, +#' backdrop.theme = "countries", +#' xlim = c(-10,5), ylim = c(35,44), +#' zcol = 1:4, +#' scales = list(draw = TRUE)) +#' +#' # Changing the default color palette and ranges: +#' spatialPlot(clim, +#' backdrop.theme = "coastline", +#' zcol = 1:4, +#' col.regions = cm.colors(27), at = seq(10,37,1)) +#' +#' # For ensemble means climatology should be called with 'by.member' set to FALSE: +#' clim <- climatology(CFS_Iberia_tas, by.member = FALSE) +#' +#' # Adding contours to the plot is direct with argument 'contour': +#' spatialPlot(clim, +#' scales = list(draw = TRUE), +#' contour = TRUE, +#' main = "tas Predictions July Ensemble Mean") +#' +#' ## Example of multigrid plotting +#' data("NCEP_Iberia_psl") +#' ## Winter data are split into monthly climatologies +#' monthly.clim.grids <- lapply(getSeason(NCEP_Iberia_psl), function(x) { +#' climatology(subsetGrid(NCEP_Iberia_psl, season = x)) +#' }) +#' ## Skip the temporal checks, as grids correspond to different time slices +#' mg <- do.call("makeMultiGrid", +#' c(monthly.clim.grids, skip.temporal.check = TRUE)) +#' ## We change the panel names +#' spatialPlot(mg, +#' backdrop.theme = "coastline", +#' names.attr = c("DEC","JAN","FEB"), +#' main = "Mean PSL climatology 1991-2010", +#' scales = list(draw = TRUE)) +#' +#' # Station data: +#' data("VALUE_Iberia_pr") +#' spatialPlot(climatology(VALUE_Iberia_pr), backdrop.theme = "countries") + + +spatialPlot <- function(grid, backdrop.theme = "none", set.min = NULL, set.max = NULL, lonCenter = NULL, ...) { + arg.list <- list(...) + bt <- match.arg(backdrop.theme, choices = c("none", "coastline", "countries")) + if (!is.null(set.min) && !is.numeric(set.min)) stop("Invalid 'set.min' value") + if (!is.null(set.min) && !is.numeric(set.max)) stop("Invalid 'set.max' value") + ## add climatology:fun attribute if getShape(grid, "time") = 1 + if (getShape(grid, "time") == 1L) attr(grid$Data, "climatology:fun") <- "undefined" + ## Change lon center + if (!is.null(lonCenter)) { + indcenter <- which(abs(grid$xyCoords$x - lonCenter) == min(abs(grid$xyCoords$x - lonCenter))) + indcenter <- abs(length(grid$xyCoords$x)/2 - indcenter) + newlon <- if (indcenter == 0) grid$xyCoords$x else c(tail(grid$xyCoords$x, -indcenter), + head(grid$xyCoords$x, indcenter)) + dimNames <- attr(grid$Data, "dimensions") + climfun <- attr(grid$Data, "climatology:fun") + ind <- match(newlon, grid$xyCoords$x) + grid$Data <- asub(grid$Data, idx = ind, + dims = which(getDim(grid) == "lon"), drop = FALSE) + attr(grid$Data, "dimensions") <- dimNames + attr(grid$Data, "climatology:fun") <- climfun + grid$xyCoords$x <- newlon + } + ## Convert to spatial object + df <- clim2sgdf(clim = grid, set.min, set.max) + ## Backdrop theme + if (bt != "none") { + uri <- switch(bt, + "coastline" = system.file("coastline.rda", package = "visualizeR"), + "countries" = system.file("countries.rda", package = "visualizeR")) + l1 <- get(load(uri)) + if (is.null(arg.list[["sp.layout"]])) { + arg.list[["sp.layout"]] <- list(l1) + } else { + arg.list[["sp.layout"]][[length(arg.list[["sp.layout"]]) + 1]] <- list(l1) + } + } + ## Default colorbar + if (is.null(arg.list[["col.regions"]])) { + jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", + "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")) + arg.list[["col.regions"]] <- jet.colors(101) + } + ## Other args + arg.list[["obj"]] <- df + arg.list[["asp"]] <- 1 + do.call("spplot", arg.list) +} + +#'@title Climatology to SpatialGridDataFrame or SpatialPointsDataFrame +#'@description Convert a climatological grid to a SpatialGridDataFrame object from package sp +#'@param clim A climatological grid, as returned by function \code{\link{climatology}} +#'@param set.min Minimum value, as passed by \code{spatialPlot} +#'@param set.max Maximum value, as passed by \code{spatialPlot} +#'@seealso \code{\link{climatology}}, \code{\link{spatialPlot}} +#'@return A \pkg{sp} object of the class \code{\link[sp]{SpatialGridDataFrame}} +#'@details This function is intended for internal usage by \code{\link{spatialPlot}}, +#'that accepts all possible arguments that can be passed to \code{\link[sp]{spplot}} for plotting. +#' However, it may be useful for advanced \pkg{sp} users in different contexts +#' (e.g. for reprojecting via \code{\link[sp]{spTransform}} etc.) +#'@keywords internal +#'@author J. Bedia +#'@export +#'@importFrom sp GridTopology SpatialGridDataFrame +#'@importFrom transformeR getShape aggregateGrid isRegular +#'@examples +#' data("CFS_Iberia_tas") +#' # Climatology is computed: +#' clim <- climatology(CFS_Iberia_tas, by.member = TRUE) +#' sgdf <- clim2sgdf(clim, NULL, NULL) +#' class(sgdf) + + + +clim2sgdf <- function(clim, set.min, set.max) { + grid <- clim + if (is.null(attr(grid[["Data"]], "climatology:fun"))) { + stop("The input grid is not a climatology: Use function 'climatology' first") + } + dimNames <- getDim(grid) + ## Multigrids are treated as realizations, previously aggregated by members if present + is.multigrid <- "var" %in% dimNames + if (is.multigrid) { + if ("member" %in% dimNames) { + mem.ind <- grep("member", dimNames) + n.mem <- getShape(grid, "member") + if (n.mem > 1) message("NOTE: The multimember mean will be displayed for each variable in the multigrid") + grid <- suppressMessages(aggregateGrid(grid, aggr.mem = list(FUN = "mean", na.rm = TRUE))) + dimNames <- getDim(grid) + } + attr(grid[["Data"]], "dimensions") <- gsub("var", "member", dimNames) + } + grid <- redim(grid, drop = FALSE) + dimNames <- getDim(grid) + mem.ind <- grep("member", dimNames) + n.mem <- getShape(grid, "member") + co <- getCoordinates(grid) + if (isRegular(grid)) co <- expand.grid(co$y, co$x)[2:1] + le <- nrow(co) + #############hemen nago! + if (isRegular(grid)) { + aux <- vapply(1:n.mem, FUN.VALUE = numeric(le), FUN = function(x) { + z <- asub(grid[["Data"]], idx = x, dims = mem.ind, drop = TRUE) + z <- unname(abind(z, along = -1L)) + attr(z, "dimensions") <- c("time", "lat", "lon") + array3Dto2Dmat(z) + }) + # Data reordering to match SpatialGrid coordinates + aux <- data.frame(aux[order(-co[,2], co[,1]), ]) + } else { + aux <- redim(grid, loc = !isRegular(grid), drop = TRUE)$Data + if (n.mem > 1) { + naind <- lapply(1:n.mem, function(i) which(!is.na(aux[i,]), arr.ind = TRUE)) + naind <- Reduce(intersect, naind) + aux <- data.frame(t(aux[,naind])) + } else { + naind <- which(!is.na(aux), arr.ind = TRUE) + aux <- data.frame(as.numeric(aux[naind])) + } + } + # Set min/max values, if provided + if (!is.null(set.max)) aux[aux > set.max] <- set.max + if (!is.null(set.min)) aux[aux < set.min] <- set.min + # Panel names + if (is.multigrid) { + vname <- attr(grid$Variable, "longname") + if (!is.null(grid$Variable$level)) { + auxstr <- paste(vname, grid$Variable$level, sep = "@") + vname <- gsub("@NA", "", auxstr) + } + vname <- gsub("\\s", "_", vname) + vname <- make.names(vname, unique = TRUE) + } else { + vname <- paste0("Member_", 1:n.mem) + } + names(aux) <- vname + # Defining grid topology ----------------- + aux.grid <- getGrid(grid) + if (!isRegular(grid)) { + df <- sp::SpatialPointsDataFrame(co[naind,], aux) + } else { + cellcentre.offset <- vapply(aux.grid, FUN = "[", 1L, FUN.VALUE = numeric(1L)) + cellsize <- vapply(c("resX", "resY"), FUN.VALUE = numeric(1L), FUN = function(x) attr(aux.grid, which = x)) + aux.grid <- getCoordinates(grid) + cells.dim <- vapply(aux.grid, FUN.VALUE = integer(1L), FUN = "length") + grd <- sp::GridTopology(c(cellcentre.offset[["x"]], cellcentre.offset[["y"]]), cellsize, c(cells.dim[["x"]], cells.dim[["y"]])) + df <- sp::SpatialGridDataFrame(grd, aux) + } + return(df) +} + + + + +#' @title Climatological map stippling +#' @description Create a points panel layout to add to \code{spatialPlot}. +#' Typically needed to stipple significant points in climatologies, or other types of coordinates +#' @param clim A climatology. This can be for instance a verification climatology as produced by \code{\link{easyVeri2grid}}. +#' This often contains p-values, but not necessarily. +#' @param condition Inequality operator to be applied considering the given \code{threshold}. +#' \code{"GT"} = greater than the value of \code{threshold}, \code{"GE"} = greater or equal, +#' \code{"LT"} = lower than, \code{"LE"} = lower or equal than. Default to \code{"LT"} (see the rationale in the next argument). +#' @param threshold Reference threshold value to specify stippling points. Default to \code{0.05}, +#' in combination with \code{condition = "LT"}, as tipically used for stippling statistically significant +#' values (p-values). +#' @param ... Further optional style arguments (see the examples). +#' @return A list with a \code{SpatialPoints} object, +#' along with optional style arguments like \code{col}, \code{pch}, \code{cex} etc., +#' to be passed to the \code{sp.layout} argument in \code{spatialPlot}. +#' @export +#' @importFrom sp SpatialPoints +#' @details The function generates a \code{"sp.points"} layout list. Further formatting arguments can be passed here. +#' For further details and examples see the help of \code{\link[sp]{spplot}}. +#' @seealso \code{\link{spatialPlot}}, to which its output is passed. +#' \code{\link{map.lines}}, for further map customizations. +#' @author J. Bedia +#' @examples +#' data("CFS_Iberia_tas") +#' p90clim <- climatology(CFS_Iberia_tas, +#' by.member = FALSE, +#' clim.fun = list("FUN" = quantile, prob = .9)) +#' spatialPlot(p90clim, backdrop.theme = "coastline", +#' main = "CFSv2 Ensemble mean Tmean 90th percentile (July 2001)") +#' +#' # We want to highlight the grid points with a 90th percentile > 25.5 degrees, +#' # on top of the Tmean model climatology: +#' pts <- map.stippling(p90clim, threshold = 15.5, condition = "GT") +#' spatialPlot(climatology(CFS_Iberia_tas), +#' backdrop.theme = "coastline", +#' sp.layout = list(pts)) +#' +#' # Some useful parameters that can be passed to the layout list: +#' pts <- map.stippling(p90clim, threshold = 15.5, condition = "GT", +#' pch = 19, # dots instead of default crosses +#' col = "black", # black dots +#' cex = .1) # point expansion factor (to make them very small) +#' spatialPlot(climatology(CFS_Iberia_tas), +#' backdrop.theme = "coastline", +#' sp.layout = list(pts)) +#' +#' # Suppose we want the stippling just in the first and fifth panels, for instance: +#' pts <- map.stippling(p90clim, threshold = 15.5, condition = "GT", +#' pch = 19, col = "black", cex = .1, +#' which = c(1, 5)) # which controls in which panel(s) the points are displayed +#' spatialPlot(climatology(CFS_Iberia_tas), +#' backdrop.theme = "coastline", +#' sp.layout = list(pts)) + + +map.stippling <- function(clim, threshold = 0.05, condition = "LT", ...) { + if (!("climatology:fun" %in% names(attributes(clim$Data)))) { + stop("Input grid was not recognized as a climatology") + } + condition <- match.arg(condition, choices = c("GT", "GE", "LT", "LE")) + ineq <- switch(condition, + "GT" = ">", + "GE" = ">=", + "LT" = "<", + "LE" = "<=") + arg.list <- list(...) + aux <- array3Dto2Dmat(clim$Data)[1, ] + ind <- eval(parse(text = paste("which(aux", ineq, "threshold)"))) + coords <- as.matrix(expand.grid(clim$xyCoords$y, clim$xyCoords$x)[2:1][ind, ]) + if (nrow(coords) == 0) stop("None of the grid points is below the specified threshold") + do.call("list", c("sp.points", SpatialPoints(coords), arg.list)) +} + + +#' @title Add lines and polygons to climatological maps +#' @description Draws user-defined lines or polygons on top of climatological maps. +#' @param lonLim A numeric vector of length 2, with minimum and maximum longitude coordinates (in grid units), +#' of the rectangle to be drawn. +#' @param latLim Same as \code{lonLim}, but for the selection of the latitudinal range. +#' @param coords Optional. 2-column numeric matrix with vertex coordinates (1 point by row). +#' Note that row order matters in order to define the drawing direction. +#' Also bear in mind that first point (row) should equal last coordinates (row) in the case of closed polygonal areas. +#' @param ... Further optional style arguments (see the examples). +#' @return A list with a \code{SpatialLines} object, +#' along with optional style arguments like \code{col}, \code{lty}, \code{lwd} etc., +#' to be passed to the \code{sp.layout} argument in \code{spatialPlot} +#' @details The function internally transforms the inputs into \code{\link[sp]{Line}} class objects, +#' so the displayed outputs are not actually polygons in a formal sense (they can not be filled, for instance). +#' The purpose of the function is just to highlight specific areas within climatological maps (typically rectangular +#' windows, but any other shapes like for instance storm tracks can be flexibly specified using \code{coords}). +#' @author J Bedia +#' @importFrom magrittr %>% +#' @importFrom sp Line Lines SpatialLines +#' @export +#' @seealso \code{\link{spatialPlot}}, to which its output is passed. +#' \code{\link{map.stippling}}, for further map customizations. +#' @examples +#' data("CFS_Iberia_tas") +#' # Define a rectangular window centered on the Iberian Peninsula +#' iberia <- map.lines(lonLim = c(-10,3.5), latLim = c(36,43)) +#' spatialPlot(climatology(CFS_Iberia_tas), backdrop.theme = "coastline", +#' sp.layout = list(iberia)) +#' +#' # Some customization options (awful, yes, but just for illustration): +#' iberia <- map.lines(lonLim = c(-10,3.5), latLim = c(36,44), +#' lwd = 3, # line width +#' col = "purple", # line color +#' lty = 2) # line type +#' spatialPlot(climatology(CFS_Iberia_tas, by.member = FALSE), backdrop.theme = "coastline", +#' sp.layout = list(iberia)) +#' +#' # Another window over the Alps +#' alps <- map.lines(lonLim = c(4,16), latLim = c(45,49), +#' lwd = 3, +#' col = "red") +#' spatialPlot(climatology(CFS_Iberia_tas, by.member = FALSE), backdrop.theme = "coastline", +#' sp.layout = list(iberia, alps)) +#' +#' \dontrun{ +#' # Adding a line (real data of a storm-track imported from a csv file) +#' # Source: http://www.europeanwindstorms.org/ +#' +#' # Requires the target server to be operative... +#' dat <- url("http://www.europeanwindstorms.org/repository/Jeanette/Jeanette_track.csv") +#' custom.coords <- read.csv(dat, header = FALSE)[ ,5:4] +#' storm <- map.lines(coords = custom.coords, +#' lwd = 3, +#' col = "red") +#' spatialPlot(climatology(CFS_Iberia_tas, by.member = FALSE), backdrop.theme = "coastline", +#' sp.layout = list(storm), # Add storm track +#' scales = list(draw = TRUE)) # Add coordinate axes +#' } + +map.lines <- function(lonLim = NULL, latLim = NULL, coords = NULL, ...) { + if (is.null(lonLim) && is.null(latLim) && is.null(coords)) stop("Undefined polygon coordinates") + if (is.null(lonLim) && !is.null(latLim) || !is.null(lonLim) && is.null(latLim)) stop("Undefined 'lonLim' or 'latLim'") + arg.list <- list(...) + spLines <- if (is.null(coords)) { + coords <- as.matrix(expand.grid(lonLim, latLim)) + coords[1:2, ] <- coords[2:1, ] + coords <- rbind(coords, coords[1,]) + Line(coords) %>% list() %>% Lines(ID = "id") %>% list() %>% SpatialLines() + } else { + lapply(1:(nrow(coords) - 1), function(x) { + Line(coords[x:(x + 1),]) + }) %>% Lines(ID = "id") %>% list() %>% SpatialLines() + } + do.call("list", c("sp.lines", spLines, arg.list)) +} diff --git a/R/spreadPlot.R b/R/spreadPlot.R index 8c92cb2..203d86d 100644 --- a/R/spreadPlot.R +++ b/R/spreadPlot.R @@ -36,6 +36,7 @@ #' @importFrom transformeR array3Dto2Dmat mat2Dto3Darray subsetGrid #' @importFrom stats filter #' @importFrom graphics lines +#' @importFrom stats median #' @import sm #' @import vioplot #' diff --git a/R/temporalPlot.R b/R/temporalPlot.R new file mode 100644 index 0000000..1fc6367 --- /dev/null +++ b/R/temporalPlot.R @@ -0,0 +1,199 @@ +# temporalPlot.R Lattice plot methods for climatological grids +# +# Copyright (C) 2017 Santander Meteorology Group (http://www.meteo.unican.es) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + + + +#' @title Lattice plot methods for climatological series +#' @description A wrapper for the lattice (trellis) plot methods for grid and station data. +#' @param ... Input grids (or station data) +#' @param aggr.spatial list containing the function and the related arguments to perform spatial +#' aggregation. Default is \code{list(FUN = mean, na.rm = TRUE)}. +#' @param cols Character with colors. +#' @param lwd Numeric for line width. +#' @param lty Numeric for line type. +#' @param missing.dates Logical. Not implemented (see Details). +#' @param show.na Logical. Implemented but under development (see Details). +#' @param xyplot.custom List of arguments as passed to xyplot. Argument \code{panel} cannot be modified, thus, +#' if specified, it will be ignored. +#' @details The function applies the \code{\link[lattice]{xyplot}} method after computing spatial aggregation +#' (parameter \code{aggr.spatial}) and member aggregation (mean and range) to the imput grids (or station data). +#' +#' In case of multimember grids, the function will internally compute the ensemble mean +#' and the range for plotting. The range is used to plot the shadow of the multimember spread. +#' +#' \strong{Controlling graphical parameters} +#' +#' Many different aspects of the plot can be controlled passing the relevant arguments to +#' \code{xyplot}. Fine control of graphical parameters for the trellis display can +#' be also controlled using \code{\link[lattice]{trellis.par.set}}. +#' +#' \strong{Future Work} +#' +#' Implement parameters missing.dates (default will be TRUE) to suppress +#' dates without values from xaxis when FALSE. +#' Implement parameter show.na (default will be FALSE) to fill with gray +#' NA values when TRUE. +#' Implement auxiliary functions to reduce dependencies and remove from +#' imports packages \code{data.table}, \code{padr} and \code{grDevices}. +#' +#' +# +#' @return A lattice plot of class \dQuote{trellis}. +#' +#' +#'@author M. Iturbide +#'@export +#'@importFrom lattice xyplot panel.grid panel.polygon panel.xyplot +#'@importFrom grDevices col2rgb rgb colors +#'@importFrom padr pad +#'@importFrom data.table rleid +#'@import latticeExtra +#' @examples +#' data("CFS_Iberia_pr") +#' data("EOBS_Iberia_pr") +#' data("VALUE_Iberia_pr") +#' # Combine grids with members (CFS) and without members (EOBS) +#' a <- subsetGrid(CFS_Iberia_tas, years = 1985:1992) +#' b <- subsetGrid(EOBS_Iberia_tas, years = 1985:1992) +#' temporalPlot("EOBS" = b, "CFS" = a, +#' xyplot.custom = list(main = "winter temperature", ylab = "Celsius")) +#' # Station and grid data can be combined, also different temporal ranges +#' v <- subsetGrid(VALUE_Iberia_tas, years = 1988:1990) +#' temporalPlot("VALUE" = v, "EOBS" = b, "CFS" = a, lwd = 0.9, +#' aggr.spatial = list(FUN = min, na.rm = TRUE), +#' xyplot.custom = list(main = "winter temperature", +#' ylab = "Celsius", ylim = c(-20, 10))) +#' # Use subsetGrid to obtain and plot a single location (no spatial aggregation) +#' a1 <- subsetGrid(a, lonLim = 2, latLim = 42) +#' b1 <- subsetGrid(b, lonLim = 2, latLim = 42) +#' +#' temporalPlot("EOBS" = b1, "CFS" = a1, +#' cols = c("green", "deeppink"), show.na = TRUE, +#' xyplot.custom = list(main = "winter temperature", ylab = "Celsius")) + + +temporalPlot <- function(..., + aggr.spatial = list(FUN = mean, na.rm = TRUE), + cols = NULL, + lwd = 1, + lty = 1, + missing.dates = TRUE, + show.na = FALSE, + xyplot.custom = list()) { + obj.list <- list(...) + if (is.null(names(obj.list))) { + nmes <- as.list(substitute(list(...)))[-1L] + names(obj.list) <- as.character(nmes) + } + obj.list <- lapply(obj.list, FUN = redim) + # spatial aggregation + aggr.spatial[["MARGIN"]] <- c(1,2) + data <- lapply(1:length(obj.list), function(x){ + aggr.spatial[["X"]] <- obj.list[[x]]$Data + do.call("apply", aggr.spatial) + }) + # extract dates + dates <- lapply(1:length(obj.list), function(x){ + as.Date(obj.list[[x]]$Dates$start) + }) + # member aggregation + mm <- lapply(data, FUN = apply, MARGIN = 2, mean) + mx <- lapply(data, FUN = apply, MARGIN = 2, max) + mn <- lapply(data, FUN = apply, MARGIN = 2, min) + # bind to data frames + df <- lapply(1:length(obj.list), function(x){ + df0 <- cbind(as.data.frame(mm[[x]]), as.data.frame(dates[[x]]), + as.data.frame(mn[[x]]), as.data.frame(mx[[x]])) + colnames(df0) <- c("Value", "Dates", "mini", "maxi") + return(df0) + }) + # complete temporal series + df <- lapply(df, FUN = pad) # uses package padr to fill dates in the data.frame. Should we implement alternative code? + # prepare inter-NA chunks for panel.polygon + df.polys <- lapply(1:length(df), function(i){ + chunkid <- rle(!is.na(df[[i]]$Value))$values + split(df[[i]], rleid(is.na(df[[i]]$Value)))[chunkid] + }) + df.polys.na <- lapply(1:length(df), function(i){ + chunkid <- rle(is.na(df[[i]]$Value))$values + split(df[[i]], rleid(is.na(df[[i]]$Value)))[chunkid] + }) + # define graphical parameters + ylim <- round(range(c(unlist(mm), unlist(mx)), na.rm = TRUE), digits = 2) + xlim <- range(do.call("c", dates)) + # xdates0 <- unique(do.call("c", dates))[order(unique(do.call("c", dates)))] + # seqval <- round(length(xdates0)/10) + # xdates <- xdates0[seq(1, length(xdates0), seqval)] + colors2 <- colors()[-c(552, 254, 26, 24)] + colors2 <- colors2[sample(1:length(colors2), size = length(colors2))] + if (is.null(cols)) cols <- c("black","red", "blue", "green", colors2) + if (length(cols) < length(obj.list)) stop("Please, add ", length(obj.list) - length(cols), " more color/s to 'cols', or keep the default option.") + if (is.null(xyplot.custom[["x"]])) xyplot.custom[["x"]] <- Value ~ Dates + if (is.null(xyplot.custom[["type"]])) xyplot.custom[["type"]] = "l" + if (is.null(xyplot.custom[["ylim"]])) xyplot.custom[["ylim"]] <- ylim + if (is.null(xyplot.custom[["xlim"]])) xyplot.custom[["xlim"]] <- xlim + if (is.null(xyplot.custom[["lwd"]])) xyplot.custom[["lwd"]] <- 2 + ylim <- xyplot.custom[["ylim"]] + xlim <- xyplot.custom[["xlim"]] + if (is.null(xyplot.custom[["scales"]])) xyplot.custom[["scales"]] <- list(x = list(at = seq(xlim[1], xlim[2],(xlim[2] - xlim[1])/10), + labels = seq(xlim[1], xlim[2],(xlim[2] - xlim[1])/10), rot = 45), + y = list(at = seq(ylim[1], ylim[2],round((ylim[2] - ylim[1])/10)), + labels = seq(ylim[1], ylim[2],round((ylim[2] - ylim[1])/10))), + cex = .6, col = "black") + if (is.null(xyplot.custom[["key"]])) xyplot.custom[["key"]] <- list(space = "right", points = list(pch = 15, + col = cols[1:length(obj.list)], + cex = .5), + text = list(names(obj.list), cex = .8)) + # crate trellis objects + xy <- lapply(1:length(df), function(i){ + col <- cols[i] + colsrgb <- do.call("rgb", as.list(c(col2rgb(col)/255, 0.15))) + xyplot.custom[["data"]] <- df[[i]] + xyplot.custom[["col"]] <- col + xyplot.custom[["panel"]] <- function(x, y, z, ...) { + for (l in 1:length(df.polys[[i]])) { + panel.polygon(na.omit(c(df.polys[[i]][[l]]$Dates, rev(df.polys[[i]][[l]]$Dates))), na.omit(c(df.polys[[i]][[l]]$mini, rev(df.polys[[i]][[l]]$maxi))), + border = NA, col = colsrgb) + } + if (show.na) { + for (l in 1:length(df.polys.na[[i]])) { + panel.polygon(na.omit(c(df.polys.na[[i]][[l]]$Dates, rev(df.polys.na[[i]][[l]]$Dates))), + na.omit(c(ylim[1], ylim[2])), + border = NA, col = "gray90") + } + } + panel.xyplot(df[[i]]$Dates, df[[i]]$Value, type = "l", lwd = lwd, lty = lty, col = col) + panel.abline(h = seq(ylim[1], ylim[2],round((ylim[2] - ylim[1])/10)), + v = seq(xlim[1], xlim[2],(xlim[2] - xlim[1])/10), + col = "gray65", lwd = 0.5, lty = 2) + } + + do.call("xyplot", xyplot.custom) + }) + # evaluate trellis objects + p0 <- lapply(1:length(obj.list), function(x){ + if (x == length(obj.list)) { + paste0("xy[[", x, "]]") + } else { + paste0("xy[[", x, "]] +") + } + }) + output <- eval(parse(text = do.call("paste", p0))) + return(output) +} + +#end diff --git a/R/tercileMap.R b/R/tercileMap.R new file mode 100644 index 0000000..f359b02 --- /dev/null +++ b/R/tercileMap.R @@ -0,0 +1,117 @@ +# tercileMap.R Tercile maps for seasonal forecasts based on C3S style +# +# Copyright (C) 2017 Santander Meteorology Group (http://www.meteo.unican.es) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + + + +#' @title Tercile map for visualization of forecast skill of seasonal climate predictions. +#' +#' @description Tercile map for the visualization of forecast skill of seasonal climate predictions. +#' This function is prepared to plot the data sets loaded from the ECOMS User Data Gateway (ECOMS-UDG). See +#' the loadeR.ECOMS R package for more details (http://meteo.unican.es/trac/wiki/udg/ecoms/RPackage). +#' +#' @param hindcast A multi-member list with the hindcast for verification. See details. +#' @param forecast A multi-member list with the forecasts. Default is NULL. +#' @param ... Optional arguments passed to \code{\link{spatialPlot}} +#' +#' @details The function applies the \code{\link{spatialPlot}} function, that in turn uses \code{lattice-methods}. +#' +#' @note The function is based on the tercile maps delivered by +#' C3S seasonal forecasts (\url{http://climate.copernicus.eu/s/charts/c3s_seasonal/}). +#' +#' Many different aspects of the plot can be controlled passing the relevant arguments to +#' \code{\link[sp]{spplot}}. +#' +#' @return A lattice plot of class \dQuote{trellis}. +#' +#' +#'@author M. Iturbide +#'@export +#'@importFrom grDevices colorRampPalette +#'@importFrom transformeR aggregateGrid subsetGrid getShape redim climatology +#' @examples +#' library(transformeR) +#' hindcast <- subsetGrid(CFS_Iberia_tas, years = 1983:2001) +#' forecast <- subsetGrid(CFS_Iberia_tas, years = 2002) +#' tercileMap(hindcast, forecast) + +tercileMap <- function(hindcast, forecast, ...) { + custom <- list(...) + message("[", Sys.time(), "] Preparing data...") + suppressMessages(hindyear <- aggregateGrid(hindcast, aggr.y = list(FUN = mean, na.rm = TRUE))) + suppressMessages(foreyear <- aggregateGrid(forecast, aggr.y = list(FUN = mean, na.rm = TRUE))) + nmem <- getShape(hindyear)["member"] + # hind <- redim(downscaleR:::flatMemberDim(hindyear), drop = TRUE) + hind <- lapply(1:nmem, function(x){ + subsetGrid(hindyear, members = x) + }) + hind <- do.call("bindGrid.time", hind) + hind <- redim(hind, drop = TRUE) + message("[", Sys.time(), "] Calculating terciles...") + a <- apply(hind$Data, MARGIN = c(2,3), FUN = quantile, probs = c(1/3, 2/3), na.rm = TRUE) + nmem <- getShape(forecast)["member"] + nlon <- getShape(forecast)["lon"] + nlat <- getShape(forecast)["lat"] + co <- expand.grid(1:nlat, 1:nlon) + indtercile <- array(dim = c(nmem, nlat, nlon)) + for (m in 1:nmem) { + formem <- subsetGrid(foreyear, members = m)$Data + for (i in 1:nrow(co)) { + indtercile[m, co[i,1],co[i,2]] <- which(order(c(formem[co[i,1],co[i,2]], a[,co[i,1],co[i,2]])) == 1) + } + } + message("[", Sys.time(), "] Calculating probabilities...") + probtercile <- array(dim = c(nlat, nlon)) + for (i in 1:nrow(co)) { + z <- table(indtercile[,co[i,1],co[i,2]]) + tercile <- which(z == max(z))[1] + probtercile[co[i,1],co[i,2]] <- round(z[tercile]/sum(z)*100, digits = 2) + if (tercile == 2) { + probtercile[co[i,1],co[i,2]] <- 0 + } else if (tercile == 1) { + probtercile[co[i,1],co[i,2]] <- (-1)*probtercile[co[i,1],co[i,2]] + } + } + suppressMessages(output <- aggregateGrid(foreyear, aggr.mem = list(FUN = mean, na.rm = TRUE))) + output$Data <- probtercile + attr(output$Data, "dimensions") <- c("lat", "lon") + jet.colors <- c("#00007F", "blue", "#007FFF", "cyan", "white", "yellow", "#FF7F00", "red", "red3") + suppressMessages( + suppressWarnings(custom[["grid"]] <- climatology(redim(output, member = FALSE))) + ) + if (is.null(custom[["backdrop.theme"]])) custom[["backdrop.theme"]] <- "coastline" + if (is.null(custom[["at"]])) custom[["at"]] <- c(-100, -70, -60, -50, -40, 40, 50, 60, 70, 100) + if (is.null(custom[["col.regions"]])) custom[["col.regions"]] <- jet.colors + if (is.null(custom[["colorkey"]])) custom[["colorkey"]] <- FALSE + if (is.null(custom[["scales"]])) custom[["scales"]] <- list(draw = TRUE, alternating = 3, cex = .6) + if (is.null(custom[["key"]])) custom[["key"]] <- list(space = "top", between = 0.3, between.columns = 0.5, + columns = 9, points = list( pch = 22, + col = "black", + fill = custom[["col.regions"]], + cex = 1.3), + text = list(c("100..70%", "70..60%", "60..50%", "50..40%", + "other", + "40..50%", "50..60%", "60..70%", "70..100%"), cex = .65)) + + message("[", Sys.time(), "] Done.") + suppressWarnings(pl <- do.call("spatialPlot", custom)) + return(pl) +} + +# end + + + diff --git a/README.md b/README.md index 41a32ec..a64cc26 100644 --- a/README.md +++ b/README.md @@ -1,16 +1,19 @@ # What is visualizeR? -`visualizeR` is an R package implementing a set of advanced visualization tools for forecast verification. It is fully integrated (yet independent) from the R climate data structures generated by the loading functions of the `climate4R` bundle, thus providing seamless integration with all steps of forecast data analysis, from data loading to the ECOMS-UDG (loadeR and loadeR.ECOMS packages, Cofiño et al. 2017), to data post-processing (transformeR package, Bedia and Iturbide 2016), downscaling and bias correction (downscaleR package, Bedia et al. 2017) and verification and visualization (visualizeR package, Frias et al. 2017). visualizeR only depends on the [transformeR](https://github.com/SantanderMetGroup/transformeR) capabilities for climata data manipulation. +`visualizeR` is an R package for climate data visualization, with special focus on ensemble forecasting and uncertainty communication. It includes functions for visualizing climatological, forecast and evaluation products, and combinations of them. Find out more about this package at the [visualizeR wiki](https://github.com/SantanderMetGroup/visualizeR/wiki).  -The recommended installation procedure is to use the `install_github` command from the devtools R package. +This package is part of the [climate4R bundle](http://www.meteo.unican.es/climate4r), formed by `loadeR`, `transformeR`, `downscaleR` and `visualizeR`. + +The recommended installation procedure is to use the `install_github` command from the devtools R package (see the [installation info](https://github.com/SantanderMetGroup/visualizeR/wiki/installation) in the wiki). ```r devtools::install_github(c("SantanderMetGroup/transformeR", "SantanderMetGroup/visualizeR")) ``` -Alternatively, you can download the sources from the [Releases Tab](https://github.com/SantanderMetGroup/visualizeR/releases) +**IMPORTANT:** Note that `transformeR` must be previously installed on your system +--- +Reference and further information: -**IMPORTANT:** Note that `transformeR` must be previously installed on your system ([installation info](https://github.com/SantanderMetGroup/downscaleR/wiki/installation)) +Frías M.D., Iturbide M., Manzanas R., Bedia J., Fernández J., Herrera S., Cofiño A.S., Gutiérrez J.M. (2018) An R package to visualize and communicate uncertainty in seasonal climate prediction. **Environmental Modelling and Software**, 99, 101-110, http://doi.org/10.1016/j.envsoft.2017.09.008 ---- -To know more about the `climate4R` bundle and its components, read the paper: Cofiño et al, The ECOMS User Data Gateway: Towards seasonal forecast data provision and research reproducibility in the era of Climate Services. Climate Services, 2017, ISSN 2405-8807, http://dx.doi.org/10.1016/j.cliser.2017.07.001. +Cofiño et al. (2018) The ECOMS User Data Gateway: Towards seasonal forecast data provision and research reproducibility in the era of Climate Services. **Climate Services**, http://dx.doi.org/10.1016/j.cliser.2017.07.001. diff --git a/inst/coastline.rda b/inst/coastline.rda new file mode 100644 index 0000000..de5aab8 Binary files /dev/null and b/inst/coastline.rda differ diff --git a/inst/countries.rda b/inst/countries.rda new file mode 100644 index 0000000..c5a6248 Binary files /dev/null and b/inst/countries.rda differ diff --git a/man/clim2sgdf.Rd b/man/clim2sgdf.Rd new file mode 100644 index 0000000..5674383 --- /dev/null +++ b/man/clim2sgdf.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatialPlot.R +\name{clim2sgdf} +\alias{clim2sgdf} +\title{Climatology to SpatialGridDataFrame or SpatialPointsDataFrame} +\usage{ +clim2sgdf(clim, set.min, set.max) +} +\arguments{ +\item{clim}{A climatological grid, as returned by function \code{\link{climatology}}} + +\item{set.min}{Minimum value, as passed by \code{spatialPlot}} + +\item{set.max}{Maximum value, as passed by \code{spatialPlot}} +} +\value{ +A \pkg{sp} object of the class \code{\link[sp]{SpatialGridDataFrame}} +} +\description{ +Convert a climatological grid to a SpatialGridDataFrame object from package sp +} +\details{ +This function is intended for internal usage by \code{\link{spatialPlot}}, +that accepts all possible arguments that can be passed to \code{\link[sp]{spplot}} for plotting. +However, it may be useful for advanced \pkg{sp} users in different contexts +(e.g. for reprojecting via \code{\link[sp]{spTransform}} etc.) +} +\examples{ +data("CFS_Iberia_tas") +# Climatology is computed: +clim <- climatology(CFS_Iberia_tas, by.member = TRUE) +sgdf <- clim2sgdf(clim, NULL, NULL) +class(sgdf) +} +\seealso{ +\code{\link{climatology}}, \code{\link{spatialPlot}} +} +\author{ +J. Bedia +} +\keyword{internal} diff --git a/man/map.lines.Rd b/man/map.lines.Rd new file mode 100644 index 0000000..3e38198 --- /dev/null +++ b/man/map.lines.Rd @@ -0,0 +1,78 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatialPlot.R +\name{map.lines} +\alias{map.lines} +\title{Add lines and polygons to climatological maps} +\usage{ +map.lines(lonLim = NULL, latLim = NULL, coords = NULL, ...) +} +\arguments{ +\item{lonLim}{A numeric vector of length 2, with minimum and maximum longitude coordinates (in grid units), +of the rectangle to be drawn.} + +\item{latLim}{Same as \code{lonLim}, but for the selection of the latitudinal range.} + +\item{coords}{Optional. 2-column numeric matrix with vertex coordinates (1 point by row). +Note that row order matters in order to define the drawing direction. + Also bear in mind that first point (row) should equal last coordinates (row) in the case of closed polygonal areas.} + +\item{...}{Further optional style arguments (see the examples).} +} +\value{ +A list with a \code{SpatialLines} object, +along with optional style arguments like \code{col}, \code{lty}, \code{lwd} etc., +to be passed to the \code{sp.layout} argument in \code{spatialPlot} +} +\description{ +Draws user-defined lines or polygons on top of climatological maps. +} +\details{ +The function internally transforms the inputs into \code{\link[sp]{Line}} class objects, + so the displayed outputs are not actually polygons in a formal sense (they can not be filled, for instance). + The purpose of the function is just to highlight specific areas within climatological maps (typically rectangular + windows, but any other shapes like for instance storm tracks can be flexibly specified using \code{coords}). +} +\examples{ +data("CFS_Iberia_tas") +# Define a rectangular window centered on the Iberian Peninsula +iberia <- map.lines(lonLim = c(-10,3.5), latLim = c(36,43)) +spatialPlot(climatology(CFS_Iberia_tas), backdrop.theme = "coastline", + sp.layout = list(iberia)) + +# Some customization options (awful, yes, but just for illustration): +iberia <- map.lines(lonLim = c(-10,3.5), latLim = c(36,44), + lwd = 3, # line width + col = "purple", # line color + lty = 2) # line type +spatialPlot(climatology(CFS_Iberia_tas, by.member = FALSE), backdrop.theme = "coastline", + sp.layout = list(iberia)) + +# Another window over the Alps +alps <- map.lines(lonLim = c(4,16), latLim = c(45,49), + lwd = 3, + col = "red") +spatialPlot(climatology(CFS_Iberia_tas, by.member = FALSE), backdrop.theme = "coastline", + sp.layout = list(iberia, alps)) + +\dontrun{ +# Adding a line (real data of a storm-track imported from a csv file) +# Source: http://www.europeanwindstorms.org/ + +# Requires the target server to be operative... +dat <- url("http://www.europeanwindstorms.org/repository/Jeanette/Jeanette_track.csv") +custom.coords <- read.csv(dat, header = FALSE)[ ,5:4] +storm <- map.lines(coords = custom.coords, + lwd = 3, + col = "red") +spatialPlot(climatology(CFS_Iberia_tas, by.member = FALSE), backdrop.theme = "coastline", + sp.layout = list(storm), # Add storm track + scales = list(draw = TRUE)) # Add coordinate axes +} +} +\seealso{ +\code{\link{spatialPlot}}, to which its output is passed. + \code{\link{map.stippling}}, for further map customizations. +} +\author{ +J Bedia +} diff --git a/man/map.stippling.Rd b/man/map.stippling.Rd new file mode 100644 index 0000000..8263a56 --- /dev/null +++ b/man/map.stippling.Rd @@ -0,0 +1,74 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatialPlot.R +\name{map.stippling} +\alias{map.stippling} +\title{Climatological map stippling} +\usage{ +map.stippling(clim, threshold = 0.05, condition = "LT", ...) +} +\arguments{ +\item{clim}{A climatology. This can be for instance a verification climatology as produced by \code{\link{easyVeri2grid}}. +This often contains p-values, but not necessarily.} + +\item{threshold}{Reference threshold value to specify stippling points. Default to \code{0.05}, +in combination with \code{condition = "LT"}, as tipically used for stippling statistically significant +values (p-values).} + +\item{condition}{Inequality operator to be applied considering the given \code{threshold}. +\code{"GT"} = greater than the value of \code{threshold}, \code{"GE"} = greater or equal, + \code{"LT"} = lower than, \code{"LE"} = lower or equal than. Default to \code{"LT"} (see the rationale in the next argument).} + +\item{...}{Further optional style arguments (see the examples).} +} +\value{ +A list with a \code{SpatialPoints} object, +along with optional style arguments like \code{col}, \code{pch}, \code{cex} etc., +to be passed to the \code{sp.layout} argument in \code{spatialPlot}. +} +\description{ +Create a points panel layout to add to \code{spatialPlot}. +Typically needed to stipple significant points in climatologies, or other types of coordinates +} +\details{ +The function generates a \code{"sp.points"} layout list. Further formatting arguments can be passed here. + For further details and examples see the help of \code{\link[sp]{spplot}}. +} +\examples{ +data("CFS_Iberia_tas") +p90clim <- climatology(CFS_Iberia_tas, + by.member = FALSE, + clim.fun = list("FUN" = quantile, prob = .9)) +spatialPlot(p90clim, backdrop.theme = "coastline", + main = "CFSv2 Ensemble mean Tmean 90th percentile (July 2001)") + +# We want to highlight the grid points with a 90th percentile > 25.5 degrees, +# on top of the Tmean model climatology: +pts <- map.stippling(p90clim, threshold = 15.5, condition = "GT") +spatialPlot(climatology(CFS_Iberia_tas), + backdrop.theme = "coastline", + sp.layout = list(pts)) + +# Some useful parameters that can be passed to the layout list: +pts <- map.stippling(p90clim, threshold = 15.5, condition = "GT", + pch = 19, # dots instead of default crosses + col = "black", # black dots + cex = .1) # point expansion factor (to make them very small) +spatialPlot(climatology(CFS_Iberia_tas), + backdrop.theme = "coastline", + sp.layout = list(pts)) + +# Suppose we want the stippling just in the first and fifth panels, for instance: +pts <- map.stippling(p90clim, threshold = 15.5, condition = "GT", + pch = 19, col = "black", cex = .1, + which = c(1, 5)) # which controls in which panel(s) the points are displayed +spatialPlot(climatology(CFS_Iberia_tas), + backdrop.theme = "coastline", + sp.layout = list(pts)) +} +\seealso{ +\code{\link{spatialPlot}}, to which its output is passed. + \code{\link{map.lines}}, for further map customizations. +} +\author{ +J. Bedia +} diff --git a/man/reliabilityCategories.Rd b/man/reliabilityCategories.Rd index 859311f..c62550d 100644 --- a/man/reliabilityCategories.Rd +++ b/man/reliabilityCategories.Rd @@ -6,8 +6,9 @@ \usage{ reliabilityCategories(hindcast, obs, regions = NULL, n.events = 3, labels = NULL, n.bins = 10, n.boot = 100, conf.level = 0.75, - diagrams = TRUE, cex0 = 0.5, cex.scale = 20, layout = c(1, n.events), - backdrop.theme = "countries", return.diagrams = FALSE) + na.rate = 0.75, diagrams = TRUE, cex0 = 0.5, cex.scale = 20, + layout = c(1, n.events), backdrop.theme = "countries", + return.diagrams = FALSE) } \arguments{ \item{hindcast}{Grid of forecast data} @@ -29,6 +30,9 @@ the relaiability is calculated. Default is NULL (See details).} \item{conf.level}{Confidence interval for the reliability line. By default \code{conf.level} = 0.75 (two sided), as in Weisheimer et al. 2014} +\item{na.rate}{Allowed proportion of NA values in each region. Regions with proportions higher than na.rate +are excuded from the analysis. Default is 0.75.} + \item{diagrams}{Logical (default = TRUE). Plotting results.} \item{cex0}{numeric (default is 0.5). Minimum size of the points shown in the reliability diagrams, i.e. size of the point diff --git a/man/skillMap.Rd b/man/skillMap.Rd new file mode 100644 index 0000000..687408c --- /dev/null +++ b/man/skillMap.Rd @@ -0,0 +1,113 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/skillMap.R +\name{skillMap} +\alias{skillMap} +\title{A wrapper of \code{spatialPlot} for the creation of verification maps for seasonal forecast systems} +\usage{ +skillMap(easyVeriGrid, stippling = list(threshold = NULL, condition = NULL), + stippling.point.options = NULL, backdrop.theme = NULL, title = NULL, + ...) +} +\arguments{ +\item{easyVeriGrid}{A climatological grid with a verification output. See details} + +\item{stippling}{A key-value list of arguments passed to \code{map.stippling}: +\code{threshold} and \code{condition}. Ignored by default, returning a map without stippling.} + +\item{stippling.point.options}{Default to \code{NULL}. Further graphical arguments passed to points +(e.g. \code{cex}, \code{pch} etc.)} + +\item{backdrop.theme}{See \code{\link{spatialPlot}}} + +\item{title}{Title of the plot} + +\item{...}{Further graphical options passed to \code{\link{spatialPlot}}.} +} +\value{ +A lattice plot of class \dQuote{trellis}. +} +\description{ +A wrapper of \code{spatialPlot} for the creation of verification maps for seasonal forecast systems. + It provides a convenient interface for \code{\link{map.stippling}} +} +\details{ +The function applies the \code{\link{spatialPlot}} function, that in turn uses \code{lattice-methods}. + +\strong{Graphical options} + +Some examples of specific map graphical options are available in the help of function \code{\link[sp]{spplot}}. + In addtion, fine-tuning of the resulting plots can be obtained using the arguments of \pkg{lattice} plots. For + an overview, see the help of function \code{\link[lattice]{xyplot}}. +} +\examples{ +\dontrun{ + +# The package 'easyVerification' will be used to calculate the RPSS: + +require(easyVerification) +require(transformeR) +# First of all, a data subset is done, considering a spatial domain centered on the North Atlantic: +tas.cfs2 <- subsetGrid(tas.cfs, lonLim = c(-100, 40), latLim = c(-5, 75)) +# The same is done with the reanalysis dataset: +tas.ncep2 <- subsetGrid(tas.ncep, lonLim = c(-100, 40), latLim = c(-5, 75)) +# In the next step, the reanalysis data are interpolated to the hindcast grid: +tas.ncep2.int <- interpGrid(tas.ncep2, new.coordinates = getGrid(tas.cfs2), method = "nearest") +# We compute the Ranked Probabiloity SKill Score using veriApply, and a cross-validation strategy: +ev <- easyVerification::veriApply(verifun = "EnsRpss", + fcst = tas.cfs2$Data, + obs = tas.ncep2.int$Data, + prob = 1:2/3, + tdim = 2, + ensdim = 1, + parallel = TRUE, + ncpus = 3, + strategy = "crossval") +# The bridging function 'easyVeri2grid' converts the object returned by 'veriApply' +# to a 'climate4R' climatological grid: + +easyVeriGrid <- easyVeri2grid(ev$skillscore, + obs.grid = tas.ncep2.int, + verifun = "EnsRpss") + +# A basic RPSS map, using the spatialPlot defaults: +skillMap(easyVeriGrid = easyVeriGrid, backdrop.theme = "coastline") + +# Stippling. Mark significant RPSS values at the 95\% ci +thresh <- ev$skillscore.sd*qnorm(0.95) + +skillMap(easyVeriGrid = easyVeriGrid, + stippling = list(threshold = th, condition = "GT"), # GT = greater than + stippling.point.options = list(pch = 19, cex = .2, col = "black"), + backdrop.theme = "coastline") + +# Further customization: For instance a more elaborated title, colorblind-friendly palette etc.: + +cb.colors <- colorRampPalette(rev(RColorBrewer::brewer.pal(11, "Spectral"))) +skillMap(easyVeriGrid = easyVeriGrid, + stippling = list(threshold = th, condition = "GT"), + stippling.point.options = list(pch = 19, cex = .2, col = "black"), + backdrop.theme = "coastline", + at = seq(-0.7, 0.7, .05), + set.max = 0.7, set.min = -0.7, + scales = list(draw = TRUE, alternating = 3), + colorkey = list(space = "bottom"), + col.regions = cb.colors(51), + title = paste("Ranked Probability Skill Score", + "Forecasting System: CFSv2 - 24 members", + "November Initializations", + "Observing System: NCEP/NCAR Reanalysis 1", + "t2m DJF 1983-2010", sep = "\\n") +) +} +} +\seealso{ +The bridging function \code{\link[transformeR]{easyVeri2grid}} from package \pkg{transformeR} allows for +the conversion of verification outputs from package \pkg{easyVerification} to the \code{climate4R} data structure +used by the function. + + Many different aspects of the plot can be controlled passing the relevant arguments to + \code{\link[sp]{spplot}}. +} +\author{ +J. Bedia +} diff --git a/man/spatialPlot.Rd b/man/spatialPlot.Rd new file mode 100644 index 0000000..4e27707 --- /dev/null +++ b/man/spatialPlot.Rd @@ -0,0 +1,136 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatialPlot.R +\name{spatialPlot} +\alias{spatialPlot} +\title{Lattice plot methods for climatological grids} +\usage{ +spatialPlot(grid, backdrop.theme = "none", set.min = NULL, set.max = NULL, + lonCenter = NULL, ...) +} +\arguments{ +\item{grid}{Input grid} + +\item{backdrop.theme}{Reference geographical lines to be added to the plot. See Details.} + +\item{set.min}{Numeric value indicating an absolute minimum value (default to \code{NULL}). All grid values below this are mapped to \code{set.min}. See details.} + +\item{set.max}{Same as \code{set.min} argument, but to force a ceiling.} + +\item{lonCenter}{Value of the longitude to be centered in the plot.} + +\item{...}{Further arguments passed to \code{spplot}} +} +\value{ +As spplot, \code{spatialPlot} returns a lattice plot of class \dQuote{trellis}. + If you fail to \dQuote{see} it, explicitly call \code{print(spatialPlot(...))}. +} +\description{ +A wrapper for the lattice (trellis) plot methods for spatial data in \code{sp::spplot} +} +\details{ +The function applies the \code{\link[sp]{spplot}} method after conversion of the climatological map(s) to a + \code{SpatialGridDataFrame}. + + The \code{set.min} and \code{set.max} options are useful in order to preserve adequate ranges for map representation, + avoiding the influence of extreme values. Note that this is different than setting a range of values with an + interval using the \code{at} argument. The latter choice, that overrides \code{set.min} and \code{set.max}, + leaves blank grid points for outlying values. + + \strong{Multigrids} + + Multigrids of climatologies can be created using \code{makeMultiGrid} + for trellis visualization of different variables, or for instance, for the comparison of + raw and corrected/downscaled scenarios side to side. In case of multimember multigrids, + the function will internally compute the ensemble mean of each variable in the multigrid + for representation (with a message). + + \strong{Backdrop theme} + + Current implemented options are \code{"none"} and \code{"coastline"}, which contains + a simplied vector theme delineating the world coastlines. Any other themes can be introduced + by the user using the \code{sp.layout} options in \code{spplot}. + + \strong{Controlling graphical parameters} + + Many different aspects of the map can be controlled passing the relevant arguments to + \code{spplot}. Fine control of graphical parameters for the trellis display can + be also controlled using \code{\link[lattice]{trellis.par.set}}. + +Some examples of specific map graphical options are available in the help of function \code{\link[sp]{spplot}}. + In addtion, fine-tuning of the resulting plots can be obtained using the arguments of \pkg{lattice} plots. For + an overview, see the help of function \code{\link[lattice]{xyplot}}. +} +\examples{ +require(transformeR) +data("CFS_Iberia_tas") +# Climatology is computed: +clim <- climatology(CFS_Iberia_tas, by.member = TRUE) +spatialPlot(clim) +# Geographical lines can be added using the argument 'backdrop.theme': +spatialPlot(clim, backdrop.theme = "coastline") +spatialPlot(clim, backdrop.theme = "coastline", center = 180) + +# Further arguments can be passed to 'spplot'... + +# ... a subset of members to be displayed, using 'zcol': +spatialPlot(clim, + backdrop.theme = "coastline", + zcol = 1:4) + +# ... regional focuses (e.g. the Iberian Peninsula). +spatialPlot(clim, + backdrop.theme = "countries", + xlim = c(-10,5), ylim = c(35,44), + zcol = 1:4, + scales = list(draw = TRUE)) + +# Changing the default color palette and ranges: +spatialPlot(clim, + backdrop.theme = "coastline", + zcol = 1:4, + col.regions = cm.colors(27), at = seq(10,37,1)) + +# For ensemble means climatology should be called with 'by.member' set to FALSE: +clim <- climatology(CFS_Iberia_tas, by.member = FALSE) + +# Adding contours to the plot is direct with argument 'contour': +spatialPlot(clim, + scales = list(draw = TRUE), + contour = TRUE, + main = "tas Predictions July Ensemble Mean") + +## Example of multigrid plotting +data("NCEP_Iberia_psl") +## Winter data are split into monthly climatologies +monthly.clim.grids <- lapply(getSeason(NCEP_Iberia_psl), function(x) { + climatology(subsetGrid(NCEP_Iberia_psl, season = x)) +}) +## Skip the temporal checks, as grids correspond to different time slices +mg <- do.call("makeMultiGrid", + c(monthly.clim.grids, skip.temporal.check = TRUE)) +## We change the panel names +spatialPlot(mg, + backdrop.theme = "coastline", + names.attr = c("DEC","JAN","FEB"), + main = "Mean PSL climatology 1991-2010", + scales = list(draw = TRUE)) + +# Station data: +data("VALUE_Iberia_pr") +spatialPlot(climatology(VALUE_Iberia_pr), backdrop.theme = "countries") +} +\references{ +\itemize{ +\item Bivand, R.S., Pebesma, E.J., Gomez-Rubio, V., 2013. Applied Spatial Data Analysis with R, 2nd ed, useR! Springer, NY. +\item For some graticulate customization examples, visit the \emph{sp Gallery}: \url{https://edzer.github.io/sp/} +} +} +\seealso{ +\code{\link{climatology}} for details on climatology calculation. + \code{\link{map.stippling}}, for adding a custom point layer on top of the map. + \code{\link{map.lines}}, to add lines and polygons to climatological maps +Also see \code{\link[sp]{spplot}} in package \pkg{sp} for further information on plotting capabilities and options +} +\author{ +J. Bedia +} diff --git a/man/temporalPlot.Rd b/man/temporalPlot.Rd new file mode 100644 index 0000000..536fc5e --- /dev/null +++ b/man/temporalPlot.Rd @@ -0,0 +1,83 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/temporalPlot.R +\name{temporalPlot} +\alias{temporalPlot} +\title{Lattice plot methods for climatological series} +\usage{ +temporalPlot(..., aggr.spatial = list(FUN = mean, na.rm = TRUE), + cols = NULL, lwd = 1, lty = 1, missing.dates = TRUE, + show.na = FALSE, xyplot.custom = list()) +} +\arguments{ +\item{...}{Input grids (or station data)} + +\item{aggr.spatial}{list containing the function and the related arguments to perform spatial +aggregation. Default is \code{list(FUN = mean, na.rm = TRUE)}.} + +\item{cols}{Character with colors.} + +\item{lwd}{Numeric for line width.} + +\item{lty}{Numeric for line type.} + +\item{missing.dates}{Logical. Not implemented (see Details).} + +\item{show.na}{Logical. Implemented but under development (see Details).} + +\item{xyplot.custom}{List of arguments as passed to xyplot. Argument \code{panel} cannot be modified, thus, +if specified, it will be ignored.} +} +\value{ +A lattice plot of class \dQuote{trellis}. +} +\description{ +A wrapper for the lattice (trellis) plot methods for grid and station data. +} +\details{ +The function applies the \code{\link[lattice]{xyplot}} method after computing spatial aggregation +(parameter \code{aggr.spatial}) and member aggregation (mean and range) to the imput grids (or station data). + +In case of multimember grids, the function will internally compute the ensemble mean +and the range for plotting. The range is used to plot the shadow of the multimember spread. + + \strong{Controlling graphical parameters} + + Many different aspects of the plot can be controlled passing the relevant arguments to + \code{xyplot}. Fine control of graphical parameters for the trellis display can + be also controlled using \code{\link[lattice]{trellis.par.set}}. + + \strong{Future Work} + + Implement parameters missing.dates (default will be TRUE) to suppress + dates without values from xaxis when FALSE. + Implement parameter show.na (default will be FALSE) to fill with gray + NA values when TRUE. + Implement auxiliary functions to reduce dependencies and remove from + imports packages \code{data.table}, \code{padr} and \code{grDevices}. +} +\examples{ +data("CFS_Iberia_pr") +data("EOBS_Iberia_pr") +data("VALUE_Iberia_pr") +# Combine grids with members (CFS) and without members (EOBS) +a <- subsetGrid(CFS_Iberia_tas, years = 1985:1992) +b <- subsetGrid(EOBS_Iberia_tas, years = 1985:1992) +temporalPlot("EOBS" = b, "CFS" = a, + xyplot.custom = list(main = "winter temperature", ylab = "Celsius")) +# Station and grid data can be combined, also different temporal ranges +v <- subsetGrid(VALUE_Iberia_tas, years = 1988:1990) +temporalPlot("VALUE" = v, "EOBS" = b, "CFS" = a, lwd = 0.9, + aggr.spatial = list(FUN = min, na.rm = TRUE), + xyplot.custom = list(main = "winter temperature", + ylab = "Celsius", ylim = c(-20, 10))) +# Use subsetGrid to obtain and plot a single location (no spatial aggregation) +a1 <- subsetGrid(a, lonLim = 2, latLim = 42) +b1 <- subsetGrid(b, lonLim = 2, latLim = 42) + +temporalPlot("EOBS" = b1, "CFS" = a1, + cols = c("green", "deeppink"), show.na = TRUE, + xyplot.custom = list(main = "winter temperature", ylab = "Celsius")) +} +\author{ +M. Iturbide +} diff --git a/man/tercileMap.Rd b/man/tercileMap.Rd new file mode 100644 index 0000000..5221170 --- /dev/null +++ b/man/tercileMap.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/tercileMap.R +\name{tercileMap} +\alias{tercileMap} +\title{Tercile map for visualization of forecast skill of seasonal climate predictions.} +\usage{ +tercileMap(hindcast, forecast, ...) +} +\arguments{ +\item{hindcast}{A multi-member list with the hindcast for verification. See details.} + +\item{forecast}{A multi-member list with the forecasts. Default is NULL.} + +\item{...}{Optional arguments passed to \code{\link{spatialPlot}}} +} +\value{ +A lattice plot of class \dQuote{trellis}. +} +\description{ +Tercile map for the visualization of forecast skill of seasonal climate predictions. + This function is prepared to plot the data sets loaded from the ECOMS User Data Gateway (ECOMS-UDG). See + the loadeR.ECOMS R package for more details (http://meteo.unican.es/trac/wiki/udg/ecoms/RPackage). +} +\details{ +The function applies the \code{\link{spatialPlot}} function, that in turn uses \code{lattice-methods}. +} +\note{ +The function is based on the tercile maps delivered by +C3S seasonal forecasts (\url{http://climate.copernicus.eu/s/charts/c3s_seasonal/}). + + Many different aspects of the plot can be controlled passing the relevant arguments to + \code{\link[sp]{spplot}}. +} +\examples{ +library(transformeR) +hindcast <- subsetGrid(CFS_Iberia_tas, years = 1983:2001) +forecast <- subsetGrid(CFS_Iberia_tas, years = 2002) +tercileMap(hindcast, forecast) +} +\author{ +M. Iturbide +}