multilandR: multilandr: Landscape Analysis at Multiple Spatial Scales

Documented in mland_plot

# Plotting classes function
.plot_classes <- function(plot, sub_int, arg_name, st_classes){

  sub_int_df <- terra::as.data.frame(sub_int, xy = T)
  colnames(sub_int_df)[3] <- "value"
  if(arg_name == "raster"){
    fills <- data.frame(class = st_classes$fill[seq(1, length(st_classes$fill), 2)],
                        fill = st_classes$fill[seq(2, length(st_classes$fill), 2)])
    fills <- fills[order(fills$class), ]
    alphas <- data.frame(class = st_classes$alpha[seq(1, length(st_classes$alpha), 2)],
                         alpha = st_classes$alpha[seq(2, length(st_classes$alpha), 2)])
    alphas <- alphas[order(alphas$class), ]
    fills_cols  <- fills[fills$class %in% levels(factor(sub_int_df$value)), "fill"]
    fills_alpha <- alphas[alphas$class %in% levels(factor(sub_int_df$value)), "alpha"]
    plot <- plot +
      ggplot2::geom_raster(data = sub_int_df,
                           ggplot2::aes(x = x, y = y, fill = factor(value), alpha = factor(value))) +
      ggplot2::scale_fill_manual(values = fills_cols, na.value = st_classes$na_value[1]) +
      ggplot2::scale_alpha_manual(values = fills_alpha, na.value = as.numeric(st_classes$na_value[2]))
  } else {
    plot <- plot +
      ggplot2::geom_raster(data = sub_int_df,
                           ggplot2::aes(x = x, y = y, fill = value), na.rm = T)
  }

  return(plot)
}

# Plotting buffers, points and labels
.plot_bpl <- function(x, plot, p, pos_buffers, points, tit, radii, st_buffers, st_points){

  # Plot buffers
  for(r in radii){
    buff <- x@buffers[pos_buffers, ]
    plot <- plot + suppressMessages(
      tidyterra::geom_spatvector(data      = buff,
                                 fill      = NA,
                                 colour    = ggplot2::alpha(st_buffers$col, st_buffers$alpha),
                                 linewidth = st_buffers$lwd,
                                 linetype  = st_buffers$lty))
  }

  # Plot points
  plot <- plot +
    suppressMessages(
      tidyterra::geom_spatvector(data = x@points[points[p], ],
                                 size = st_points$size,
                                 shape = st_points$shape,
                                 fill = st_points$fill,
                                 colour =  st_points$col,
                                 alpha =  st_points$alpha))

  # Plot points labels
  plot <- plot +
    ggplot2::ggtitle(tit) +
    ggplot2::theme(plot.title   = ggplot2::element_text(size = ggplot2::rel(0.9)),
                   panel.border = ggplot2::element_blank())

  # Fixes visualization
  plot <- plot +
    #ggplot2::coord_equal() +
    ggplot2::theme(panel.background = ggplot2::element_rect(fill = "white"),
                   legend.position = "none",
                   plot.title = ggplot2::element_text(hjust = 0.5))

  # Remove annoying axes
  plot <- plot +
    ggplot2::labs(x = "", y ="") +
    ggplot2::theme(axis.ticks.y = ggplot2::element_blank(), axis.text.y = ggplot2::element_blank(),
                   axis.ticks.x = ggplot2::element_blank(), axis.text.x = ggplot2::element_blank(),
                   axis.ticks.length = ggplot2::unit(0, "mm"))

  return(plot)
}

get_ext_range <- function(x, ext_range){
  if(is.null(ext_range)){
    ext_range <- c(terra::minmax(x)[1, 1], terra::minmax(x)[2, 1])
  } else {
    if(!all(is.na(terra::minmax(x)[, 1]))){
      if(terra::minmax(x)[1, 1] < ext_range[1]) ext_range[1] <- terra::minmax(x)[1, 1]
      if(terra::minmax(x)[2, 1] > ext_range[2]) ext_range[2] <- terra::minmax(x)[2, 1]
    }
  }
  ext_range
}

