R/ggplot.R

In inlabru-org/inlabru: Bayesian Latent Gaussian Modelling using INLA and Extensions

#' Plot a map using extent of a spatial object
#'
#' Uses `ggmap::get_map()` to query map services like Google Maps for a region centered around
#' the spatial object provided. Then calls `ggmap()` to plot the map.
#'
#' This function requires the `ggmap` package.
#'
#' @aliases gmap
#' @name gmap
#' @export
#' @param data A Spatial* object.
#' @param ... Arguments passed on to get_map().
#' @return a ggplot object
#'
#' @examples
#' \dontrun{
#' if (requireNamespace("ggmap", quietly = TRUE) &&
#'   require("ggplot2", quietly = TRUE)) {
#'   # Load the Gorilla data
#'   data(gorillas, package = "inlabru")
#'
#'   # Create a base map centred around the nests and plot the boundary as well
#'   # as the nests
#'   gmap(gorillas$nests, maptype = "satellite") +
#'     gm(gorillas$boundary) +
#'     gm(gorillas$nests, color = "white", size = 0.5)
#' }
#' }
gmap <- function(data, ...) {
  data <- fm_transform(data, crs = fm_crs("longlat_globe"))
  df <- cbind(coordinates(data), data@data)

  # Figure out a sensible bounding box (range of data plus 30%)
  extend <- 2.5
  crange <- apply(coordinates(data), MARGIN = 2, range)
  lonlim <- (extend * (crange[, 1] - mean(crange[, 1]))) + mean(crange[, 1])
  latlim <- (extend * (crange[, 2] - mean(crange[, 2]))) + mean(crange[, 2])

  # Create map
  requireNamespace("ggmap")
  myMap <-
    ggmap::get_map(
      c(lonlim[1], latlim[1], lonlim[2], latlim[2]),
      ...
    )

  # Return map
  ggmap::ggmap(myMap)
}

#' ggplot2 geomes for inlabru related objects
#'
#' @description
#' gg is a generic function for generating geomes from various kinds of spatial objects, e.g. Spatial* data,
#' meshes, Raster objects and inla/inlabru predictions. The function invokes particular methods which depend
#' on the [class] of the first argument.
#'
#' @name gg
#' @export
#' @param data an object for which to generate a geom.
#' @param ... Arguments passed on to the geom method.
#' @return The form of the value returned by gg depends on the class of its argument. See the documentation of the particular methods for details of what is produced by that method.
#'
#' @family geomes for inla and inlabru predictions
#' @family geomes for spatial data
#' @family geomes for meshes
#' @family geomes for Raster data
#'
#' @examples
#' if (require("ggplot2", quietly = TRUE)) {
#'   # Load Gorilla data
#'
#'   data(gorillas, package = "inlabru")
#'
#'   # Invoke ggplot and add geomes for the Gorilla nests and the survey boundary
#'
#'   ggplot() +
#'     gg(gorillas$boundary) +
#'     gg(gorillas$nests)
#' }
gg <- function(data, ...) {
  UseMethod("gg")
}


#' ggplot geom for spatial data
#'
#' gm is a wrapper for the [gg] method. It will take the first argument and transform its
#' coordinate system to latitude and longitude. Thereafter, [gg] is called using the transformed
#' data and the arguments provided via `...`. gm is intended to replace gg whenever the
#' data is supposed to be plotted over a spatial map generated by [gmap], which only works
#' if the coordinate system is latitude/longitude.
#'
#' @aliases gm
#' @name gm
#' @export
#' @param data an object for which to generate a geom.
#' @param ... Arguments passed on to [gg()].
#' @return The form of the value returned by gm depends on the class of its argument. See the documentation of the particular methods for details of what is produced by that method.
#' @family geomes for inla and inlabru predictions
#' @family geomes for spatial data
#' @family geomes for meshes
#' @family geomes for Raster data
#'
#' @examples
#' \dontrun{
#' if (require("ggplot2", quietly = TRUE)) {
#'   # Load the Gorilla data
#'   data(gorillas, package = "inlabru")
#'
#'   # Create a base map centered around the nests and plot the boundary as well as the nests
#'   gmap(gorillas$nests, maptype = "satellite") +
#'     gm(gorillas$boundary) +
#'     gm(gorillas$nests, color = "white", size = 0.5)
#' }
#' }
#'
gm <- function(data, ...) {
  gg(data, crs = CRS("+proj=longlat"), ...)
}


