R/geofacet-plot.R

In bbsBayes: Hierarchical Bayesian Analysis of North American BBS Data

#' Generate a geofacet plot of population trajectories by province/state
#'
#' \code{geofacet_plot} allows you to generate a faceted plot of population trajectories
#'   for each strata by province/state. Given a model stratified by "state", "bbs_cws", or "bbs_usgs"
#'   and indices generated by \code{generate_strata_indices} or \code{generate_regional_indices}, this
#'   function will generate a faceted plot showing the population trajectories. All geofacet plots have
#'   one facet per state/province, so if strata-level indices from the "bbs_cws" or "bbs_usgs" are given,
#'   the function plots multiple trajectories (one for each of the relevant strata) within each facet.
#'
#' @param indices_list Dataframe of strata or state/province indices produced by
#'   \code{generate_strata_indices} or \code{generate_regional_indices}
#' @param select logical flag to indicate if the strata_level data need to be selected out of an indices_list object that includes stratum, national, or other region-types. Default is FALSE
#' @param stratify_by How were the data stratified?
#' @param multiple Logical, if TRUE, multiple strata-level trajectories are plotted within each prov/state facet
#' @param trends Optional dataframe of matching strata or state/province trends produced by
#'   \code{generate_strata_trends} or \code{generate_regional_trends}. If included trajectories are coloured based on the same
#'   colour scale used in \code{generate_map}
#' @param slope Optional Logical, if dataframe of trends is included, colours in the plot are based on slope trends, Default = FALSE
#' @param add_observed_means Should the facet plots include points indicating the observed mean counts. Defaults to FALSE.  Note: scale of observed means and annual indices may not match due to imbalanced sampling among strata
#' @param species Species name to be added onto the plot
#' @param ci_width quantile to define the width of the plotted credible interval. Defaults to 0.95, lower = 0.025 and upper = 0.975
#' @param col_viridis Logical flag to use "viridis" colour-blind friendly palette. Default is FALSE
#'
#'
#'
#' @return ggplot object
#'
#' @importFrom geofacet facet_geo
#' @importFrom ggplot2 ggplot theme element_blank element_line
#' labs geom_line geom_ribbon aes element_text element_rect margin
#' @importFrom grDevices grey
#' @importFrom ggrepel geom_text_repel
#'
#' @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  stratum indices
#' indices <- generate_indices(jags_mod = jags_mod,
#'                             jags_data = jags_data,
#'                             regions = c("stratum"))
#'
#' # Now make the geofacet plot.
#' gp <- geofacet_plot(indices_list = indices,
#'                     stratify_by = "bbs_cws",
#'                     species = "Pacific Wren",
#'                     multiple = TRUE)
#'
#' # There is an unfortunate conflict between geofacet function in the geofacet package
#' # and the S3 +.gg method in other ggplot-extension-packages like ggmcmc
#' # The geofacet_plot function may fail with the following error message:
#' #  Error: Don't know how to add e2 to a plot
#' # If this happens, you can fix the problem by following these steps
#' #   1 - save your model output
#' #   2 - restart your R-session
#' #   3 - reload the bbsBayes package (do not re-load the other conflicting package, e.g., ggmcmc)
#'
#' @export
#'