#' Plots landscapes from 'MultiLand' objects
#'
#' Returns multiple plots for each landscape generated from each point and buffer, with their radii and classes,
#' defined by the user through a 'MultiLand' object (generated by [mland()]).
#' Aesthetic parameters of plots can be customized.
#'
#' @param points Numeric or character vector of points to be plotted. See Details.
#' @param radii Numeric vector of radii to be plotted.
#' @param st_points List of aesthetic arguments for points plotting:
#'   \code{shape} for points shape, \code{size} for points size, \code{col} for
#'   points border color, \code{fill} for points fill color and \code{alpha} for point transparency.
#' @param st_buffers List of aesthetic arguments for buffers plotting:
#'   \code{lty} for buffers linetype, \code{lwd} for buffers linewidth,
#'   \code{col} for buffers border color and \code{alpha} for border transparency.
#' @param st_classes List of aesthetic arguments for classes plotting:
#'   \code{palette}, for classes color palette, \code{fill} a vector of fill colors for classes, \code{alpha}, a vector of
#'   alpha values for classes, and \code{na_value} for the color of NA values. See Details.
#' @param st_ext Character vector of length 2, depicting the color for the minimum and maximum values
#' of the raster defined in `ext_raster`.
#' @param x An object of class 'MultiLand' generated with [mland()].
#' @param raster,ext_raster Numeric. The rasterlayer to be plotted. Only one rasterlayer can be
#' plotted at the same time, either defined in `raster` or `ext_raster`.
#' @param title One of the following: "id" to plot titles as each point id (default), or "sitename" to
#' plot titles as each pre-defined point name in `x`. See Details.
#' @param ncol,nrow Number of columns and rows wherein individual plots will be arranged.
#'
#' @details
#' If argument `points` is a character vector,
#' [mland_plot()] will assume that the 'MultiLand' object inputted in argument `x` was created with
#' `site_ref = TRUE`. This is, there is a column/attribute in points layer data with the names for
#' each distinct point. Therefore, the inputted values in argument `points` will be taken as these
#' identification names. Otherwise, if a numeric vector is inputted, these values
#' will be taken as the automatically generated point ids (created when running [mland()]).
#'
#' If `title = "sitename"`, the title of individual plots will be the names of each point. For this,
#' the names of the points in `x` must had been defined when the object was created with [mland()]
#' (i.e. `x@site_ref = TRUE`). Otherwise, the argument will be ignored and the titles will be the
#' ids of the points.
#'
#' A pre-defined palette can be chosen to differentiate classes inside `palette = "palette_name"`,
#' inside the list defined in `st_classes`. Any palette from [hcl.pals()] can be chosen. Otherwise,
#' the user can define specific colors for each class, inside `fill`. This must be a vector built
#' with concatenated pair of values, the first value being the class (or class name, if defined
#' during `x` generation), and the second value the color (either the name of the color or the hex
#' code of the color). For example, in the case the rasterlayer has four unique values: (1, 2, 3 and 4), a plausible color definition
#' could be the following:
#'
#' \preformatted{
#'  list(c(1, "green", 2, "red", 3, "black", 4, "yellow"))
#' }
#'
#' @return Multiple plots (in a unique plotting device) of landscapes around
#'   defined points, radii and classes of a MultiLand object.
#'
#' @export
#'
#' @examples
#' \dontrun{
#' # Loads a 'MultiLand' object
#' ernesdesign <- system.file("extdata", "ernesdesign.zip", package = "multilandr")
#' ernesdesign <- load_mland(ernesdesign)
#'
#' # Plots all points and radii
#' mland_plot(ernesdesign)
#'
#' # Plots points 1 to 3 and only radius 3000 m
#' mland_plot(ernesdesign, points = 1:3, radii = 3000)
#'
#' # Plot with pre-defined colors, and specifying other arguments
#' cols <- c(1, "forestgreen",
#'           2, "darkolivegreen2",
#'           3, "firebrick3",
#'           4, "goldenrod1",
#'           5, "deepskyblue3",
#'           6, "black")
#'
#' mland_plot(ernesdesign, points = 9:11, radii = c(1000, 2000, 3000),
#'            title = "sitename", nrow = 1,
#'            st_points = list(shape = 9),
#'            st_buffers = list(lty = "dashed"),
#'            st_classes = list(fill = cols))
#'
#' # Plot a unique landscape by calling it with its name
#' mland_plot(ernesdesign, points = "Peje", title = "sitename",
#'            st_points = list(shape = 15, col = "red"),
#'            st_classes = list(palette = "Hawaii"))
#'
#' # Plot extra rasterlaer
#' mland_plot(ernesdesign, radii = 3000, ext_raster = 1, title = "sitename")
#'
#' # Plot extra rasterlater with customized colors
#' mland_plot(ernesdesign, radii = 3000, ext_raster = 1, title = "sitename",
#'            st_ext = c("blue", "red"))
#' }
mland_plot <- function(x, raster = NULL, points = NULL, radii = NULL, ext_raster = NULL,
                       title = "id", ncol = NULL, nrow = NULL,
                       st_points  = list(shape = 21, size = 2, col = "black", fill = "white",
                                         alpha = 1),
                       st_buffers = list(lty = 1, lwd = 1, col = "black", alpha = 0.6),
                       st_classes = list(palette = "Spectral", fill = NULL, alpha = NULL,
                                         na_value = c("white", 1)),
                       st_ext = c("chartreuse", "firebrick1")){

  # Check arguments
  if(!is(x, "MultiLand")) stop("- argument 'x' must be an object of class 'MultiLand'.")
  environment(.mland_plot_check_args) <- environment()
  chk <- .mland_plot_check_args()
  if(length(chk[[1]]) > 0)
    for(w in 1:length(chk[[1]])){
      warning(strwrap(chk[[1]], prefix = "\n", initial = ""), call. = FALSE)
    }
  if(length(chk[[2]]) > 0){
    errors <- chk[[2]]
    stop(strwrap(errors, prefix = "\n", initial = "\n"))
  } else {
    objs <- names(chk)
    for(i in 3:length(chk)){ assign(objs[i], chk[[i]]) }
  }

  df_reference <- x@l_ref

  # if points and/or radii are null, take all of points and radii defined in x
  if(is.null(points)){
    points <- 1:length(x@points)
  } else { points <- as.numeric(points) }
  if(is.null(radii)) radii <- x@radii

  df_reference <- df_reference[df_reference$point_id %in% points &
                                 df_reference$radius %in% radii, ]

  points <- sort(unique(df_reference$point_id))
  radii  <- sort(unique(df_reference$radius))
  total_points <- length(points)

  # Check grid size
  if(length(nrow*ncol) > 0){
    if(nrow*ncol < total_points){
      ncol <- nrow <- NULL
      warning(strwrap("- the size of the grid define through arguments ncol and nrow is smaller than
                      the number of points to be plotted. Default NULL for both ncol and nrow was
                      taken.", prefix = "\n", initial = ""), call. = FALSE)
    }
  }

  # Asks if it is okay to plot so many plots
  if(total_points > 100){
    ask <- askYesNo("You are attempting to plot more than one hundred plots. Are you sure?")
    if(is.na(ask) | !ask) stop("Operation cancelled")
  }

  plots      <- vector("list", total_points + 1)
  max_radius <- max(radii)
  ext_range <- NULL
  # Plot local landscape of each point
  for(p in 1:total_points){
    pp <- points[p]
    pos <- df_reference[df_reference$point_id == pp & df_reference$radius == max_radius, "row_id"]
    pos_buffers <- df_reference[df_reference$point_id == pp & df_reference$radius %in% radii,
                                "row_id"]

    if(arg_name == "ext_raster"){
      t_slot <- x@landscapes$ext_rasters
    } else {
      t_slot <- x@landscapes$lsm_rasters
    }

    if(!x@onthefly){
      sub_int <- t_slot[[raster]][[pos]]
    } else {
      clip <- suppressWarnings(tryCatch(terra::crop(t_slot[[raster]],
                                                     terra::ext(x@buffers[pos, ])), error = c))
      if(!is.list(clip)){
        sub_int <- terra::mask(clip, x@buffers[pos, ])
      } else {
        empty_raster <- terra::rast(nrows = 1, ncols = 1, crs = terra::crs(x@buffers[1, ]), vals = NA)
        sub_int <- empty_raster
      }
    }

    if(arg_name == "ext_raster") ext_range <- get_ext_range(sub_int, ext_range)

    # Generate plot
    plots[[p]] <-  ggplot2::ggplot()

    # Plot all classes of each point
    plots[[p]] <- .plot_classes(plots[[p]], sub_int, arg_name, st_classes)

    # Plot buffers, point, label and extras
    if(title == "sitename"){
      tit <- unique(df_reference[df_reference$point_id == points[p], "site"])
    } else {
      tit <- points[p]
    }
    plots[[p]] <- .plot_bpl(x, plots[[p]], p, pos_buffers, points, tit, radii, st_buffers,
                            st_points)
  }

  # Plot legend: extract legend from artificial plot
  if(arg_name == "ext_raster"){
    for(i in 1:(length(plots)-1)){
      plots[[i]] <- plots[[i]] + ggplot2::scale_fill_gradient(low = st_ext[1], high = st_ext[2],
                                                     na.value = "transparent", limits = ext_range)
    }
    df <- data.frame(value = ext_range[1]:ext_range[2])
    legend_name <- x@layer_names[[2]][x@layer_names[[2]]$rasterlayer == raster, "name"]
    if(is.na(legend_name)) legend_name <- "value"
    art_plot <- ggplot2::ggplot(data = df, ggplot2::aes(x = 0, y = value, fill = value)) +
      ggplot2::geom_point() +
      ggplot2::scale_fill_gradient(low = st_ext[1], high = st_ext[2],
                                  na.value = "transparent", limits = ext_range,
                                  name = legend_name)
    legend   <- .g_legend(art_plot)
  } else {
    if(all(is.na(x@classes$classname))){
      cl_names <- x@classes[x@classes$rasterlayer == raster, "class"]
    } else {
      cl_names <- x@classes[x@classes$rasterlayer == raster, "classname"]
    }
    art_pos  <- rep(0, length(cl_names))
    df <- data.frame(art_pos, classes = cl_names)
    df$classes <- factor(df$classes, levels = cl_names)
    legend_name <- x@layer_names[[1]][x@layer_names[[1]]$rasterlayer == raster, "name"]
    if(is.na(legend_name)) legend_name <- "classes"
    art_plot <- ggplot2::ggplot(data = df, ggplot2::aes(art_pos, fill = classes, alpha = classes)) +
      ggplot2::geom_bar() +
      ggplot2::scale_fill_manual(values = st_classes$fill[seq(2, length(st_classes$fill), 2)],
                                 name = legend_name) +
      ggplot2::scale_alpha_manual(values = st_classes$alpha[seq(2, length(st_classes$alpha), 2)],
                                  name = legend_name)
    legend   <- .g_legend(art_plot)
  }

  # Only plot legend if actually are any intersections
  if(!is.null(legend)){
    plots[[length(plots)]] <- legend
  } else {
    plots[[length(plots)]] <- ggplot2::ggplot() +
      ggplot2::annotate("text", x = 4, y = 25, size = 4, label = "No intersections to be plotted") +
      ggplot2::labs(x = "", y ="") +
      ggplot2::theme(axis.ticks.y = ggplot2::element_blank(), axis.text.y = ggplot2::element_blank(),
                     axis.ticks.x = ggplot2::element_blank(), axis.text.x = ggplot2::element_blank(),
                     panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
                     panel.background = ggplot2::element_rect(fill = "white"))
  }

  pps <- gridExtra::grid.arrange(grobs = plots, nrow = nrow, ncol = ncol)
  invisible(pps)
}