In bmcnellis/RSFIA: FIA Tree Mortality Analysis

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This vignette is for creating the FIA/mortality dataset to be used in model building. The data is in FIA_mortality_with_explanatory.

#stop('Dont knit this.')

# Methods subtitle: 'Forest Inventory Analysis data'
rm(list = ls())
library(RSFIA)
library(ClelandEcoregions)
library(sp)
library(raster)
library(rgdal)
orig_dir <- getwd()

if (Sys.info()['sysname'] == 'Linux') {
  fia_dir <- '/media/bem/Q_STORAGE/extdata/FIA/data'
  data_dir <- '/home/bem/Documents/docs_linux/R_workspace/RSFIA/inst/extdata'
  fd <- NULL # TODO fix this
  cond_dir <- '/media/bem/Q_STORAGE/extdata/FIA/data'
} else {
  fia_dir <- 'C:/Users/Brandon/Documents/docs/PHD/FIA/data'
  data_dir <- "C:/Users/Brandon/Documents/docs/R_workspace/RSFIA/inst/extdata"
  # Fire perim directory:
  fd <- 'C:/Users/Brandon/Documents/docs/PHD/RSFIA_project/GeoMAC'
  cond_dir <- 'C:/Users/Brandon/Documents/docs/PHD/FIA/data'
}
# Flags:
download_soil <- F
drop_fire <- F
drop_low_ftype <- F
bin_multiple <- F
drop_nonstocked <- F

# Section 1: Data Imports
data(Cleland_meta_df)
# Divisions used here are Temperate Desert, Warm Continent, Medit., Temp Steppe,
# and Hot Continental:
divis <- unique(Cleland_meta_df$division)[c(1:5, 7)]
provs <- unique(Cleland_meta_df$province_code[which(Cleland_meta_df$division %in% divis)])
# ImportPlots requires the full code to work with ECOSUBCD:
subsecs <- Cleland_meta_df$subsection_code[which(Cleland_meta_df$province_code %in% provs)]
subsecs <- unique(subsecs)
plots <- ImportPlots(ecoregion = subsecs, dir = fia_dir, remeas = 4,
                            status = c(1, 2))
unique_key <- GenerateUniquePlot(plots)
colnames(unique_key)[1] <- 'CN'
plots <- dplyr::full_join(plots, unique_key, by = 'CN')
bad_plots <- CheckUniquePlot(plots, plot_var = 'unique_plot_key')
if (length(bad_plots) > 0) {
  plots <- plots[-which(plots$CN %in% bad_plots), ]
}
bad_plots <- CheckUniquePlot(plots, plot_var = 'unique_plot_key')
if (length(bad_plots) > 0) stop('Didnt catch all bad plots!')
if (length(which(plots$ECOSUBCD == 'Water')) > 0) {
  plots <- plots[-which(plots$ECOSUBCD == 'Water'), ]
}
# Drop the stray plots from 2006/2016:
plots <- plots[-which(plots$INVYR == 2016), ]
plots <- plots[-which(plots$INVYR == 2006), ]
# Add Cleland tags:
Cleland_section <- Cleland_meta_df$section_code[match(plots$ECOSUBCD, 
                                                      Cleland_meta_df$subsection_code)]
Cleland_province <- Cleland_meta_df$province_code[match(plots$ECOSUBCD, 
                                                        Cleland_meta_df$subsection_code)]
plots <- data.frame(plots, Cleland_province, Cleland_section, stringsAsFactors = F)

conds <- FilterMortByCondition(dir = cond_dir, data = plots)

trees <- ImportTrees(plots = plots, dir = fia_dir)
trees$DIA[which(is.na(trees$DIA))] <- trees$DIACALC[which(is.na(trees$DIA))]
trees$DIA_cm <- trees$DIA * 2.54
if (length(which(is.na(trees$DIA))) > 0) {
  trees <- trees[-which(is.na(trees$DIA)), ]
}

# Calculate mortality:
mort_df <- MakeMortalityRates(trees = trees, plots = plots, debug = F)
tree_CNs_used <- mort_df[[2]]
mort_df <- mort_df[[1]]
trees <- trees[which(trees$CN %in% tree_CNs_used), ]

# Take out the plots which broke the mort_rate calculation
if (sum(is.na(mort_df$mort_rate)) > 0) {
  mort_df <- mort_df[-which(is.na(mort_df$mort_rate)), ]
}
# Drop the small-province plots:
excl_provs <- names(which(table(mort_df$Cleland_province) < 500))
mort_df <- mort_df[-which(mort_df$Cleland_province %in% excl_provs), ]
excl_sects <- names(which(table(mort_df$Cleland_section) < 100))
mort_df <- mort_df[-which(mort_df$Cleland_section %in% excl_sects), ]

