R/plotsample.R

In inlabru: Bayesian Latent Gaussian Modelling using INLA and Extensions

# @title Make rectangular \code{SpatialPolygonsDataFrame}.
#
# @description
# Makes a rectangular \code{SpatialPolygonsDataFrame} using its bottom left
# coordinates, its width and its height.
#
# @param start Bottom left coordinates (x,y).
# @param width Width of base of the rectangle.
# @param height Height of the rectangle.
#
# @return A \code{SpatialPolygonsDataFrame}.
#
# @export

makepoly <- function(start, width, height) {
  poly <- matrix(
    c(
      start[1], start[2],
      start[1], start[2] + height,
      start[1] + width, start[2] + height,
      start[1] + width, start[2],
      start[1], start[2]
    ),
    ncol = 2, byrow = TRUE
  )
  return(sp::Polygon(poly))
}

#' @title Create a plot sample.
#'
#' @description
#' Creates a plot sample on a regular grid with a random start location.
#'
#' @param spdf A `SpatialPointsDataFrame` defining the points that are to be
#'   sampled by the plot sample.
#' @param boundary A `SpatialPolygonsDataFrame` defining the survey boundary
#'   within which the  points occur.
#' @param x.ppn The proportion of the x=axis that is to be included in the
#'   plots.
#' @param y.ppn The proportion of the y=axis that is to be included in the
#'   plots.
#' @param nx The number of plots in the x-dimension.
#' @param ny The number of plots in the y-dimension.
#'
#' @return A list with three components:
#'  \describe{
#'    \item{`plots`}{A `SpatialPolygonsDataFrame` object containing the plots
#'      that were sampled.}
#'    \item{`dets`}{A `SpatialPointsDataFrame` object containing the locations
#'      of the points within the plots.}
#'    \item{`counts`}{A `dataframe` containing the following columns
#'    \describe{
#'    \item{`x`}{The x-coordinates of the centres of the plots within the
#'      boundary.}
#'    \item{`y`}{The y-coordinates of the centres of the plots within the
#'      boundary.}
#'    \item{`n`}{The numbers of points in each plot.}
#'    \item{`area`}{The areas of the plots within the boundary}
#'    }}
#'  }.
#'
#' @examples
#' \donttest{
#' # Some features require the raster package
#' if (bru_safe_sp() &&
#'   require("sp") &&
#'   require("raster", quietly = TRUE) &&
#'   require("ggplot2", quietly = TRUE) &&
#'   require("terra", quietly = TRUE) &&
#'   require("sf", quietly = TRUE)) {
#'   gorillas <- gorillas_sp()
#'   plotpts <- plotsample(gorillas$nests, gorillas$boundary,
#'     x.ppn = 0.4, y.ppn = 0.4, nx = 5, ny = 5
#'   )
#'   ggplot() +
#'     gg(plotpts$plots) +
#'     gg(plotpts$dets, pch = "+", cex = 2) +
#'     gg(gorillas$boundary)
#' }
#' }
#'
#' @export
#'
plotsample <- function(spdf,
                       boundary,
                       x.ppn = 0.25,
                       y.ppn = 0.25,
                       nx = 5,
                       ny = 5) {
  if (x.ppn <= 0 ||
    x.ppn >= 1) {
    stop("'x.ppn' must greater than 0 and less than 1")
  }
  if (y.ppn <= 0 ||
    y.ppn >= 1) {
    stop("'y.ppn' must greater than 0 and less than 1")
  }

  srange <- raster::extent(boundary)
  xrange <- srange[1:2]
  yrange <- srange[3:4]
  nxtot <- round(nx / x.ppn)
  nytot <- round(ny / y.ppn)
  width <- diff(xrange) / nxtot
  height <- diff(yrange) / nytot
  dx <- diff(xrange) / nx
  dy <- diff(yrange) / ny
  startx <- runif(
    1,
    xrange[1] - 0.99999 * width,
    xrange[1] + dx - 0.99999 * width
  )
  starty <- runif(
    1,
    yrange[1] - 0.99999 * height,
    yrange[1] + dy - 0.99999 * height
  )
  xs <- startx + (0:nx) * dx
  ys <- starty + (0:ny) * dy

  nxs <- length(xs)
  nys <- length(ys)
  starts <- data.frame(x = rep(xs, nys), y = rep(ys, rep(nxs, nys)))
  nplots <- dim(starts)[1]

  polys <- vector("list", nplots)
  for (i in 1:nplots) {
    polys[[i]] <- sp::Polygons(
      list(makepoly(
        as.numeric(starts[i, ]),
        width,
        height
      )),
      i
    )
  }
  plots <- sp::SpatialPolygons(polys, proj4string = fm_CRS(spdf))
  # remove bits of plot outside boundary
  plots <- raster::intersect(boundary, plots)
  # extract only those nests inside the polygons (neat!)
  dets <- spdf[plots, ]

  return(list(plots = plots, dets = dets))
}