#' Geom for matrix
#'
#' Creates a tile geom for plotting a matrix
#'
#' Requires the `ggplot2` package.
#'
#' @name gg.matrix
#' @export
#' @param data A `matrix` object.
#' @param mapping a set of aesthetic mappings created by `aes`. These are passed on to `geom_tile`.
#' @param ... Arguments passed on to `geom_tile`.
#' @return A `geom_tile` with reversed y scale.
#' @family geomes for inla and inlabru predictions
#'
#' @importFrom rlang .data
#' @examples
#' if (require("ggplot2", quietly = TRUE)) {
#'   A <- matrix(runif(100), nrow = 10)
#'   ggplot() +
#'     gg(A)
#' }
gg.matrix <- function(data, mapping = NULL, ...) {
  requireNamespace("ggplot2")
  A <- as.matrix(data)
  grd <- expand.grid(row = seq_len(nrow(A)), column = seq_len(ncol(A)))
  df <- data.frame(value = as.vector(A), grd)

  dmap <- ggplot2::aes(x = .data$column, y = .data$row, fill = .data$value)

  if (!is.null(mapping)) {
    dmap <- modifyList(dmap, mapping)
  }

  ggp <- c(
    ggplot2::geom_tile(dmap, data = df, ...),
    ggplot2::scale_y_reverse()
  )
}


#' Geom for data.frame
#'
#' This geom constructor will simply call [gg.bru_prediction()] for the data provided.
#'
#' Requires the `ggplot2` package.
#'
#' @name gg.data.frame
#' @export
#' @param ... Arguments passed on to [gg.bru_prediction()].
#' @return Concatenation of a `geom_line` value and optionally a `geom_ribbon` value.
#' @family geomes for inla and inlabru predictions
#' @example inst/examples/gg.prediction.R

gg.data.frame <- function(...) {
  gg.bru_prediction(...)
}


#' Geom for predictions
#'
#' @description
#'
#' This geom serves to visualize `prediction` objects which usually results from a call to
#' [predict.bru()]. Predictions objects provide summary statistics (mean, median, sd, ...) for
#' one or more random variables. For single variables (or if requested so by setting `bar = TRUE`),
#' a boxplot-style geom is constructed to show the statistics. For multivariate predictions the
#' mean of each variable (y-axis) is plotted agains the row number of the varriable in the prediction
#' data frame (x-axis) using `geom_line`. In addition, a `geom_ribbon` is used to show
#' the confidence interval.
#'
#' Note: `gg.bru_prediction` also understands the format of INLA-style posterior summaries, e.g.
#' `fit$summary.fixed` for an inla object `fit`
#'
#' Requires the `ggplot2` package.
#'
#' @name gg.bru_prediction
#' @export
#' @importFrom utils modifyList
#' @param data A prediction object, usually the result of a [predict.bru()] call.
#' @param mapping a set of aesthetic mappings created by `aes`. These are passed on to `geom_line`.
#' @param ribbon If TRUE, plot a ribbon around the line based on the smalles and largest
#'   quantiles present in the data, found by matching names starting with `q` and
#'   followed by a numerical value.  `inla()`-style `numeric+"quant"` names are converted
#'   to inlabru style before matching.
#' @param alpha The ribbons numeric alpha (transparency) level in `[0,1]`.
#' @param bar If TRUE plot boxplot-style summary for each variable.
#' @param \dots Arguments passed on to `geom_line`.
#' @return Concatenation of a `geom_line` value and optionally a `geom_ribbon` value.
#' @family geomes for inla and inlabru predictions
#' @example inst/examples/gg.prediction.R

