R/DiamFunctions.R

In ryanmismith/inventoryfunctions: Forestry Inventory Functions

#' A variety of plot level diameter metrics for even and uneven aged stands
#'
#'
#' This function defaults to plot level QMD. To change metric enter QMD = FALSE and
#' METRICOFINTEREST = TRUE. Only a single metric may be calculated and output at one time.
#' The function can calculate either QMD, the weighted mean, the
#' mean diameter of largest 100 trees, the mean diameter of each species, or the
#' mean diameter of a single species with all other trees being lumped into a single
#' diameter group - all at the plot level.
#' This lst metric can be useful when looking at relative importance of a single
#' species or for any species specific analysis.
#'
#' All of these functions are meant to provide
#' different tools for plot level diameter metrics.
#'
#'
#'@param Stand Unique Stand ID for each stand in your dataset.
#'Can be set to 1 if only one stand is of interest.
#'@param Plot Unique Plot ID for each Plot. Use Unique.ID function to ensure
#'plot numbers are not recycled between stands.
#'@param Tree Unique Tree ID.
#'@param DBH Diameter at Breast Height (cm)
#'@param EXPF Expansion factor of each tree. See EXP.F function.
#'@param SPP FVS Species code for each tree. Defaults to NA. Is only necessary for
#'Species Specific Diameter and Interest SPP outputs.
#'@param QMD Plot level Quadratic Mean Diameter. Default Output.
#'Defaults to TRUE and must be set as FALSE to use another metric.
#'@param TrueMean True weighted mean diameter for each plot. Defaults to FALSE.
#'@param PlotLargest Mean value of 100 largest diameter trees in plot. Defaults to FALSE.
#'@param AllSpeciesDiam Provides weighted mean of diameter for each species in plot.
#'requires entry of SPP diameter. Defaults to FALSE.
#'@param SpecificSPP Obtain the weighted mean of a specific species' diameters and the
#'weighted mean of all the other species combined. Requires a InterestSPP be named.
#'Defaults to FALSE.
#'@param InterestSPP The SPP that will be singled out for analysis using the
#'SpecificSPP option in this function.
#'
#'@seealso [inventoryfunctions::Unique.ID]
#'@seealso [inventoryfunctions::EXP.F]
#'
#'@family Plot Level Functions
#'
#'@details This function uses the dplyr package to select, arrange, and mutate the data.
#'
#' This function defaults to plot level QMD. To change metric enter QMD = FALSE and
#' METRICOFINTEREST = TRUE.
#'
#' Only a single metric may be calculated and output at one time.
#'
#' Your data should include unique IDs specific to each plot. If you have not done this, it can be done using the
#' Unique.ID function included in this package.
#'
#'@return The function returns a vector of length n with the descriptive metric chosen by the user.
#'@author Ryan Smith and Brian Smith
#'@examples
#'
#'\dontrun{
#'  ###### RUNNING THE SCRIPT ######
#'  ###### REQUIRES FULL VECTORS #####
#'  ###### AS SEEN HERE ######
#'  Stand <- c(1,1,1,1,rep(2,5))
#'  Plot <- c(1,1,1,1,3,3,4,4,4)
#'  Tree <- c(1:9)
#'  DBH <- seq(10, 18, 1)
#'  EXPF <- c(13,22,16,12,7,12,34,7,11)
#'  SPP <- c("WP", "WP", rep("RO", 7))
#'
#'  DBHFuns(Stand, Plot, Tree, DBH, EXPF, SPP) # provides QMD
#'  DBHFuns(Stand, Plot, Tree, DBH, EXPF, SPP, QMD = FALSE, TrueMean = TRUE)
#'  DBHFuns(Stand, Plot, Tree, DBH, EXPF, SPP, QMD = FALSE, PlotLargest = TRUE)
#'  DBHFuns(Stand, Plot, Tree, DBH, EXPF, SPP, QMD = FALSE, AllSpeciesDiam = TRUE)
#'  DBHFuns(Stand, Plot, Tree, DBH, EXPF, SPP, QMD = FALSE, SpecificSPP = TRUE, InterestSPP = "WP")
#'
#' }
#'
#'
#'@export

