R/generate-map.R

In BrandonEdwards/bbsBayes: Hierarchical Bayesian Analysis of North American BBS Data

#' Generate a map of trends by strata.
#'
#' \code{generate_map} allows you to generate a colour-coded map of species trends
#'   for each strata. Given trends generated by \code{generate_strata_trends}, this
#'   function will shade in each stratum based on the percent
#'   change in that stratum.
#'
#' @param trend Dataframe of strata trends produced by
#'   \code{generate_strata_trends} or \code{generate_regional_trends(..., regions = "stratum")}
#' @param select logical flag to indicate if the stratum data need to be selected out of an trends object that includes continental, national, or other region-types. Default is FALSE
#' @param stratify_by How were the data stratified?
#' @param slope Logical, if TRUE, maps values of the alternative trend metric if slope = TRUE was used in \code{generate_strata_trends}, the slope of a log-linear regression through the annual indices. Default FALSE.
#' @param species Text, optional species name to add plot title. if left blank "" no title is added
#' @param col_viridis Logical flag to use "viridis" colour-blind friendly palette. Default is FALSE
#'
#' @return spplot object
#'
#' @importFrom sf read_sf
#' @importFrom dplyr left_join
#' @importFrom ggplot2 geom_sf
#'
#' @examples
#' # Toy example with Pacific Wren sample data
#' # First, stratify the sample data
#' strat_data <- stratify(by = "bbs_cws", sample_data = TRUE)
#'
#' # Prepare the stratified data for use in a JAGS model.
#' jags_data <- prepare_jags_data(strat_data = strat_data,
#'                                species_to_run = "Pacific Wren",
#'                                model = "firstdiff",
#'                                min_year = 2009,
#'                                max_year = 2018)
#'
#' # Now run a JAGS model.
#' jags_mod <- run_model(jags_data = jags_data,
#'                       n_adapt = 0,
#'                       n_burnin = 0,
#'                       n_iter = 10,
#'                       n_thin = 1)
#'
#' # Generate the continental and stratum indices
#' indices <- generate_indices(jags_mod = jags_mod,
#'                             jags_data = jags_data)
#'
#' # Now, generate the trends
#' trends <- generate_trends(indices = indices)
#'
#' # Generate the map. Set select = TRUE because we are passing a
#' #    dataframe of trends of more than just the stratum regions
#' map <- generate_map(trend = trends,
#'                     stratify_by = "bbs_cws",
#'                     select = TRUE,
#'                     species = "Pacific Wren")
#' @export
#'

generate_map <- function(trend = NULL,
                         select = FALSE,
                         stratify_by = NULL,
                         slope = FALSE,
                         species = "",
                         col_viridis = FALSE)
{
  # Silly things to remove "visible binding" note in check
  Trend <- NULL
  rm(Trend)
  Tplot <- NULL
  rm(Tplot)

  if(select){
    trend = trend[which(trend$Region_type == "stratum"),]
  }
  if (is.null(stratify_by))
  {
    stop("No stratification specified."); return(NULL)
  }

  if(isFALSE(is.element(stratify_by, c("state", "bcr", "latlong", "bbs_cws", "bbs_usgs"))))
  {
    stop("Invalid stratification specified, choose one of state, bcr, latlong, bbs_cws, or bbs_usgs"); return(NULL)
  }

  if (is.null(trend))
  {
    stop("Argument trend is empty, or stratification of model does not match stratify_by argument"); return(NULL)
  }


  fyr = min(trend$Start_year)
  lyr = min(trend$End_year)

  map <- sf::read_sf(dsn = system.file("maps",
                                       package = "bbsBayes"),
                     layer = maps[[stratify_by]],
                     quiet = TRUE)


  breaks <- c(-7, -4, -2, -1, -0.5, 0.5, 1, 2, 4, 7)
  labls = c(paste0("< ",breaks[1]),paste0(breaks[-c(length(breaks))],":", breaks[-c(1)]),paste0("> ",breaks[length(breaks)]))
  labls = paste0(labls, " %")

  if(slope){
    trend$Tplot <- as.numeric(as.character(trend$Slope_Trend))
  }else{
    trend$Tplot <- as.numeric(as.character(trend$Trend))
  }
  trend$Tplot <- cut(trend$Tplot,breaks = c(-Inf, breaks, Inf),labels = labls)

  trend$ST_12 = trend$Region
  map = dplyr::left_join(x = map,y = trend,by = "ST_12")

  if(species != ""){
    ptit = paste(species,"trends",fyr,"-",lyr)
  }else{
    ptit = ""
  }

  if (col_viridis)
  {
    map_palette <- c("#fde725", "#dce319", "#b8de29", "#95d840", "#73d055", "#55c667",
                     "#238a8d", "#2d708e", "#39568c", "#453781", "#481567")
  }else
  {
    map_palette <- c("#a50026", "#d73027", "#f46d43", "#fdae61", "#fee090", "#ffffbf",
                     "#e0f3f8", "#abd9e9", "#74add1", "#4575b4", "#313695")
  }

  names(map_palette) <- labls

  mp.plot = ggplot2::ggplot()+
    ggplot2::geom_sf(data = map,ggplot2::aes(fill = Tplot),colour = grey(0.4),size = 0.1)+
    ggplot2::theme_minimal()+
    ggplot2::ylab("")+
    ggplot2::xlab("")+
    ggplot2::labs(title = ptit)+
    ggplot2::theme(legend.position = "right", line = ggplot2::element_line(size = 0.4),
                   rect = ggplot2::element_rect(size = 0.1),
          axis.text = ggplot2::element_blank(),
          axis.line = ggplot2::element_blank())+
    ggplot2::scale_colour_manual(values = map_palette, aesthetics = c("fill"),
                        guide = ggplot2::guide_legend(reverse=TRUE),
                        name = paste0("Trend\n",fyr,"-",lyr))


  return(mp.plot)
}