gg.bru_prediction <- function(data, mapping = NULL, ribbon = TRUE, alpha = NULL, bar = FALSE, ...) {
  if (inherits(data, c("Spatial", "SpatRaster", "RasterLayer", "sf"))) {
    return(NextMethod())
  }

  if (is.null(alpha)) {
    alpha <- 0.3
  }

  requireNamespace("ggplot2")
  # Convert from old and inla style names
  if (("median" %in% names(data)) &&
    !("0.5quant" %in% names(data)) &&
    !("q0.5" %in% names(data))) {
    names(data)[names(data) == "median"] <- "q0.5"
  }
  new_quant_names <- list(
    q0.025 = "0.025quant",
    q0.5 = "0.5quant",
    q0.975 = "0.975quant"
  )

  for (quant_name in names(new_quant_names[new_quant_names %in% names(data)])) {
    names(data)[names(data) == new_quant_names[[quant_name]]] <- quant_name
  }
  # Find quantile levels
  quant_names <- names(data)[grepl("^q[01]\\.?[0-9]*$", names(data))]
  if (length(quant_names) > 0) {
    quant_probs <- as.numeric(sub("^q", "", quant_names))
    quant_names <- quant_names[order(quant_probs)]
    quant_probs <- quant_probs[order(quant_probs)]
    lqname <- quant_names[1]
    uqname <- quant_names[length(quant_names)]
  }

  if (bar || (nrow(data) == 1)) {
    sz <- 10 # bar width
    med_sz <- 4 # median marker size

    data <- cbind(
      data.frame(data),
      data.frame(
        variable = rownames(data),
        summary = data$mean,
        sdmax = data$mean + data$sd,
        sdmin = data$mean - data$sd
      )
    )

    geom <- c(
      ggplot2::geom_point(
        data = data,
        mapping = ggplot2::aes(
          x = .data$variable,
          y = .data$summary,
          color = .data$variable
        ),
        shape = 95, size = 0
      )
    )
    if (length(quant_names) > 0) {
      geom <- c(
        geom,
        ggplot2::geom_segment(
          data = data,
          mapping = ggplot2::aes(
            y = .data[[lqname]],
            yend = .data[[uqname]],
            x = .data$variable,
            xend = .data$variable,
            color = .data$variable
          ),
          linewidth = sz
        )
      )
    }

    # Mean and median
    geom <- c(
      geom,
      ggplot2::geom_point(
        data = data,
        mapping = ggplot2::aes(
          x = .data$variable,
          y = .data$mean
        ),
        color = "black", shape = 20, size = med_sz
      )
    )
    if ("q0.5" %in% quant_names) {
      geom <- c(
        geom,
        ggplot2::geom_point(
          data = data,
          mapping = ggplot2::aes(
            x = .data$variable,
            y = .data$q0.5
          ),
          color = "black", shape = 3, size = sz * 2 / 3
        )
      )
    }
    geom <- c(
      geom,
      ggplot2::coord_flip()
    )
  } else {
    if ("pdf" %in% names(data)) {
      y.str <- "pdf"
      ribbon <- FALSE
    } else if ("mean" %in% names(data)) {
      y.str <- "mean"
    } else if ("q0.5" %in% names(data)) {
      y.str <- "q0.5"
    } else if ("median" %in% names(data)) {
      y.str <- "q0.5"
    } else {
      stop("Prediction has neither mean nor median/q0.5 or pdf as column. Don't know what to plot.")
    }

    line.map <- ggplot2::aes(
      x = .data[[names(data)[1]]],
      y = .data[[y.str]]
    )

    if (!is.null(mapping)) {
      line.map <- utils::modifyList(line.map, mapping)
    }

    geom <- ggplot2::geom_line(data = data, line.map, ...)

    if (ribbon && length(quant_names) > 0) {
      ribbon.map <- ggplot2::aes(
        x = .data[[names(data)[1]]],
        ymin = .data[[lqname]],
        ymax = .data[[uqname]]
      )
      # Use line color for ribbon filling
      if ("colour" %in% names(line.map)) {
        ribbon.map <- modifyList(
          ribbon.map,
          ggplot2::aes(fill = .data[[line.map[["colour"]]]])
        )
      }
      geom <- c(geom, ggplot2::geom_ribbon(data = data, ribbon.map, alpha = alpha))
    }
  }
  geom
}

#' @rdname gg.bru_prediction
#' @export
gg.prediction <- function(data, ...) {
  gg.bru_prediction(data = data, ...)
}

#' Geom for SpatialPoints objects
#'
#' This function coerces the `SpatialPoints` into a `data.frame` and uses `geom_point`
#' to plot the points. Requires the `ggplot2` package.
#'
#' @aliases gg.SpatialPoints
#' @name gg.SpatialPoints
#' @export
#' @param data A SpatialPoints object.
#' @param mapping Aesthetic mappings created by `aes` used to update the default
#'                mapping. The default mapping is
#'                `ggplot2::aes(x = .data[[coordnames(data)[1]]],
#'                 y = .data[[coordnames(data)[2]]])`.
#' @param crs A [CRS] object defining the coordinate system to project the data to before plotting.
#' @param ... Arguments passed on to `geom_point`.
#' @return A `geom_point` return value
#' @family geomes for spatial data
#'
#' @example inst/examples/gg.spatial.R

gg.SpatialPoints <- function(data, mapping = NULL, crs = NULL, ...) {
  requireNamespace("ggplot2")
  if (!is.null(crs)) {
    data <- fm_transform(data, crs)
  }

  df <- as.data.frame(data)
  cnames <- coordnames(data)
  if (is.null(cnames)) {
    cnames <- c("x", "y")
  }

  dmap <- ggplot2::aes(
    x = .data[[cnames[1]]],
    y = .data[[cnames[2]]]
  )

  if (!is.null(mapping)) {
    dmap <- modifyList(dmap, mapping)
  }

  ggplot2::geom_point(data = df, mapping = dmap, ...)
}

#' Geom helper for sf objects
#'
#' This function uses `geom_sf()`, unless overridden by the geom argument.
#' Requires the `ggplot2` package.
#'
#' @name gg.sf
#' @export
#' @param data An `sf` object.
#' @param mapping Default mapping is `ggplot2::aes(geometry = ...)`,
#' where the geometry name is obtained from `attr(data, "sf_column")`.
#' This is merged with the user supplied mapping.
#' @param geom Either "sf" (default) or "tile". For "tile", uses
#' `geom_tile(..., stat = "sf_coordinates")`, intended for converting point data
#' to grid tiles with the `fill` aesthetic, which is by default set to the first
#' data column.
#' @param ... Arguments passed on to `geom_sf` or `geom_tile`.
#' @return A ggplot return value
#' @family geomes for spatial data

