RSFIA: FIA Tree Mortality Analysis

#' @title Summarizes FIA_mortality_with_explanatory
#'
#' @description Provides summary statistics for the FIA mortality
#'
#' @param ftype_min Minimim percentage of forest type to display with print
#'
#' @examples
#' ForestSummary(ftype_min = 4)
#' @export
ForestSummary <- function(ftype_min = 3) {
  # To use with Rmarkdown: make the output a list
  plot_df <- FIA_mortality_with_explanatory
  tree_df <- FIA_mortality_TREE_level
  sect_tbl <- MakePlotCountTable(by = 'section', csv = F)
  sect_tbl <- sect_tbl[order(sect_tbl$`N plots`, decreasing = T), ]
  prov_counts <- aggregate(sect_tbl$`N plots`, by = list(sect_tbl$Province), FUN = sum)
  prov_counts <- prov_counts[order(prov_counts$x, decreasing = T), ]
  state_tbl <- table(KeyStateCode(plot_df$STATECD, 7.2))
  plots_missing_ftype <- table(is.na(plot_df$forest_type))['TRUE']
  swoods1 <- table(KeyForestGroup(plot_df$forest_type))
  swoods2 <- table(KeyForestGroup(unique(plot_df$forest_type)))
  for_tbl <- table(plot_df$forest_type)
  plt_cnt_tbl <- MakePlotCountTable(by = 'section', csv = F)

  cat('\n\n')
  message('Summary statistics:')
  cat('Total plots:', sum(sect_tbl$`N plots`))
  cat('\nN Provinces/Sections:',
      length(unique(sect_tbl$Province)),
      '/',
      length(unique(sect_tbl$Section)))
  ORWACA <- round(state_tbl[c('OR', 'WA', 'CA')] / sum(state_tbl) * 100, 2)
  cat('\nOR/WA/CA %:', ORWACA[1], ',', ORWACA[2], ',', ORWACA[3])
  cat('\nOR/WA/CA total %:', sum(ORWACA))
  mpn <- round(prov_counts$x[1] / sum(prov_counts$x) * 100, 2)
  mpp <- paste0('(', mpn, '%)')
  cat('\n1-Most plots in province:', prov_counts$x[1], 'in', prov_counts$Group.1[1], mpp)
  sw <- 1
  if (mpn < 50) {
    sw <- sw + 1
    mpn2 <- round(prov_counts$x[2] / sum(prov_counts$x) * 100, 2)
    mpp2 <- paste0('(', mpn2, '%)')
    cat('\n2-Most plots in province:', prov_counts$x[2], 'in', prov_counts$Group.1[2], mpp2)
  }
  if ((mpn + mpn2) < 75) {
    sw <- sw + 1
    mpn3 <- round(prov_counts$x[3] / sum(prov_counts$x) * 100, 2)
    mpp3 <- paste0('(', mpn3, '%)')
    cat('\n3-Most plots in province:', prov_counts$x[3], 'in', prov_counts$Group.1[3], mpp3)
  }
  sw0 <- switch(sw, mpn, mpn + mpn2, mpn + mpn2 + mpn3)
  sw0 <- paste0(sw0, ' %')
  cat('\nSum-% of', sw, '-most plots in province:', sw0)
  lpp <- paste0('(',
                round(prov_counts$x[nrow(prov_counts)] / sum(prov_counts$x) * 100, 2),
                '%)')
  cat('\nLeast plots in province:',
      prov_counts$x[nrow(prov_counts)],
      'in',
      prov_counts$Group.1[nrow(prov_counts)],
      lpp)
  cat('\nMean/min/max N plots per section:',
      'Mean @', round(mean(sect_tbl$`N plots`)), '\nmin-max:',
      round(min(sect_tbl$`N plots`)), '@', sect_tbl$Section[which.min(sect_tbl$`N plots`)],
      round(max(sect_tbl$`N plots`)), '@', sect_tbl$Section[which.max(sect_tbl$`N plots`)]
  )
  cat('\nTotal forest types:', length(unique(plot_df$forest_type)))
  cat('\nTotal N species:', length(unique(tree_df$SPCD)))
  ncon <- round(swoods2['conifer'] / sum(swoods2) * 100, 2)
  cat('\nN Conifer types:', swoods2['conifer'])
  cat('\n% Conifer types:', ncon, '%')
  con_plot_perc <- round(swoods1['conifer'] / sum(swoods1) * 100, 2)
  cat('\n% Conifer plots:', con_plot_perc)
  cat('\nN plots with no forest type:', plots_missing_ftype)
  perc_missing_ftype <- round(plots_missing_ftype / nrow(plot_df) * 100, 2)
  cat('\n% plots with no forest type:', perc_missing_ftype)

  cat('\n\n')
  message('Plot counts by state:')
  for (i in names(state_tbl)) {
    if (i == names(state_tbl)[1]) {
      aa <- ''
    } else {
      aa <- '\n'
    }
    pst <- paste0('(', round(as.numeric(state_tbl[i]) / sum(state_tbl) * 100, 1), '%)')
    cat(aa, i, ':', as.numeric(state_tbl[i]), pst)
  }

  cat('\n\n')
  message('Plot counts by forest type:')
  p1f <- numeric()
  for (i in names(for_tbl)) {
    if (i == names(for_tbl)[1]) {
      bb <- ''
    } else {
      bb <- '\n'
    }
    pftp <- round(as.numeric(for_tbl[i]) / sum(for_tbl) * 100, 1)
    if (pftp > ftype_min) {
      p1f <- append(p1f, pftp)
      names(p1f)[length(p1f)] <- i
    }
    pft <- paste0('(', pftp, '%)')
    cat(bb, i, ':', as.numeric(for_tbl[i]), pft)
  }
  cat('\n\n')
  message('Most abundant forest types:')
  message('(cutoff =' , ftype_min, '%)')
  p1f <- as.data.frame(p1f)
  colnames(p1f) <- c('')
  print(p1f)

  ftype_dom <- for_tbl[which(names(for_tbl) %in% row.names(p1f))]
  dom_perc <- round(sum(ftype_dom) / sum(for_tbl) * 100, 2)
  cat('\nPercent of all plots represented by dominant ftypes:', dom_perc)
  uncom_names <- names(for_tbl)[-which(names(for_tbl) %in% names(ftype_dom))]
  cat('\nNumber of ftypes in uncommon:', length(uncom_names))
  mean_dom_mort <- mean(plot_df[which(plot_df$forest_type %in% names(ftype_dom)), 'mort_rate'])
  mean_dom_mort <- round(mean_dom_mort * 100, 2)
  mdb_df <- plot_df[which(plot_df$mort_outlier == F), c('mort_rate', 'forest_type')]
  mean_dom_back <- mean(mdb_df[which(mdb_df$forest_type %in% names(ftype_dom)), 'mort_rate'])
  mean_dom_back <- round(mean_dom_back * 100, 2)
  mean_sub_back <- mean(mdb_df[which(mdb_df$forest_type %in% uncom_names), 'mort_rate'])
  mean_sub_back <- round(mean_sub_back * 100, 2)
  mdo_df <- plot_df[which(plot_df$mort_outlier == T), c('mort_rate', 'forest_type')]
  mean_dom_out <- mean(mdo_df[which(mdo_df$forest_type %in% names(ftype_dom)), 'mort_rate'])
  mean_dom_out <- round(mean_dom_out * 100, 2)
  mean_sub_out <- mean(mdo_df[which(mdo_df$forest_type %in% uncom_names), 'mort_rate'])
  mean_sub_out <- round(mean_sub_out * 100, 2)
  cat('\nMean mort rate, dominant ftypes, full mort:', mean_dom_mort)
  cat('\nMean mort rate, dominant ftypes, background:', mean_dom_back)
  cat('\nMean mort rate, dominant ftypes, outlier:', mean_dom_out)
  cat('\nMean mort rate, subdom ftypes, background:', mean_sub_back)
  cat('\nMean mort rate, subdom ftypes, outlier:', mean_sub_out)


  mean_sub_mort <- mean(plot_df[which(plot_df$forest_type %in% uncom_names), 'mort_rate'])
  mean_sub_mort <- round(mean_sub_mort * 100, 2)
  cat('\nMean mort rate, subdom ftypes:', mean_sub_mort)


  cat('\n')
  message('Mortality rates:')
  for (i in c('All Mort', 'Minus Log', 'Minus Log/Fire', 'Minus Log/Fire/Insect')) {
    cat(i, 'by section:\n',
        'mean:', round(mean(sect_tbl[[i]]) * 100, 2), '\nmin/max:',
        min(sect_tbl[[i]]) * 100, '@', sect_tbl$Section[which.min(sect_tbl[[i]])],
        max(sect_tbl[[i]]) * 100, '@', sect_tbl$Section[which.max(sect_tbl[[i]])],
        '\n'
    )
  }
  mean_out <- round(mean(plot_df$mort_rate[plot_df$mort_outlier]) * 100, 2)
  mean_bkg <- round(mean(plot_df$mort_rate[!plot_df$mort_outlier]) * 100, 2)
  cat('\nMean overall background/outlier:\n', mean_bkg, '/', mean_out)
  for (i in c('All Mort, Background', 'All Mort, Outlier')) {
    tg <- tolower(strsplit(i, ',')[[1]][2])
    il <- round(plt_cnt_tbl[which.min(plt_cnt_tbl[[i]]), i] * 100, 2)
    is <- plt_cnt_tbl[which.min(plt_cnt_tbl[[i]]), 'Section']
    ip <- plt_cnt_tbl[which.min(plt_cnt_tbl[[i]]), 'Province']
    cat('\n', paste0('Lowest', tg), 'rate/sect/prov:\n',
        paste(il, is, ip, sep = ' / '))
    hl <- round(plt_cnt_tbl[which.max(plt_cnt_tbl[[i]]), i] * 100, 2)
    hs <- plt_cnt_tbl[which.max(plt_cnt_tbl[[i]]), 'Section']
    hp <- plt_cnt_tbl[which.max(plt_cnt_tbl[[i]]), 'Province']
    cat('\n', paste0('Highest', tg), 'rate/sect/prov:\n',
        paste(hl, hs, hp, sep = ' / '))
  }

  cat('\nMean mortality rate by province:\n\n')
  pagg <- aggregate(plot_df$mort_rate, by = list(plot_df$province), FUN = mean)
  pagg[, 2] <- round(pagg[, 2] * 100, 2)
  colnames(pagg) <- c('Province', 'Overall Mortality')
  print(pagg)

  browser()
  cat('\n')
  invisible()
}