R/map.R

In bioRad: Biological Analysis and Visualization of Weather Radar Data

#' Map a plan position indicator (`ppi`) on a map
#'
#' Plots a plan position indicator (`ppi`) on a base layer
#'
#' @param x A `ppi` object.
#' @param map Basemap to use, one of `rosm::osm.types()`
#' @param param Character. Scan parameter to plot, e.g. `DBZH` or `VRADH`. See
#'   [summary.param()] for commonly available parameters.
#' @param alpha Numeric. Transparency of the data, value between 0 and 1.
#' @param xlim Numeric vector of length 2. Range of x values (degrees longitude)
#'   to plot.
#' @param ylim Numeric vector of length 2. Range of y values (degrees latitude)
#'   to plot.
#' @param zlim Numeric vector of length 2. The range of values to plot.
#' @param ratio Numeric. Aspect ratio between x and y scale, by default
#' \eqn{1/cos(latitude radar * pi/180)}.
#' @param radar_size Numeric. Size of the symbol indicating the radar position.
#' @param radar_color Character. Color of the symbol indicating the radar
#'   position.
#' @param n_color Numeric. Number of colors (>=1) to use in the palette.
#' @param palette Character vector. Hexadecimal color values defining the plot
#'   color scale, e.g. output from [viridisLite::viridis()].
#' @param ... Arguments passed to [ggplot2::ggplot()].
#' @importFrom methods as
#' @return A ggplot object
#' @details
#' Available scan parameters for mapping can by printed to screen by
#' `summary(x)`. Commonly available parameters are:
#' * `DBZH`, `DBZ`: (Logged) reflectivity factor (dBZ)
#' * `TH`, `T`: (Logged) uncorrected reflectivity factor (dBZ)
#' * `VRADH`, `VRAD`: Radial velocity (m/s). Radial velocities towards the radar
#'   are negative, while radial velocities away from the radar are positive
#' * `RHOHV`: Correlation coefficient (unitless) Correlation between vertically
#'   polarized and horizontally polarized reflectivity factor
#' * `PHIDP`: Differential phase (degrees)
#' * `ZDR`: (Logged) differential reflectivity (dB)
#' The scan parameters are named according to the OPERA data information
#' model (ODIM), see Table 16 in the
#' [ODIM specification](https://github.com/adokter/vol2bird/blob/master/doc/OPERA2014_O4_ODIM_H5-v2.2.pdf).
#'
#' @export
#'
#' @seealso
#' * [project_as_ppi()]
#'
#' @examples
#' # Project a scan as a ppi
#' ppi <- project_as_ppi(example_scan)
#' \donttest{
#' # Choose a basemap
#' basemap <- rosm::osm.types()[1]
#'
#' # Map the radial velocity of the ppi onto the basemap
#' map(ppi, map = basemap, param = "VRADH")
#'
#' # Extend the plotting range of velocities, from -50 to 50 m/s
#' map(ppi, map = basemap, param = "VRADH", zlim = c(-50, 50))
#'
#' # Map the reflectivity
#' map(ppi, map = basemap, param = "DBZH")
#'
#' # Change the color palette to Viridis colors
#' map(ppi, map = basemap, param = "DBZH", palette = viridis::viridis(100), zlim=c(-10,10))
#'
#' # Give the data more transparency
#' map(ppi, map = basemap, param = "DBZH", alpha = 0.3)
#'
#' # Change the appearance of the symbol indicating the radar location
#' map(ppi, map = basemap, radar_size = 5, radar_color = "blue")
#'
#' # Crop the map
#' map(ppi, map = basemap, xlim = c(12.4, 13.2), ylim = c(56, 56.5))
#' }
map <- function(x, ...) {
  UseMethod("map", x)
}