gg.sf <- function(data, mapping = NULL, ..., geom = "sf") {
  requireNamespace("ggplot2")

  sf_column <- attr(data, "sf_column")
  geom <- match.arg(geom, c("sf", "tile"))

  if (identical(geom, "sf")) {
    dmap <- ggplot2::aes(
      geometry = .data[[sf_column]]
    )
    if (!is.null(mapping)) {
      dmap <- modifyList(dmap, mapping)
    }
    result <- ggplot2::geom_sf(data = data, mapping = dmap, ...)
  } else {
    data_column <- setdiff(names(data), sf_column)
    if (length(data_column) > 0) {
      data_column <- data_column[1]
      dmap <- ggplot2::aes(
        geometry = .data[[sf_column]],
        fill = .data[[data_column]]
      )
    } else {
      dmap <- ggplot2::aes(
        geometry = .data[[sf_column]]
      )
    }
    if (!is.null(mapping)) {
      dmap <- modifyList(dmap, mapping)
    }
    result <-
      c(
        ggplot2::geom_tile(
          data = data,
          mapping = dmap,
          stat = "sf_coordinates",
          ...
        ),
        if (fmesher::fm_crs_is_null(fm_crs(data))) {
          ggplot2::coord_sf(default = TRUE)
        } else {
          ggplot2::coord_sf(default = TRUE, crs = fm_crs(data))
        }
      )
  }

  result
}

#' Geom for SpatialLines objects
#'
#' @description
#'
#' Extracts start and end points of the lines and calls `geom_segment` to plot
#' lines between them. Requires the `ggplot2` package.
#'
#' @aliases gg.SpatialLines
#' @name gg.SpatialLines
#' @export
#' @param data A `SpatialLines` or `SpatialLinesDataFrame` object.
#' @param mapping Aesthetic mappings created by `ggplot2::aes` or `ggplot2::aes_`
#'                used to update the default
#'                mapping. The default mapping is
#'                `ggplot2::aes(x = .data[[coordnames(data)[1]]],
#'                     y = .data[[coordnames(data)[2]]],
#'                     xend = .data[[paste0("end.", coordnames(data)[1])]],
#'                     yend = .data[[paste0("end.", coordnames(data)[2])]])`.
#' @param crs A `CRS` object defining the coordinate system to project the data
#'   to before plotting.
#' @param ... Arguments passed on to `ggplot2::geom_segment`.
#' @return A `geom_segment`` return value.
#' @family geomes for spatial data
#' @example inst/examples/gg.spatial.R

gg.SpatialLines <- function(data, mapping = NULL, crs = NULL, ...) {
  requireNamespace("ggplot2")
  if (!is.null(crs)) {
    data <- fm_transform(data, crs)
  }

  qq <- coordinates(data)
  cnames <- coordnames(data)
  if (is.null(cnames)) {
    cnames <- c("x", "y")
  }
  sp <- do.call(rbind, lapply(
    qq,
    function(k) do.call(rbind, lapply(k, function(x) x[1:(nrow(x) - 1), ]))
  ))
  ep <- do.call(rbind, lapply(
    qq,
    function(k) do.call(rbind, lapply(k, function(x) x[2:(nrow(x)), ]))
  ))
  colnames(sp) <- cnames
  colnames(ep) <- paste0("end.", cnames)
  if (!inherits(data, "SpatialLinesDataFrame")) {
    df <- data.frame(cbind(sp, ep))
  } else {
    df <- data.frame(
      cbind(sp, ep),
      data@data[
        rep(
          seq_len(NROW(data)),
          times = vapply(
            qq,
            FUN = function(x) {
              NROW(x[[1]]) - 1
            },
            0
          )
        ), ,
        drop = FALSE
      ]
    )
  }
  dmap <- ggplot2::aes(
    x = .data[[cnames[1]]],
    y = .data[[cnames[2]]],
    xend = .data[[paste0("end.", cnames[1])]],
    yend = .data[[paste0("end.", cnames[2])]]
  )

  if (!is.null(mapping)) {
    dmap <- modifyList(dmap, mapping)
  }

  ggplot2::geom_segment(data = df, mapping = dmap, ...)
}

#' Geom for SpatialPolygons objects
#'
#' Uses the `ggplot2::fortify()` function to turn the `SpatialPolygons` objects into a
#' `data.frame`. Then
#' calls `geom_polygon` to plot the polygons. Requires the `ggplot2` package.
#'
#'
#' @aliases gg.SpatialPolygons
#' @name gg.SpatialPolygons
#' @export
#' @param data A `SpatialPolygons` or `SpatialPolygonsDataFrame` object.
#' @param mapping Aesthetic mappings created by `aes` used to update the default
#'                mapping.
#' @param crs A `CRS` object defining the coordinate system to project the data to before plotting.
#' @param ... Arguments passed on to `geom_sf`.
#' Unless specified by the user,
#' the argument `alpha = 0.2` (alpha level for polygon filling) is added.
#' @return A `geom_sf` object.
#' @details Up to version `2.10.0`, the `ggpolypath` package was used to ensure
#' proper plotting, since the `ggplot2::geom_polygon` function doesn't always
#' handle geometries with holes properly. After `2.10.0`, the object is converted
#' to `sf` format and passed on to [gg.sf()] instead, as `ggplot2` version `3.4.4`
#' deprecated the intenrally used `ggplot2::fortify()` method for
#' `SpatialPolygons/DataFrame` objects.
#' @family geomes for spatial data
#' @example inst/examples/gg.spatial.R

