|
| 1 | +#' Add rerddap::griddap() data to a plotdap map |
| 2 | +#' |
| 3 | +#' \code{add_griddap} adds the data from an 'rerddap::griddap() call to |
| 4 | +#' a 'plotdap' map |
| 5 | +#' @param plot a \link{plotdap} object. |
| 6 | +#' @param grid a \link{griddap} object. |
| 7 | +#' @param var a formula defining a variable, or function of variables to visualize. |
| 8 | +#' @param fill either a character string of length 1 matching a name in \\link[rerddap]{colors} |
| 9 | +#' or a vector of color codes. This defines the colorscale used to encode values |
| 10 | +#' of \code{var}. |
| 11 | +#' @param maxpixels integer > 0. Maximum number of cells to use for the plot. |
| 12 | +#' If maxpixels < ncell(x), sampleRegular is used before plotting. |
| 13 | +#' If gridded=TRUE maxpixels may be ignored to get a larger sample |
| 14 | +#' @param time how to resolve multiple time frames. Choose one of the following: |
| 15 | +#' \itemize{ |
| 16 | +#' \item A function to apply to each observation at a particular location |
| 17 | +#' (\link{mean} is the default). |
| 18 | +#' \item A character string (of length 1) matching a time value. |
| 19 | +#' } |
| 20 | +#' @param animate whether to animate over the \code{time} variable (if it exists). |
| 21 | +#' Currently only implemented for \code{method='ggplot2'} and requires the |
| 22 | +#' gganimate package. |
| 23 | +#' @param cumulative - if animation should be cumulative -default FALSE |
| 24 | +#' @param ... arguments passed along to \code{geom_sf()} |
| 25 | +#' (if \code{method='ggplot2'}, otherwise ignored). |
| 26 | +#' @return A plotdap object |
| 27 | +#' @export |
| 28 | +#' @rdname add_griddap |
| 29 | +#' @examples |
| 30 | +#' |
| 31 | +#' # base plotting tends to be faster, |
| 32 | +#' # but is less extensible plotdap("base") |
| 33 | +#' |
| 34 | +#' # actual datasets in data folder to meet execution timings |
| 35 | +#' |
| 36 | +#'\donttest{ |
| 37 | +#' murSST <- griddap( |
| 38 | +#' 'jplMURSST41', latitude = c(35, 40), longitude = c(-125, -1205), |
| 39 | +#' time = c('last', 'last'), fields = 'analysed_sst' |
| 40 | +#' ) |
| 41 | +#' |
| 42 | +#' QMwind <- griddap( |
| 43 | +#' 'erdQMwindmday', time = c('2016-11-16', '2017-01-16'), |
| 44 | +#' latitude = c(30, 50), longitude = c(210, 240), |
| 45 | +#' fields = 'x_wind' |
| 46 | +#' ) |
| 47 | +#' |
| 48 | +#' p <- plotdap(crs = "+proj=robin") |
| 49 | +#' add_griddap(p, murSST, ~analysed_sst) |
| 50 | +#' |
| 51 | +#' p <- plotdap(mapTitle = "Average wind over time") |
| 52 | +#' add_griddap(p, QMwind, ~x_wind) |
| 53 | +#' |
| 54 | +#'} |
| 55 | +#' |
| 56 | +# |
| 57 | +#' p <- plotdap("base", crs = "+proj=robin") |
| 58 | +#' p <- add_griddap(p, murSST, ~analysed_sst) |
| 59 | +#' |
| 60 | +#' # layer tables on top of grids |
| 61 | +#' require(magrittr) |
| 62 | +#' p <- plotdap("base") %>% |
| 63 | +#' add_griddap(murSST, ~analysed_sst) %>% |
| 64 | +#' add_tabledap(sardines, ~subsample_count) |
| 65 | +#' |
| 66 | +#' # multiple time periods |
| 67 | +#' p <- plotdap("base", mapTitle = "Average wind over time") |
| 68 | +#' p <- add_griddap(p, QMwind, ~x_wind) |
| 69 | +#' |
| 70 | + |
| 71 | +add_griddap <- function(plot, grid, var, fill = "viridis", |
| 72 | + maxpixels = 10000, time = mean, animate = FALSE, |
| 73 | + cumulative = FALSE, ...) { |
| 74 | + if (!is.grid(grid)) |
| 75 | + stop("The `grid` argument must be a `griddap()` object", call. = FALSE) |
| 76 | + if (!lazyeval::is_formula(var)) |
| 77 | + stop("The `var` argument must be a formula", call. = FALSE) |
| 78 | + if (!is.function(time) && !is.character(time)) |
| 79 | + stop("The `time` argument must be a function or a character string", |
| 80 | + call. = FALSE) |
| 81 | + |
| 82 | + # create raster object from filename; |
| 83 | + # otherwise create a sensible raster from data |
| 84 | + r <- get_raster(grid, var) |
| 85 | + |
| 86 | + # checks for naming and numeric lat/lon |
| 87 | + latlon_is_valid(r) |
| 88 | + # adjust to ensure everthing is on standard lat/lon scale |
| 89 | + r <- latlon_adjust(r) |
| 90 | + |
| 91 | + # if necessary, reduce a RasterBrick to a RasterLayer |
| 92 | + # http://gis.stackexchange.com/questions/82390/summarize-values-from-a-raster-brick-by-latitude-bands-in-r |
| 93 | + if (body(time) == 'UseMethod("mean")') { |
| 94 | + time <- function(x) mean(x, na.rm = TRUE) |
| 95 | + } |
| 96 | + if (raster::nlayers(r) > 1) { |
| 97 | + if (is.function(time)) { |
| 98 | + r <- raster::calc(r, time) |
| 99 | + } else { |
| 100 | + nm <- make.names(time) |
| 101 | + if (!nm %in% names(r)) { |
| 102 | + warning( |
| 103 | + "The `time` argument doesn't match any of time values.\n", |
| 104 | + sprintf( |
| 105 | + "Valid options include: '%s'", |
| 106 | + paste(unique(grid$data$time), collapse = "', '") |
| 107 | + ), |
| 108 | + call. = FALSE |
| 109 | + ) |
| 110 | + } |
| 111 | + r <- r[[nm]] |
| 112 | + } |
| 113 | + |
| 114 | + if (raster::nlayers(r) > 1 && !animate) { |
| 115 | + stop( |
| 116 | + "The `time` argument hasn't reduced the raster down to a single layer.\n", |
| 117 | + "Either set `animate=TRUE` or provide a suitable value to `time`.", |
| 118 | + call. = FALSE |
| 119 | + ) |
| 120 | + } |
| 121 | + } |
| 122 | + |
| 123 | + # simplify raster, if necessary |
| 124 | + n <- raster::ncell(r) |
| 125 | + if (n > maxpixels) { |
| 126 | + message("grid object contains more than ", maxpixels, " pixels") |
| 127 | + message("increase `maxpixels` for a finer resolution") |
| 128 | + rnew <- raster::raster( |
| 129 | + nrow = floor(raster::nrow(r) * sqrt(maxpixels / n)), |
| 130 | + ncol = floor(raster::ncol(r) * sqrt(maxpixels / n)), |
| 131 | + crs = raster::crs(r), |
| 132 | + ext = raster::extent(r) |
| 133 | + ) |
| 134 | + if (inherits(r, "RasterBrick")) { |
| 135 | + for (i in seq_len(raster::nlayers(r))) { |
| 136 | + r[[i]] <- raster::resample(r[[i]], rnew, method = 'bilinear') |
| 137 | + } |
| 138 | + } else { |
| 139 | + r <- raster::resample(r, rnew, method = 'bilinear') |
| 140 | + } |
| 141 | + } |
| 142 | + |
| 143 | + # assumes we apply sf::st_crs() to plot on initiation |
| 144 | + if (inherits(plot$crs, "crs")) { |
| 145 | + crs_string <- plot$crs$proj4string |
| 146 | + if (!is.na(plot$crs$epsg)) { |
| 147 | + epsg_string <- paste0("+init=epsg:", plot$crs$epsg) |
| 148 | + crs_string <- paste(epsg_string, crs_string) |
| 149 | + } |
| 150 | + #r <- raster::projectRaster(r, crs = plot$crs$proj4string) |
| 151 | + r <- raster::projectRaster(r, crs = crs_string) |
| 152 | + } |
| 153 | + |
| 154 | + # color scale |
| 155 | + cols <- if (length(fill) == 1) rerddap::colors[[fill]] else fill |
| 156 | + |
| 157 | + if (is_ggplotdap(plot)) { |
| 158 | + # TODO: not the most efficient approach, but it will have to do for now |
| 159 | + # https://twitter.com/hadleywickham/status/841763265344487424 |
| 160 | + s <- sf::st_as_sf(raster::rasterToPolygons(r)) |
| 161 | + vars <- setdiff(names(s), "geometry") |
| 162 | + sg <- sf::st_as_sf(tidyr::gather_(s, "variable", "value", vars)) |
| 163 | + if (animate) { |
| 164 | + try_gganimate() |
| 165 | + plot$animate <- TRUE |
| 166 | + plot$nper <- length(sg) |
| 167 | + plot$ggplot <- plot$ggplot + |
| 168 | + gganimate::transition_manual(variable, cumulative = cumulative) + |
| 169 | + ggplot2::labs(title = "{current_frame}") |
| 170 | + } |
| 171 | + |
| 172 | + return( |
| 173 | + add_ggplot( |
| 174 | + plot, |
| 175 | + geom_sf(data = sg, |
| 176 | + mapping = aes_string(fill = "value", colour = "value"), ...), |
| 177 | + scale_fill_gradientn(name = lazyeval::f_text(var), colors = cols), |
| 178 | + scale_colour_gradientn(colors = cols), |
| 179 | + guides(colour = FALSE) |
| 180 | + ) |
| 181 | + ) |
| 182 | + } |
| 183 | + |
| 184 | + if (animate) { |
| 185 | + warning( |
| 186 | + "Animations are currently only implemented for `method='ggplot2'`", |
| 187 | + call. = FALSE |
| 188 | + ) |
| 189 | + } |
| 190 | + |
| 191 | + # TODO: more props! |
| 192 | + grid <- structure( |
| 193 | + r, props = list( |
| 194 | + name = lazyeval::f_text(var), |
| 195 | + values = raster::values(r), |
| 196 | + color = cols |
| 197 | + ) |
| 198 | + ) |
| 199 | + |
| 200 | + # Throw a warning if the grid extent overlaps with another grid? |
| 201 | + plot$layers <- c( |
| 202 | + plot$layers, list(grid) |
| 203 | + ) |
| 204 | + |
| 205 | + plot |
| 206 | +} |
0 commit comments