R/PlotFullMortMap.R

In bmcnellis/RSFIA: FIA Tree Mortality Analysis

#' @title Plots all mortality in a map format
#'
#' @description Used for Figure 1 in results
#'
#' @examples
#' PlotFullMort()
#' @export
PlotFullMortMap <- function() {
  require(ggplot2)
  # Sanity checks:
  if (!('spData' %in% installed.packages()[, 1])) {
    stop('Need suggested package spData to plot the mortality map.')
  }

  # Obj prep for both graphs:
  mort_df <- PrepPlotEnvir(clear = F)
  full_mort <- mort_df[, c('LON', 'LAT', 'mort_rate')]
  state_overlay <- spData::us_states[which(spData::us_states$REGION == 'West'), ]

  # Obj prep for point graph:
  amin <- 0.1
  highcol <- 'black'
  lowcol.hex <- as.hexmode(round(col2rgb(highcol) * amin + 255 * (1 - amin)))
  lowcol <- paste0("#",   sep = "",
                   paste(format(lowcol.hex, width = 2), collapse = ""))
  highcol <- '#000000'

  # Obj prep for section graph:
  sect_mort <- MakePlotCountTable(by = 'section', csv = F)
  sect_mort <- sect_mort[, c('Section', 'All Mort')]
  sects <- Cleland2007_eco_map[[2]]
  sects$MAP_UNIT_S <- KeyClelandCode(sects$MAP_UNIT_S, lvl = 'section')
  sects <- sects[which(sects$MAP_UNIT_S %in% sect_mort$Section), ]
  vals <- sect_mort$`All Mort`[match(sects$MAP_UNIT_S, sect_mort$Section)]
  a <- 0
  while(any(duplicated(vals))) {
    a <- a + 1
    vals[which(duplicated(vals))] <- vals[which(duplicated(vals))] + 1e-9
    if (a > 100) stop('broke shit')
  }
  row.names(sects) <- as.character(vals)
  sects <- ggplot2::fortify(sects)
  sects$id <- as.numeric(sects$id)
  sects$id <- round(sects$id, 2)

  # Make graphs:
  out_plot0 <- ggplot(data = full_mort) + theme_bw() +
    # State outline:
    geom_sf(data = state_overlay, fill = 'white') +

    # Points, need 2x for the legend to work:
    geom_point(aes(x = LON, y = LAT, alpha = mort_rate), size = 1.0, shape = 16) +
    geom_point(aes(x = LON, y = LAT, colour = mort_rate), alpha = 0) +
    # Legend:
    scale_colour_continuous(high = highcol, low = lowcol) +
    guides(alpha = F) +
    labs(colour = '') +
    xlab('') + ylab('') +
    annotate('text', x = -123, y = 32.5, label = '(a)', size = 12)

  out_plot1 <- ggplot() +
    geom_sf(aes_(), data = state_overlay) +
    geom_polygon(data = sects, mapping = aes(x = long, y = lat, group = group, fill = id),
                 colour = 'black', size = 1) +
    scale_fill_gradient2(low = lowcol, mid = 'grey30', high = highcol, midpoint = 0.15) +
    #scale_fill_gradient(low = lowcol, high = highcol) +
    guides(fill = guide_legend('')) +
    xlab('') + ylab('') +
    annotate('text', x = -123, y = 32.5, label = '(b)', size = 12)

  # Two-panel plot:
  Multiplot(out_plot0, out_plot1)
  invisible(list(out_plot0, out_plot1))
}