gg.SpatialPolygons <- function(data, mapping = NULL, crs = NULL, ...) {
  data <- sf::st_as_sf(data)
  if (!is.null(crs)) {
    data <- fm_transform(data, crs, passthrough = TRUE)
  }
  arg <- list(...)
  if ("alpha" %in% names(arg)) {
    gg(data, mapping = mapping, ...)
  } else {
    gg(data, mapping = mapping, ..., alpha = 0.2)
  }
}

#' Geom for SpatialGridDataFrame objects
#'
#' Coerces input `SpatialGridDataFrame` to `SpatialPixelsDataFrame` and calls
#' [gg.SpatialPixelsDataFrame()] to plot it.
#' Requires the `ggplot2` package.
#'
#' @aliases gg.SpatialGridDataFrame
#' @name gg.SpatialGridDataFrame
#' @export
#' @param data A SpatialGridDataFrame object.
#' @param ... Arguments passed on to [gg.SpatialPixelsDataFrame()].
#' @return A `geom_tile` value.
#' @family geomes for spatial data
#' @example inst/examples/gg.spatial.R

gg.SpatialGridDataFrame <- function(data, ...) {
  data <- as(data, "SpatialPixelsDataFrame")
  gg(data, ...)
}

#' Geom for SpatialPixelsDataFrame objects
#'
#' Coerces input SpatialPixelsDataFrame to data.frame and uses `geom_tile` to plot it.
#' Requires the `ggplot2` package.
#'
#' @aliases gg.SpatialPixelsDataFrame
#' @name gg.SpatialPixelsDataFrame
#' @export
#' @param data A SpatialPixelsDataFrame object.
#' @param mapping Aesthetic mappings created by `aes` used to update the default
#'                mapping. The default mapping is `ggplot2::aes(x = .data[[coordnames(data)[1]]],
#'                y = .data[[coordnames(data)[2]]], fill = .data[[names(data)[[1]]]])`.
#' @param crs A [CRS] object defining the coordinate system to project the data to before plotting.
#' @param mask A SpatialPolygon defining the region that is plotted.
#' @param ... Arguments passed on to `geom_tile`.
#' @return A `geom_tile` return value.
#' @family geomes for spatial data
#' @example inst/examples/gg.spatial.R

gg.SpatialPixelsDataFrame <- function(data,
                                      mapping = NULL,
                                      crs = NULL,
                                      mask = NULL, ...) {
  if (!is.null(crs)) {
    data <- fm_transform(data, crs)
  }
  if (!is.null(mask)) {
    data <- data[as.vector(!is.na(over(data, mask))), ]
  }

  df <- as.data.frame(data)

  if (length(names(data)) == 0) {
    dmap <- ggplot2::aes(
      x = .data[[coordnames(data)[1]]],
      y = .data[[coordnames(data)[2]]]
    )
  } else {
    dmap <- ggplot2::aes(
      x = .data[[coordnames(data)[1]]],
      y = .data[[coordnames(data)[2]]],
      fill = .data[[names(data)[[1]]]]
    )
  }

  if (!is.null(mapping)) {
    dmap <- modifyList(dmap, mapping)
  }
  gm <- ggplot2::geom_tile(data = df, mapping = dmap, ...)

  gm
}


#' Geom for SpatialPixels objects
#'
#' Uses `geom_point` to plot the pixel centers.
#' Requires the `ggplot2` package.
#'
#' @aliases gg.SpatialPixels
#' @name gg.SpatialPixels
#' @export
#' @param data A [SpatialPixels] object.
#' @param ... Arguments passed on to `geom_tile`.
#' @return A `geom_tile` return value.
#' @family geomes for spatial data
#' @examples
#' if (require("ggplot2", quietly = TRUE) &&
#'   bru_safe_sp()) {
#'   # Load Gorilla data
#'
#'   data(gorillas, package = "inlabru")
#'
#'   # Turn elevation covariate into SpatialPixels
#'   pxl <- sp::SpatialPixels(sp::SpatialPoints(gorillas$gcov$elevation))
#'
#'   # Plot the pixel centers
#'   ggplot() +
#'     gg(pxl, size = 0.1)
#' }
gg.SpatialPixels <- function(data, ...) {
  gg(SpatialPoints(data), ...)
}


