R/CountTrees.R

In bmcnellis/RSFIA: FIA Tree Mortality Analysis

#' @title Creates tree counts by species for a given plot
#' @description
#'
#' @export
CountTrees <- function(PLT_CNs, trees, cum_BA = 1, drop_unknown = T,
                       use_metric = T, yr_inc = 10) {
  message('Counting trees...')
  require(RSFIA)
  data('FIA_table_metadata')
  stopifnot(
    'PLT_CN' %in% colnames(trees),
    all(cum_BA > 0, cum_BA <= 1),
    length(unique(PLT_CNs)) == length(PLT_CNs)
  )
  SP_full <- FIA_table_metadata$APPENDIX_F
  SP_names <- SP_full$Scientific.Name[match(trees$SPCD, SP_full$SPCD)]
  if (use_metric) {
    trees$DIA <- round(trees$DIA * 2.54, 2)
    trees$PREVDIA <- round(trees$PREVDIA * 2.54, 2)
  }
  BA <- round(pi * (trees$DIA / 2) ^ 2, 1)
  trees <- data.frame(trees, SP_names, BA, stringsAsFactors = F)
  if (drop_unknown) {
    drop_ids <- c('Tree unknown', 'Tree evergreen', 'Tree broadleaf')
    if (any(drop_ids %in% trees$SP_names)) {
      drop_rows <- which(trees$SP_names %in% drop_ids)
      trees <- trees[-drop_rows, ]
      cat('\nDropped', length(drop_rows), 'unidentified trees from dataset.\n\n')
    }
  }
  i_cnt <- 0
  out_df <- data.frame(matrix(ncol = 1, nrow = length(PLT_CNs)))
  out_df[, 1] <- PLT_CNs
  colnames(out_df) <- c('PLT_CN')
  for (i in PLT_CNs) {
    i_cnt <- i_cnt + 1
    tree_sub <- trees[which(trees$PLT_CN == i), ]
    n_subp <- length(table(tree_sub$SUBP))
    if (cum_BA < 1) {
      # CountTrees specific inputs:
      nmax <- 10000
      in_vec <- tree_sub$BA
      top <- cum_BA
      # Generalized:
      cut <- sum(in_vec, na.rm = T) * top
      c <- numeric()
      j <- 0
      w_vec <- in_vec
      while (sum(in_vec[c], na.rm = T) < cut) {
        j <- j + 1
        m <- which.max(w_vec)
        c <- append(c, m)
        w_vec <- w_vec[-m]
        if (length(w_vec) < 1) stop()
        if (length(c) >= length(in_vec)) stop()
        if (j > nmax) stop('While loop overrun')
      }
      tree_sub <- tree_sub[c[-length(c)], ]
    }
    for (j in tree_sub$SP_names) {
      j_sub <- tree_sub[which(tree_sub$SP_names == j), ]
      j_tag <- paste(unlist(strsplit(j, ' ')), collapse = '_')

      n_trees <- as.numeric(table(tree_sub$SP_names)[j])
      n_tag <- paste0('n_', j_tag)
      out_df[which(PLT_CNs == i), n_tag] <- n_trees

      p_trees <- n_trees / nrow(tree_sub)
      p_tag <- paste0('percent_', j_tag)
      out_df[which(PLT_CNs == i), p_tag] <- p_trees

      d_trees <- n_trees / RSFIA::SubsToAcres(n_subp, r = 24, hectare = use_metric)
      a_tag <- paste0(ifelse(use_metric, 'hectares', 'acres'), '_')
      d_tag <- paste0('density_', a_tag, j_tag)
      out_df[which(PLT_CNs == i), d_tag] <- round(d_trees, 3)

      sp_BA <- mean(j_sub$BA, na.rm = T)
      BA_tag <- paste0('mean_BA_', j_tag)
      out_df[which(PLT_CNs == i), BA_tag] <- round(sp_BA, 2)

      sp_dia <- mean(j_sub$DIA, na.rm = T)
      dia_tag <- paste0('mean_diameter_', j_tag)
      out_df[which(PLT_CNs == i), dia_tag] <- round(sp_dia, 2)

      sp_ht <- mean(j_sub$HT, na.rm = T)
      ht_tag <- paste0('mean_height_', j_tag)
      out_df[which(PLT_CNs == i), ht_tag] <- round(sp_ht, 2)

      sp_age <- mean(j_sub$TOTAGE, na.rm = T)
      sp_age <- ifelse(is.nan(sp_age), NA, sp_age)
      age_tag <- paste0('mean_age_', j_tag)
      out_df[which(PLT_CNs == i), age_tag] <- round(sp_age, 2)

      sp_inc <- round(mean(j_sub$DIA - j_sub$PREVDIA, na.rm = T), 2)
      inc_tag <- paste0('mean_growth_inc_', j_tag)
      out_df[which(PLT_CNs == i), inc_tag] <- sp_inc

      # RGR:
      ln_sp_inc <- mean(log(j_sub$DIA) - log(j_sub$PREVDIA), na.rm = T)
      sp_RGR <- round(ln_sp_inc / yr_inc, 4)
      rgr_tag <- paste0('mean_RGR_', j_tag)
      out_df[which(PLT_CNs == i), rgr_tag] <- sp_RGR

      if (nrow(out_df) > length(unique(PLT_CNs))) stop()
    } # end j
    bmcUtils::LoopStatus(i_cnt, length(PLT_CNs), digits = 1)
  } # end i
  return(out_df)
}