R/spacetime.R

In reconhub/epimaps: Mapping Tools for Epidemics

#' Spatio-temporal maps from incidence objects
#'
#' This function implements similar features to the S4 generic \code{stplot}
#' from the \code{spacetime} package for \code{incidence} objects (from the
#' \code{incidence} package), which are stratified by locations.
#'
#' @export
#'
#' @importFrom methods as
#'
#' @author Isobel Blake (\email{isobel.blake@@imperial.ac.uk})
#' and Thibaut Jombart (\email{thibautjombart@@gmail.com})
#'
#' @seealso
#' \itemize{
#'  \item \code{\link[raster]{getData}} to fetch shapefiles freely available from GDAM
#'  \item the incidence package documented at: \url{http://www.repidemicsconsortium.org/incidence/}
#' }
#'
#' @examples
#'
#' \dontrun{
#' if (require("sp") &&
#'     require("raster") &&
#'     require("outbreaks") &&
#'     require("incidence")) {
#'
#' ## fetch shapefiles for different admin levels
#' adm2data <- raster::getData('GADM', country = 'SLE', level = 2)
#' adm1data <- raster::getData('GADM', country = 'SLE', level = 1)
#' adm0data <- raster::getData('GADM', country = 'SLE', level = 0)
#'
#' ## plot data
#' plot(adm2data, border = "grey")
#' plot(adm1data, add = TRUE)
#' plot(adm0data, add = TRUE, lwd = 2)
#' text(sp::coordinates(adm2data), labels = adm2data$NAME_2, cex = 0.6)
#'
#' ## make fake data
#' dat <- ebola_sim$linelist
#'
#' ## randomly assign district location
#' p <- c(0.15,0.1,0.01,0.01,0.02,0,0,0.05,0.01,0,0.05,0.1,0.25,0.25)
#' dat$adm2 <- sample(adm2data$NAME_2, size = nrow(dat),
#'                    replace = TRUE, prob = p)
#'
#' ## make an incidence object
#' x <- incidence::incidence(dat$date_of_onset, 30, groups = dat$adm2)
#'
#' ## make plot
#' res <- spacetime(x, adm2data)
#' res
#'
#' ## extra data stored in attributes
#' attributes(res)
#' }
#' }
#'
#' @return
#' The function returns a plot with extra attributes:
#' \itemize{
#'  \item \code{incience_mat}: an incidence matrix
#'  \item \code{stfdf_obj}: an stfdf object
#' }
#'
#' @param obj An \code{incidence} objects, as generated by the code{incidence}
#'   package, stratified by geographic units matching that of a provided
#'   shapefile.
#'
#' @param sf A shapefile stored as a \code{SpatialPolygonsDataFrame}
#'   object from the \code{sp} package.
#'
#' @param field A character string indicating the field containing the name of
#'   the spatial units to be used for plotting - the same as the ones used for
#'   incidence stratification.
#'
#' @param type A character string indicating the type of graphics to produce: a
#'   map ("map"), or a heatmap ("heatmap").
#'
#' @param breaks A numeric vector indicating break points for the color scale.
#'
#' @param pal A color palette to be used for plotting.
#'
#' @param ... Further arguments to be passed to \code{stplot}.
#'
spacetime <- function(obj, sf, type = c("map", "heatmap"),
                      breaks = NULL, field = "NAME_2",
                      pal = viridis::viridis, ...) {

  type <- match.arg(type)

  if (!inherits(sf, "SpatialPolygonsDataFrame")) {
    msg <- "'sf' must be a 'SpatialPolygonsDataFrame'"
    stop(msg)
  }

  if (!field %in% names(sf)) {
    msg <- paste("column", field, "is not in the sf")
    stop(msg)
  }

  units <- sf[[field]]

  missing_strat <- colnames(obj$counts)[!colnames(obj$counts) %in% units]
  if (length(missing_strat) > 0) {
    msg <- paste("The following stratification units are missing from sf:",
                 paste(missing_strat, collapse = ",")
                 )
    stop(msg)
  }

  ## We need to build a matrix of case counts which contains information for all
  ## polygons of the sf being used. In pratice, this means adding a bunch
  ## of 'zeros' for polygons for which there are no cases reported.

  matrix_inc <- t(obj$counts)
  colnames(matrix_inc) <- seq_len(ncol(matrix_inc))


  ## add districts with no incidence

  extra_units <- units[which(!(units %in% colnames(obj$counts)))]
  matrix_inc_ex <- matrix(0L,
                          nrow = length(extra_units),
                          ncol = ncol(matrix_inc))
  rownames(matrix_inc_ex) <- extra_units
  colnames(matrix_inc_ex) <- colnames(matrix_inc)

  all_matrix <- t(merge(t(matrix_inc),
                        t(matrix_inc_ex),
                        by = "row.names", all = TRUE)[,-1])

  ## order rows by order of sf
  all_matrix <- all_matrix[units, ]


  ## Now that the matrix of incidence has been defined for all spatial units, we
  ## convert the shapefile to a space-time obect (STFDF).

  sf_spatialpoly <- as(sf, "SpatialPolygons")
  xts_obj <- xts::xts(1:length(obj$dates), obj$dates)
  df <- data.frame(spt = c(all_matrix), # same as as.vector ;)
                   row.names = 1:(nrow(all_matrix) * ncol(all_matrix)))

  stfdf_obj <- spacetime::STFDF(sf_spatialpoly,
                                xts_obj,
                                df)


  ## make plots

  if (is.null(breaks)) {
    breaks <- pretty(c(0, max(obj$counts), 5))
    breaks <- c(0, 1, breaks[breaks > 1])
  }

  my_col <- c("white", pal(length(breaks) - 1))

  if (type == "map") {

    out_plot <- spacetime::stplot(stfdf_obj, obj$dates, at = breaks,
                                  col.regions = my_col,
                                  ...)
    # +
    ## latticeExtra::layer(sp.polygons(adm1data,lwd=1.2)) + latticeExtra::layer(sp.polygons(adm0data,lwd=2))

  } else if (type == "heatmap") {
    out_plot <- spacetime::stplot(stfdf_obj, obj$dates, at = breaks,
                                  col.regions = my_col, mode = "xt", ...)
  }


  ## generate output: only the plot, and other potentially useful data are
  ## attached as attributes

  out <- out_plot
  attr(out, "incidence_mat") <- all_matrix
  attr(out, "stfdf_obj") <- stfdf_obj

  return(out)
}