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 zoomin Numeric. Maps to [ggspatial::annotation_map_tile()]
#' @param cachedir Character. Maps to [ggspatial::annotation_map_tile()], defaults to
#' `tools::R_user_dir("bioRad")`
#' @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{
#' if (all(sapply(c("ggspatial","prettymapr", "rosm"), requireNamespace, quietly = TRUE))) {
#' # 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, zoomin = -2, cachedir = tools::R_user_dir("bioRad"), ...) {

  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."
    )
  }

  # validate xlim and ylim
  if(!missing(xlim)) assertthat::assert_that(is.numeric(xlim), length(xlim) == 2, msg = "xlim must be a numeric vector of length 2")
  if(!missing(ylim)) assertthat::assert_that(is.numeric(ylim), length(ylim) == 2, msg = "ylim must be a numeric vector of length 2")

  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", "rosm"),'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
  wgs84 <- sp::CRS(SRS_string = sf::st_crs(4326)$wkt)
  epsg3857 <- sp::CRS(SRS_string = sf::st_crs(3857)$wkt)

  mybbox <-
    sp::spTransform(
      sp::SpatialPoints(t(data@bbox),
                        proj4string = data@proj4string
      ),
      epsg3857
    )

  mybbox.wgs <-
    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 <- as.data.frame(sp::spTransform(methods::as(data,"SpatialPointsDataFrame"), epsg3857))

  # 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
                               ))

  # Initialize projected limits as NULL
  projected_xlim <- NULL
  projected_ylim <- NULL

  # Check if xlim is provided; if so, generate projected xlim
  if (!missing(xlim)) {
    # Create a temporary bounding box for x limits only (use full latitude range for the bbox)
    bbox_latlon_x <- sf::st_bbox(c(xmin = xlim[1], xmax = xlim[2], ymin = -90, ymax = 90), crs = 4326)
    bbox_projected_x <- sf::st_transform(sf::st_as_sfc(bbox_latlon_x), crs = 3857)
    projected_xlim <- sf::st_bbox(bbox_projected_x)[c("xmin", "xmax")]
  }

  # Check if ylim is provided; if so, generate projected ylim
  if (!missing(ylim)) {
    # Create a temporary bounding box for y limits only (use full longitude range for the bbox)
    bbox_latlon_y <- sf::st_bbox(c(xmin = -180, xmax = 180, ymin = ylim[1], ymax = ylim[2]), crs = 4326)
    bbox_projected_y <- sf::st_transform(sf::st_as_sfc(bbox_latlon_y), crs = 3857)
    projected_ylim <- sf::st_bbox(bbox_projected_y)[c("ymin", "ymax")]
  }

  mymap <-
    ggplot2::ggplot() +
    ggspatial::annotation_map_tile(type = map, zoomin = zoomin, cachedir = cachedir) +
    ggplot2::coord_sf(xlim = projected_xlim, ylim = projected_ylim, expand=FALSE) +
    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")

  mymap
}