RSFIA: FIA Tree Mortality Analysis

#' @title RSFIA Utility functions
#'
#' @description
#'
#' Function families are:
#'
#' generic_utilities
#' FIA_processing
#' plot_wrappers
#' package_utilities
#' soil_data_utilities
#' daymet_utilities
#'
#' @name RSFIA_utils
NULL
#' @describeIn RSFIA_utils Prep function for other plot functions
#' @family package_utilities
#' @export
PrepPlotEnvir <- function(ret = 'plots', clear = T) {
  require(RSFIA)
  require(ClelandEcoregions)
  require(sf)
  require(ggplot2)
  require(grid)
  require(RColorBrewer)
  #require(scales)
  #require(wesanderson)
  if (!('spData' %in% installed.packages())) stop('Need spData installed!')
  if (clear) {
    if (!(names(dev.cur()) == 'null device')) dev.off()
  }
  if (ret == 'plots') {
    mort_df <- FIA_mortality_with_explanatory
    mort_df <- mort_df[order(mort_df$PLT_CN), ]
    return(mort_df)
  } else {
    return(invisible())
  }
}
#' @describeIn RSFIA_utils Summarizes metadata from a tree file
#' @family FIA_processing
#' @export
PullTreeMeta <- function(df, out_cols) {
  #out_cols <- c('PLT_CN', 'INVYR', 'STATECD', 'COUNTYCD', 'PLOT')
  if (class(df) != 'data.frame') stop('Input must be dataframe')
  if (!(all(out_cols %in% colnames(df)))) stop('Missing columns')
  out <- data.frame(matrix(ncol = length(out_cols), nrow = 1))
  colnames(out) <- out_cols
  for (i in out_cols) {
    if (!(i %in% colnames(df))) next
    icol <- df[[i]]
    if (length(unique(icol)) != 1) {
      stop(paste('Meta column', i, 'has multiple values'))
    }
    out[[i]] <- unique(icol)
  }
  return(out)
}
#' @describeIn RSFIA_utils Summarizes metadata from a plot file
#' @family FIA_processing
#' @export
PullPlotMeta <- function(CNs, df, out_cols, db_ver = 7.2) {
  if (!(all(out_cols %in% colnames(df)))) stop('Missing columns')
  if (!db_ver == 7.2) stop('Requires FIADB version 7.2')
  out_df <- data.frame(matrix(ncol = length(out_cols) + 1, nrow = length(CNs)))
  colnames(out_df) <- c('PLT_CN', out_cols)
  for (i in 1:length(CNs)) {
    for (j in out_cols) {
      out_df[i, j] <- df[which(df$PLT_CN == CNs[i]), j]
    }
  }
  return(out_df)
}
#' @describeIn RSFIA_utils Pulls lat/longs
#' @family FIA_processing
#' @export
PullLatLongs <- function(plots) {
  if (!('LAT' %in% colnames(plots)) && ('LON' %in% colnames(plots))) {
    stop('Missing LAT and/or LON from plots dataframe, cant reference LAT/LON')
  }
  LAT <- vector('numeric', nrow(plots))
  LON <- vector('numeric', nrow(plots))
  for (i in 1:nrow(plots)) {
    LAT[i] <- plots$LAT[which(plots$CN == plots$CN[i])]
    LON[i] <- plots$LON[which(plots$CN == plots$CN[i])]
  }
  out_df <- data.frame(plots$CN, LAT, LON)
  colnames(out_df) <- c('PLT_CN', 'LAT', 'LON')
  return(out_df)
}
#' @describeIn RSFIA_utils Pulls a plot-level variable
#' @family FIA_processing
#' @export
PullPlotVars <- function(plots, vars) {
  out_df <- data.frame(matrix(ncol = length(vars) + 1, nrow = nrow(plots)))
  out_df[, 1] <- plots$CN
  cnt <- 1
  for (i in vars) {
    cnt <- cnt + 1
    if (!(i %in% colnames(plots))) stop('Missing var from plots df')
    out_df[, cnt] <- plots[[i]]
  }
  colnames(out_df) <- c('PLT_CN', vars)
  return(out_df)
}
#' @describeIn RSFIA_utils Converts the RDDISTCD to minimum distance to road in meters
#' @family FIA_processing
#' @export
RoadCodeConvert <- function(vec, db_ver = 7.2, na_fix = T) {
  if (db_ver != 7.2) stop('Unsupported db version')
  min_dist <- vector('numeric', length(vec))
  fix_val <- round(mean(vec))
  cnt <- 0
  for (i in vec) {
    cnt <- cnt + 1
    if (is.na(i)) {
      if (na_fix) {
        min_dist[cnt] <- fix_val
      } else {
        min_dist[cnt] <- NA
      }
      next
    }
    if (db_ver == 7.2) {
      min_dist[cnt] <- switch(
        i, 0, 101 * 0.3048, 301 * 0.3048, 501 * 0.3048, 1001 * 0.3048,
        0.5 * 5280 * 0.3048, 1 * 5280 * 0.3048, 3 * 5280 * 0.3048, 5 * 5280
      )
    }
  }
  return(min_dist)
}
#' @describeIn RSFIA_utils Converts WATERCD to categorical variable
#' @family FIA_processing
#' @export
WaterCodeConvert <- function(vec, db_ver = 7.2, na_fix = T) {
  if (db_ver != 7.2) stop('Unsupported db version')
  water_cat <- vector('character', length(vec))
  cnt <- 0
  # Fixes the vec code for 'other'
  vec <- ifelse(vec == 9, 6, vec)
  for (i in vec) {
    cnt <- cnt + 1
    if (is.na(i)) {
      if (na_fix) {
        water_cat[cnt] <- 0
      } else {
        water_cat[cnt] <- NA
      }
      next
    }
    if (db_ver == 7.2) {
      water_cat[cnt] <- switch(
        # As written, canals/ditches are 'temp'
        i + 1, 'none', 'perm', 'perm', 'temp', 'temp', 'flood', 'temp'
      )
    }
  }
  return(water_cat)
}
#' @describeIn RSFIA_utils Brute-force correlations with given explanatory variables
#' @family generic_utilities
#' @export
IterativeStatistics <- function(df, vars, type = 'adj_r2') {
  if (is.data.frame(df) == FALSE) stop('Need input data.frame')
  if (!(vars %in% colnames(df))) stop('Vars not in dataframe')
  accept_types <- c('adj_r2', 'signif')
  if (!(type %in% accept_types)) {
    message('Type invalid, accepted types are:')
    print(accept_types)
    stop()
  }
  ex <- colnames(df)[which(!(colnames(df) %in% vars))]
  out_df <- data.frame(matrix(ncol = length(vars), nrow = length(ex)))
  icnt <- 0
  jcnt <- 0
  for (i in vars) {
    icnt <- icnt + 1
    for (j in ex) {
      jcnt <- jcnt + 1
      fm <- as.formula(paste(i, j, sep = ' ~ '))
      mod <- lm(fm, df)
      if (type == 'adj_r2') {
        out_df[jcnt, icnt] <- summary(mod)$adj.r.squared
      }
      if (type == 'signif') {
        if (all(nrow(summary(mod)$coefficients) == 2,
                ncol(summary(mod)$coefficients) == 4)) {
          out_df[jcnt, icnt] <- summary(mod)$coefficients[2, 4]
        } else {
          out_df[jcnt, icnt] <- NA
        }
      }
    }
  }
  colnames(out_df) <- paste(vars, type, sep = '_')
  row.names(out_df) <- ex
  return(out_df)
}
#' @describeIn RSFIA_utils Numerical soil layer integration using trapezoidal rule
#' @family soil_data_utilities
#' @export
TrapezoidalIntegrateSoil <- function(d1, d2, ref) {
  # d1/d2 are the min/max depths of desired range, respectively, while
  # ref is the reference values at 0, 5, 15, 30, 60, 100, 200 cm respectively
  d <- c(0, 5, 15, 30, 60, 100, 200)
  nl <- which(d %in% c(d1:d2))
  z <- 0
  for (k in 1:(length(nl) - 1)) {
    dd <- d[nl[k + 1]] - d[nl[k]]
    val <- ref[nl[k]] + ref[nl[k + 1]]
    z <- z + (dd * val)
  }
  l <- 0.5 * (1 / (d2 - d1))
  return(round(z * l, 2))
}
#' @describeIn RSFIA_utils Displays Daymet variable info
#' @family daymet_utilities
#' @export
ShowDaymetMeta <- function() {
  raw_import <- c('year', 'yday', 'dayl..s.', 'prcp..mm.day.', 'srad..W.m.2.',
                  'swe..kg.m.2.', 'tmax..deg.c.', 'tmin..deg.c.', 'vp..Pa.')
  fixed_name <- c('year', 'doy', 'day_length', 'precipitation', 'solar_radiation',
                  'snow_water_equivalent', 'max_temp', 'min_temp', 'water_vapor_pressure')
  units_short <- c(NA, NA, 's', 'mm day^{-1}', 'W m^{-2}',
                   'kg m^{-2}', 'deg C', 'deg C', 'Pa')
  units_long <- c(NA, NA, 'seconds', 'mm / day', 'Watts / square meter',
                  'kilograms / square meter', 'degrees Celcius', 'degrees Celcius',
                  'Pascals')
  out_df <- data.frame(raw_import, fixed_name, units_short, units_long,
                       stringsAsFactors = F)
  return(out_df)
}
#' @describeIn RSFIA_utils Find min/max quarterly value
#' @family daymet_utilities
#' @export
MinMaxByQuarter <- function(x, n = 90, fun = sum) {
  roll <- zoo::rollapply(x, n, fun)
  min_val <- min(roll, na.rm = T)
  max_val <- max(roll, na.rm = T)
  min_loc <- round(median(which(roll == min_val), na.rm = T))
  max_loc <- round(median(which(roll == max_val), na.rm = T))
  data.frame(min_val, min_loc, max_val, max_loc)
}
#' @describeIn RSFIA_utils Pastes a unique factor variable to other columns in input df
#' @family generic_utilities
#' @export
PasteUniqueFactor <- function(df, uniq = 1) {
  # Unique identifies the column# of the paste variable
  fac <- df[, uniq]
  if (length(unique(fac)) > 1) stop('Uniq col has >1 value')
  out_df <- df[, -uniq]
  new_cols <- paste(colnames(df)[uniq], unique(fac), colnames(out_df), sep = '_')
  colnames(out_df) <- new_cols
  return(out_df)
}
#' @describeIn RSFIA_utils Combines identical variables, e.g. with weather data
#' @family package_utilities
CollapseIdenticalVars <- function(df) {
  # Might not work, the issue this was supposed to fix was fixed upstream instead
  out_vals <- data.frame(matrix(ncol = 0, nrow = nrow(df)))
  out_df <- df
  kill_cols <- numeric()
  for (i in 1:ncol(df)) {
    chk_df <- df[, -i]
    i_eq <- colSums(chk_df != df[, i], na.rm = T)
    dups <- which(i_eq == 0)
    if (length(dups) == 0) {
      next
    }
    kill_cols <- append(kill_cols, which(colnames(df) %in% names(dups)))
    list_names <- as.data.frame(strsplit(names(dups), '_'), stringsAsFactors = F)
    matching <- rowSums(list_names[, 1] != list_names)
    new_name <- list_names[!as.logical(matching), 1]
    new_name <- paste(new_name[which(!(new_name %in% 'year'))], collapse = '_')
    new_val <- df[, i]
    out_vals <- data.frame(out_vals, new_val)
    colnames(out_vals)[ncol(out_vals)] <- new_name
    out_df <- out_df[, -kill_cols]
  }
  out_df <- data.frame(out_df, out_vals)
  return(out_df)
}
#' @describeIn RSFIA_utils Lists available mortality variables
#' @family package_utilities
#' @export
ListMortVars <- function() {
  if (length(grep('Windows', sessionInfo()$running)) < 1) {
    stop('\nNot sure this will work under not-windows!\n
         Check for FIA_mortality_variables_list.csv manually in RSFIA/metadata')
  }
  pkg_loc <- find.package('RSFIA')
  chk_file <- paste(pkg_loc, 'metadata/FIA_mortality_variable_list.csv', sep = '/')
  out_df <- read.csv(chk_file, stringsAsFactors = F)
  return(out_df)
}
#' @describeIn RSFIA_utils Print r2 for RF object
#' @family generic_utilities
#' @export
r2_rf <- function(rf) {
  require(randomForest, quietly = T)
  out_pred <- predict(rf)
  y_in <- rf$y
  r2 <- round(summary(lm(y_in ~ out_pred))$r.squared, 3)
  return(r2)
}
#' @describeIn RSFIA_utils Adds Landsat/Daymet data by lat/long
#' @family remote_sensing
#' @export
CombineLandsatDaymet <- function(df, ...) {
  args <- list(...)
  for (i in 1:length(args)) {
    i_arg <- args[[i]]
    colnames(i_arg)[1:2] <- c('LAT', 'LON')
    df <- dplyr::left_join(df, i_arg, by = c('LAT', 'LON'))
  }
  return(df)
}
#' @describeIn RSFIA_utils Converts subplot size, # subplots to acres
#' @family FIA_processing
#' @export
SubsToAcres <- function(n, r = 24, hectare = F) {
  # Radius is in feet
  a <- n * pi * r ^ 2
  # a is in square feet, 43560 per acre
  ac <- a / 43560
  if (hectare) {
    return(ac / 2.47105)
  } else {
    return(ac)
  }
}
#' @describeIn RSFIA_utils Provides a keymatch function for AGENTCD
#' @family FIA_processing
#' @export
KeyAgentCode <- function(db_ver) {
  if (db_ver == 7.2) {
    out_fun <- function(x) {
      if (length(x) == 0) return(character())
      c <- character(length(x))
      for (i in 1:length(x)) {
        if (is.na(x[i])) {
          c[i] <- 'Missing'
        } else if (x[i] == 0) {
          c[i] <- 'No Mortality'
        } else {
          c[i] <- switch(x[i] / 10,
                         'Insect', 'Disease', 'Fire', 'Animal',
                         'Weather', 'Vegetation', 'Other', 'Harvest')
        }
      }
      return(c)
    }
    return(out_fun)
  } else {
    stop('only db 7.2 currently supported')
  }
}
#' @describeIn RSFIA_utils Provides a keymatch for SPCD
#' @family FIA_processing
#' @export
KeySpeciesCode <- function(db_ver, by = 'species') {
  if (db_ver == 7.2) {
    out_fun <- function(x) {
      require(RSFIA)
      sp_key <- FIA_table_metadata$APPENDIX_F # F is species codes
      if (length(x) == 0) return(character())
      c <- sp_key$Scientific.Name[match(x, sp_key$SPCD)]
      if (by == 'genera') {
        browser()
      }
      return(c)
    }
  }
  return(out_fun)
}
#' @describeIn RSFIA_utils Returns STATECD metadata
#' @family FIA_processing
#' @export
MakeStateCodeKey <- function(ver) {
  ret_df <- data.frame(state.abb, rep(NA, length(state.abb)), stringsAsFactors = F)
  if (ver == 7.2) {
    colnames(ret_df) <- c('state_abb', 'STATECD')
    # From main lineup: missing is 3, 7, 11, 14, 43, 52
    # Last in main lineup is Wyoming at 56
    cd <- 1:56
    cd <- cd[-c(3, 7, 11, 14, 43, 52)]
    ret_df[1:50, 2] <- cd
    # Territories: Dist. Columbia, Samoa, Micronesia, Guam, Marshall,
    # Mariana, Palau, PR, USVI.
    # 11 is the District of Columbia
    abb <- c('DC', 'AS', 'FM', 'GU', 'MH', 'MP', 'PW', 'PR', 'VI')
    cd <- c(11, 60, 64, 66, 68, 69, 70, 72, 78)
    ret_df[51:59, 1] <- abb
    ret_df[51:59, 2] <- cd
  } else {
    stop('Unsupported FIADB Version')
  }
  return(ret_df)
}
#' @describeIn RSFIA_utils Keys STATECD
#' @family FIA_processing
#' @export
KeyStateCode <- function(x, db_ver) {
  out <- character(length = length(x))
  if (db_ver == 7.2) {
    key <- RSFIA::MakeStateCodeKey(ver = 7.2)
    out <- key$state_abb[match(x, key$STATECD)]
  }
  return(out)
}
#' @describeIn RSFIA_utils Subtraction operator that doesnt fail on NAs
#' @family generic_utilities
#' @export
`%NA-%` <- function(e1, e2) {
  r <- e1
  for (i in e2) {
    if (is.na(i)) next
    r <- r - i
  }
  r
}
#' @describeIn RSFIA_utils Replaces species codes with names
#' @family FIA_processing
#' @export
ReplaceSPCD <- function(x, by = 'species', db_ver = 7.2) {
  key <- KeySpeciesCode(db_ver, by = by)
  cd_indx <- grep('SPCD', x)
  code <- sapply(strsplit(x[cd_indx], '_'), function(z) {
    suppressWarnings(na.omit(as.numeric(z)))
  })
  sp <- key(code)
  sp_list <- sapply(as.list(sp), function(z) strsplit(z, ' '))
  sp_list <- sapply(sp_list, function(z) paste(z, collapse = '_'))
  in_names <- x[cd_indx]
  out_names <- character(length(in_names))
  for (i in 1:length(in_names)) {
    out_names[i] <- sub(paste0('SPCD_', code[i]), sp_list[i], in_names[i])
  }
  x[cd_indx] <- out_names
  return(x)
}
#' @describeIn RSFIA_utils Keys the forest type codes
#' @family FIA_processing
#' @export
KeyForestType <- function(x, db_ver = 7.2) {
  key_df <- FIA_table_metadata[['APPENDIX_D']]
  out_vec <- key_df$Forest.type...type.group[match(x, key_df$Code)]
  return(out_vec)
}
#' @describeIn RSFIA_utils Multi-panel ggplot wrapper
#' @family generic_utilities
#' @export
Multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  # BEM note: this is pretty much verbatim copied from the R cookbook
  # Multiple plot function
  #
  # ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
  # - cols:   Number of columns in layout
  # - layout: A matrix specifying the layout. If present, 'cols' is ignored.
  #
  # If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
  # then plot 1 will go in the upper left, 2 will go in the upper right, and
  # 3 will go all the way across the bottom.
  require(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }

  if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
#' @describeIn RSFIA_utils Keys OWNCD
#' @family FIA_processing
#' @export
KeyOwnerCode <- function(x, ret = 'logical', db_ver = 7.2) {
  if (db_ver == 7.2) {
    # Key: 11-13 USFS, 21 NPS, 22 BLM, 23 FWS, 24 DOD/DOE, 25 Other Fed
    # 31 State, 32 Local, 33 Other Public, 46 Private/NA
    if (ret == 'logical') {
      is_public <- ifelse(x < 46, TRUE, FALSE)
      out_vec <- is_public
    } else {
      stop('not implemented')
    }
  } else {
    stop('unsupported FIA database version')
  }
  return(out_vec)
}
#' @describeIn RSFIA_utils Aggregation function for calculating dominant AGENTCD
#' @family FIA_processing
#' @export
AggregateDominantAgent <- function(in_df, cname, join_df = NULL, FUN0, join = F,
                                   scale = NULL) {
  x <- in_df$percent_AGENTCD
  y <- in_df$dominant_AGENTCD
  sx <- ifelse(length(scale) == 0, nrow(in_df), scale)
  z_df <- aggregate(x, by = list(y), FUN = FUN0)
  z_df[, 2] <- round(z_df[, 2] / sx, 3)
  cat('\nSum of aggregation categories:', sum(z_df[, 2]), '\n')
  if (sum(z_df[, 2]) == 0) stop('scaling failed')
  colnames(z_df) <- c('Dominant_AGENTCD', cname)
  if (cname == '%_Total') {
    z_df$perc_as_char <- as.character(z_df[[cname]])
  }
  if (join) {
    df0 <- dplyr::left_join(join_df, z_df, by = 'Dominant_AGENTCD')
  } else {
    df0 <- z_df
  }
  return(df0)
}
#' @describeIn RSFIA_utils Wrapper for stacked bar plot
#' @family plot_wrappers
#' @export
PlotStacked <- function(stacked_df, by = 'Section', bar_fill = NULL, llab, yylab) {
  stacked_bar <- ggplot(stacked_df) + theme_bw() +
    geom_bar(aes_string(by, fill = bar_fill)) +
    xlab('') + ylab(yylab) +
    guides(fill = guide_legend(title = llab)) +
    scale_fill_brewer(palette = 'Paired') +
    coord_flip()
  return(stacked_bar)
}
#' @describeIn RSFIA_utils Standard error
#' @family generic_utilities
#' @export
se <- function(x) {
  y <- sd(x, na.rm = T) / sqrt(length(x))
  y
}
#' @describeIn RSFIA_utils Barplot wrapper
#' @family plot_wrappers
#' @export
PlotBarPlot <- function(df, xx = 1, yy = 2, SE = 3, ylab) {
  xn <- colnames(df)[xx]
  yn <- colnames(df)[yy]
  y0 <- df[, yy]
  df[, 'SE_min'] <- y0 - df[, SE]
  df[, 'SE_max'] <- y0 + df[, SE]
  mx <- max(df[, 'SE_max']) * 1.1
  mn <- min(df[, 'SE_min']) * 0.9
  p0 <- ggplot(data = df) + theme_bw() +
    geom_bar(aes_string(x = xn, y = yn), stat = 'identity') +
    geom_errorbar(aes_string(
      x = xn,
      ymin = 'SE_min',
      ymax = 'SE_max'), width = 0.2) +
    #scale_y_continuous(limits = c(0, 100)) # breaks if limit[1] > 0
    coord_cartesian(ylim = c(mn, mx)) +
    xlab('') + ylab(ylab)
  print(p0)
  invisible()
}
#' @describeIn RSFIA_utils Prints loop progress
#' @family generic_utilities
#' @export
LoopStatus <- function(from, to, big_inc = 1000, digits = 1) {
  if (to > big_inc) {
    if (from %% (big_inc / 10) == 0) cat('[]')
    if (from %% big_inc != 0) {
      return(invisible())
    }
  }
  prog <- from / to * 100
  prog <- round(prog, digits)
  cat(' ', prog, '%\n')
  invisible()
}
#' @describeIn RSFIA_utils Bins common forest types
#' @family package_utilities
#' @export
BinUncommonForestType <- function(mort_df = NULL) {
  if (length(mort_df) == 0) {
    mort_df <- FIA_mortality_with_explanatory
  }
  uc_tbl <- table(mort_df[which(mort_df$is_common_ftype == F), 'forest_type'])
  uc0 <- KeyForestGroup(names(uc_tbl))
  cf <- which(mort_df$forest_type %in% names(uc_tbl)[which(uc0 == 'conifer')])
  hf <- which(mort_df$forest_type %in% names(uc_tbl)[which(uc0 == 'hwoods')])
  mort_df[cf, 'forest_type'] <- 'Other conifers'
  mort_df[hf, 'forest_type'] <- 'Other hardwoods'
  mort_df[which(is.na(mort_df$forest_type)), 'forest_type'] <- 'Non-stocked'
  # Returns modified mort_df as well as a list of the uncommon forest types + counts
  m0 <- list(mort_df, uc_tbl)
  return(m0)
}
#' @describeIn RSFIA_utils Subsets/excludes dominant forest types
#' @family package_utilities
#' @export
SubDominantForestType <- function(mort_df = NULL, sub = T, n_plot_cutoff = 400) {
  if (length(mort_df) == 0) {
    mort_df <- FIA_mortality_with_explanatory
  }
  ftype_tbl <- table(mort_df$forest_type)
  dom_ftype <- names(ftype_tbl)[which(ftype_tbl > n_plot_cutoff)]
  if (sub) {
    mort0 <- mort_df[which(mort_df$forest_type %in% dom_ftype), ]
  } else {
    mort0 <- mort_df[-which(mort_df$forest_type %in% dom_ftype), ]
  }
  return(mort0)
}
#' @describeIn RSFIA_utils Keys hardwood/softwood from forest type
#' @family FIA_processing
#' @export
KeyForestGroup <- function(x, db_ver = 7.2) {
  x0 <- character(length(x))
  if (db_ver == 7.2) {
    conifers <- c(
      'Alaska-yellow-cedar', 'Blue spruce', 'Foxtail pine / bristlecone pine',
      'Gray pine', 'Incense-cedar', 'Knobcone pine', 'Miscellaneous western softwoods',
      'Port-Orford-cedar', 'Redwood', 'Sitka spruce', 'Sugar pine', 'Western white pine',
      'California mixed conifer', 'Douglas-fir', 'Engelmann spruce',
      'Engelmann spruce / subalpine fir', 'Grand fir', 'Jeffrey pine', 'Juniper woodland',
      'Mountain hemlock', 'Noble fir', 'Pacific silver fir', 'Ponderosa pine',
      'Pinyon / juniper woodland', 'Red fir', 'Rocky Mountain juniper',
      'Subalpine fir', 'Western juniper', 'Western hemlock', 'White fir',
      'Whitebark pine', 'Western larch', 'Lodgepole pine', 'Western redcedar',
      'Limber pine'
      )
    hwoods <- c(
      'California laurel', 'California white oak (valley oak)', 'Coast live oak',
      'Cottonwood', 'Cottonwood / willow', 'Giant chinkapin', 'Oregon ash', 'Paper birch',
      'Bigleaf maple', 'Blue oak', 'California black oak', 'Canyon live oak',
      'Cercocarpus (mountain brush) woodland', 'Deciduous oak woodland',
      'Interior live oak', 'Intermountain maple woodland', 'Oregon white oak',
      'Other hardwoods', 'Pacific madrone', 'Red alder', 'Tanoak', 'Aspen'
    )
    if (all(x %in% c(NA, conifers, hwoods)) == F) {
      mis <- which((x %in% c(NA, conifers, hwoods)) == F)
      cat('Not in KeyForestGroup key:', unique(x[mis]), '\n')
      stop('KeyForestGroup failed, missing forest groups need to be added by hand')
    }
    x0[which(x %in% conifers)] <- 'conifer'
    x0[which(x %in% hwoods)] <- 'hwoods'
    x0[which(x0 == '')] <- NA
  } else {
    stop('unsupported db_ver')
  }
  return(x0)
}
#' @describeIn RSFIA_utils Aggregates forest_type by dominant Cleland_province
#' @family package_utilities
#' @export
AggForTypeByProv <- function(mort_df = NULL) {
  if (is.null(mort_df)) {
    mort_df <- FIA_mortality_with_explanatory
  }
  m0 <- mort_df
  for (i in unique(mort_df$forest_type)) {
    i_sub <- mort_df[which(mort_df$forest_type == i), ]
    wdom <- which.max(table(i_sub$Cleland_province))
    np <- names(table(i_sub$Cleland_province))[wdom]
    m0[which(m0$forest_type == i), 'Cleland_province'] <- np
  }
  return(m0)
}
#' @describeIn RSFIA_utils Transpose/mirror an uneven matrix
#' @family generic_utilities
#' @export
TransposeUnevenMatrix <- function(x, fill = F) {
  stopifnot(is.matrix(x))
  cx <- colnames(x)
  rx <- row.names(x)
  stopifnot(!any(is.null(cx), is.null(rx)))
  c0 <- which(!(cx %in% rx))
  r0 <- which(!(rx %in% cx))

  x <- rbind(rep(NA, ncol(x)), x)
  row.names(x)[c0] <- cx[c0]
  x <- cbind(x, rep(NA, nrow(x)))
  colnames(x)[r0 + 1] <- rx[r0]
  y <- t(x)

  if (fill) {
    z <- ifelse(is.na(x), y, x)
    diag(z) <- NA
  } else {
    z <- y
  }
  z
}
#' @describeIn RSFIA_utils Group a logical matrix by row/col names
#' @family generic_utilities
#' @export
GroupLogicalMatrix <- function(x) {
  stopifnot(is.matrix(x))
  cx <- colnames(x)
  rx <- row.names(x)
  stopifnot(!any(is.null(cx), is.null(rx)))
  stopifnot(all(cx == rx))
  letters_ext <- c(letters, LETTERS)
  colnames(x) <- letters_ext[1:ncol(x)]
  y <- x
  y[which(y == F)] <- NA
  y[, ] <- as.character(y[, ])
  y[which(is.na(y))] <- ''
  for (i in 1:ncol(x)) {
    y[which(x[, i] == T), i] <- colnames(y)[i]
  }
  z <- apply(y, 1, function(x) paste(x, collapse = ''))
  #zl <- strsplit(z, '')
  #ll <- which(!letters_ext %in% unlist(strsplit(z, '')))[1]
  #pl <- letters_ext[1:ncol(x)]
  #for (i in pl[1:length(zl)]) {
  #  ml <- lapply(zl, function(x) i %in% x)
  #  if (all(unlist(ml))) print('what')
  #}
  z
}
#' @describeIn RSFIA_utils Generates summary table for binned mortality using IQR
#' @family generic_utilities
#' @export
CalcCategoricalMortByIQR <- function(x, group) {
  stopifnot(length(x) == length(group))
  in0 <- data.frame(x, group, stringsAsFactors = F)
  iqr <- aggregate(x, by = list(group), FUN = function(z) {
    y <- IQR(z) * 1.5
  })
  colnames(iqr) <- c('group', 'outlier_val')
  out <- dplyr::left_join(in0, iqr, by = 'group')
  out$is_outlier <- ifelse(out$x > out$outlier_val, TRUE, FALSE)
  colnames(out) <- c('mort_val', 'group', 'outlier_val', 'high_mort_logical')

  return(out)
}
#' @describeIn RSFIA_utils Generates mort cause cross-table
#' @family package_utilities
#' @export
ComorbidAgentTable <- function(mort_df = NULL,
                               only_multiple = F, only_mort = F) {
  if (only_multiple && only_mort) {
    warning('Multiple category is by definition only mort, ignoring only_mort')
  }
  if (is.null(mort_df)) {
    mort_df <- FIA_mortality_with_explanatory
  }
  if (only_mort) {
    mort_df <- mort_df[which(mort_df$mort_rate > 0), ]
  }
  if (only_multiple) {
    mort_df <- mort_df[which(mort_df$dominant_AGENTCD == 'Multiple'), ]
  }
  pos_agent <- c('Animal', 'Disease', 'Fire', 'Harvest', 'Insect',
                 'Other', 'Vegetation', 'Weather')
  suffix <- paste0('_ndead')
  m0 <- diag(length(pos_agent))
  colnames(m0) <- pos_agent
  row.names(m0) <- pos_agent
  for (i in 1:length(pos_agent)) {
    ia <- paste0(pos_agent[i], suffix)
    #ja <- pos_agent[-i]
    ja <- pos_agent
    for (j in 1:length(ja)) {
      if (j <= i) {
        m0[j, i] <- NA
        next
      }
      jb <- paste0(ja[j], suffix)
      idf <- mort_df[, c(ia, jb)]
      idf[, ] <- ifelse(idf[, ] > 0, T, F)
      m0[j, i] <- round(sum(rowSums(idf) > 1) / nrow(idf) * 100, 3)
    }
  }
  diag(m0) <- NA

  return(m0)
}
#' @describeIn RSFIA_utils Statistical test for tree vs plot AGENTCD tallies
#' @family package_utilities
AgentTypeDiffTest <- function() {
  stop()
  tree_tbl0 <- table(FIA_mortality_TREE_level[, 'AGENTCD'])
  tree_tbl1 <- table(FIA_mortality_TREE_level[, 'AGENTCD'])
  mort_tbl0 <- table(FIA_mortality_with_explanatory[, 'dominant_AGENTCD'])

  filtTree <- function(tree_tbl) {
    names(tree_tbl) <- KeyAgentCode(7.2)(as.numeric(names(tree_tbl)))
    tree_tbl <- tree_tbl[order(names(tree_tbl))]
    tree_tbl
  }

  filtMort <- function(mort_tbl) {
    mort_tbl <- mort_tbl[-which(names(mort_tbl) == 'No Mortality')]
    mort_tbl
  }


  calcTbl <- function(tbl1, tbl2, cnm = c('n_tree', 'n_plot')) {
    tbl0 <- matrix(ncol = 2, nrow = length(tbl1))
    tbl0[, 1] <- tbl1
    tbl0[, 2] <- tbl2
    row.names(tbl0) <- names(tbl1)
    colnames(tbl0) <- cnm
    tbl0
  }

  ftbl <- calcTbl(filtTree(tree_tbl0), filtMort(mort_tbl0))
  btbl <- calcTbl(filtTree(tree_tbl1), filtMort(mort_tbl1))
  otbl <- calcTbl(filtTree(tree_tbl2), filtMort(mort_tbl2))

  tp_df <- tp_df[order(tp_df[, 1]), ]

  type <- c(rep('plot', nrow(pp_df)), rep('tree', nrow(tp_df)))
  df0 <- data.frame(matrix(ncol = 1 + ncol(pp_df), nrow = length(type)),
                    stringsAsFactors = F)
  df0[, 1] <- type
  df0[c(1:nrow(pp_df)), c(2:ncol(df0))] <- pp_df[, ]
  df0[c((nrow(pp_df) + 1):(nrow(df0))), c(2:ncol(df0))] <- tp_df[, ]
  colnames(df0) <- c('type', 'AGENTCD', colnames(pp_df)[2:ncol(pp_df)])

  for (i in colnames(df0)[3:ncol(df0)]) {
    df1 <- df0[, c('type', 'AGENTCD', i)]
    browser()
  }
}
#' @describeIn RSFIA_utils Keys damage codes
#' @family package_utilities
#' @export
KeyDamageAgentCode <- function(x, db_ver = 7.2) {
  if (db_ver == 7.2) {
    stopifnot(all(range(x, na.rm = T)[1] >= 0, range(x, na.rm = T)[2] <= 1e5))
    y <- character(length(x))
    y[which(x == 0)] <- 'No Damage'
    y[intersect(which(x > 0), which(x < 11000))] <- 'General Insects'
    y[intersect(which(x >= 11000), which(x < 12000))] <- 'Bark Beetles'
    y[intersect(which(x >= 12000), which(x < 13000))] <- 'Defoliators'
    y[intersect(which(x >= 13000), which(x < 14000))] <- 'Chewing Insects'
    y[intersect(which(x >= 14000), which(x < 15000))] <- 'Sucking Insects'
    y[intersect(which(x >= 15000), which(x < 16000))] <- 'Boring Insects'
    y[intersect(which(x >= 19000), which(x < 21000))] <- 'General Diseases'
    y[intersect(which(x >= 21000), which(x < 22000))] <- 'Root/Butt Diseases'
    y[intersect(which(x >= 22000), which(x < 22500))] <- 'Cankers'
    y[which(x == 22500)] <- 'Stem Decay Fungi'
    y[intersect(which(x >= 23000), which(x < 24000))] <- 'Parasitic Plants'
    y[intersect(which(x >= 24000), which(x < 25000))] <- 'General Dieback'
    y[intersect(which(x >= 25000), which(x < 26000))] <- 'Leaf Diseases'
    y[intersect(which(x >= 26000), which(x < 27000))] <- 'Stem Rust'
    y[intersect(which(x >= 27000), which(x < 30000))] <- 'Broom Rust'
    y[intersect(which(x >= 30000), which(x < 41000))] <- 'Fire'
    y[intersect(which(x >= 41000), which(x < 42000))] <- 'Animals, Wild'
    y[intersect(which(x >= 42000), which(x < 50000))] <- 'Animals, Domestic'
    y[intersect(which(x >= 50000), which(x < 60000))] <- 'Abiotic'
    y[intersect(which(x >= 60000), which(x < 70000))] <- 'Competition'
    y[intersect(which(x >= 70000), which(x < 80000))] <- 'Humans'
    y[intersect(which(x >= 90000), which(x < 99000))] <- 'Other'
    y[which(x == 99000)] <- 'Unknown'

  } else {
    stop('unsupported db_ver')
  }
  y <- ifelse(y == '', NA, y)
  return(y)
}