geofacet_plot <- function(indices_list = NULL,
                          select = FALSE,
                          stratify_by = NULL,
                          ci_width = 0.95,
                          multiple = FALSE,
                          trends = NULL,
                          slope = FALSE,
                          add_observed_means = FALSE,
                          species = "",
                          col_viridis = FALSE)
{
  # Annoying thing to get rid of check notes
  Year <- NULL; rm(Year)
  Index <- NULL; rm(Index)
  Region <- NULL; rm(Region)
  lci <- NULL; rm(lci)
  uci <- NULL; rm(uci)
  Trendcat <- NULL; rm(Trendcat)
  lbl <- NULL; rm(lbl)
  obs_mean <- NULL; rm(obs_mean)

  if (is.null(stratify_by))
  {
    stop("Argument stratify_by is empty."); return(NULL)
  }

  if (is.null(indices_list))
  {
    stop("Argument indices_list is empty. Please supply indices from generate_strata_indices() or generate_regional_indices()"); return(NULL)
  }

  alpha_ribbon = 0.5

  facets <- utils::read.csv(system.file("geofacet-grids", strata[[stratify_by]], package = "bbsBayes"),stringsAsFactors = FALSE)

  indices = indices_list$data_summary
  if(select){
    if(multiple){
    indices <- indices[which(indices$Region_type == "stratum"),]
    }else{
      indices <- indices[which(indices$Region_type == "prov_state"),]
    }
    }

  lq = (1-ci_width)/2
  uq = ci_width+lq
  lqc = paste0("Index_q_",lq)
  uqc = paste0("Index_q_",uq)

  if (any(grepl(pattern = lqc,x = names(indices))) == FALSE | any(grepl(pattern = uqc,x = names(indices))) == FALSE )
  {
    stop("Desired quantiles are not included in the indices_list object. Re-run generate_x_indices() funtion with desired quantiles."); return(NULL)
  }

  indices$lci = indices[,lqc]
  indices$uci = indices[,uqc]


  mny = min(indices$Year)
  mxy = max(indices$Year)
  yys = as.integer(seq(mny,mxy,length.out = 3))


  uplim = max(indices$Index)

  if(multiple){

     #
     # indices$code = indices$Region
     #

    region_names <- utils::read.csv(system.file("composite-regions", strata[[stratify_by]], package = "bbsBayes"),stringsAsFactors = FALSE)

    indices = merge(indices,region_names[,c("prov_state","region","bcr")],by.x = "Region",by.y = "region")
    indices$code = indices$prov_state

    map_palette <- c("#39568c")

    if(!is.null(trends)){

      trends = merge(trends,region_names[,c("prov_state","region","bcr")],by.x = "Region",by.y = "region")
      trends$code = trends$prov_state

      breaks <- c(-7, -4, -2, -1, -0.5, 0.5, 1, 2, 4, 7)

      if(slope){
        trends$Trend <- trends$Slope_Trend
      }
      trends$Trendcat <- cut(trends$Trend, breaks = c(-Inf, breaks, Inf),ordered_result = TRUE)

      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) <- levels(trends$Trendcat)
      indices <- merge(indices,trends[,c("Region","Trend","Trendcat")],by = "Region")
      indices = indices[order(indices$Region,indices$Year),]

      trlabs = indices[which(indices$Year == yys[2]),]
      trlabs$lbl = paste0(signif(round(trlabs$Trend,1),2)," BCR",trlabs$bcr)


    }else{

      indices$Trendcat = factor(rep("#313695",nrow(indices),levels = map_palette))
      trlabs = indices[which(indices$Year == yys[2]),]
      trlabs$lbl = paste0("BCR",trlabs$bcr)

       }


    ptraj <- ggplot2::ggplot(data = indices) +
      ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
                     panel.grid.minor = ggplot2::element_blank(),
                     panel.background = ggplot2::element_blank(),
                     axis.text.x = ggplot2::element_text(colour = grey(0.2),size = 5,angle = 90),
                     axis.text.y = ggplot2::element_text(colour = grey(0.2),size = 5),
                     strip.background = ggplot2::element_rect(fill = grDevices::grey(0.97)),#strcol #, colour = grey(0.9), size = NULL, linetype = NULL, color = NULL, inherit.blank = FALSE
                     strip.text = ggplot2::element_text(size = 6,margin = ggplot2::margin()),#
                     legend.position = "none") +
      ggplot2::labs(title = paste(species,"trajectories within Provinces and States"), x = "", y = "Annual indices") +

      ggplot2::geom_line(data = indices, ggplot2::aes(x = Year, y = Index,group = Region),colour = grDevices::grey(0.6)) +
      #        ggplot2::geom_line(data = indices, ggplot2::aes(x = Year, y = Index,colour = Trendcat)) +
      ggplot2::geom_ribbon(data = indices, ggplot2::aes(x = Year, ymin = lci, ymax = uci,fill = Trendcat,group = Region), alpha = alpha_ribbon)+
      ggplot2::scale_x_continuous(breaks = yys)+
      ggplot2::coord_cartesian(ylim = c(0,uplim))+
      ggplot2::scale_colour_manual(values = map_palette, aesthetics = c("colour","fill"))+
      ggrepel::geom_text_repel(data = trlabs, mapping = ggplot2::aes(x = Year,y = uci,label = lbl,group = Region),colour = grDevices::grey(0.6), size = 2,nudge_y = 0.2*uplim,segment.alpha = 0.1)


    if(add_observed_means){
      ptraj <- ptraj+ggplot2::geom_point(data = indices,ggplot2::aes(x = Year,y = obs_mean,group = Region),colour = grDevices::grey(0.6),size = 0.5, alpha = alpha_ribbon)
    }


   outplot <- suppressMessages(ptraj+geofacet::facet_geo(facets = ~ code,grid = facets,label = "code"))
   # messages above suppressed so user does not receive invitation to submit the geofacet





  }else{

    indices$code = indices$Region

    map_palette <- c("#313695")

    if(!is.null(trends)){

      trends$code = trends$Region
      breaks <- c(-7, -4, -2, -1, -0.5, 0.5, 1, 2, 4, 7)

      if(slope){
        trends$Trend <- trends$Slope_Trend
      }
      trends$Trendcat <- cut(trends$Trend, breaks = c(-Inf, breaks, Inf),ordered_result = TRUE)

      map_palette <- c("#a50026", "#d73027", "#f46d43", "#fdae61", "#fee090", "#ffffbf",
                       "#e0f3f8", "#abd9e9", "#74add1", "#4575b4", "#313695")
      names(map_palette) <- levels(trends$Trendcat)
      indices <- merge(indices,trends[,c("Region","Trend","Trendcat")],by = "Region")
      trlabs = indices[which(indices$Year == yys[2]),]
      trlabs$lbl = paste(signif(round(trlabs$Trend,1),2),"%/yr")


       }else{

      indices$Trendcat = factor(rep("#313695",nrow(indices),levels = map_palette))

      }


      ptraj <- ggplot2::ggplot(data = indices) +
        ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
                       panel.grid.minor = ggplot2::element_blank(),
                       panel.background = ggplot2::element_blank(),
                       axis.text.x = ggplot2::element_text(colour = grey(0.2),size = 5,angle = 90),
                       axis.text.y = ggplot2::element_text(colour = grey(0.2),size = 5),
                       strip.background = ggplot2::element_rect(fill = grDevices::grey(0.97)),#strcol #, colour = grey(0.9), size = NULL, linetype = NULL, color = NULL, inherit.blank = FALSE
                       strip.text = ggplot2::element_text(size = 6, margin = ggplot2::margin()),#
                       legend.position = "none") +
        ggplot2::labs(title = paste(species,"trajectories within Provinces and States"), x = "", y = "Annual indices") +
        ggplot2::geom_line(data = indices, ggplot2::aes(x = Year, y = Index),colour = grDevices::grey(0.6)) +
