#' 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.