#' Geom wrapper for SpatRaster objects
#'
#' Convenience wrapper function for `tidyterra::geom_spatraster()`.
#' Requires the `ggplot2` and `tidyterra` packages.
#'
#' @aliases gg.SpatRaster
#' @name gg.SpatRaster
#' @export
#' @param data A SpatRaster object.
#' @param ... Arguments passed on to `geom_spatraster`.
#' @return The output from `geom_spatraster.
#' @family geomes for spatial data

gg.SpatRaster <- function(data, ...) {
  if (!requireNamespace("tidyterra", quietly = TRUE)) {
    stop(paste0(
      "gg.SpatRaster requires the 'tidyterra' package, ",
      "which is not installed.\n",
      "Please install it and try again!"
    ))
  }
  if (".mask" %in% names(data)) {
    data <-
      terra::mask(
        data,
        mask = data$.mask,
        maskvalues = FALSE,
        updatevalue = NA
      )
  }
  tidyterra::geom_spatraster(data = data, ...)
}





#' Geom for inla.mesh objects
#'
#' @description
#'
#' This function extracts the graph of an inla.mesh object and uses `geom_line` to visualize
#' the graph's edges. Alternatively, if the `color` argument is provided, interpolates the colors
#' across for a set of SpatialPixels covering the mesh area and calls [gg.SpatialPixelsDataFrame()]
#' to plot the interpolation.
#' Requires the `ggplot2` package.
#'
#' @aliases gg.fm_mesh_2d
#' @name gg.fm_mesh_2d
#' @export

#' @param data An `fm_mesh_2d` object.
#' @param color A vector of scalar values to fill the mesh with colors.
#' The length of the vector mus correspond to the number of mesh vertices.
#' The alternative name `colour` is also recognised.
#' @param alpha A vector of scalar values setting the alpha value of the colors provided.
#' @param edge.color Color of the regular mesh edges.
#' @param edge.linewidth Line width for the regular mesh edges. Default 0.25
#' @param interior If TRUE, plot the interior boundaries of the mesh.
#' @param int.color Color used to plot the interior constraint edges.
#' @param int.linewidth Line width for the interior constraint edges. Default 0.5
#' @param exterior If TRUE, plot the exterior boundaries of the mesh.
#' @param ext.color Color used to plot the exterior boundary edges.
#' @param ext.linewidth Line width for the exterior boundary edges. Default 1
#' @param crs A CRS object supported by [fm_transform()] defining the coordinate
#' system to project the mesh to before plotting.
#' @param nx Number of pixels in x direction (when plotting using the color parameter).
#' @param ny Number of pixels in y direction (when plotting using the color parameter).
#' @param mask A SpatialPolygon defining the region that is plotted.
#' @param ... ignored arguments (S3 generic compatibility).
#' @return `geom_line` return values or, if the color argument is used, the
#'   values of [gg.SpatialPixelsDataFrame()].
#' @family geomes for meshes
#' @example inst/examples/gg.inla.mesh.R

gg.fm_mesh_2d <- function(data,
                          color = NULL,
                          alpha = NULL,
                          edge.color = "grey",
                          edge.linewidth = 0.25,
                          interior = TRUE,
                          int.color = "blue",
                          int.linewidth = 0.5,
                          exterior = TRUE,
                          ext.color = "black",
                          ext.linewidth = 1,
                          crs = NULL,
                          mask = NULL,
                          nx = 500, ny = 500,
                          ...) {
  requireNamespace("INLA")
  requireNamespace("ggplot2")
  if (is.null(color) && ("colour" %in% names(list(...)))) {
    color <- list(...)[["colour"]]
  }
  if (!is.null(color)) {
    if (inherits(mask, "Spatial")) {
      # For backwards compatibility with old plotting code, use Spatial
      # format, and mask only in the gg.Spatial* method.
      px <- fm_pixels(data, dims = c(nx, ny), format = "sp")
      A <- fm_evaluator(data, px)$proj$A
      px$color <- as.vector(A %*% color)
      if (!is.null(alpha)) {
        px$alpha <- as.vector(A %*% alpha)
        gg <- gg(px, mask = mask, alpha = "alpha")
      } else {
        gg <- gg(px, mask = mask)
      }
    } else {
      px <- fm_pixels(data, dims = c(nx, ny), mask = mask, format = "sf")
      proj <- fm_evaluator(data, px)
      px$color <- fm_evaluate(proj, field = color)
      if (!is.null(alpha)) {
        if (length(alpha) == 1) {
          px$alpha <- alpha
        } else {
          px$alpha <- fm_evaluate(proj, field = alpha)
        }
        gg <- gg(px,
          ggplot2::aes(fill = .data[["color"]]),
          alpha = px[["alpha"]],
          geom = "tile"
        )
      } else {
        gg <- gg(px,
          ggplot2::aes(fill = .data[["color"]]),
          geom = "tile"
        )
      }
    }

    return(gg)
  }

  if (data$manifold == "S2") {
    stop("Geom not implemented for spherical meshes (manifold = S2)")
  }
  if (!is.null(crs)) {
    data <- fm_transform(data, crs = crs)
  }

  df <- rbind(
    data.frame(a = data$loc[data$graph$tv[, 1], c(1, 2)], b = data$loc[data$graph$tv[, 2], c(1, 2)]),
    data.frame(a = data$loc[data$graph$tv[, 2], c(1, 2)], b = data$loc[data$graph$tv[, 3], c(1, 2)]),
    data.frame(a = data$loc[data$graph$tv[, 1], c(1, 2)], b = data$loc[data$graph$tv[, 3], c(1, 2)])
  )

  colnames(df) <- c("x", "y", "xend", "yend")
  mp <- ggplot2::aes(x = .data$x, y = .data$y, xend = .data$xend, yend = .data$yend)
  msh <- ggplot2::geom_segment(data = df, mapping = mp, color = edge.color, linewidth = edge.linewidth)

  # Outer boundary
  if (exterior) {
    df <- data.frame(
      data$loc[data$segm$bnd$idx[, 1], 1:2],
      data$loc[data$segm$bnd$idx[, 2], 1:2]
    )
    colnames(df) <- c("x", "y", "xend", "yend")
    bnd <- ggplot2::geom_segment(data = df, mapping = mp, color = ext.color, linewidth = ext.linewidth)
  } else {
    bnd <- NULL
  }

  if (interior) {
    # Interior boundary
    df <- data.frame(
      data$loc[data$segm$int$idx[, 1], 1:2],
      data$loc[data$segm$int$idx[, 2], 1:2]
    )
    colnames(df) <- c("x", "y", "xend", "yend")
    if (nrow(df) == 0) {
      int <- NULL
    } else {
      int <- ggplot2::geom_segment(data = df, mapping = mp, color = int.color, linewidth = int.linewidth)
    }
  } else {
    int <- NULL
  }

  # Return combined geomes
  c(msh, bnd, int)
}