#        ggplot2::geom_line(data = indices, ggplot2::aes(x = Year, y = Index,colour = Trendcat)) +
        ggplot2::geom_ribbon(data = indices, ggplot2::aes(x = Year, ymin = lci, ymax = uci,fill = Trendcat), alpha = alpha_ribbon)+
        ggplot2::scale_x_continuous(breaks = yys)+
        ggplot2::coord_cartesian(ylim = c(0,uplim))+
        ggplot2::scale_colour_manual(values = map_palette, aesthetics = c("colour","fill"))


       if(add_observed_means){
        ptraj <- ptraj+ggplot2::geom_point(data = indices,ggplot2::aes(x = Year,y = obs_mean),colour = grDevices::grey(0.6),size = 0.5, alpha = alpha_ribbon)
      }

      if(!is.null(trends)){
        ptraj <- ptraj+ggrepel::geom_text_repel(data = trlabs, mapping = ggplot2::aes(x = Year,y = uci,label = lbl),colour = grDevices::grey(0.6), size = 2,nudge_y = 0.2*uplim,segment.alpha = 0.1)
      }

      outplot <- suppressMessages(ptraj+geofacet::facet_geo(facets = ~ code,grid = facets,label = "code"))
      # messages above suppressed so user does not receive invitation to submit the geofacet




  }

     return(outplot)

}#end function