# Add condition variables:
colnames(conds)[1] <- c('PREV_PLT_CN')
mort_df <- dplyr::left_join(mort_df, conds, by = 'PREV_PLT_CN')

# Add a logical variable for 'common forest types'
#n_bef <- nrow(mort_df)
for_typ_counts <- as.data.frame(table(mort_df$FORTYPCD))
for_typ_counts <- for_typ_counts[which(for_typ_counts$Freq >= 30), ]
# This 'mort_df <-' line drops all the dead fire trees...
ftype_indx <- which(mort_df$FORTYPCD %in% for_typ_counts$Var1)
is_common_ftype <- logical(nrow(mort_df))
is_common_ftype[ftype_indx] <- TRUE
mort_df <- cbind(mort_df, is_common_ftype)

bef <- nrow(mort_df)
if (drop_fire) {
  # Drop fire perimeters:
  debug_fire <- DropFirePerimeters(mort_df, fire_dir = fd, yrs = c(2000, 2015))
  #mort_df <- DropFirePerimeters(mort_df, fire_dir = fd, yrs = c(2000, 2015))
}
aft <- nrow(mort_df)
message('# plots dropped in fire perimeters:')
cat(bef - aft, '\n')

# Fix trees from dropped plots:
trees <- trees[which(trees$PLT_CN %in% mort_df$PLT_CN), ]

soil_grid_vars <- read.csv(list.files()[grep('soilgrids', list.files())])
soil_grid_vars <- soil_grid_vars[, -c(1, 2)]
soil_grid_vars[, c(3:ncol(soil_grid_vars))] <- lapply(
  soil_grid_vars[, c(3:ncol(soil_grid_vars))], function(x) ifelse(x < 0, NA, x))
colnames(soil_grid_vars)[1:2] <- c('LAT', 'LON')
mort_df <- dplyr::left_join(mort_df, soil_grid_vars, by = c('LAT', 'LON'))
mort_df$BDTICM <- round(mort_df$BDTICM / 1000, 3) # scales depth-of-bedrock to m
mort_df$BLDFIE <- round(mort_df$BLDFIE / 1000, 3) # scales bulk density
if (drop_nonstocked) {
  bef <- nrow(mort_df)
  if (length(which(mort_df$FORTYPCD == 999) > 0)) {
    mort_df <- mort_df[-which(mort_df$FORTYPCD == 999), ]
  }
  cat('Dropped', bef - nrow(mort_df), 'non-stocked plots.\n')
}
mort_df$forest_type <- KeyForestType(mort_df$FORTYPCD) # Keys forest type
wth <- FIA_weather_by_year_bin
mort_df <- dplyr::left_join(mort_df, wth, by = c('LAT', 'LON'))

# Mort outliers:
mort_outlier <- CalcCategoricalMortByIQR(
  mort_df$mort_rate, group = mort_df$Cleland_section)[, 4]
mort_df <- data.frame(mort_df, mort_outlier, stringsAsFactors = F)

# Dominant mort causes:
if (bin_multiple) {
  mc0 <- 75
} else {
  mc0 <- 0
}
dom_mort <- CalcDominantMortCause(tree_data = trees, multiple_cutoff = mc0)
b_ndead <- mort_df$n_dead
mort_df <- dplyr::left_join(mort_df, dom_mort, by = c('PLT_CN', 'n_dead'))
stopifnot(identical(b_ndead, mort_df$n_dead))

# Pub/priv logical:
pub_priv <- PullOwnerCode(CNs = mort_df$PLT_CN, dir = fia_dir, quietly = T)
pub_priv <- pub_priv[, c('PLT_CN', 'is_public')]
mort_df <- dplyr::left_join(mort_df, pub_priv, by = 'PLT_CN')

# Cleland by name instead of code:
mort_df$section <- KeyClelandCode(mort_df$Cleland_section, lvl = 'section')
mort_df$province <- KeyClelandCode(mort_df$Cleland_province, lvl = 'province')

# Add fire perimeter logical:
mort_df <- DropFirePerimeters(mort_df, fire_dir = fd, as_logical = T)

mort_df <- mort_df[, order(colnames(mort_df))]
FIA_mortality_with_explanatory <- mort_df
FIA_mortality_TREE_level <- trees
setwd(orig_dir)
message('Sampled ecoregions, with plot counts:')
print(table(mort_df$Cleland_province))
message('Included columns:')
print(colnames(mort_df))
message('Year range and plot sample size, first census:')
print(table(mort_df$prev_year))
message('Year range and plot sample size, second census:')
print(table(mort_df$INVYR))
message('Total number of plots:')
print(nrow(mort_df))
message('Total number of trees:')
print(nrow(trees))
message('Date analysis run:')
print(Sys.time())
library(devtools)
use_data(FIA_mortality_TREE_level, overwrite = T)
use_data(FIA_mortality_with_explanatory, overwrite = T)