#' Geom for fm_mesh_1d objects
#'
#' This function generates a `geom_point` object showing the knots (vertices)
#' of a 1D mesh.
#' Requires the `ggplot2` package.
#'
#' @aliases gg.fm_mesh_1d gg.inla.mesh.1d
#' @name gg.fm_mesh_1d
#' @export
#' @param data An inla.mesh.1d object.
#' @param mapping aesthetic mappings created by `aes`. These are passed on to `geom_point`.
#' @param y Single or vector numeric defining the y-coordinates of the mesh knots to plot.
#' @param shape Shape of the knot markers.
#' @param ... parameters passed on to `geom_point`.
#' @return An object generated by `geom_point`.
#' @family geomes for meshes
#'
#' @examples
#' \donttest{
#' if (require("fmesher", quietly = TRUE) &&
#'   require("ggplot2", quietly = TRUE)) {
#'   # Create a 1D mesh
#'
#'   mesh <- fm_mesh_1d(seq(0, 10, by = 0.5))
#'
#'   # Plot it
#'
#'   ggplot() +
#'     gg(mesh)
#'
#'   # Plot it using a different shape and size for the mesh nodes
#'
#'   ggplot() +
#'     gg(mesh, shape = "|", size = 5)
#' }
#' }
#'
gg.fm_mesh_1d <- function(data, mapping = ggplot2::aes(.data[["x"]], .data[["y"]]),
                          y = 0, shape = 4, ...) {
  df <- data.frame(x = data$loc, y = y)
  ggplot2::geom_point(data = df, mapping = mapping, shape = shape, ...)
}

#' @describeIn gg.fm_mesh_2d Alias for `gg.fm_mesh_2d`, supporting `inla.mesh`
#' objects.
#' @export
gg.inla.mesh <- gg.fm_mesh_2d

#' @describeIn gg.fm_mesh_1d Alias for `gg.fm_mesh_1d`, supporting
#' `inla.mesh.1d` objects.
#' @export
gg.inla.mesh.1d <- gg.fm_mesh_1d



#' Geom for RasterLayer objects
#'
#' This function takes a RasterLayer object, converts it into a
#' `SpatialPixelsDataFrame` and uses `geom_tile` to plot the data.
#'
#' This function requires the `raster` and `ggplot2` packages.
#'
#' @aliases gg.RasterLayer
#' @name gg.RasterLayer
#' @export
#' @param data A RasterLayer object.
#' @param mapping aesthetic mappings created by `aes`. These are passed on to `geom_tile`.
#' @param ... Arguments passed on to `geom_tile`.
#' @return An object returned by `geom_tile`
#' @family geomes for Raster data
#'
#' @examples
#' \dontrun{
#' # Some features require the raster and spatstat.data packages.
#' if (require("spatstat.data", quietly = TRUE) &&
#'   require("raster", quietly = TRUE) &&
#'   require("ggplot2", quietly = TRUE)) {
#'   # Load Gorilla data
#'   data("gorillas", package = "spatstat.data")
#'
#'   # Convert elevation covariate to RasterLayer
#'
#'   elev <- as(gorillas.extra$elevation, "RasterLayer")
#'
#'   # Plot the elevation
#'
#'   ggplot() +
#'     gg(elev)
#' }
#' }
gg.RasterLayer <- function(data, mapping = ggplot2::aes(x = .data[["x"]], y = .data[["y"]], fill = .data[["layer"]]), ...) {
  requireNamespace("ggplot2")
  requireNamespace("raster")
  spdf <- as(data, "SpatialPixelsDataFrame")
  df <- as.data.frame(spdf)
  # head(r.df)
  # g <- ggplot(r.df, ggplot2::aes(x=x, y=y)) + geom_tile(ggplot2::aes(fill = layer)) + coord_equal()
  ggplot2::geom_tile(data = df, mapping = mapping, ...)
}