#' @title Convert a plot sample of points into one of counts.
#'
#' @description
#' Converts a plot sample with locations of each point within each plot, into a
#' plot sample with only the count within each plot.
#'
#' @param plots A `SpatialPolygonsDataFrame` object containing the plots that
#'   were sampled.
#' @param dets A `SpatialPointsDataFrame` object containing the locations of the
#'   points within the plots.
#'
#' @return A `SpatialPolygonsDataFrame` with counts in each plot contained in
#'   slot `@data$n`.
#'
#' @examples
#' \donttest{
#' # Some features require the raster package
#' if (bru_safe_sp() &&
#'   require("sp") &&
#'   require("raster", quietly = TRUE) &&
#'   require("ggplot2", quietly = TRUE) &&
#'   require("terra", quietly = TRUE) &&
#'   require("sf", quietly = TRUE)) {
#'   gorillas <- gorillas_sp()
#'   plotpts <- plotsample(gorillas$nests, gorillas$boundary,
#'     x.ppn = 0.4, y.ppn = 0.4, nx = 5, ny = 5
#'   )
#'   p1 <- ggplot() +
#'     gg(plotpts$plots) +
#'     gg(plotpts$dets) +
#'     gg(gorillas$boundary)
#'   countdata <- point2count(plotpts$plots, plotpts$dets)
#'   x <- sp::coordinates(countdata)[, 1]
#'   y <- sp::coordinates(countdata)[, 2]
#'   count <- countdata@data$n
#'   p2 <- ggplot() +
#'     gg(gorillas$boundary) +
#'     gg(plotpts$plots) +
#'     geom_text(aes(label = count, x = x, y = y))
#'   multiplot(p1, p2, cols = 2)
#' }
#' }
#'
#' @export
#'
point2count <- function(plots, dets) {
  np <- length(plots)
  x <- y <- plotarea <- count <- numeric(length = np)
  for (i in 1:np) {
    polylist <- list(sp::Polygons(list(plots@polygons[[i]]@Polygons[[1]]), 1))
    spoly <- sp::SpatialPolygons(
      polylist,
      proj4string = sp::CRS(as.character(sp::proj4string(plots)))
    )
    count[i] <- dim(dets[spoly, ])[1]

    plotarea[i] <- plots@polygons[[i]]@area
    xs <- unique(plots@polygons[[i]]@Polygons[[1]]@coords[, 1])
    ys <- unique(plots@polygons[[i]]@Polygons[[1]]@coords[, 2])
    x[i] <- min(xs) + abs(diff(range(xs)) / 2)
    y[i] <- min(ys) + abs(diff(range(ys)) / 2)
  }

  # make a data frame of it
  countdf <- data.frame(n = count, area = plotarea, x = x, y = y)
  # make SpatialPointsDataFrame of it
  plotcounts <- sp::SpatialPointsDataFrame(
    coords = data.frame(x = x, y = y),
    data = data.frame(n = count, area = plotarea),
    proj4string = sp::CRS(as.character(sp::proj4string(plots)))
  )
  return(plotcounts)
}