#' @describeIn map Plot a `ppi` object on a map.
#' @export
map.ppi <- function(x, map="cartolight", param, alpha = 0.7, xlim, ylim, zlim = c(-20, 20),
                    ratio, radar_size = 3, radar_color = "#202020", n_color = 1000,
                    palette = NA, ...) {

  stopifnot(inherits(x, "ppi"))

  if (methods::hasArg("quantity")) stop("unknown function argument 'quantity`. Did you mean `param`?")

  if (missing(param)) {
    if ("DBZH" %in% names(x$data)) {
      param <- "DBZH"
    } else {
      param <- names(x$data)[1]
    }
  } else if (!is.character(param)) {
    stop(
      "'param' should be a character string with a valid ",
      "scan parameter name."
    )
  }
  if (missing(zlim)) {
    zlim <- get_zlim(param, zlim)
  }
  if (!(param %in% names(x$data))) {
    stop(paste("no scan parameter '", param, "' in this ppi", sep = ""))
  }
  if (is.na(raster::crs(x$data))) {
    stop("Not a projected ppi, map() expects a ppi generated by project_as_ppi() with argument project=TRUE")
  }

  # check that suggested dependencies are presetn
  rlang::check_installed(c("ggspatial","prettymapr"),'to map ppi\'s')

  assert_valid_basemap <- function(basemap) {
    valid_types <- rosm::osm.types()
    if (!(basemap %in% valid_types)) {
      stop(paste("Invalid basemap argument. Must be one of:", paste(valid_types, collapse=", ")))
    }
  }

  assert_valid_basemap(map)

  # set color scales and palettes
  if (!assertthat::are_equal(palette, NA)) {
    if(!(is.character(palette) && length(palette) > 1)) stop("palette should be a character vector with hex color values")
    # apply transparancy
    palette <- ggplot2::alpha(palette,alpha)
    n_color = length(palette)
    colorscale <- color_palette_to_scale_colour(param, zlim, palette, na.value = "transparent")
  }
  else{
    palette <- color_palette(param = param, n_color = n_color, alpha = alpha)
    colorscale <- color_scale(param, zlim)
  }

  # extract the scan parameter
  data <- do.call(function(y) x$data[y], list(param))
  wgs84 <- sp::CRS("+proj=longlat +datum=WGS84")
  epsg3857 <- sp::CRS("+init=epsg:3857") # this is the google mercator projection
  mybbox <- suppressWarnings(
    sp::spTransform(
      sp::SpatialPoints(t(data@bbox),
        proj4string = data@proj4string
      ),
      sp::CRS("+init=epsg:3857")
    )
  )
  mybbox.wgs <- suppressWarnings(
    sp::spTransform(
      sp::SpatialPoints(t(data@bbox),
        proj4string = data@proj4string
      ),
      wgs84
    )
  )
  e <- raster::extent(mybbox.wgs)
  r <- raster::raster(raster::extent(mybbox),
    ncol = data@grid@cells.dim[1] * .9,
    nrow = data@grid@cells.dim[2] * .9, crs = sp::CRS(sp::proj4string(mybbox))
  )

  # convert to google earth mercator projection
  data <- suppressWarnings(
    as.data.frame(sp::spTransform(as(data,"SpatialPointsDataFrame"), sp::CRS("+init=epsg:3857")))
  )
  # bring z-values within plotting range
  index <- which(data$z < zlim[1])
  if (length(index) > 0) {
    data[index, ]$z <- zlim[1]
  }
  index <- which(data$z > zlim[2])
  if (length(index) > 0) {
    data[index, ]$z <- zlim[2]
  }

  # rasterize data on the raster, if there is valid data
  if(FALSE %in% is.na(data[,1])){
    r <- raster::rasterize(data[, 2:3], r, data[, 1])
  } else{
    raster::values(r) <- NA
  }

  # function to convert values to hex color strings
  col_func <- function(value, lim) {
    output <- rep(0, length(value))
    output <- round((value - lim[1]) / (lim[2] - lim[1]) * n_color)
    output[output > n_color] <- n_color
    output[output < 1] <- 1
    return(palette[output])
  }

  # convert data values to hex color string values.

  rdf <- raster::as.data.frame(r, xy = TRUE)
  names(rdf) <- c("x", "y", "fill")
  rdf$fill <- col_func(r@data@values, zlim)

  # these declarations prevent generation of NOTE "no visible binding for
  # global variable" during package Check
  fill <- lon <- lat <- X <- Y <- y <- z <- NA
  # extract unique radar locations

  # dummy is a hack to be able to include the ggplot2 color scale,
  # radarpoint is the actual plotting of radar positions.

  radar_location = data.frame(lon = x$geo$lon, lat = x$geo$lat) %>%
    sf::st_as_sf(coords = c("lon", "lat"))  %>%
    sf::st_set_crs(4326)  %>%
    sf::st_transform(3857)

  radar_df <- data.frame(sf::st_coordinates(radar_location))

  dummy <- ggplot2::geom_point(ggplot2::aes(x = lon, y = lat, colour = z),
                               size = 0,
                               data = data.frame(
                                 lon = radar_df$X,
                                 lat = radar_df$Y,
                                 z = 0
                               ))

  mymap <- suppressWarnings(
    ggplot2::ggplot() +
      ggspatial::annotation_map_tile(type = map) +
      ggplot2::geom_raster(data = rdf, ggplot2::aes(x = x, y = y, fill = fill), na.rm = TRUE, interpolate = FALSE) +
      ggplot2::scale_fill_identity(na.value = "transparent") +
      ggplot2::geom_point(data = radar_df, ggplot2::aes(x = X, y = Y),
                 shape = 21,
                 fill = "transparent",
                 colour = radar_color,
                 size = radar_size,
                 stroke = 1.5,
                 show.legend = FALSE) +
      dummy +
      colorscale +
      ggplot2::labs(x="lon", y="lat")
  )

  suppressWarnings(mymap)
}