#' Plot method for posterior marginals estimated by bru
#'
#' @description
#'
#' [bru()] uses `INLA::inla()` to fit models. The latter estimates the posterior densities of
#' all random effects in the model. This function serves to access and plot the posterior
#' densities in a convenient way.
#'
#' Requires the `ggplot2` package.
#'
#' @method plot bru
#' @export
#' @param x a fitted [bru()] model.
#' @param ... A character naming the effect to plot, e.g. "Intercept". For random
#' effects, adding `index = ...` plots the density for a single component of the
#' latent model.
#' @return an object of class `gg`
#'
#'
#' @examples
#' \dontrun{
#' if (require("ggplot2", quietly = TRUE)) {
#'   # Generate some data and fit a simple model
#'   input.df <- data.frame(x = cos(1:10))
#'   input.df <- within(input.df, y <- 5 + 2 * cos(1:10) + rnorm(10, mean = 0, sd = 0.1))
#'   fit <- bru(y ~ x, family = "gaussian", data = input.df)
#'   summary(fit)
#'
#'   # Plot the posterior of the model's x-effect
#'   plot(fit, "x")
#' }
#' }
#'
plot.bru <- function(x, ...) {
  plotmarginal.inla(x, ...)
}

#' Plot prediction using ggplot2
#'
#' Generates a base ggplot2 using `ggplot()` and adds a geom for input `x` using [gg].
#'
#' Requires the `ggplot2` package.
#'
#' @name plot.bru_prediction
#' @param x a prediction object.
#' @param y Ignored argument but required for S3 compatibility.
#' @param ... Arguments passed on to [gg.prediction()].
#' @return an object of class `gg`
#' @example inst/examples/gg.prediction.R
#' @export
#' @method plot bru_prediction

plot.bru_prediction <- function(x, y = NULL, ...) {
  requireNamespace("ggplot2")
  ggplot2::ggplot() +
    gg(x, ...)
}

#' @rdname plot.bru_prediction
#' @export
#' @method plot prediction
plot.prediction <- function(x, y = NULL, ...) {
  requireNamespace("ggplot2")
  ggplot2::ggplot() +
    gg(x, ...)
}



#' @title Multiple ggplots on a page.
#'
#' @description
#' Renders multiple ggplots on a single page.
#'
#' @param ... Comma-separated `ggplot` objects.
#' @param plotlist A list of `ggplot` objects - an alternative to the comma-separated argument above.
#' @param cols Number of columns of plots on the page.
#' @param layout A matrix specifying the layout. If present, 'cols' is ignored. If the layout is
#' something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE), then plot 1 will go in the upper left,
#' 2 will go in the upper right, and 3 will go all the way across the bottom.
#'
#' @author David L. Borchers \email{dlb@@st-andrews.ac.uk}
#'
#' @source
#' <http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/>
#'
#' @examples
#' if (require("ggplot2", quietly = TRUE)) {
#'   df <- data.frame(x = 1:10, y = 1:10, z = 11:20)
#'   pl1 <- ggplot(data = df) +
#'     geom_line(mapping = aes(x, y), color = "red")
#'   pl2 <- ggplot(data = df) +
#'     geom_line(mapping = aes(x, z), color = "blue")
#'   multiplot(pl1, pl2, cols = 2)
#' }
#' @export
#
multiplot <- function(..., plotlist = NULL, cols = 1, layout = NULL) {
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots <- length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots / cols)),
      ncol = cols, nrow = ceiling(numPlots / cols)
    )
  }

  if (numPlots == 1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid::grid.newpage()
    grid::pushViewport(grid::viewport(layout = grid::grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = grid::viewport(
        layout.pos.row = matchidx$row,
        layout.pos.col = matchidx$col
      ))
    }
  }
}



# Default color palette for plots using \link{gg}
#
# @aliases bru.pal
# @name bru.pal
# @export
# @param
# @return Color values as returned by \code{RColorBrewer::brewer.pal(9, "YlOrRd")}

bru.pal <- function() {
  # library(RColorBrewer)
  # brewer.pal(9, "YlOrRd")
  pal <- c("#FFFFCC", "#FFEDA0", "#FED976", "#FEB24C", "#FD8D3C", "#FC4E2A", "#E31A1C", "#BD0026", "#800026")
}