R/diversityNew.R
In hunter-stanke/rFIA_TX: Duplicate rFIA for TX
Defines functions
diversity_backup
diversity
diversityStarter
Documented in
diversity
diversityStarter <- function(x,
db,
grpBy_quo = NULL,
polys = NULL,
returnSpatial = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
stateVar = TPA_UNADJ,
grpVar = SPCD,
treeDomain = NULL,
areaDomain = NULL,
byPlot = FALSE,
totals = FALSE,
nCores = 1,
remote,
mr){
reqTables <- c('PLOT', 'TREE', 'COND', 'POP_PLOT_STRATUM_ASSGN', 'POP_ESTN_UNIT', 'POP_EVAL',
'POP_STRATUM', 'POP_EVAL_TYP', 'POP_EVAL_GRP')
if (remote){
## Store the original parameters here
params <- db
## Read in one state at a time
db <- readFIA(dir = db$dir, common = db$common,
tables = reqTables, states = x, ## x is the vector of state names
nCores = nCores)
## If a clip was specified, run it now
if ('mostRecent' %in% names(params)){
db <- clipFIA(db, mostRecent = params$mostRecent,
mask = params$mask, matchEval = params$matchEval,
evalid = params$evalid, designCD = params$designCD,
nCores = nCores)
}
} else {
## Really only want the required tables
db <- db[names(db) %in% reqTables]
}
## Handle the modifier if it was given
db$TREE <- db$TREE %>%
mutate(state = !!stateVar,
grp = !!grpVar)
## Need a plotCN, and a new ID
db$PLOT <- db$PLOT %>% mutate(PLT_CN = CN,
pltID = paste(UNITCD, STATECD, COUNTYCD, PLOT, sep = '_'))
## don't have to change original code
#grpBy_quo <- enquo(grpBy)
# Probably cheating, but it works
if (quo_name(grpBy_quo) != 'NULL'){
## Have to join tables to run select with this object type
plt_quo <- filter(db$PLOT, !is.na(PLT_CN))
## We want a unique error message here to tell us when columns are not present in data
d_quo <- tryCatch(
error = function(cnd) {
return(0)
},
plt_quo[10,] %>% # Just the first row
left_join(select(db$COND, PLT_CN, names(db$COND)[names(db$COND) %in% names(db$PLOT) == FALSE]), by = 'PLT_CN') %>%
inner_join(select(db$TREE, PLT_CN, names(db$TREE)[names(db$TREE) %in% c(names(db$PLOT), names(db$COND)) == FALSE]), by = 'PLT_CN') %>%
select(!!grpBy_quo)
)
# If column doesnt exist, just returns 0, not a dataframe
if (is.null(nrow(d_quo))){
grpName <- quo_name(grpBy_quo)
stop(paste('Columns', grpName, 'not found in PLOT, TREE, or COND tables. Did you accidentally quote the variables names? e.g. use grpBy = ECOSUBCD (correct) instead of grpBy = "ECOSUBCD". ', collapse = ', '))
} else {
# Convert to character
grpBy <- names(d_quo)
}
} else {
grpBy <- NULL
}
reqTables <- c('PLOT', 'TREE', 'COND', 'POP_PLOT_STRATUM_ASSGN', 'POP_ESTN_UNIT', 'POP_EVAL',
'POP_STRATUM', 'POP_EVAL_TYP', 'POP_EVAL_GRP')
if (!is.null(polys) & first(class(polys)) %in% c('sf', 'SpatialPolygons', 'SpatialPolygonsDataFrame') == FALSE){
stop('polys must be spatial polygons object of class sp or sf. ')
}
if (landType %in% c('timber', 'forest') == FALSE){
stop('landType must be one of: "forest" or "timber".')
}
if (any(reqTables %in% names(db) == FALSE)){
missT <- reqTables[reqTables %in% names(db) == FALSE]
stop(paste('Tables', paste (as.character(missT), collapse = ', '), 'not found in object db.'))
}
if (str_to_upper(method) %in% c('TI', 'SMA', 'LMA', 'EMA', 'ANNUAL') == FALSE) {
warning(paste('Method', method, 'unknown. Defaulting to Temporally Indifferent (TI).'))
}
# I like a unique ID for a plot through time
if (byPlot) {grpBy <- c('pltID', grpBy)}
# Save original grpBy for pretty return with spatial objects
grpByOrig <- grpBy
## IF the object was clipped
if ('prev' %in% names(db$PLOT)){
## Only want the current plots, no grm
db$PLOT <- filter(db$PLOT, prev == 0)
}
### AREAL SUMMARY PREP
if(!is.null(polys)) {
# # Convert polygons to an sf object
# polys <- polys %>%
# as('sf')%>%
# mutate_if(is.factor,
# as.character)
# ## A unique ID
# polys$polyID <- 1:nrow(polys)
#
# # Add shapefile names to grpBy
grpBy = c(grpBy, 'polyID')
## Make plot data spatial, projected same as polygon layer
pltSF <- select(db$PLOT, c('LON', 'LAT', pltID)) %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
distinct(pltID, .keep_all = TRUE)
coordinates(pltSF) <- ~LON+LAT
proj4string(pltSF) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
pltSF <- as(pltSF, 'sf') %>%
st_transform(crs = st_crs(polys))
## Split up polys
polyList <- split(polys, as.factor(polys$polyID))
suppressWarnings({suppressMessages({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
cl <- makeCluster(nCores)
clusterEvalQ(cl, {
library(dplyr)
library(stringr)
library(rFIA)
})
out <- parLapply(cl, X = names(polyList), fun = areal_par, pltSF, polyList)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(polyList), FUN = areal_par, pltSF, polyList, mc.cores = nCores)
}
})})
pltSF <- bind_rows(out)
# A warning
if (length(unique(pltSF$pltID)) < 1){
stop('No plots in db overlap with polys.')
}
## Add polygon names to PLOT
db$PLOT <- db$PLOT %>%
left_join(select(pltSF, polyID, pltID), by = 'pltID')
# Test if any polygons cross state boundaries w/ different recent inventory years (continued w/in loop)
if ('mostRecent' %in% names(db) & length(unique(db$POP_EVAL$STATECD)) > 1){
mergeYears <- pltSF %>%
right_join(select(db$PLOT, PLT_CN, pltID), by = 'pltID') %>%
inner_join(select(db$POP_PLOT_STRATUM_ASSGN, c('PLT_CN', 'EVALID', 'STATECD')), by = 'PLT_CN') %>%
inner_join(select(db$POP_EVAL, c('EVALID', 'END_INVYR')), by = 'EVALID') %>%
group_by(polyID) %>%
summarize(maxYear = max(END_INVYR, na.rm = TRUE))
}
## TO RETURN SPATIAL PLOTS
}
if (byPlot & returnSpatial){
grpBy <- c(grpBy, 'LON', 'LAT')
} # END AREAL
## Build domain indicator function which is 1 if observation meets criteria, and 0 otherwise
# Land type domain indicator
if (tolower(landType) == 'forest'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1, 1, 0)
} else if (tolower(landType) == 'timber'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1 & db$COND$SITECLCD %in% c(1, 2, 3, 4, 5, 6) & db$COND$RESERVCD == 0, 1, 0)
}
# Tree Type domain indicator
if (tolower(treeType) == 'live'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 1, 1, 0)
} else if (tolower(treeType) == 'dead'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 2 & db$TREE$STANDING_DEAD_CD == 1, 1, 0)
} else if (tolower(treeType) == 'gs'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 1 & db$TREE$DIA >= 5 & db$TREE$TREECLCD == 2, 1, 0)
} else if (tolower(treeType) == 'all'){
db$TREE$typeD <- 1
}
# update spatial domain indicator
if(!is.null(polys)){
db$PLOT$sp <- ifelse(db$PLOT$pltID %in% pltSF$pltID, 1, 0)
} else {
db$PLOT$sp <- 1
}
# User defined domain indicator for area (ex. specific forest type)
pcEval <- left_join(db$PLOT, select(db$COND, -c('STATECD', 'UNITCD', 'COUNTYCD', 'INVYR', 'PLOT')), by = 'PLT_CN')
#areaDomain <- substitute(areaDomain)
pcEval$aD <- rlang::eval_tidy(areaDomain, pcEval) ## LOGICAL, THIS IS THE DOMAIN INDICATOR
if(!is.null(pcEval$aD)) pcEval$aD[is.na(pcEval$aD)] <- 0 # Make NAs 0s. Causes bugs otherwise
if(is.null(pcEval$aD)) pcEval$aD <- 1 # IF NULL IS GIVEN, THEN ALL VALUES TRUE
pcEval$aD <- as.numeric(pcEval$aD)
db$COND <- left_join(db$COND, select(pcEval, c('PLT_CN', 'CONDID', 'aD')), by = c('PLT_CN', 'CONDID')) %>%
mutate(aD_c = aD)
aD_p <- pcEval %>%
group_by(PLT_CN) %>%
summarize(aD_p = as.numeric(any(aD > 0)))
db$PLOT <- left_join(db$PLOT, aD_p, by = 'PLT_CN')
rm(pcEval)
# Same as above for tree (ex. trees > 20 ft tall)
#treeDomain <- substitute(treeDomain)
tD <- rlang::eval_tidy(treeDomain, db$TREE) ## LOGICAL, THIS IS THE DOMAIN INDICATOR
if(!is.null(tD)) tD[is.na(tD)] <- 0 # Make NAs 0s. Causes bugs otherwise
if(is.null(tD)) tD <- 1 # IF NULL IS GIVEN, THEN ALL VALUES TRUE
db$TREE$tD <- as.numeric(tD)
### Snag the EVALIDs that are needed
db$POP_EVAL<- db$POP_EVAL %>%
select('CN', 'END_INVYR', 'EVALID', 'ESTN_METHOD', STATECD) %>%
inner_join(select(db$POP_EVAL_TYP, c('EVAL_CN', 'EVAL_TYP')), by = c('CN' = 'EVAL_CN')) %>%
filter(EVAL_TYP == 'EXPVOL' | EVAL_TYP == 'EXPCURR') %>%
filter(!is.na(END_INVYR) & !is.na(EVALID) & END_INVYR >= 2003) %>%
distinct(END_INVYR, EVALID, .keep_all = TRUE)# %>%
#group_by(END_INVYR) %>%
#summarise(id = list(EVALID)
## If a most-recent subset, make sure that we don't get two reporting years in
## western states
if (mr) {
db$POP_EVAL <- db$POP_EVAL %>%
group_by(EVAL_TYP, STATECD) %>%
filter(END_INVYR == max(END_INVYR, na.rm = TRUE)) %>%
ungroup()
}
## Cut STATECD
db$POP_EVAL <- select(db$POP_EVAL, -c(STATECD))
### The population tables
pops <- select(db$POP_EVAL, c('EVALID', 'ESTN_METHOD', 'CN', 'END_INVYR', 'EVAL_TYP')) %>%
rename(EVAL_CN = CN) %>%
left_join(select(db$POP_ESTN_UNIT, c('CN', 'EVAL_CN', 'AREA_USED', 'P1PNTCNT_EU')), by = c('EVAL_CN')) %>%
rename(ESTN_UNIT_CN = CN) %>%
left_join(select(db$POP_STRATUM, c('ESTN_UNIT_CN', 'EXPNS', 'P2POINTCNT', 'CN', 'P1POINTCNT', 'ADJ_FACTOR_SUBP', 'ADJ_FACTOR_MICR', "ADJ_FACTOR_MACR")), by = c('ESTN_UNIT_CN')) %>%
rename(STRATUM_CN = CN) %>%
left_join(select(db$POP_PLOT_STRATUM_ASSGN, c('STRATUM_CN', 'PLT_CN', 'INVYR', 'STATECD')), by = 'STRATUM_CN') %>%
mutate_if(is.factor,
as.character)
### Which estimator to use?
if (str_to_upper(method) %in% c('ANNUAL')){
## Want to use the year where plots are measured, no repeats
## Breaking this up into pre and post reporting becuase
## Estimation units get weird on us otherwise
popOrig <- pops
pops <- pops %>%
group_by(STATECD) %>%
filter(END_INVYR == INVYR) %>%
ungroup()
prePops <- popOrig %>%
group_by(STATECD) %>%
filter(INVYR < min(END_INVYR, na.rm = TRUE)) %>%
distinct(PLT_CN, .keep_all = TRUE) %>%
ungroup()
pops <- bind_rows(pops, prePops) %>%
mutate(YEAR = INVYR)
} else { # Otherwise temporally indifferent
pops <- rename(pops, YEAR = END_INVYR)
}
## P2POINTCNT column is NOT consistent for annnual estimates, plots
## within individual strata and est units are related to different INVYRs
p2_INVYR <- pops %>%
group_by(ESTN_UNIT_CN, STRATUM_CN, INVYR) %>%
summarize(P2POINTCNT_INVYR = length(unique(PLT_CN)))
## Want a count of p2 points / eu, gets screwed up with grouping below
p2eu_INVYR <- p2_INVYR %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, INVYR, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN, INVYR) %>%
summarize(p2eu_INVYR = sum(P2POINTCNT_INVYR, na.rm = TRUE))
p2eu <- pops %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(p2eu = sum(P2POINTCNT, na.rm = TRUE))
## Rejoin
pops <- pops %>%
left_join(p2_INVYR, by = c('ESTN_UNIT_CN', 'STRATUM_CN', 'INVYR')) %>%
left_join(p2eu_INVYR, by = c('ESTN_UNIT_CN', 'INVYR')) %>%
left_join(p2eu, by = 'ESTN_UNIT_CN')
## Recode a few of the estimation methods to make things easier below
pops$ESTN_METHOD = recode(.x = pops$ESTN_METHOD,
`Post-Stratification` = 'strat',
`Stratified random sampling` = 'strat',
`Double sampling for stratification` = 'double',
`Simple random sampling` = 'simple',
`Subsampling units of unequal size` = 'simple')
## Break into size classes
if (bySizeClass){
grpBy <- c(grpBy, 'sizeClass')
grpByOrig <- c(grpByOrig, 'sizeClass')
db$TREE$sizeClass <- makeClasses(db$TREE$DIA, interval = 2, numLabs = TRUE)
db$TREE <- db$TREE[!is.na(db$TREE$sizeClass),]
}
## Only the necessary plots for EVAL of interest
db$PLOT <- filter(db$PLOT, PLT_CN %in% pops$PLT_CN)
## Which grpByNames are in which table? Helps us subset below
grpP <- names(db$PLOT)[names(db$PLOT) %in% grpBy]
grpC <- names(db$COND)[names(db$COND) %in% grpBy & names(db$COND) %in% grpP == FALSE]
grpT <- names(db$TREE)[names(db$TREE) %in% grpBy & names(db$TREE) %in% c(grpP, grpC) == FALSE]
db$PLOT <- select(db$PLOT, c('PLT_CN', 'STATECD', 'COUNTYCD', 'MACRO_BREAKPOINT_DIA', 'INVYR', 'MEASYEAR', 'PLOT_STATUS_CD', all_of(grpP), 'aD_p', 'sp'))
db$COND <- select(db$COND, c('PLT_CN', 'CONDPROP_UNADJ', 'PROP_BASIS', 'COND_STATUS_CD', 'CONDID', all_of(grpC), 'aD_c', 'landD'))
db$TREE <- db$TREE %>%
mutate(TRE_CN = CN) %>%
select(c('PLT_CN', 'TRE_CN', 'CONDID', 'DIA', 'grp', 'state', 'SUBP', 'TREE', all_of(grpT), 'tD', 'typeD'))
## Merging state and county codes
plts <- split(db$PLOT, as.factor(paste(db$PLOT$COUNTYCD, db$PLOT$STATECD, sep = '_')))
#plts <- split(db$PLOT, as.factor(db$PLOT$STATECD))
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
cl <- makeCluster(nCores)
clusterEvalQ(cl, {
library(dplyr)
library(stringr)
library(rFIA)
})
out <- parLapply(cl, X = names(plts), fun = divHelper1, plts, db, grpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = divHelper1, plts, db, grpBy, byPlot, mc.cores = nCores)
}
})
if (byPlot){
## back to dataframes
out <- unlist(out, recursive = FALSE)
tOut <- bind_rows(out[names(out) == 't'])
## Make it spatial
if (returnSpatial){
tOut <- tOut %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
st_as_sf(coords = c('LON', 'LAT'),
crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
}
out <- list(tEst = tOut, grpBy = grpBy, grpByOrig = grpByOrig)
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
t <- bind_rows(out[names(out) == 't'])
full <- bind_rows(out[names(out) == 'full'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
## Splitting up by STATECD and groups of 25 ESTN_UNIT_CNs
#estunit <- distinct(pops, ESTN_UNIT_CN) #%>%
#mutate(estID)
#estID <- seq(1, nrow(estunit), 50)
#estunit$estID <- rep_len(estID, length.out = nrow(estunit))
#pops <- pops %>%
# left_join(estunit, by = 'ESTN_UNIT_CN') #%>%
#mutate(estBreaks = )
#popState <- split(pops, as.factor(pops$estID))
popState <- split(pops, as.factor(pops$STATECD))
#
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
## Use the same cluster as above
# cl <- makeCluster(nCores)
# clusterEvalQ(cl, {
# library(dplyr)
# library(stringr)
# library(rFIA)
# })
out <- parLapply(cl, X = names(popState), fun = divHelper2, popState, t, full, grpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = divHelper2, popState, t, full, grpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
full <- bind_rows(out[names(out) == 'full'])
##### ----------------- MOVING AVERAGES
if (str_to_upper(method) %in% c("SMA", 'EMA', 'LMA')){
### ---- SIMPLE MOVING AVERAGE
if (str_to_upper(method) == 'SMA'){
## Assuming a uniform weighting scheme
wgts <- pops %>%
group_by(ESTN_UNIT_CN) %>%
summarize(wgt = 1 / length(unique(INVYR)))
#aEst <- left_join(aEst, wgts, by = 'ESTN_UNIT_CN')
tEst <- left_join(tEst, wgts, by = 'ESTN_UNIT_CN')
#### ----- Linear MOVING AVERAGE
} else if (str_to_upper(method) == 'LMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
## Want a number of INVYRs per EU
neu <- wgts %>%
group_by(ESTN_UNIT_CN) %>%
summarize(n = sum(rank, na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = rank / n) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
#### ----- EXPONENTIAL MOVING AVERAGE
} else if (str_to_upper(method) == 'EMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
if (length(lambda) < 2){
## Want sum of weighitng functions
neu <- wgts %>%
mutate(l = lambda) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
} else {
grpBy <- c('lambda', grpBy)
#aGrpBy <- c('lambda', aGrpBy)
## Duplicate weights for each level of lambda
yrWgts <- list()
for (i in 1:length(unique(lambda))) {
yrWgts[[i]] <- mutate(wgts, lambda = lambda[i])
}
wgts <- bind_rows(yrWgts)
## Want sum of weighitng functions
neu <- wgts %>%
group_by(lambda, ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = c('lambda', 'ESTN_UNIT_CN')) %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(lambda, ESTN_UNIT_CN, INVYR, wgt)
}
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
}
### Applying the weights
tEst <- tEst %>%
mutate_at(vars(aEst:sEst), ~(.*wgt)) %>%
mutate_at(vars(aVar:cvEst_s), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
}
out <- list(tEst = tEst, full = full, grpBy = grpBy, grpByOrig = grpByOrig)
}
return(out)
}
#' @export
diversity <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
stateVar = TPA_UNADJ,
grpVar = SPCD,
treeDomain = NULL,
areaDomain = NULL,
byPlot = FALSE,
totals = FALSE,
variance = FALSE,
nCores = 1) {
## don't have to change original code
grpBy_quo <- rlang::enquo(grpBy)
areaDomain <- rlang::enquo(areaDomain)
treeDomain <- rlang::enquo(treeDomain)
stateVar <- rlang::enquo(stateVar)
grpVar <- rlang::enquo(grpVar)
### Is DB remote?
remote <- ifelse(class(db) == 'Remote.FIA.Database', 1, 0)
if (remote){
iter <- db$states
## In memory
} else {
## Some warnings
if (class(db) != "FIA.Database"){
stop('db must be of class "FIA.Database". Use readFIA() to load your FIA data.')
}
## an iterator for remote
iter <- 1
}
## Check for a most recent subset
if (remote){
if ('mostRecent' %in% names(db)){
mr = db$mostRecent # logical
} else {
mr = FALSE
}
## In-memory
} else {
if ('mostRecent' %in% names(db)){
mr = TRUE
} else {
mr = FALSE
}
}
### AREAL SUMMARY PREP
if(!is.null(polys)) {
# Convert polygons to an sf object
polys <- polys %>%
as('sf')%>%
mutate_if(is.factor,
as.character)
## A unique ID
polys$polyID <- 1:nrow(polys)
}
## Run the main portion
out <- lapply(X = iter, FUN = diversityStarter, db,
grpBy_quo = grpBy_quo, polys, returnSpatial,
bySizeClass,
landType, treeType, method,
lambda, stateVar, grpVar,
treeDomain, areaDomain,
byPlot, totals, nCores, remote, mr)
## Bring the results back
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
full <- bind_rows(out[names(out) == 'full'])
grpBy <- out[names(out) == 'grpBy'][[1]]
grpByOrig <- out[names(out) == 'grpByOrig'][[1]]
if (byPlot){
tOut <- tEst
} else {
suppressMessages({suppressWarnings({
## If a clip was specified, handle the reporting years
if (mr){
## If a most recent subset, ignore differences in reporting years across states
## instead combine most recent information from each state
# ID mr years by group
maxyearsT <- tEst %>%
select(grpBy) %>%
group_by(.dots = grpBy[!c(grpBy %in% 'YEAR')]) %>%
summarise(YEAR = max(YEAR, na.rm = TRUE))
# Combine estimates
tEst <- tEst %>%
ungroup() %>%
select(-c(YEAR)) %>%
left_join(maxyearsT, by = grpBy[!c(grpBy %in% 'YEAR')])
}
})})
##--------------------- TOTALS and RATIOS
suppressWarnings({
tOut <- tEst %>%
group_by(.dots = grpBy) %>%
summarize_all(sum,na.rm = TRUE) %>%
mutate(H_a = hEst / aEst,
Eh_a = ehEst / aEst,
S_a = sEst / aEst,
AREA_TOTAL = aEst,
## Ratio variance
haVar = (1/AREA_TOTAL^2) * (hVar + (H_a^2 * aVar) - 2 * H_a * cvEst_h),
ehaVar = (1/AREA_TOTAL^2) * (ehVar + (Eh_a^2 * aVar) - (2 * Eh_a * cvEst_eh)),
saVar = (1/AREA_TOTAL^2) * (sVar + (S_a^2 * aVar) - 2 * S_a * cvEst_s),
H_a_VAR = haVar,
Eh_a_VAR = ehaVar,
S_a_VAR = saVar,
## RATIO SE
H_a_SE = sqrt(haVar) / H_a * 100,
Eh_a_SE = sqrt(ehaVar) / Eh_a * 100,
S_a_SE = sqrt(saVar) / S_a * 100,
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL * 100,
AREA_TOTAL_VAR = aVar,
nPlots = plotIn_AREA)
})
### Up a few spatial scales
full <- full %>%
group_by(.dots = grpBy[!c(grpBy %in% 'lambda')]) %>%
summarize(H_g = divIndex(grp, state, index = 'H'),
Eh_g = divIndex(grp, state, index = 'Eh'),
S_g = divIndex(grp, state, index = 'S'))
tOut <- tOut %>%
left_join(full, by = grpBy[!c(grpBy %in% 'lambda')]) %>%
mutate(H_b = H_g - H_a,
Eh_b = Eh_g - Eh_a,
S_b = S_g - S_a)
if (totals) {
if (variance){
tOut <- tOut %>%
select(grpBy, H_a, H_b, H_g, Eh_a, Eh_b, Eh_g, S_a, S_b, S_g, AREA_TOTAL, H_a_VAR, Eh_a_VAR, S_a_VAR, AREA_TOTAL_VAR, nPlots, N)
} else {
tOut <- tOut %>%
select(grpBy, H_a, H_b, H_g, Eh_a, Eh_b, Eh_g, S_a, S_b, S_g, AREA_TOTAL, H_a_SE, Eh_a_SE, S_a_SE, AREA_TOTAL_SE, nPlots)
}
} else {
if (variance){
tOut <- tOut %>%
select(grpBy, H_a, H_b, H_g, Eh_a, Eh_b, Eh_g, S_a, S_b, S_g, H_a_VAR, Eh_a_VAR, S_a_VAR, nPlots, N)
} else {
tOut <- tOut %>%
select(grpBy, H_a, H_b, H_g, Eh_a, Eh_b, Eh_g, S_a, S_b, S_g, H_a_SE, Eh_a_SE, S_a_SE, nPlots)
}
}
# Snag the names
tNames <- names(tOut)[names(tOut) %in% grpBy == FALSE]
}
## Pretty output
tOut <- tOut %>%
ungroup() %>%
mutate_if(is.factor, as.character) %>%
drop_na(grpBy) %>%
arrange(YEAR) %>%
as_tibble()
# Return a spatial object
if (!is.null(polys) & byPlot == FALSE) {
## NO IMPLICIT NA
nospGrp <- unique(grpBy[grpBy %in% c('SPCD', 'SYMBOL', 'COMMON_NAME', 'SCIENTIFIC_NAME') == FALSE])
nospSym <- syms(nospGrp)
tOut <- complete(tOut, !!!nospSym)
## If species, we don't want unique combos of variables related to same species
## but we do want NAs in polys where species are present
if (length(nospGrp) < length(grpBy)){
spGrp <- unique(grpBy[grpBy %in% c('SPCD', 'SYMBOL', 'COMMON_NAME', 'SCIENTIFIC_NAME')])
spSym <- syms(spGrp)
tOut <- complete(tOut, nesting(!!!nospSym))
}
suppressMessages({suppressWarnings({tOut <- left_join(tOut, polys) %>%
select(c('YEAR', grpByOrig, tNames, names(polys))) %>%
filter(!is.na(polyID) & !is.na(nPlots))})})
## Makes it horrible to work with as a dataframe
if (returnSpatial == FALSE) tOut <- select(tOut, -c(geometry))
} else if (!is.null(polys) & byPlot){
polys <- as.data.frame(polys)
tOut <- left_join(tOut, select(polys, -c(geometry)), by = 'polyID')
}
## For spatial plots
if (returnSpatial & byPlot) grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
## Above converts to tibble
if (returnSpatial) tOut <- st_sf(tOut)
# ## remove any duplicates in byPlot (artifact of END_INYR loop)
if (byPlot) tOut <- unique(tOut)
return(tOut)
}
diversity_backup <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
stateVar = TPA_UNADJ,
grpVar = SPCD,
treeDomain = NULL,
areaDomain = NULL,
byPlot = FALSE,
totals = FALSE,
nCores = 1) {
## Need a plotCN, and a new ID
db$PLOT <- db$PLOT %>% mutate(PLT_CN = CN,
pltID = paste(UNITCD, STATECD, COUNTYCD, PLOT, sep = '_'))
## Converting names given in grpBy to character vector (NSE to standard)
## don't have to change original code
grpBy_quo <- enquo(grpBy)
# Probably cheating, but it works
if (quo_name(grpBy_quo) != 'NULL'){
## Have to join tables to run select with this object type
plt_quo <- filter(db$PLOT, !is.na(PLT_CN))
## We want a unique error message here to tell us when columns are not present in data
d_quo <- tryCatch(
error = function(cnd) {
return(0)
},
plt_quo[10,] %>% # Just the first row
left_join(select(db$COND, PLT_CN, names(db$COND)[names(db$COND) %in% names(db$PLOT) == FALSE]), by = 'PLT_CN') %>%
inner_join(select(db$TREE, PLT_CN, names(db$TREE)[names(db$TREE) %in% c(names(db$PLOT), names(db$COND)) == FALSE]), by = 'PLT_CN') %>%
select(!!grpBy_quo)
)
# If column doesnt exist, just returns 0, not a dataframe
if (is.null(nrow(d_quo))){
grpName <- quo_name(grpBy_quo)
stop(paste('Columns', grpName, 'not found in PLOT, TREE, or COND tables. Did you accidentally quote the variables names? e.g. use grpBy = ECOSUBCD (correct) instead of grpBy = "ECOSUBCD". ', collapse = ', '))
} else {
# Convert to character
grpBy <- names(d_quo)
}
}
reqTables <- c('PLOT', 'TREE', 'COND', 'POP_PLOT_STRATUM_ASSGN', 'POP_ESTN_UNIT', 'POP_EVAL',
'POP_STRATUM', 'POP_EVAL_TYP', 'POP_EVAL_GRP')
## Some warnings
if (class(db) != "FIA.Database"){
stop('db must be of class "FIA.Database". Use readFIA() to load your FIA data.')
}
if (!is.null(polys) & first(class(polys)) %in% c('sf', 'SpatialPolygons', 'SpatialPolygonsDataFrame') == FALSE){
stop('polys must be spatial polygons object of class sp or sf. ')
}
if (landType %in% c('timber', 'forest') == FALSE){
stop('landType must be one of: "forest" or "timber".')
}
if (any(reqTables %in% names(db) == FALSE)){
missT <- reqTables[reqTables %in% names(db) == FALSE]
stop(paste('Tables', paste (as.character(missT), collapse = ', '), 'not found in object db.'))
}
if (str_to_upper(method) %in% c('TI', 'SMA', 'LMA', 'EMA', 'ANNUAL') == FALSE) {
warning(paste('Method', method, 'unknown. Defaulting to Temporally Indifferent (TI).'))
}
##### CHECK TO SEE IF STATE VARIABLE AND GROUP VARIABLE EXIST IN TREE
# Need to quote all variables for NSE
stateVar = enquo(stateVar)
grpVar = enquo(grpVar)
## If a modifier was given to a variable, handle it (ex. y = TPA_PERC / 100)
db$TREE <- db$TREE %>%
mutate(state = !!stateVar,
grp = !!grpVar)
# if (quo_name(stateVar) %in% names(db$TREE) == FALSE){
# stop(paste('Column',quo_name(stateVar), 'not found in TREE table. Did you accidentally quote the variables name? e.g. use stateVar = TPA_UNADJ (correct) instead of grpBy = "TPA_UNADJ". ', collapse = ', '))
# }
# if (quo_name(grpVar) %in% names(db$TREE) == FALSE){
# stop(paste('Column',quo_name(grpVar), 'not found in TREE table. Did you accidentally quote the variables name? e.g. use stateVar = TPA_UNADJ (correct) instead of grpBy = "TPA_UNADJ". ', collapse = ', '))
# }
# I like a unique ID for a plot through time
if (byPlot) {grpBy <- c('pltID', grpBy)}
# Save original grpBy for pretty return with spatial objects
grpByOrig <- grpBy
## IF the object was clipped
if ('prev' %in% names(db$PLOT)){
## Only want the current plots, no grm
db$PLOT <- filter(db$PLOT, prev == 0)
}
### DEAL WITH TEXAS
if (any(db$POP_EVAL$STATECD %in% 48)){
## Will require manual updates, fix your shit texas
txIDS <- db$POP_EVAL %>%
filter(STATECD %in% 48) %>%
filter(END_INVYR < 2017) %>%
filter(END_INVYR > 2006) %>%
## Removing any inventory that references east or west, sorry
filter(str_detect(str_to_upper(EVAL_DESCR), 'EAST', negate = TRUE) &
str_detect(str_to_upper(EVAL_DESCR), 'WEST', negate = TRUE))
db$POP_EVAL <- bind_rows(filter(db$POP_EVAL, !(STATECD %in% 48)), txIDS)
}
### AREAL SUMMARY PREP
if(!is.null(polys)) {
# Convert polygons to an sf object
polys <- polys %>%
as('sf')%>%
mutate_if(is.factor,
as.character)
## A unique ID
polys$polyID <- 1:nrow(polys)
# Add shapefile names to grpBy
#grpBy = c(names(polys)[str_detect(names(polys), 'geometry') == FALSE], 'polyID', grpBy)
grpBy = c(grpBy, 'polyID')
## Make plot data spatial, projected same as polygon layer
pltSF <- select(db$PLOT, c('LON', 'LAT', pltID)) %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
distinct(pltID, .keep_all = TRUE)
coordinates(pltSF) <- ~LON+LAT
proj4string(pltSF) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
pltSF <- as(pltSF, 'sf') %>%
st_transform(crs = st_crs(polys))
## Split up polys
polyList <- split(polys, as.factor(polys$polyID))
suppressWarnings({suppressMessages({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
cl <- makeCluster(nCores)
clusterEvalQ(cl, {
library(dplyr)
library(stringr)
library(rFIA)
})
out <- parLapply(cl, X = names(polyList), fun = areal_par, pltSF, polyList)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(polyList), FUN = areal_par, pltSF, polyList, mc.cores = nCores)
}
})})
pltSF <- bind_rows(out)
# A warning
if (length(unique(pltSF$pltID)) < 1){
stop('No plots in db overlap with polys.')
}
## Add polygon names to PLOT
db$PLOT <- db$PLOT %>%
left_join(select(pltSF, polyID, pltID), by = 'pltID')
# Test if any polygons cross state boundaries w/ different recent inventory years (continued w/in loop)
if ('mostRecent' %in% names(db) & length(unique(db$POP_EVAL$STATECD)) > 1){
mergeYears <- pltSF %>%
right_join(select(db$PLOT, PLT_CN, pltID), by = 'pltID') %>%
inner_join(select(db$POP_PLOT_STRATUM_ASSGN, c('PLT_CN', 'EVALID', 'STATECD')), by = 'PLT_CN') %>%
inner_join(select(db$POP_EVAL, c('EVALID', 'END_INVYR')), by = 'EVALID') %>%
group_by(polyID) %>%
summarize(maxYear = max(END_INVYR, na.rm = TRUE))
}
## TO RETURN SPATIAL PLOTS
}
if (byPlot & returnSpatial){
grpBy <- c(grpBy, 'LON', 'LAT')
} # END AREAL
## Build domain indicator function which is 1 if observation meets criteria, and 0 otherwise
# Land type domain indicator
if (tolower(landType) == 'forest'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1, 1, 0)
} else if (tolower(landType) == 'timber'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1 & db$COND$SITECLCD %in% c(1, 2, 3, 4, 5, 6) & db$COND$RESERVCD == 0, 1, 0)
}
# Tree Type domain indicator
if (tolower(treeType) == 'live'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 1, 1, 0)
} else if (tolower(treeType) == 'dead'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 2 & db$TREE$STANDING_DEAD_CD == 1, 1, 0)
} else if (tolower(treeType) == 'gs'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 1 & db$TREE$DIA >= 5 & db$TREE$TREECLCD == 2, 1, 0)
} else if (tolower(treeType) == 'all'){
db$TREE$typeD <- 1
}
# update spatial domain indicator
if(!is.null(polys)){
db$PLOT$sp <- ifelse(db$PLOT$pltID %in% pltSF$pltID, 1, 0)
} else {
db$PLOT$sp <- 1
}
# User defined domain indicator for area (ex. specific forest type)
pcEval <- left_join(db$PLOT, select(db$COND, -c('STATECD', 'UNITCD', 'COUNTYCD', 'INVYR', 'PLOT')), by = 'PLT_CN')
areaDomain <- substitute(areaDomain)
pcEval$aD <- eval(areaDomain, pcEval) ## LOGICAL, THIS IS THE DOMAIN INDICATOR
if(!is.null(pcEval$aD)) pcEval$aD[is.na(pcEval$aD)] <- 0 # Make NAs 0s. Causes bugs otherwise
if(is.null(pcEval$aD)) pcEval$aD <- 1 # IF NULL IS GIVEN, THEN ALL VALUES TRUE
pcEval$aD <- as.numeric(pcEval$aD)
db$COND <- left_join(db$COND, select(pcEval, c('PLT_CN', 'CONDID', 'aD')), by = c('PLT_CN', 'CONDID')) %>%
mutate(aD_c = aD)
aD_p <- pcEval %>%
group_by(PLT_CN) %>%
summarize(aD_p = as.numeric(any(aD > 0)))
db$PLOT <- left_join(db$PLOT, aD_p, by = 'PLT_CN')
rm(pcEval)
# Same as above for tree (ex. trees > 20 ft tall)
treeDomain <- substitute(treeDomain)
tD <- eval(treeDomain, db$TREE) ## LOGICAL, THIS IS THE DOMAIN INDICATOR
if(!is.null(tD)) tD[is.na(tD)] <- 0 # Make NAs 0s. Causes bugs otherwise
if(is.null(tD)) tD <- 1 # IF NULL IS GIVEN, THEN ALL VALUES TRUE
db$TREE$tD <- as.numeric(tD)
### Snag the EVALIDs that are needed
db$POP_EVAL<- db$POP_EVAL %>%
select('CN', 'END_INVYR', 'EVALID', 'ESTN_METHOD') %>%
inner_join(select(db$POP_EVAL_TYP, c('EVAL_CN', 'EVAL_TYP')), by = c('CN' = 'EVAL_CN')) %>%
filter(EVAL_TYP == 'EXPVOL' | EVAL_TYP == 'EXPCURR') %>%
filter(!is.na(END_INVYR) & !is.na(EVALID) & END_INVYR >= 2003) %>%
distinct(END_INVYR, EVALID, .keep_all = TRUE)# %>%
#group_by(END_INVYR) %>%
#summarise(id = list(EVALID)
## Make an annual panel ID, associated with an INVYR
### The population tables
pops <- select(db$POP_EVAL, c('EVALID', 'ESTN_METHOD', 'CN', 'END_INVYR')) %>%
rename(EVAL_CN = CN) %>%
left_join(select(db$POP_ESTN_UNIT, c('CN', 'EVAL_CN', 'AREA_USED', 'P1PNTCNT_EU')), by = c('EVAL_CN')) %>%
rename(ESTN_UNIT_CN = CN) %>%
left_join(select(db$POP_STRATUM, c('ESTN_UNIT_CN', 'EXPNS', 'P2POINTCNT', 'CN', 'P1POINTCNT', 'ADJ_FACTOR_SUBP', 'ADJ_FACTOR_MICR', "ADJ_FACTOR_MACR")), by = c('ESTN_UNIT_CN')) %>%
rename(STRATUM_CN = CN) %>%
left_join(select(db$POP_PLOT_STRATUM_ASSGN, c('STRATUM_CN', 'PLT_CN', 'INVYR', 'STATECD')), by = 'STRATUM_CN') %>%
ungroup() %>%
mutate_if(is.factor,
as.character)
### Which estimator to use?
if (str_to_upper(method) %in% c('ANNUAL')){
## Want to use the year where plots are measured, no repeats
## Breaking this up into pre and post reporting becuase
## Estimation units get weird on us otherwise
popOrig <- pops
pops <- pops %>%
group_by(STATECD) %>%
filter(END_INVYR == INVYR) %>%
ungroup()
prePops <- popOrig %>%
group_by(STATECD) %>%
filter(INVYR < min(END_INVYR, na.rm = TRUE)) %>%
distinct(PLT_CN, .keep_all = TRUE) %>%
ungroup()
pops <- bind_rows(pops, prePops) %>%
mutate(YEAR = INVYR)
} else { # Otherwise temporally indifferent
pops <- rename(pops, YEAR = END_INVYR)
}
## P2POINTCNT column is NOT consistent for annnual estimates, plots
## within individual strata and est units are related to different INVYRs
p2_INVYR <- pops %>%
group_by(ESTN_UNIT_CN, STRATUM_CN, INVYR) %>%
summarize(P2POINTCNT_INVYR = length(unique(PLT_CN)))
## Want a count of p2 points / eu, gets screwed up with grouping below
p2eu_INVYR <- p2_INVYR %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, INVYR, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN, INVYR) %>%
summarize(p2eu_INVYR = sum(P2POINTCNT_INVYR, na.rm = TRUE))
p2eu <- pops %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(p2eu = sum(P2POINTCNT, na.rm = TRUE))
## Rejoin
pops <- pops %>%
left_join(p2_INVYR, by = c('ESTN_UNIT_CN', 'STRATUM_CN', 'INVYR')) %>%
left_join(p2eu_INVYR, by = c('ESTN_UNIT_CN', 'INVYR')) %>%
left_join(p2eu, by = 'ESTN_UNIT_CN')
## Recode a few of the estimation methods to make things easier below
pops$ESTN_METHOD = recode(.x = pops$ESTN_METHOD,
`Post-Stratification` = 'strat',
`Stratified random sampling` = 'strat',
`Double sampling for stratification` = 'double',
`Simple random sampling` = 'simple',
`Subsampling units of unequal size` = 'simple')
## Break into size classes
if (bySizeClass){
grpBy <- c(grpBy, 'sizeClass')
grpByOrig <- c(grpByOrig, 'sizeClass')
db$TREE$sizeClass <- makeClasses(db$TREE$DIA, interval = 2, numLabs = TRUE)
db$TREE <- db$TREE[!is.na(db$TREE$sizeClass),]
}
## Only the necessary plots for EVAL of interest
db$PLOT <- filter(db$PLOT, PLT_CN %in% pops$PLT_CN)
## Merging state and county codes
plts <- split(db$PLOT, as.factor(paste(db$PLOT$COUNTYCD, db$PLOT$STATECD, sep = '_')))
#plts <- split(db$PLOT, as.factor(db$PLOT$STATECD))
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
cl <- makeCluster(nCores)
clusterEvalQ(cl, {
library(dplyr)
library(stringr)
library(rFIA)
})
out <- parLapply(cl, X = names(plts), fun = divHelper1, plts, db, grpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = divHelper1, plts, db, grpBy, byPlot, mc.cores = nCores)
}
})
if (byPlot){
## back to dataframes
out <- unlist(out, recursive = FALSE)
tOut <- bind_rows(out[names(out) == 't'])
## Make it spatial
if (returnSpatial){
tOut <- tOut %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
st_as_sf(coords = c('LON', 'LAT'),
crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
}
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
## Splitting up by ESTN_UNIT
popState <- split(pops, as.factor(pops$STATECD))
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
## Use the same cluster as above
# cl <- makeCluster(nCores)
# clusterEvalQ(cl, {
# library(dplyr)
# library(stringr)
# library(rFIA)
# })
out <- parLapply(cl, X = names(popState), fun = divHelper2, popState, t, grpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = divHelper2, popState, t, grpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
##### ----------------- MOVING AVERAGES
if (str_to_upper(method) %in% c("SMA", 'EMA', 'LMA')){
### ---- SIMPLE MOVING AVERAGE
if (str_to_upper(method) == 'SMA'){
## Assuming a uniform weighting scheme
wgts <- pops %>%
group_by(ESTN_UNIT_CN) %>%
summarize(wgt = 1 / length(unique(INVYR)))
#aEst <- left_join(aEst, wgts, by = 'ESTN_UNIT_CN')
tEst <- left_join(tEst, wgts, by = 'ESTN_UNIT_CN')
#### ----- Linear MOVING AVERAGE
} else if (str_to_upper(method) == 'LMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
## Want a number of INVYRs per EU
neu <- wgts %>%
group_by(ESTN_UNIT_CN) %>%
summarize(n = sum(rank, na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = rank / n) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
#### ----- EXPONENTIAL MOVING AVERAGE
} else if (str_to_upper(method) == 'EMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
if (length(lambda) < 2){
## Want sum of weighitng functions
neu <- wgts %>%
mutate(l = lambda) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
} else {
grpBy <- c('lambda', grpBy)
#aGrpBy <- c('lambda', aGrpBy)
## Duplicate weights for each level of lambda
yrWgts <- list()
for (i in 1:length(unique(lambda))) {
yrWgts[[i]] <- mutate(wgts, lambda = lambda[i])
}
wgts <- bind_rows(yrWgts)
## Want sum of weighitng functions
neu <- wgts %>%
group_by(lambda, ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = c('lambda', 'ESTN_UNIT_CN')) %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(lambda, ESTN_UNIT_CN, INVYR, wgt)
}
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
}
### Applying the weights
tEst <- tEst %>%
mutate_at(vars(aEst:sEst), ~(.*wgt)) %>%
mutate_at(vars(aVar:cvEst_s), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
}
##--------------------- TOTALS and RATIOS
suppressWarnings({
tOut <- tEst %>%
group_by(.dots = grpBy) %>%
summarize_all(sum,na.rm = TRUE) %>%
mutate(H_a = hEst / aEst,
Eh_a = ehEst / aEst,
S_a = sEst / aEst,
AREA_TOTAL = aEst,
## Ratio variance
haVar = (1/AREA_TOTAL^2) * (hVar + (H_a^2 * aVar) - 2 * H_a * cvEst_h),
ehaVar = (1/AREA_TOTAL^2) * (ehVar + (Eh_a^2 * aVar) - (2 * Eh_a * cvEst_eh)),
saVar = (1/AREA_TOTAL^2) * (sVar + (S_a^2 * aVar) - 2 * S_a * cvEst_s),
## RATIO SE
H_a_SE = sqrt(haVar) / H_a * 100,
Eh_a_SE = sqrt(ehaVar) / Eh_a * 100,
S_a_SE = sqrt(saVar) / S_a * 100,
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL * 100,
nPlots = plotIn_AREA)
})
### UP a few spatial scales
## Which grpByNames are in which table? Helps us subset below
grpP <- names(db$PLOT)[names(db$PLOT) %in% grpBy]
grpC <- names(db$COND)[names(db$COND) %in% grpBy & names(db$COND) %in% grpP == FALSE]
grpT <- names(db$TREE)[names(db$TREE) %in% grpBy & names(db$TREE) %in% c(grpP, grpC) == FALSE]
### REALLY DON'T LIKE WORKING WITH FULL TABLES, WORK ON A FIX
### Only joining tables necessary to produce plot level estimates, adjusted for non-response
fullArea <- select(db$PLOT, c('PLT_CN', 'STATECD', 'MACRO_BREAKPOINT_DIA', 'INVYR', 'MEASYEAR', 'PLOT_STATUS_CD', grpP, 'aD_p', 'sp')) %>%
left_join(select(db$COND, c('PLT_CN', 'CONDPROP_UNADJ', 'PROP_BASIS', 'COND_STATUS_CD', 'CONDID', grpC, 'aD_c', 'landD')), by = c('PLT_CN')) %>%
left_join(select(db$TREE, c('PLT_CN', 'CONDID', 'DIA', 'grp', 'state', 'SUBP', 'TREE', grpT, 'tD', 'typeD')), by = c('PLT_CN', 'CONDID')) %>%
left_join(select(pops, PLT_CN, YEAR), by = 'PLT_CN') %>%
mutate(tDI = landD * aD_p * aD_c * tD * typeD * sp) %>%
group_by(.dots = grpBy) %>%
summarize(H_g = divIndex(grp, state * tDI, index = 'H'),
Eh_g = divIndex(grp, state * tDI, index = 'Eh'),
S_g = divIndex(grp, state * tDI, index = 'S'))
tOut <- left_join(tOut, fullArea, by = grpBy) %>%
mutate(H_b = H_g - H_a,
Eh_b = Eh_g - Eh_a,
S_b = S_g - S_a)
if (totals) {
tOut <- tOut %>%
select(grpBy, H_a, Eh_a, S_a, H_b, Eh_b, S_b, H_g, Eh_g, S_g, AREA_TOTAL, H_a_SE, Eh_a_SE, S_a_SE, AREA_TOTAL_SE, nPlots)
} else {
tOut <- tOut %>%
select(grpBy, H_a, Eh_a, S_a, H_b, Eh_b, S_b, H_g, Eh_g, S_g, H_a_SE, Eh_a_SE, S_a_SE, nPlots)
}
# Snag the names
tNames <- names(tOut)[names(tOut) %in% grpBy == FALSE]
}
## Pretty output
tOut <- tOut %>%
ungroup() %>%
mutate_if(is.factor, as.character) %>%
drop_na(grpBy) %>%
arrange(YEAR) %>%
as_tibble()
# Return a spatial object
if (!is.null(polys) & byPlot == FALSE) {
## NO IMPLICIT NA
nospGrp <- unique(grpBy[grpBy %in% c('SPCD', 'SYMBOL', 'COMMON_NAME', 'SCIENTIFIC_NAME') == FALSE])
nospSym <- syms(nospGrp)
tOut <- complete(tOut, !!!nospSym)
## If species, we don't want unique combos of variables related to same species
## but we do want NAs in polys where species are present
if (length(nospGrp) < length(grpBy)){
spGrp <- unique(grpBy[grpBy %in% c('SPCD', 'SYMBOL', 'COMMON_NAME', 'SCIENTIFIC_NAME')])
spSym <- syms(spGrp)
tOut <- complete(tOut, nesting(!!!nospSym))
}
suppressMessages({suppressWarnings({tOut <- left_join(tOut, polys) %>%
select(c('YEAR', grpByOrig, tNames, names(polys))) %>%
filter(!is.na(polyID))})})
## Makes it horrible to work with as a dataframe
if (returnSpatial == FALSE) tOut <- select(tOut, -c(geometry))
} else if (!is.null(polys) & byPlot){
polys <- as.data.frame(polys)
tOut <- left_join(tOut, select(polys, -c(geometry)), by = 'polyID')
}
## For spatial plots
if (returnSpatial & byPlot) grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
## Above converts to tibble
if (returnSpatial) tOut <- st_sf(tOut)
# ## remove any duplicates in byPlot (artifact of END_INYR loop)
if (byPlot) tOut <- unique(tOut)
return(tOut)
}
hunter-stanke/rFIA_TX documentation built on Jan. 1, 2021, 3:22 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions diversity_backup diversity diversityStarter
Documented in diversity
diversityStarter <- function(x,
db,
grpBy_quo = NULL,
polys = NULL,
returnSpatial = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
stateVar = TPA_UNADJ,
grpVar = SPCD,
treeDomain = NULL,
areaDomain = NULL,
byPlot = FALSE,
totals = FALSE,
nCores = 1,
remote,
mr){
reqTables <- c('PLOT', 'TREE', 'COND', 'POP_PLOT_STRATUM_ASSGN', 'POP_ESTN_UNIT', 'POP_EVAL',
'POP_STRATUM', 'POP_EVAL_TYP', 'POP_EVAL_GRP')
if (remote){
## Store the original parameters here
params <- db
## Read in one state at a time
db <- readFIA(dir = db$dir, common = db$common,
tables = reqTables, states = x, ## x is the vector of state names
nCores = nCores)
## If a clip was specified, run it now
if ('mostRecent' %in% names(params)){
db <- clipFIA(db, mostRecent = params$mostRecent,
mask = params$mask, matchEval = params$matchEval,
evalid = params$evalid, designCD = params$designCD,
nCores = nCores)
}
} else {
## Really only want the required tables
db <- db[names(db) %in% reqTables]
}
## Handle the modifier if it was given
db$TREE <- db$TREE %>%
mutate(state = !!stateVar,
grp = !!grpVar)
## Need a plotCN, and a new ID
db$PLOT <- db$PLOT %>% mutate(PLT_CN = CN,
pltID = paste(UNITCD, STATECD, COUNTYCD, PLOT, sep = '_'))
## don't have to change original code
#grpBy_quo <- enquo(grpBy)
# Probably cheating, but it works
if (quo_name(grpBy_quo) != 'NULL'){
## Have to join tables to run select with this object type
plt_quo <- filter(db$PLOT, !is.na(PLT_CN))
## We want a unique error message here to tell us when columns are not present in data
d_quo <- tryCatch(
error = function(cnd) {
return(0)
},
plt_quo[10,] %>% # Just the first row
left_join(select(db$COND, PLT_CN, names(db$COND)[names(db$COND) %in% names(db$PLOT) == FALSE]), by = 'PLT_CN') %>%
inner_join(select(db$TREE, PLT_CN, names(db$TREE)[names(db$TREE) %in% c(names(db$PLOT), names(db$COND)) == FALSE]), by = 'PLT_CN') %>%
select(!!grpBy_quo)
)
# If column doesnt exist, just returns 0, not a dataframe
if (is.null(nrow(d_quo))){
grpName <- quo_name(grpBy_quo)
stop(paste('Columns', grpName, 'not found in PLOT, TREE, or COND tables. Did you accidentally quote the variables names? e.g. use grpBy = ECOSUBCD (correct) instead of grpBy = "ECOSUBCD". ', collapse = ', '))
} else {
# Convert to character
grpBy <- names(d_quo)
}
} else {
grpBy <- NULL
}
reqTables <- c('PLOT', 'TREE', 'COND', 'POP_PLOT_STRATUM_ASSGN', 'POP_ESTN_UNIT', 'POP_EVAL',
'POP_STRATUM', 'POP_EVAL_TYP', 'POP_EVAL_GRP')
if (!is.null(polys) & first(class(polys)) %in% c('sf', 'SpatialPolygons', 'SpatialPolygonsDataFrame') == FALSE){
stop('polys must be spatial polygons object of class sp or sf. ')
}
if (landType %in% c('timber', 'forest') == FALSE){
stop('landType must be one of: "forest" or "timber".')
}
if (any(reqTables %in% names(db) == FALSE)){
missT <- reqTables[reqTables %in% names(db) == FALSE]
stop(paste('Tables', paste (as.character(missT), collapse = ', '), 'not found in object db.'))
}
if (str_to_upper(method) %in% c('TI', 'SMA', 'LMA', 'EMA', 'ANNUAL') == FALSE) {
warning(paste('Method', method, 'unknown. Defaulting to Temporally Indifferent (TI).'))
}
# I like a unique ID for a plot through time
if (byPlot) {grpBy <- c('pltID', grpBy)}
# Save original grpBy for pretty return with spatial objects
grpByOrig <- grpBy
## IF the object was clipped
if ('prev' %in% names(db$PLOT)){
## Only want the current plots, no grm
db$PLOT <- filter(db$PLOT, prev == 0)
}
### AREAL SUMMARY PREP
if(!is.null(polys)) {
# # Convert polygons to an sf object
# polys <- polys %>%
# as('sf')%>%
# mutate_if(is.factor,
# as.character)
# ## A unique ID
# polys$polyID <- 1:nrow(polys)
#
# # Add shapefile names to grpBy
grpBy = c(grpBy, 'polyID')
## Make plot data spatial, projected same as polygon layer
pltSF <- select(db$PLOT, c('LON', 'LAT', pltID)) %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
distinct(pltID, .keep_all = TRUE)
coordinates(pltSF) <- ~LON+LAT
proj4string(pltSF) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
pltSF <- as(pltSF, 'sf') %>%
st_transform(crs = st_crs(polys))
## Split up polys
polyList <- split(polys, as.factor(polys$polyID))
suppressWarnings({suppressMessages({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
cl <- makeCluster(nCores)
clusterEvalQ(cl, {
library(dplyr)
library(stringr)
library(rFIA)
})
out <- parLapply(cl, X = names(polyList), fun = areal_par, pltSF, polyList)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(polyList), FUN = areal_par, pltSF, polyList, mc.cores = nCores)
}
})})
pltSF <- bind_rows(out)
# A warning
if (length(unique(pltSF$pltID)) < 1){
stop('No plots in db overlap with polys.')
}
## Add polygon names to PLOT
db$PLOT <- db$PLOT %>%
left_join(select(pltSF, polyID, pltID), by = 'pltID')
# Test if any polygons cross state boundaries w/ different recent inventory years (continued w/in loop)
if ('mostRecent' %in% names(db) & length(unique(db$POP_EVAL$STATECD)) > 1){
mergeYears <- pltSF %>%
right_join(select(db$PLOT, PLT_CN, pltID), by = 'pltID') %>%
inner_join(select(db$POP_PLOT_STRATUM_ASSGN, c('PLT_CN', 'EVALID', 'STATECD')), by = 'PLT_CN') %>%
inner_join(select(db$POP_EVAL, c('EVALID', 'END_INVYR')), by = 'EVALID') %>%
group_by(polyID) %>%
summarize(maxYear = max(END_INVYR, na.rm = TRUE))
}
## TO RETURN SPATIAL PLOTS
}
if (byPlot & returnSpatial){
grpBy <- c(grpBy, 'LON', 'LAT')
} # END AREAL
## Build domain indicator function which is 1 if observation meets criteria, and 0 otherwise
# Land type domain indicator
if (tolower(landType) == 'forest'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1, 1, 0)
} else if (tolower(landType) == 'timber'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1 & db$COND$SITECLCD %in% c(1, 2, 3, 4, 5, 6) & db$COND$RESERVCD == 0, 1, 0)
}
# Tree Type domain indicator
if (tolower(treeType) == 'live'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 1, 1, 0)
} else if (tolower(treeType) == 'dead'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 2 & db$TREE$STANDING_DEAD_CD == 1, 1, 0)
} else if (tolower(treeType) == 'gs'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 1 & db$TREE$DIA >= 5 & db$TREE$TREECLCD == 2, 1, 0)
} else if (tolower(treeType) == 'all'){
db$TREE$typeD <- 1
}
# update spatial domain indicator
if(!is.null(polys)){
db$PLOT$sp <- ifelse(db$PLOT$pltID %in% pltSF$pltID, 1, 0)
} else {
db$PLOT$sp <- 1
}
# User defined domain indicator for area (ex. specific forest type)
pcEval <- left_join(db$PLOT, select(db$COND, -c('STATECD', 'UNITCD', 'COUNTYCD', 'INVYR', 'PLOT')), by = 'PLT_CN')
#areaDomain <- substitute(areaDomain)
pcEval$aD <- rlang::eval_tidy(areaDomain, pcEval) ## LOGICAL, THIS IS THE DOMAIN INDICATOR
if(!is.null(pcEval$aD)) pcEval$aD[is.na(pcEval$aD)] <- 0 # Make NAs 0s. Causes bugs otherwise
if(is.null(pcEval$aD)) pcEval$aD <- 1 # IF NULL IS GIVEN, THEN ALL VALUES TRUE
pcEval$aD <- as.numeric(pcEval$aD)
db$COND <- left_join(db$COND, select(pcEval, c('PLT_CN', 'CONDID', 'aD')), by = c('PLT_CN', 'CONDID')) %>%
mutate(aD_c = aD)
aD_p <- pcEval %>%
group_by(PLT_CN) %>%
summarize(aD_p = as.numeric(any(aD > 0)))
db$PLOT <- left_join(db$PLOT, aD_p, by = 'PLT_CN')
rm(pcEval)
# Same as above for tree (ex. trees > 20 ft tall)
#treeDomain <- substitute(treeDomain)
tD <- rlang::eval_tidy(treeDomain, db$TREE) ## LOGICAL, THIS IS THE DOMAIN INDICATOR
if(!is.null(tD)) tD[is.na(tD)] <- 0 # Make NAs 0s. Causes bugs otherwise
if(is.null(tD)) tD <- 1 # IF NULL IS GIVEN, THEN ALL VALUES TRUE
db$TREE$tD <- as.numeric(tD)
### Snag the EVALIDs that are needed
db$POP_EVAL<- db$POP_EVAL %>%
select('CN', 'END_INVYR', 'EVALID', 'ESTN_METHOD', STATECD) %>%
inner_join(select(db$POP_EVAL_TYP, c('EVAL_CN', 'EVAL_TYP')), by = c('CN' = 'EVAL_CN')) %>%
filter(EVAL_TYP == 'EXPVOL' | EVAL_TYP == 'EXPCURR') %>%
filter(!is.na(END_INVYR) & !is.na(EVALID) & END_INVYR >= 2003) %>%
distinct(END_INVYR, EVALID, .keep_all = TRUE)# %>%
#group_by(END_INVYR) %>%
#summarise(id = list(EVALID)
## If a most-recent subset, make sure that we don't get two reporting years in
## western states
if (mr) {
db$POP_EVAL <- db$POP_EVAL %>%
group_by(EVAL_TYP, STATECD) %>%
filter(END_INVYR == max(END_INVYR, na.rm = TRUE)) %>%
ungroup()
}
## Cut STATECD
db$POP_EVAL <- select(db$POP_EVAL, -c(STATECD))
### The population tables
pops <- select(db$POP_EVAL, c('EVALID', 'ESTN_METHOD', 'CN', 'END_INVYR', 'EVAL_TYP')) %>%
rename(EVAL_CN = CN) %>%
left_join(select(db$POP_ESTN_UNIT, c('CN', 'EVAL_CN', 'AREA_USED', 'P1PNTCNT_EU')), by = c('EVAL_CN')) %>%
rename(ESTN_UNIT_CN = CN) %>%
left_join(select(db$POP_STRATUM, c('ESTN_UNIT_CN', 'EXPNS', 'P2POINTCNT', 'CN', 'P1POINTCNT', 'ADJ_FACTOR_SUBP', 'ADJ_FACTOR_MICR', "ADJ_FACTOR_MACR")), by = c('ESTN_UNIT_CN')) %>%
rename(STRATUM_CN = CN) %>%
left_join(select(db$POP_PLOT_STRATUM_ASSGN, c('STRATUM_CN', 'PLT_CN', 'INVYR', 'STATECD')), by = 'STRATUM_CN') %>%
mutate_if(is.factor,
as.character)
### Which estimator to use?
if (str_to_upper(method) %in% c('ANNUAL')){
## Want to use the year where plots are measured, no repeats
## Breaking this up into pre and post reporting becuase
## Estimation units get weird on us otherwise
popOrig <- pops
pops <- pops %>%
group_by(STATECD) %>%
filter(END_INVYR == INVYR) %>%
ungroup()
prePops <- popOrig %>%
group_by(STATECD) %>%
filter(INVYR < min(END_INVYR, na.rm = TRUE)) %>%
distinct(PLT_CN, .keep_all = TRUE) %>%
ungroup()
pops <- bind_rows(pops, prePops) %>%
mutate(YEAR = INVYR)
} else { # Otherwise temporally indifferent
pops <- rename(pops, YEAR = END_INVYR)
}
## P2POINTCNT column is NOT consistent for annnual estimates, plots
## within individual strata and est units are related to different INVYRs
p2_INVYR <- pops %>%
group_by(ESTN_UNIT_CN, STRATUM_CN, INVYR) %>%
summarize(P2POINTCNT_INVYR = length(unique(PLT_CN)))
## Want a count of p2 points / eu, gets screwed up with grouping below
p2eu_INVYR <- p2_INVYR %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, INVYR, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN, INVYR) %>%
summarize(p2eu_INVYR = sum(P2POINTCNT_INVYR, na.rm = TRUE))
p2eu <- pops %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(p2eu = sum(P2POINTCNT, na.rm = TRUE))
## Rejoin
pops <- pops %>%
left_join(p2_INVYR, by = c('ESTN_UNIT_CN', 'STRATUM_CN', 'INVYR')) %>%
left_join(p2eu_INVYR, by = c('ESTN_UNIT_CN', 'INVYR')) %>%
left_join(p2eu, by = 'ESTN_UNIT_CN')
## Recode a few of the estimation methods to make things easier below
pops$ESTN_METHOD = recode(.x = pops$ESTN_METHOD,
`Post-Stratification` = 'strat',
`Stratified random sampling` = 'strat',
`Double sampling for stratification` = 'double',
`Simple random sampling` = 'simple',
`Subsampling units of unequal size` = 'simple')
## Break into size classes
if (bySizeClass){
grpBy <- c(grpBy, 'sizeClass')
grpByOrig <- c(grpByOrig, 'sizeClass')
db$TREE$sizeClass <- makeClasses(db$TREE$DIA, interval = 2, numLabs = TRUE)
db$TREE <- db$TREE[!is.na(db$TREE$sizeClass),]
}
## Only the necessary plots for EVAL of interest
db$PLOT <- filter(db$PLOT, PLT_CN %in% pops$PLT_CN)
## Which grpByNames are in which table? Helps us subset below
grpP <- names(db$PLOT)[names(db$PLOT) %in% grpBy]
grpC <- names(db$COND)[names(db$COND) %in% grpBy & names(db$COND) %in% grpP == FALSE]
grpT <- names(db$TREE)[names(db$TREE) %in% grpBy & names(db$TREE) %in% c(grpP, grpC) == FALSE]
db$PLOT <- select(db$PLOT, c('PLT_CN', 'STATECD', 'COUNTYCD', 'MACRO_BREAKPOINT_DIA', 'INVYR', 'MEASYEAR', 'PLOT_STATUS_CD', all_of(grpP), 'aD_p', 'sp'))
db$COND <- select(db$COND, c('PLT_CN', 'CONDPROP_UNADJ', 'PROP_BASIS', 'COND_STATUS_CD', 'CONDID', all_of(grpC), 'aD_c', 'landD'))
db$TREE <- db$TREE %>%
mutate(TRE_CN = CN) %>%
select(c('PLT_CN', 'TRE_CN', 'CONDID', 'DIA', 'grp', 'state', 'SUBP', 'TREE', all_of(grpT), 'tD', 'typeD'))
## Merging state and county codes
plts <- split(db$PLOT, as.factor(paste(db$PLOT$COUNTYCD, db$PLOT$STATECD, sep = '_')))
#plts <- split(db$PLOT, as.factor(db$PLOT$STATECD))
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
cl <- makeCluster(nCores)
clusterEvalQ(cl, {
library(dplyr)
library(stringr)
library(rFIA)
})
out <- parLapply(cl, X = names(plts), fun = divHelper1, plts, db, grpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = divHelper1, plts, db, grpBy, byPlot, mc.cores = nCores)
}
})
if (byPlot){
## back to dataframes
out <- unlist(out, recursive = FALSE)
tOut <- bind_rows(out[names(out) == 't'])
## Make it spatial
if (returnSpatial){
tOut <- tOut %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
st_as_sf(coords = c('LON', 'LAT'),
crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
}
out <- list(tEst = tOut, grpBy = grpBy, grpByOrig = grpByOrig)
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
t <- bind_rows(out[names(out) == 't'])
full <- bind_rows(out[names(out) == 'full'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
## Splitting up by STATECD and groups of 25 ESTN_UNIT_CNs
#estunit <- distinct(pops, ESTN_UNIT_CN) #%>%
#mutate(estID)
#estID <- seq(1, nrow(estunit), 50)
#estunit$estID <- rep_len(estID, length.out = nrow(estunit))
#pops <- pops %>%
# left_join(estunit, by = 'ESTN_UNIT_CN') #%>%
#mutate(estBreaks = )
#popState <- split(pops, as.factor(pops$estID))
popState <- split(pops, as.factor(pops$STATECD))
#
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
## Use the same cluster as above
# cl <- makeCluster(nCores)
# clusterEvalQ(cl, {
# library(dplyr)
# library(stringr)
# library(rFIA)
# })
out <- parLapply(cl, X = names(popState), fun = divHelper2, popState, t, full, grpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = divHelper2, popState, t, full, grpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
full <- bind_rows(out[names(out) == 'full'])
##### ----------------- MOVING AVERAGES
if (str_to_upper(method) %in% c("SMA", 'EMA', 'LMA')){
### ---- SIMPLE MOVING AVERAGE
if (str_to_upper(method) == 'SMA'){
## Assuming a uniform weighting scheme
wgts <- pops %>%
group_by(ESTN_UNIT_CN) %>%
summarize(wgt = 1 / length(unique(INVYR)))
#aEst <- left_join(aEst, wgts, by = 'ESTN_UNIT_CN')
tEst <- left_join(tEst, wgts, by = 'ESTN_UNIT_CN')
#### ----- Linear MOVING AVERAGE
} else if (str_to_upper(method) == 'LMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
## Want a number of INVYRs per EU
neu <- wgts %>%
group_by(ESTN_UNIT_CN) %>%
summarize(n = sum(rank, na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = rank / n) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
#### ----- EXPONENTIAL MOVING AVERAGE
} else if (str_to_upper(method) == 'EMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
if (length(lambda) < 2){
## Want sum of weighitng functions
neu <- wgts %>%
mutate(l = lambda) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
} else {
grpBy <- c('lambda', grpBy)
#aGrpBy <- c('lambda', aGrpBy)
## Duplicate weights for each level of lambda
yrWgts <- list()
for (i in 1:length(unique(lambda))) {
yrWgts[[i]] <- mutate(wgts, lambda = lambda[i])
}
wgts <- bind_rows(yrWgts)
## Want sum of weighitng functions
neu <- wgts %>%
group_by(lambda, ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = c('lambda', 'ESTN_UNIT_CN')) %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(lambda, ESTN_UNIT_CN, INVYR, wgt)
}
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
}
### Applying the weights
tEst <- tEst %>%
mutate_at(vars(aEst:sEst), ~(.*wgt)) %>%
mutate_at(vars(aVar:cvEst_s), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
}
out <- list(tEst = tEst, full = full, grpBy = grpBy, grpByOrig = grpByOrig)
}
return(out)
}
#' @export
diversity <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
stateVar = TPA_UNADJ,
grpVar = SPCD,
treeDomain = NULL,
areaDomain = NULL,
byPlot = FALSE,
totals = FALSE,
variance = FALSE,
nCores = 1) {
## don't have to change original code
grpBy_quo <- rlang::enquo(grpBy)
areaDomain <- rlang::enquo(areaDomain)
treeDomain <- rlang::enquo(treeDomain)
stateVar <- rlang::enquo(stateVar)
grpVar <- rlang::enquo(grpVar)
### Is DB remote?
remote <- ifelse(class(db) == 'Remote.FIA.Database', 1, 0)
if (remote){
iter <- db$states
## In memory
} else {
## Some warnings
if (class(db) != "FIA.Database"){
stop('db must be of class "FIA.Database". Use readFIA() to load your FIA data.')
}
## an iterator for remote
iter <- 1
}
## Check for a most recent subset
if (remote){
if ('mostRecent' %in% names(db)){
mr = db$mostRecent # logical
} else {
mr = FALSE
}
## In-memory
} else {
if ('mostRecent' %in% names(db)){
mr = TRUE
} else {
mr = FALSE
}
}
### AREAL SUMMARY PREP
if(!is.null(polys)) {
# Convert polygons to an sf object
polys <- polys %>%
as('sf')%>%
mutate_if(is.factor,
as.character)
## A unique ID
polys$polyID <- 1:nrow(polys)
}
## Run the main portion
out <- lapply(X = iter, FUN = diversityStarter, db,
grpBy_quo = grpBy_quo, polys, returnSpatial,
bySizeClass,
landType, treeType, method,
lambda, stateVar, grpVar,
treeDomain, areaDomain,
byPlot, totals, nCores, remote, mr)
## Bring the results back
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
full <- bind_rows(out[names(out) == 'full'])
grpBy <- out[names(out) == 'grpBy'][[1]]
grpByOrig <- out[names(out) == 'grpByOrig'][[1]]
if (byPlot){
tOut <- tEst
} else {
suppressMessages({suppressWarnings({
## If a clip was specified, handle the reporting years
if (mr){
## If a most recent subset, ignore differences in reporting years across states
## instead combine most recent information from each state
# ID mr years by group
maxyearsT <- tEst %>%
select(grpBy) %>%
group_by(.dots = grpBy[!c(grpBy %in% 'YEAR')]) %>%
summarise(YEAR = max(YEAR, na.rm = TRUE))
# Combine estimates
tEst <- tEst %>%
ungroup() %>%
select(-c(YEAR)) %>%
left_join(maxyearsT, by = grpBy[!c(grpBy %in% 'YEAR')])
}
})})
##--------------------- TOTALS and RATIOS
suppressWarnings({
tOut <- tEst %>%
group_by(.dots = grpBy) %>%
summarize_all(sum,na.rm = TRUE) %>%
mutate(H_a = hEst / aEst,
Eh_a = ehEst / aEst,
S_a = sEst / aEst,
AREA_TOTAL = aEst,
## Ratio variance
haVar = (1/AREA_TOTAL^2) * (hVar + (H_a^2 * aVar) - 2 * H_a * cvEst_h),
ehaVar = (1/AREA_TOTAL^2) * (ehVar + (Eh_a^2 * aVar) - (2 * Eh_a * cvEst_eh)),
saVar = (1/AREA_TOTAL^2) * (sVar + (S_a^2 * aVar) - 2 * S_a * cvEst_s),
H_a_VAR = haVar,
Eh_a_VAR = ehaVar,
S_a_VAR = saVar,
## RATIO SE
H_a_SE = sqrt(haVar) / H_a * 100,
Eh_a_SE = sqrt(ehaVar) / Eh_a * 100,
S_a_SE = sqrt(saVar) / S_a * 100,
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL * 100,
AREA_TOTAL_VAR = aVar,
nPlots = plotIn_AREA)
})
### Up a few spatial scales
full <- full %>%
group_by(.dots = grpBy[!c(grpBy %in% 'lambda')]) %>%
summarize(H_g = divIndex(grp, state, index = 'H'),
Eh_g = divIndex(grp, state, index = 'Eh'),
S_g = divIndex(grp, state, index = 'S'))
tOut <- tOut %>%
left_join(full, by = grpBy[!c(grpBy %in% 'lambda')]) %>%
mutate(H_b = H_g - H_a,
Eh_b = Eh_g - Eh_a,
S_b = S_g - S_a)
if (totals) {
if (variance){
tOut <- tOut %>%
select(grpBy, H_a, H_b, H_g, Eh_a, Eh_b, Eh_g, S_a, S_b, S_g, AREA_TOTAL, H_a_VAR, Eh_a_VAR, S_a_VAR, AREA_TOTAL_VAR, nPlots, N)
} else {
tOut <- tOut %>%
select(grpBy, H_a, H_b, H_g, Eh_a, Eh_b, Eh_g, S_a, S_b, S_g, AREA_TOTAL, H_a_SE, Eh_a_SE, S_a_SE, AREA_TOTAL_SE, nPlots)
}
} else {
if (variance){
tOut <- tOut %>%
select(grpBy, H_a, H_b, H_g, Eh_a, Eh_b, Eh_g, S_a, S_b, S_g, H_a_VAR, Eh_a_VAR, S_a_VAR, nPlots, N)
} else {
tOut <- tOut %>%
select(grpBy, H_a, H_b, H_g, Eh_a, Eh_b, Eh_g, S_a, S_b, S_g, H_a_SE, Eh_a_SE, S_a_SE, nPlots)
}
}
# Snag the names
tNames <- names(tOut)[names(tOut) %in% grpBy == FALSE]
}
## Pretty output
tOut <- tOut %>%
ungroup() %>%
mutate_if(is.factor, as.character) %>%
drop_na(grpBy) %>%
arrange(YEAR) %>%
as_tibble()
# Return a spatial object
if (!is.null(polys) & byPlot == FALSE) {
## NO IMPLICIT NA
nospGrp <- unique(grpBy[grpBy %in% c('SPCD', 'SYMBOL', 'COMMON_NAME', 'SCIENTIFIC_NAME') == FALSE])
nospSym <- syms(nospGrp)
tOut <- complete(tOut, !!!nospSym)
## If species, we don't want unique combos of variables related to same species
## but we do want NAs in polys where species are present
if (length(nospGrp) < length(grpBy)){
spGrp <- unique(grpBy[grpBy %in% c('SPCD', 'SYMBOL', 'COMMON_NAME', 'SCIENTIFIC_NAME')])
spSym <- syms(spGrp)
tOut <- complete(tOut, nesting(!!!nospSym))
}
suppressMessages({suppressWarnings({tOut <- left_join(tOut, polys) %>%
select(c('YEAR', grpByOrig, tNames, names(polys))) %>%
filter(!is.na(polyID) & !is.na(nPlots))})})
## Makes it horrible to work with as a dataframe
if (returnSpatial == FALSE) tOut <- select(tOut, -c(geometry))
} else if (!is.null(polys) & byPlot){
polys <- as.data.frame(polys)
tOut <- left_join(tOut, select(polys, -c(geometry)), by = 'polyID')
}
## For spatial plots
if (returnSpatial & byPlot) grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
## Above converts to tibble
if (returnSpatial) tOut <- st_sf(tOut)
# ## remove any duplicates in byPlot (artifact of END_INYR loop)
if (byPlot) tOut <- unique(tOut)
return(tOut)
}
diversity_backup <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
stateVar = TPA_UNADJ,
grpVar = SPCD,
treeDomain = NULL,
areaDomain = NULL,
byPlot = FALSE,
totals = FALSE,
nCores = 1) {
## Need a plotCN, and a new ID
db$PLOT <- db$PLOT %>% mutate(PLT_CN = CN,
pltID = paste(UNITCD, STATECD, COUNTYCD, PLOT, sep = '_'))
## Converting names given in grpBy to character vector (NSE to standard)
## don't have to change original code
grpBy_quo <- enquo(grpBy)
# Probably cheating, but it works
if (quo_name(grpBy_quo) != 'NULL'){
## Have to join tables to run select with this object type
plt_quo <- filter(db$PLOT, !is.na(PLT_CN))
## We want a unique error message here to tell us when columns are not present in data
d_quo <- tryCatch(
error = function(cnd) {
return(0)
},
plt_quo[10,] %>% # Just the first row
left_join(select(db$COND, PLT_CN, names(db$COND)[names(db$COND) %in% names(db$PLOT) == FALSE]), by = 'PLT_CN') %>%
inner_join(select(db$TREE, PLT_CN, names(db$TREE)[names(db$TREE) %in% c(names(db$PLOT), names(db$COND)) == FALSE]), by = 'PLT_CN') %>%
select(!!grpBy_quo)
)
# If column doesnt exist, just returns 0, not a dataframe
if (is.null(nrow(d_quo))){
grpName <- quo_name(grpBy_quo)
stop(paste('Columns', grpName, 'not found in PLOT, TREE, or COND tables. Did you accidentally quote the variables names? e.g. use grpBy = ECOSUBCD (correct) instead of grpBy = "ECOSUBCD". ', collapse = ', '))
} else {
# Convert to character
grpBy <- names(d_quo)
}
}
reqTables <- c('PLOT', 'TREE', 'COND', 'POP_PLOT_STRATUM_ASSGN', 'POP_ESTN_UNIT', 'POP_EVAL',
'POP_STRATUM', 'POP_EVAL_TYP', 'POP_EVAL_GRP')
## Some warnings
if (class(db) != "FIA.Database"){
stop('db must be of class "FIA.Database". Use readFIA() to load your FIA data.')
}
if (!is.null(polys) & first(class(polys)) %in% c('sf', 'SpatialPolygons', 'SpatialPolygonsDataFrame') == FALSE){
stop('polys must be spatial polygons object of class sp or sf. ')
}
if (landType %in% c('timber', 'forest') == FALSE){
stop('landType must be one of: "forest" or "timber".')
}
if (any(reqTables %in% names(db) == FALSE)){
missT <- reqTables[reqTables %in% names(db) == FALSE]
stop(paste('Tables', paste (as.character(missT), collapse = ', '), 'not found in object db.'))
}
if (str_to_upper(method) %in% c('TI', 'SMA', 'LMA', 'EMA', 'ANNUAL') == FALSE) {
warning(paste('Method', method, 'unknown. Defaulting to Temporally Indifferent (TI).'))
}
##### CHECK TO SEE IF STATE VARIABLE AND GROUP VARIABLE EXIST IN TREE
# Need to quote all variables for NSE
stateVar = enquo(stateVar)
grpVar = enquo(grpVar)
## If a modifier was given to a variable, handle it (ex. y = TPA_PERC / 100)
db$TREE <- db$TREE %>%
mutate(state = !!stateVar,
grp = !!grpVar)
# if (quo_name(stateVar) %in% names(db$TREE) == FALSE){
# stop(paste('Column',quo_name(stateVar), 'not found in TREE table. Did you accidentally quote the variables name? e.g. use stateVar = TPA_UNADJ (correct) instead of grpBy = "TPA_UNADJ". ', collapse = ', '))
# }
# if (quo_name(grpVar) %in% names(db$TREE) == FALSE){
# stop(paste('Column',quo_name(grpVar), 'not found in TREE table. Did you accidentally quote the variables name? e.g. use stateVar = TPA_UNADJ (correct) instead of grpBy = "TPA_UNADJ". ', collapse = ', '))
# }
# I like a unique ID for a plot through time
if (byPlot) {grpBy <- c('pltID', grpBy)}
# Save original grpBy for pretty return with spatial objects
grpByOrig <- grpBy
## IF the object was clipped
if ('prev' %in% names(db$PLOT)){
## Only want the current plots, no grm
db$PLOT <- filter(db$PLOT, prev == 0)
}
### DEAL WITH TEXAS
if (any(db$POP_EVAL$STATECD %in% 48)){
## Will require manual updates, fix your shit texas
txIDS <- db$POP_EVAL %>%
filter(STATECD %in% 48) %>%
filter(END_INVYR < 2017) %>%
filter(END_INVYR > 2006) %>%
## Removing any inventory that references east or west, sorry
filter(str_detect(str_to_upper(EVAL_DESCR), 'EAST', negate = TRUE) &
str_detect(str_to_upper(EVAL_DESCR), 'WEST', negate = TRUE))
db$POP_EVAL <- bind_rows(filter(db$POP_EVAL, !(STATECD %in% 48)), txIDS)
}
### AREAL SUMMARY PREP
if(!is.null(polys)) {
# Convert polygons to an sf object
polys <- polys %>%
as('sf')%>%
mutate_if(is.factor,
as.character)
## A unique ID
polys$polyID <- 1:nrow(polys)
# Add shapefile names to grpBy
#grpBy = c(names(polys)[str_detect(names(polys), 'geometry') == FALSE], 'polyID', grpBy)
grpBy = c(grpBy, 'polyID')
## Make plot data spatial, projected same as polygon layer
pltSF <- select(db$PLOT, c('LON', 'LAT', pltID)) %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
distinct(pltID, .keep_all = TRUE)
coordinates(pltSF) <- ~LON+LAT
proj4string(pltSF) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
pltSF <- as(pltSF, 'sf') %>%
st_transform(crs = st_crs(polys))
## Split up polys
polyList <- split(polys, as.factor(polys$polyID))
suppressWarnings({suppressMessages({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
cl <- makeCluster(nCores)
clusterEvalQ(cl, {
library(dplyr)
library(stringr)
library(rFIA)
})
out <- parLapply(cl, X = names(polyList), fun = areal_par, pltSF, polyList)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(polyList), FUN = areal_par, pltSF, polyList, mc.cores = nCores)
}
})})
pltSF <- bind_rows(out)
# A warning
if (length(unique(pltSF$pltID)) < 1){
stop('No plots in db overlap with polys.')
}
## Add polygon names to PLOT
db$PLOT <- db$PLOT %>%
left_join(select(pltSF, polyID, pltID), by = 'pltID')
# Test if any polygons cross state boundaries w/ different recent inventory years (continued w/in loop)
if ('mostRecent' %in% names(db) & length(unique(db$POP_EVAL$STATECD)) > 1){
mergeYears <- pltSF %>%
right_join(select(db$PLOT, PLT_CN, pltID), by = 'pltID') %>%
inner_join(select(db$POP_PLOT_STRATUM_ASSGN, c('PLT_CN', 'EVALID', 'STATECD')), by = 'PLT_CN') %>%
inner_join(select(db$POP_EVAL, c('EVALID', 'END_INVYR')), by = 'EVALID') %>%
group_by(polyID) %>%
summarize(maxYear = max(END_INVYR, na.rm = TRUE))
}
## TO RETURN SPATIAL PLOTS
}
if (byPlot & returnSpatial){
grpBy <- c(grpBy, 'LON', 'LAT')
} # END AREAL
## Build domain indicator function which is 1 if observation meets criteria, and 0 otherwise
# Land type domain indicator
if (tolower(landType) == 'forest'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1, 1, 0)
} else if (tolower(landType) == 'timber'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1 & db$COND$SITECLCD %in% c(1, 2, 3, 4, 5, 6) & db$COND$RESERVCD == 0, 1, 0)
}
# Tree Type domain indicator
if (tolower(treeType) == 'live'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 1, 1, 0)
} else if (tolower(treeType) == 'dead'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 2 & db$TREE$STANDING_DEAD_CD == 1, 1, 0)
} else if (tolower(treeType) == 'gs'){
db$TREE$typeD <- ifelse(db$TREE$STATUSCD == 1 & db$TREE$DIA >= 5 & db$TREE$TREECLCD == 2, 1, 0)
} else if (tolower(treeType) == 'all'){
db$TREE$typeD <- 1
}
# update spatial domain indicator
if(!is.null(polys)){
db$PLOT$sp <- ifelse(db$PLOT$pltID %in% pltSF$pltID, 1, 0)
} else {
db$PLOT$sp <- 1
}
# User defined domain indicator for area (ex. specific forest type)
pcEval <- left_join(db$PLOT, select(db$COND, -c('STATECD', 'UNITCD', 'COUNTYCD', 'INVYR', 'PLOT')), by = 'PLT_CN')
areaDomain <- substitute(areaDomain)
pcEval$aD <- eval(areaDomain, pcEval) ## LOGICAL, THIS IS THE DOMAIN INDICATOR
if(!is.null(pcEval$aD)) pcEval$aD[is.na(pcEval$aD)] <- 0 # Make NAs 0s. Causes bugs otherwise
if(is.null(pcEval$aD)) pcEval$aD <- 1 # IF NULL IS GIVEN, THEN ALL VALUES TRUE
pcEval$aD <- as.numeric(pcEval$aD)
db$COND <- left_join(db$COND, select(pcEval, c('PLT_CN', 'CONDID', 'aD')), by = c('PLT_CN', 'CONDID')) %>%
mutate(aD_c = aD)
aD_p <- pcEval %>%
group_by(PLT_CN) %>%
summarize(aD_p = as.numeric(any(aD > 0)))
db$PLOT <- left_join(db$PLOT, aD_p, by = 'PLT_CN')
rm(pcEval)
# Same as above for tree (ex. trees > 20 ft tall)
treeDomain <- substitute(treeDomain)
tD <- eval(treeDomain, db$TREE) ## LOGICAL, THIS IS THE DOMAIN INDICATOR
if(!is.null(tD)) tD[is.na(tD)] <- 0 # Make NAs 0s. Causes bugs otherwise
if(is.null(tD)) tD <- 1 # IF NULL IS GIVEN, THEN ALL VALUES TRUE
db$TREE$tD <- as.numeric(tD)
### Snag the EVALIDs that are needed
db$POP_EVAL<- db$POP_EVAL %>%
select('CN', 'END_INVYR', 'EVALID', 'ESTN_METHOD') %>%
inner_join(select(db$POP_EVAL_TYP, c('EVAL_CN', 'EVAL_TYP')), by = c('CN' = 'EVAL_CN')) %>%
filter(EVAL_TYP == 'EXPVOL' | EVAL_TYP == 'EXPCURR') %>%
filter(!is.na(END_INVYR) & !is.na(EVALID) & END_INVYR >= 2003) %>%
distinct(END_INVYR, EVALID, .keep_all = TRUE)# %>%
#group_by(END_INVYR) %>%
#summarise(id = list(EVALID)
## Make an annual panel ID, associated with an INVYR
### The population tables
pops <- select(db$POP_EVAL, c('EVALID', 'ESTN_METHOD', 'CN', 'END_INVYR')) %>%
rename(EVAL_CN = CN) %>%
left_join(select(db$POP_ESTN_UNIT, c('CN', 'EVAL_CN', 'AREA_USED', 'P1PNTCNT_EU')), by = c('EVAL_CN')) %>%
rename(ESTN_UNIT_CN = CN) %>%
left_join(select(db$POP_STRATUM, c('ESTN_UNIT_CN', 'EXPNS', 'P2POINTCNT', 'CN', 'P1POINTCNT', 'ADJ_FACTOR_SUBP', 'ADJ_FACTOR_MICR', "ADJ_FACTOR_MACR")), by = c('ESTN_UNIT_CN')) %>%
rename(STRATUM_CN = CN) %>%
left_join(select(db$POP_PLOT_STRATUM_ASSGN, c('STRATUM_CN', 'PLT_CN', 'INVYR', 'STATECD')), by = 'STRATUM_CN') %>%
ungroup() %>%
mutate_if(is.factor,
as.character)
### Which estimator to use?
if (str_to_upper(method) %in% c('ANNUAL')){
## Want to use the year where plots are measured, no repeats
## Breaking this up into pre and post reporting becuase
## Estimation units get weird on us otherwise
popOrig <- pops
pops <- pops %>%
group_by(STATECD) %>%
filter(END_INVYR == INVYR) %>%
ungroup()
prePops <- popOrig %>%
group_by(STATECD) %>%
filter(INVYR < min(END_INVYR, na.rm = TRUE)) %>%
distinct(PLT_CN, .keep_all = TRUE) %>%
ungroup()
pops <- bind_rows(pops, prePops) %>%
mutate(YEAR = INVYR)
} else { # Otherwise temporally indifferent
pops <- rename(pops, YEAR = END_INVYR)
}
## P2POINTCNT column is NOT consistent for annnual estimates, plots
## within individual strata and est units are related to different INVYRs
p2_INVYR <- pops %>%
group_by(ESTN_UNIT_CN, STRATUM_CN, INVYR) %>%
summarize(P2POINTCNT_INVYR = length(unique(PLT_CN)))
## Want a count of p2 points / eu, gets screwed up with grouping below
p2eu_INVYR <- p2_INVYR %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, INVYR, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN, INVYR) %>%
summarize(p2eu_INVYR = sum(P2POINTCNT_INVYR, na.rm = TRUE))
p2eu <- pops %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(p2eu = sum(P2POINTCNT, na.rm = TRUE))
## Rejoin
pops <- pops %>%
left_join(p2_INVYR, by = c('ESTN_UNIT_CN', 'STRATUM_CN', 'INVYR')) %>%
left_join(p2eu_INVYR, by = c('ESTN_UNIT_CN', 'INVYR')) %>%
left_join(p2eu, by = 'ESTN_UNIT_CN')
## Recode a few of the estimation methods to make things easier below
pops$ESTN_METHOD = recode(.x = pops$ESTN_METHOD,
`Post-Stratification` = 'strat',
`Stratified random sampling` = 'strat',
`Double sampling for stratification` = 'double',
`Simple random sampling` = 'simple',
`Subsampling units of unequal size` = 'simple')
## Break into size classes
if (bySizeClass){
grpBy <- c(grpBy, 'sizeClass')
grpByOrig <- c(grpByOrig, 'sizeClass')
db$TREE$sizeClass <- makeClasses(db$TREE$DIA, interval = 2, numLabs = TRUE)
db$TREE <- db$TREE[!is.na(db$TREE$sizeClass),]
}
## Only the necessary plots for EVAL of interest
db$PLOT <- filter(db$PLOT, PLT_CN %in% pops$PLT_CN)
## Merging state and county codes
plts <- split(db$PLOT, as.factor(paste(db$PLOT$COUNTYCD, db$PLOT$STATECD, sep = '_')))
#plts <- split(db$PLOT, as.factor(db$PLOT$STATECD))
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
cl <- makeCluster(nCores)
clusterEvalQ(cl, {
library(dplyr)
library(stringr)
library(rFIA)
})
out <- parLapply(cl, X = names(plts), fun = divHelper1, plts, db, grpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = divHelper1, plts, db, grpBy, byPlot, mc.cores = nCores)
}
})
if (byPlot){
## back to dataframes
out <- unlist(out, recursive = FALSE)
tOut <- bind_rows(out[names(out) == 't'])
## Make it spatial
if (returnSpatial){
tOut <- tOut %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
st_as_sf(coords = c('LON', 'LAT'),
crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
}
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
## Splitting up by ESTN_UNIT
popState <- split(pops, as.factor(pops$STATECD))
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
## Use the same cluster as above
# cl <- makeCluster(nCores)
# clusterEvalQ(cl, {
# library(dplyr)
# library(stringr)
# library(rFIA)
# })
out <- parLapply(cl, X = names(popState), fun = divHelper2, popState, t, grpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = divHelper2, popState, t, grpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
##### ----------------- MOVING AVERAGES
if (str_to_upper(method) %in% c("SMA", 'EMA', 'LMA')){
### ---- SIMPLE MOVING AVERAGE
if (str_to_upper(method) == 'SMA'){
## Assuming a uniform weighting scheme
wgts <- pops %>%
group_by(ESTN_UNIT_CN) %>%
summarize(wgt = 1 / length(unique(INVYR)))
#aEst <- left_join(aEst, wgts, by = 'ESTN_UNIT_CN')
tEst <- left_join(tEst, wgts, by = 'ESTN_UNIT_CN')
#### ----- Linear MOVING AVERAGE
} else if (str_to_upper(method) == 'LMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
## Want a number of INVYRs per EU
neu <- wgts %>%
group_by(ESTN_UNIT_CN) %>%
summarize(n = sum(rank, na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = rank / n) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
#### ----- EXPONENTIAL MOVING AVERAGE
} else if (str_to_upper(method) == 'EMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
if (length(lambda) < 2){
## Want sum of weighitng functions
neu <- wgts %>%
mutate(l = lambda) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
} else {
grpBy <- c('lambda', grpBy)
#aGrpBy <- c('lambda', aGrpBy)
## Duplicate weights for each level of lambda
yrWgts <- list()
for (i in 1:length(unique(lambda))) {
yrWgts[[i]] <- mutate(wgts, lambda = lambda[i])
}
wgts <- bind_rows(yrWgts)
## Want sum of weighitng functions
neu <- wgts %>%
group_by(lambda, ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = c('lambda', 'ESTN_UNIT_CN')) %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(lambda, ESTN_UNIT_CN, INVYR, wgt)
}
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
}
### Applying the weights
tEst <- tEst %>%
mutate_at(vars(aEst:sEst), ~(.*wgt)) %>%
mutate_at(vars(aVar:cvEst_s), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
}
##--------------------- TOTALS and RATIOS
suppressWarnings({
tOut <- tEst %>%
group_by(.dots = grpBy) %>%
summarize_all(sum,na.rm = TRUE) %>%
mutate(H_a = hEst / aEst,
Eh_a = ehEst / aEst,
S_a = sEst / aEst,
AREA_TOTAL = aEst,
## Ratio variance
haVar = (1/AREA_TOTAL^2) * (hVar + (H_a^2 * aVar) - 2 * H_a * cvEst_h),
ehaVar = (1/AREA_TOTAL^2) * (ehVar + (Eh_a^2 * aVar) - (2 * Eh_a * cvEst_eh)),
saVar = (1/AREA_TOTAL^2) * (sVar + (S_a^2 * aVar) - 2 * S_a * cvEst_s),
## RATIO SE
H_a_SE = sqrt(haVar) / H_a * 100,
Eh_a_SE = sqrt(ehaVar) / Eh_a * 100,
S_a_SE = sqrt(saVar) / S_a * 100,
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL * 100,
nPlots = plotIn_AREA)
})
### UP a few spatial scales
## Which grpByNames are in which table? Helps us subset below
grpP <- names(db$PLOT)[names(db$PLOT) %in% grpBy]
grpC <- names(db$COND)[names(db$COND) %in% grpBy & names(db$COND) %in% grpP == FALSE]
grpT <- names(db$TREE)[names(db$TREE) %in% grpBy & names(db$TREE) %in% c(grpP, grpC) == FALSE]
### REALLY DON'T LIKE WORKING WITH FULL TABLES, WORK ON A FIX
### Only joining tables necessary to produce plot level estimates, adjusted for non-response
fullArea <- select(db$PLOT, c('PLT_CN', 'STATECD', 'MACRO_BREAKPOINT_DIA', 'INVYR', 'MEASYEAR', 'PLOT_STATUS_CD', grpP, 'aD_p', 'sp')) %>%
left_join(select(db$COND, c('PLT_CN', 'CONDPROP_UNADJ', 'PROP_BASIS', 'COND_STATUS_CD', 'CONDID', grpC, 'aD_c', 'landD')), by = c('PLT_CN')) %>%
left_join(select(db$TREE, c('PLT_CN', 'CONDID', 'DIA', 'grp', 'state', 'SUBP', 'TREE', grpT, 'tD', 'typeD')), by = c('PLT_CN', 'CONDID')) %>%
left_join(select(pops, PLT_CN, YEAR), by = 'PLT_CN') %>%
mutate(tDI = landD * aD_p * aD_c * tD * typeD * sp) %>%
group_by(.dots = grpBy) %>%
summarize(H_g = divIndex(grp, state * tDI, index = 'H'),
Eh_g = divIndex(grp, state * tDI, index = 'Eh'),
S_g = divIndex(grp, state * tDI, index = 'S'))
tOut <- left_join(tOut, fullArea, by = grpBy) %>%
mutate(H_b = H_g - H_a,
Eh_b = Eh_g - Eh_a,
S_b = S_g - S_a)
if (totals) {
tOut <- tOut %>%
select(grpBy, H_a, Eh_a, S_a, H_b, Eh_b, S_b, H_g, Eh_g, S_g, AREA_TOTAL, H_a_SE, Eh_a_SE, S_a_SE, AREA_TOTAL_SE, nPlots)
} else {
tOut <- tOut %>%
select(grpBy, H_a, Eh_a, S_a, H_b, Eh_b, S_b, H_g, Eh_g, S_g, H_a_SE, Eh_a_SE, S_a_SE, nPlots)
}
# Snag the names
tNames <- names(tOut)[names(tOut) %in% grpBy == FALSE]
}
## Pretty output
tOut <- tOut %>%
ungroup() %>%
mutate_if(is.factor, as.character) %>%
drop_na(grpBy) %>%
arrange(YEAR) %>%
as_tibble()
# Return a spatial object
if (!is.null(polys) & byPlot == FALSE) {
## NO IMPLICIT NA
nospGrp <- unique(grpBy[grpBy %in% c('SPCD', 'SYMBOL', 'COMMON_NAME', 'SCIENTIFIC_NAME') == FALSE])
nospSym <- syms(nospGrp)
tOut <- complete(tOut, !!!nospSym)
## If species, we don't want unique combos of variables related to same species
## but we do want NAs in polys where species are present
if (length(nospGrp) < length(grpBy)){
spGrp <- unique(grpBy[grpBy %in% c('SPCD', 'SYMBOL', 'COMMON_NAME', 'SCIENTIFIC_NAME')])
spSym <- syms(spGrp)
tOut <- complete(tOut, nesting(!!!nospSym))
}
suppressMessages({suppressWarnings({tOut <- left_join(tOut, polys) %>%
select(c('YEAR', grpByOrig, tNames, names(polys))) %>%
filter(!is.na(polyID))})})
## Makes it horrible to work with as a dataframe
if (returnSpatial == FALSE) tOut <- select(tOut, -c(geometry))
} else if (!is.null(polys) & byPlot){
polys <- as.data.frame(polys)
tOut <- left_join(tOut, select(polys, -c(geometry)), by = 'polyID')
}
## For spatial plots
if (returnSpatial & byPlot) grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
## Above converts to tibble
if (returnSpatial) tOut <- st_sf(tOut)
# ## remove any duplicates in byPlot (artifact of END_INYR loop)
if (byPlot) tOut <- unique(tOut)
return(tOut)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.