R/MakePlotCountTable.R

In bmcnellis/RSFIA: FIA Tree Mortality Analysis

#' @title Create summary tables of plot information
#'
#' @description
#'
#' Summary statistics of dataset, used in 1st section of results
#'
#' @param by one of 'section' or 'forest_type'
#' @param csv logical flag, write the tables to a CSV in current wd?
#'
#' @examples
#' tbl1 <- MakePlotCountTable(by = 'section', csv = F)
#' print(tbl1)
#' tbl2 <- MakePlotCountTable(by = 'forest_type', csv = F)
#' print(tbl2)
#' @export
MakePlotCountTable <- function(by = 'section', csv = F) {
  stopifnot(by %in% c('section', 'forest_type'))
  trees <- FIA_mortality_TREE_level
  mort <- FIA_mortality_with_explanatory
  #wth <- FIA_weather_by_year_bin
  #mort <- dplyr::left_join(mort, wth, by = c('LAT', 'LON'))
  ntree <- aggregate(trees$CN, by = list(trees$PLT_CN),
                     FUN = length)
  colnames(ntree) <- c('PLT_CN', 'tree_count')
  tree_df <- dplyr::left_join(mort, ntree, by = 'PLT_CN')
  if (by == 'section') {
    by_var0 <- tree_df$Cleland_section
    #by_var_name0 <- 'section'
    by_var <- mort$Cleland_section
    var_name <- 'section'
    table_var_name <- 'Section'
  } else {
    by_var0 <- tree_df$forest_type
    #by_var_name0 <- 'forest_type'
    by_var <- mort$forest_type
    var_name <- 'forest_type'
    table_var_name <- 'Forest Type'
  }
  tree_counts <- aggregate(tree_df$tree_count, by = list(by_var0),
                           FUN = function(x) sum(x, na.rm = T))
  colnames(tree_counts) <- c(var_name, 'tree_count')

  plt_tbl <- table(by_var)
  plot_counts <- data.frame(matrix(ncol = 0, nrow = length(plt_tbl)), stringsAsFactors = F)
  plot_counts[[var_name]] <- names(plt_tbl)
  plot_counts$plot_count <- as.numeric(plt_tbl)

  out_table <- dplyr::left_join(plot_counts, tree_counts, by = var_name)

  if (by == 'section') {
    provs <- ScaleUpClelandCode(out_table[, 1], in_lvl = 'section', out_lvl = 'province')
    provs <- KeyClelandCode(provs, lvl = 'province')
    out_table <- data.frame(provs, out_table, stringsAsFactors = F)
  }

  # MAT:
  mat <- aggregate(mort$MAT_10, by = list(by_var),
                   FUN = function(x) round(mean(x, na.rm = T), 2))
  colnames(mat) <- c(var_name, 'MAT')
  out_table <- dplyr::left_join(out_table, mat, by = var_name)

  # MAP:
  map <- aggregate(mort$MAP_10, by = list(by_var),
                   FUN = function(x) round(mean(x, na.rm = T), 0))
  colnames(map) <- c(var_name, 'MAP')
  out_table <- dplyr::left_join(out_table, map, by = var_name)

  # ELEV, needs the conversion:
  ele <- aggregate(mort$ELEV, by = list(by_var),
                   FUN = function(x) mean(x, na.rm = T))
  ele[, 2] <- round(ele[, 2] / 3.28084, 0)
  colnames(ele) <- c(var_name, 'Mean Elevation')
  out_table <- dplyr::left_join(out_table, ele, by = var_name)

  # Mean mortality, overall ('mort_rate'):
  mcol <- c('mort_rate', 'non_harv_mort_rate', 'non_harv_fire_mort_rate',
                 'non_harv_fire_beet_mort_rate')
  mcol_long <- c('All Mort', 'Minus Log', 'Minus Log/Fire', 'Minus Log/Fire/Insect')
  coln <- character()
  aa <- 0
  for (i in mcol) {
    aa <- aa + 1
    for (j in c(NA, T, F)) {
      if (is.na(j)) {
        jt <- ''
        mort0 <- mort
        by_var0 <- by_var
      } else if (j) {
        jt <- ', Outlier'
        mort0 <- mort[which(mort$mort_outlier == T), ]
        by_var0 <- by_var[which(mort$mort_outlier == T)]
      } else if (!j) {
        jt <- ', Background'
        mort0 <- mort[which(mort$mort_outlier == F), ]
        by_var0 <- by_var[which(mort$mort_outlier == F)]
      }
      mv <- aggregate(mort0[[i]], by = list(by_var0),
                      FUN = function(x) mean(x, na.rm = T))
      mv[, 2] <- round(mv[, 2], 4)
      i_col <- paste0(mcol_long[aa], jt)
      coln <- append(coln, i_col)
      colnames(mv) <- c(var_name, i_col)
      out_table <- dplyr::left_join(out_table, mv, by = var_name)
    }
  }

  # Put it together:
  out_cols <- c(table_var_name,
                'N plots', 'N trees', 'MAT', 'MAP', 'Mean Elevation',
                coln)
  if (by == 'section') {
    out_cols <- c('Province', out_cols)
    out_table[, 2] <- KeyClelandCode(out_table[, 2], lvl = var_name)
    out_table <- out_table[order(out_table[, 1], out_table[, 2]), ]
  } else {
    out_table <- out_table[order(out_table[, 1]), ]
  }
  colnames(out_table) <- out_cols

  # Output, just requires out_table:
  if (csv) {
    message('Writing CSV')
    fl_name <- 'plot_count_table'
    fl_inc <- 0
    cnt <- 0
    fl_chk <- paste0(fl_name, as.character(fl_inc), '.csv')
    while (fl_chk %in% list.files()) {
      fl_inc <- fl_inc + 1
      fl_chk <- paste0(fl_name, as.character(fl_inc), '.csv')
      if (cnt > 100) stop('whiles fucked')
    }
    fl_name <- paste0(fl_name, fl_inc, '.csv')
    write.csv(out_table, file = fl_name)
  }
  return(out_table)
}