R/geos-plot.R

In geos: Open Source Geometry Engine ('GEOS') R API

#' Plot GEOS geometries
#'
#' @param x A [GEOS geometry vector][as_geos_geometry]
#' @inheritParams wk::wk_plot
#' @param simplify A relative tolerance to use for simplification of
#'   geometries. Use 0 to disable simplification; use a higher number
#'   to make simplification coarser.
#' @param crop Use `TRUE` to crop the input to the extent of the plot.
#'
#' @return The input, invisibly
#' @export
#'
#' @examples
#' plot(as_geos_geometry("LINESTRING (0 0, 1 1)"))
#' plot(as_geos_geometry("POINT (0.5 0.4)"), add = TRUE)
#'
plot.geos_geometry <- function(x, ..., asp = 1, bbox = NULL, xlab = "", ylab = "",
                               rule = "evenodd", add = FALSE,
                               simplify = 1, crop = TRUE) {
  # this is too hard without vctrs (already in Suggests)
  if (!requireNamespace("vctrs", quietly = TRUE)) {
    stop("Package 'vctrs' is required for plot.geos_geometry()", call. = FALSE) # nocov
  }

  if (!add) {
    plot_bbox <- if (is.null(bbox)) wk::wk_bbox(x) else bbox
    wk::wk_plot(wk::wkb(), asp = asp, bbox = plot_bbox, xlab = xlab)
  }

  if (length(x) == 0) {
    return(x)
  }

  # estimate resolution for simplifying
  usr <- graphics::par("usr")
  usr_x <- usr[1:2]
  usr_y <- usr[3:4]
  device_x <- graphics::grconvertX(usr_x, to = "device")
  device_y <- graphics::grconvertY(usr_y, to = "device")

  # Use resolution of 0.05 at the device level, scale to usr coords.
  # This rarely results in simplification that will be noticed by the user.
  scale_x <- diff(device_x) / diff(usr_x)
  scale_y <- diff(device_y) / diff(usr_y)
  scale <- min(abs(scale_x), abs(scale_y))
  resolution_usr <- 0.05 / scale

  x_plot <- x

  # we do some subsetting of x_plot so we need to keep the dots aligned
  dots <- list(...)
  dots_is_vector <- vapply(dots, vctrs::vec_is, logical(1)) &
    vapply(dots, function(x) !identical(length(x), 1L), logical(1))
  dots_scalar <- dots[!dots_is_vector]
  dots_vector <- dots[dots_is_vector]

  # if adding, we only need features that touch the plot bbox
  plot_area <- wk::rct(usr_x[1], usr_y[1], usr_x[2], usr_y[2], crs = wk::wk_crs(x_plot))
  ignore_bbox <- is.null(bbox) && !add

  if (!ignore_bbox) {
    x_tree <- geos_strtree(x_plot)
    x_touching <- geos_strtree_query(x_tree, plot_area)[[1]]
    x_plot <- x_plot[x_touching]
    dots_vector <- lapply(dots_vector, vctrs::vec_slice, x_touching)
  }

  if (simplify > 0) {
    x_plot <- geos_simplify(x_plot, resolution_usr * simplify)
  }

  if (!ignore_bbox && crop) {
    # give the crop bbox 5% so that new borders get generated outside
    x_mid <- mean(usr_x)
    y_mid <- mean(usr_y)
    usr_x_expanded <- x_mid + ((usr_x - x_mid) * 1.05)
    usr_y_expanded <- y_mid + ((usr_y - y_mid) * 1.05)

    crop_area <- wk::rct(
      usr_x_expanded[1], usr_y_expanded[1],
      usr_x_expanded[2], usr_y_expanded[2],
      crs = wk::wk_crs(x_plot)
    )

    x_plot <- geos_clip_by_rect(x_plot, crop_area)
  }

  # skip empties generated by simplification or clip by rect
  x_is_empty <- geos_is_empty(x_plot)
  x_plot <- x_plot[!x_is_empty]
  dots_vector <- lapply(dots_vector, vctrs::vec_slice, !x_is_empty)

  # pass plotting on to wk_plot()
  if (length(x_plot) > 0) {
    do.call(wk::wk_plot, c(list(x_plot), dots_scalar, dots_vector, list(add = TRUE)))
  }

  invisible(x)
}