DBHFuns <- function(Stand, Plot, Tree, DBH, EXPF = NA, SPP = NA,
                          QMD = TRUE, TrueMean = FALSE,
                          PlotLargest = FALSE, AllSpeciesDiam = FALSE,
                          SpecificSPP = FALSE,
                          InterestSPP = NA
                          ){


  if(is.na(EXPF[1]) == TRUE){
    stop("Please Enter Individual Tree's Expansion Factor as a vector of length Tree")
  } else {

   SPP <- SPP
   trees <- tidyr::tibble(Stand, Plot, Tree, DBH, EXPF, SPP)
   trees <- trees %>% dplyr::mutate(
     dplyr::across(
       .cols = c(Stand, Plot, SPP),
       .fns = as.factor
     )
   )
   trees <- trees %>% dplyr::mutate(
     ID = Unique.ID(Stand, Plot),
     ID = as.factor(ID)
   )
  }
   ############################
   ########### QMD ############
   ############################

   if(QMD == TRUE){

     trees <- trees %>%
       dplyr::mutate(                                     # Get Squared Diameters
       diamsq = DBH^2
     )

     trees <- trees %>% dplyr::group_by(Stand, Plot) %>%         # Calculate QMD using Curtis
       dplyr::mutate(
        count = length(unique(Tree)),
        sumdiamsq = sum(diamsq),
        QMD = sqrt(count/sum((1/diamsq)))
        )

     x <- trees$QMD

   }

   ############################
   ######### True Mean ########
   ############################

   else if(TrueMean == TRUE) {                           # Speaks for itself


     trees <- trees %>% dplyr::group_by(Stand, Plot) %>%
       dplyr::mutate(
         average = weighted.mean(DBH, EXPF),
       )

     x <- trees$average

   }

   ############################
   ######100 Largest ##########
   ############################

   else if(PlotLargest == TRUE) {


     Temp <- trees                                       # Make Dataframe
     Temp <- Temp[order(-DBH),]                          # Order it for Cumsum function
     Temp$csum <- ave(Temp$EXPF, Temp$ID, FUN=cumsum)    # Cumsum each ID (Plot or Stand)
     Temp$X <- NA
     Temp$Counts <- NA
     Temp$leftover <- NA
     Temp$Total <- NA
     Temp$Result <- NA
     for (i in 1:length(Temp$csum)){
       if(Temp$csum[i] <= 100){
         Temp$X[i] <- Temp$DBH[i] * Temp$EXPF[i]          # Get a total sum of DBHs for trees where cumsum <= 100
         Temp$Counts[i] <- Temp$csum[i]                   # Identify which trees are included in that sum
       } else {
         Temp$X[i] <- 0                                   # Identify which trees are not included in that sum
         Temp$Counts[i] <- 0
       }
     }

     Temp$MaxCum <- ave(Temp$csum, Temp$ID, FUN = function(x) 100 - max(x))                    # 100-MaxCum == the number of trees in plots where cumsum < 100
     Temp$remainder <- ave(Temp$Counts, Temp$ID, FUN = function(x) 100-max(x, na.rm = TRUE))   # Identify number of trees needed for X to get to 100
     Temp$max <- ave(Temp$Counts, Temp$ID, FUN = function(x) which.max(x))                     # Index for max tree if first tree has EXPF > 100
     Temp$index <- ave(Temp$Counts, Temp$ID, FUN = function(x) which.max(x) + 1)               # Index of the first tree not included in X
     Temp <- Temp %>%                                                                          # Identify the height of the remainder trees
       dplyr::group_by(ID) %>%
       dplyr::arrange(desc(DBH), .by_group = TRUE) %>%
       dplyr::mutate(
         mindbh = DBH[index]
       )
     Temp <- Temp %>%                                                                          # Identify the top tree when expf of first tree > 100
       dplyr::group_by(ID) %>%
       dplyr::arrange(desc(DBH), .by_group = TRUE) %>%
       dplyr::mutate(
         maxdbh = DBH[max]
       )


     for(i in 1:length(Temp$mindbh)){                                      # Get rid of NAs for plots with less than 100 stems
       if(is.na(Temp$mindbh[i]) == TRUE){
         Temp$mindbh[i] <- 0
       } else {
         Temp$mindbh[i] <- Temp$mindbh[i]
       }
     }

     for (i in 1:length(Temp$remainder)){                                    # Create column with total combined DBH of remainder trees
       if(Temp$remainder[i] < 100){
         Temp$leftover[i] <- Temp$remainder[i] * Temp$mindbh[i]
       } else {
         Temp$leftover[i] <- Temp$remainder[i] * Temp$maxdbh[i]                 # Sum of 100*MaxTree if Tallest tree has EXPF > 100
       }
     }

     Temp$Y <- ave(Temp$X, Temp$ID, FUN = function(x) sum(x, na.rm = TRUE))  #Combined DBH of trees in X (cumsum <= 100)

     for (i in 1:length(Temp$Y)){
       Temp$Total[i] <- Temp$leftover[i] + Temp$Y[i]                         # Create column with combined DBH of 100 largest trees
     }

     for (i in 1:length(Temp$Total)){
       if(Temp$MaxCum[i] > 0){
         Temp$Result[i] <- Temp$Total[i]/(100-Temp$MaxCum[i])                # If there are less than 100 trees divide by the # of trees.
       } else {
         Temp$Result[i] <- Temp$Total[i]/100                                 # If there are > 100 trees, combined dbh of 100 tallest trees / 100.
       }
     }


     as.numeric(round(Temp$Result, 2))                                                 # Round return.

     x <- Temp$Result                                                     # Result returned


   }

   ####################################
   ######Average Diam By SPP ##########
   ####################################

   else if (AllSpeciesDiam == TRUE) {
     if(is.na(trees$SPP[1]) == TRUE){
       stop("Error: Please provide Tree's SPP Input in vector of length Tree")
     } else {
     trees <- trees %>% dplyr::group_by(Stand, Plot, SPP) %>%
       dplyr::mutate(
         SPPMean = weighted.mean(x = DBH, w = EXPF)
       )

     x <- trees$SPPMean
     }



   ##########################################################
   ######Compare Specific Species to Rest of Stand ##########
   #########################################################

   } else if(SpecificSPP == TRUE) {


     if (is.na(InterestSPP) == TRUE) {
       stop("Error: Please Provide SPP of interest input (InterestSPP = SPP)")

     }
    if(is.na(trees$SPP[1]) == TRUE) {
      stop("Error: Please Provide SPP input")
  }
    if(is.na(trees$SPP[1]) == TRUE | is.na(trees$SPP[length(trees$SPP) - 1]) == TRUE) {

      stop("Error: Please Provide SPP input")
    } else {

      trees$group <- NA


            for (i in 1:length(trees$SPP)){

      for (i in 1:length(trees$SPP)){

        if(trees$SPP[i] == InterestSPP){
          trees$group[i] <- "group1"
        } else {
          trees$group[i] <- "group2"
        }
      }

      trees$group <- as.factor(trees$group)
      trees <- trees %>% dplyr::group_by(Stand, Plot, group) %>%
        dplyr::mutate(
          average = weighted.mean(DBH, EXPF)
        )

      x <- trees$average

    }
   }


   }
  return(x)
}