R/tpaNew.R
In hunter-stanke/rFIA_TX: Duplicate rFIA for TX
Defines functions
tpa_backup
tpa_test
tpaStarter_test
tpa
tpaStarter
Documented in
tpa
tpaStarter <- function(x,
db,
grpBy_quo = NULL,
polys = NULL,
returnSpatial = FALSE,
bySpecies = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
treeDomain = NULL,
areaDomain = NULL,
totals = FALSE,
byPlot = 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]
}
## 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)
library(sf)
})
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')
## 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') %>%
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')
## Add species to groups
if (bySpecies) {
db$TREE <- db$TREE %>%
left_join(select(intData$REF_SPECIES_2018, c('SPCD','COMMON_NAME', 'GENUS', 'SPECIES')), by = 'SPCD') %>%
mutate(SCIENTIFIC_NAME = paste(GENUS, SPECIES, sep = ' ')) %>%
mutate_if(is.factor,
as.character)
grpBy <- c(grpBy, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
grpByOrig <- c(grpByOrig, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
}
## 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),]
}
# Seperate area grouping names, (ex. TPA red oak in total land area of ingham county, rather than only area where red oak occurs)
if (!is.null(polys)){
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND) | grpBy %in% names(pltSF)])
} else {
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND)])
}
## Only the necessary plots for EVAL of interest
db$PLOT <- filter(db$PLOT, PLT_CN %in% pops$PLT_CN)
## Narrow up the tables to the necessary variables, reduces memory load
## sent out the cores
## 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]
### Only joining tables necessary to produce plot level estimates, adjusted for non-response
db$PLOT <- select(db$PLOT, c('PLT_CN', 'STATECD', 'MACRO_BREAKPOINT_DIA', 'INVYR', 'MEASYEAR', 'PLOT_STATUS_CD', grpP, 'aD_p', 'sp'))
db$COND <- select(db$COND, c('PLT_CN', 'CONDPROP_UNADJ', 'PROP_BASIS', 'COND_STATUS_CD', 'CONDID', grpC, 'aD_c', 'landD')) %>%
filter(PLT_CN %in% db$PLOT$PLT_CN)
db$TREE <- select(db$TREE, c('PLT_CN', 'CONDID', 'DIA', 'SPCD', 'TPA_UNADJ', 'SUBP', 'TREE', grpT, 'tD', 'typeD')) %>%
filter(PLT_CN %in% db$PLOT$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 = tpaHelper1, plts,
db[names(db) %in% c('COND', 'TREE')],
grpBy, aGrpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = tpaHelper1, plts,
db[names(db) %in% c('COND', 'TREE')],
grpBy, aGrpBy, 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, aGrpBy = aGrpBy, grpByOrig = grpByOrig)
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
a <- bind_rows(out[names(out) == 'a'])
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
aGrpBy <- c('YEAR', aGrpBy)
## 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 = tpaHelper2, popState, a, t, grpBy, aGrpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = tpaHelper2, popState, a, t, grpBy, aGrpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
aEst <- bind_rows(out[names(out) == 'aEst'])
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
# Area
aEst <- aEst %>%
mutate_at(vars(aEst), ~(.*wgt)) %>%
mutate_at(vars(aVar), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = aGrpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
tEst <- tEst %>%
mutate_at(vars(tEst:bTEst), ~(.*wgt)) %>%
mutate_at(vars(tVar:cvEst_bT), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(tEst:cvEst_bT), sum, na.rm = TRUE)
}
out <- list(tEst = tEst, aEst = aEst, grpBy = grpBy, aGrpBy = aGrpBy, grpByOrig = grpByOrig)
}
return(out)
}
#' @export
tpa <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
bySpecies = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
treeDomain = NULL,
areaDomain = NULL,
totals = FALSE,
variance = FALSE,
byPlot = FALSE,
nCores = 1) {
## don't have to change original code
grpBy_quo <- rlang::enquo(grpBy)
areaDomain <- rlang::enquo(areaDomain)
treeDomain <- rlang::enquo(treeDomain)
### 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 = tpaStarter, db,
grpBy_quo = grpBy_quo, polys, returnSpatial,
bySpecies, bySizeClass,
landType, treeType, method,
lambda, treeDomain, areaDomain,
totals, byPlot, nCores, remote, mr)
## Bring the results back
out <- unlist(out, recursive = FALSE)
aEst <- bind_rows(out[names(out) == 'aEst'])
tEst <- bind_rows(out[names(out) == 'tEst'])
grpBy <- out[names(out) == 'grpBy'][[1]]
aGrpBy <- out[names(out) == 'aGrpBy'][[1]]
grpByOrig <- out[names(out) == 'grpByOrig'][[1]]
if (byPlot){
tOut <- tEst
## Population estimates
} 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))
maxyearsA <- aEst %>%
select(aGrpBy) %>%
group_by(.dots = aGrpBy[!c(aGrpBy %in% 'YEAR')]) %>%
summarise(YEAR = max(YEAR, na.rm = TRUE))
# Combine estimates
aEst <- aEst %>%
ungroup() %>%
select(-c(YEAR)) %>%
left_join(maxyearsA, by = aGrpBy[!c(aGrpBy %in% 'YEAR')])
tEst <- tEst %>%
ungroup() %>%
select(-c(YEAR)) %>%
left_join(maxyearsT, by = grpBy[!c(grpBy %in% 'YEAR')])
}
})})
##--------------------- TOTALS and RATIOS
# Area
aTotal <- aEst %>%
group_by(.dots = aGrpBy) %>%
summarize(AREA_TOTAL = sum(aEst, na.rm = TRUE),
aVar = sum(aVar, na.rm = TRUE),
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL * 100,
AREA_TOTAL_VAR = aVar,
nPlots_AREA = sum(plotIn_AREA, na.rm = TRUE))
# # Tree
# tTotal <- tEst %>%
# group_by(.dots = grpBy) %>%
# #left_join(aTotal, by = c(aGrpBy)) %>%
# summarize(TREE_TOTAL = sum(tEst, na.rm = TRUE),
# BA_TOTAL = sum(bEst, na.rm = TRUE),
# ## Variances
# treeVar = sum(tVar, na.rm = TRUE),
# baVar = sum(bVar, na.rm = TRUE),
# #aVar = first(aVar),
# cvT = sum(cvEst_t, na.rm = TRUE),
# cvB = sum(cvEst_b, na.rm = TRUE),
# ## Sampling Errors
# TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
# BA_SE = sqrt(baVar) / BA_TOTAL * 100,
# nPlots_TREE = sum(plotIn_TREE, na.rm = TRUE)) #%>%
tTotal <- tEst %>%
group_by(.dots = grpBy) %>%
summarize_all(sum, na.rm = TRUE) %>%
mutate(TREE_TOTAL = tEst,
BA_TOTAL = bEst,
## Variances
treeVar = tVar,
baVar = bVar,
#aVar = first(aVar),
cvT = cvEst_t,
cvB = cvEst_b,
## Sampling Errors
TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
TREE_VAR = treeVar,
BA_SE = sqrt(baVar) / BA_TOTAL * 100,
BA_VAR = baVar,
N = sum(N, na.rm = TRUE),
nPlots_TREE = plotIn_TREE) %>%
select(grpBy, TREE_TOTAL, BA_TOTAL, treeVar, baVar, cvT, cvB, TREE_SE, BA_SE, TREE_VAR, BA_VAR,
nPlots_TREE, N)
## IF using polys, we treat each zone as a unique population
if (!is.null(polys)){
propGrp <- c('polyID', grpBy)
} else {
propGrp <- 'YEAR'
}
## Hand the proportions
tpTotal <- tEst %>%
group_by(.dots = unique(propGrp)) %>%
summarize(TREE_TOTAL_full = sum(tTEst, na.rm = TRUE), ## Need to sum this
BA_TOTAL_full = sum(bTEst, na.rm = TRUE), ## Need to sum this
tTVar = sum(tTVar, na.rm = TRUE),
bTVar = sum(bTVar, na.rm = TRUE),
cvTT = sum(cvEst_tT, na.rm = TRUE),
cvBT = sum(cvEst_bT, na.rm = TRUE))
suppressWarnings({
## Bring them together
tTotal <- tTotal %>%
left_join(aTotal, by = aGrpBy) %>%
left_join(tpTotal, by = unique(propGrp)) %>%
mutate(TPA = TREE_TOTAL / AREA_TOTAL,
BAA = BA_TOTAL / AREA_TOTAL,
tpaVar = (1/AREA_TOTAL^2) * (treeVar + (TPA^2 * aVar) - 2 * TPA * cvT),
baaVar = (1/AREA_TOTAL^2) * (baVar + (BAA^2 * aVar) - (2 * BAA * cvB)),
TPA_SE = sqrt(tpaVar) / TPA * 100,
BAA_SE = sqrt(baaVar) / BAA * 100,
TPA_VAR = tpaVar,
BAA_VAR = baaVar,
TPA_PERC = TREE_TOTAL / (TREE_TOTAL_full) * 100,
BAA_PERC = BA_TOTAL / (BA_TOTAL_full) * 100,
tpVar = (1/TREE_TOTAL_full^2) * (treeVar + (TPA_PERC^2 * tTVar) - 2 * TPA_PERC * cvTT),
bpVar = (1/BA_TOTAL_full^2) * (baVar + (BAA_PERC^2 * bTVar) - (2 * BAA_PERC * cvBT)),
TPA_PERC_SE = sqrt(tpVar) / TPA_PERC * 100,
BAA_PERC_SE = sqrt(bpVar) / BAA_PERC * 100,
TPA_PERC_VAR = tpVar,
BAA_PERC_VAR = bpVar)
})
if (totals) {
if (variance){
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TREE_TOTAL, BA_TOTAL, AREA_TOTAL, TPA_VAR, BAA_VAR,
TPA_PERC_VAR, BAA_PERC_VAR, TREE_VAR, BA_VAR, AREA_TOTAL_VAR, nPlots_TREE, nPlots_AREA, N)
} else {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TREE_TOTAL, BA_TOTAL, AREA_TOTAL, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, TREE_SE, BA_SE, AREA_TOTAL_SE, nPlots_TREE, nPlots_AREA)
}
} else {
if (variance) {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TPA_VAR, BAA_VAR,
TPA_PERC_VAR, BAA_PERC_VAR, nPlots_TREE, nPlots_AREA, N)
} else {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, nPlots_TREE, nPlots_AREA)
}
}
# 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_AREA))})})
## 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)
}
tpaStarter_test <- function(x,
db,
grpBy_quo = NULL,
polys = NULL,
returnSpatial = FALSE,
bySpecies = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
treeDomain = NULL,
areaDomain = NULL,
totals = FALSE,
byPlot = FALSE,
nCores = 1,
remote){
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]
}
## 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)
}
### 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(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') %>%
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', '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') %>%
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')
## Add species to groups
if (bySpecies) {
db$TREE <- db$TREE %>%
left_join(select(intData$REF_SPECIES_2018, c('SPCD','COMMON_NAME', 'GENUS', 'SPECIES')), by = 'SPCD') %>%
mutate(SCIENTIFIC_NAME = paste(GENUS, SPECIES, sep = ' ')) %>%
mutate_if(is.factor,
as.character)
grpBy <- c(grpBy, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
grpByOrig <- c(grpByOrig, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
}
## 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),]
}
# Seperate area grouping names, (ex. TPA red oak in total land area of ingham county, rather than only area where red oak occurs)
if (!is.null(polys)){
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND) | grpBy %in% names(pltSF)])
} else {
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND)])
}
## Only the necessary plots for EVAL of interest
db$PLOT <- filter(db$PLOT, PLT_CN %in% pops$PLT_CN)
## Narrow up the tables to the necessary variables, reduces memory load
## sent out the cores
## 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]
### Only joining tables necessary to produce plot level estimates, adjusted for non-response
db$PLOT <- select(db$PLOT, c('PLT_CN', 'STATECD', 'MACRO_BREAKPOINT_DIA', 'INVYR', 'MEASYEAR', 'PLOT_STATUS_CD', grpP, 'aD_p', 'sp'))
db$COND <- select(db$COND, c('PLT_CN', 'CONDPROP_UNADJ', 'PROP_BASIS', 'COND_STATUS_CD', 'CONDID', grpC, 'aD_c', 'landD'))
db$TREE <- select(db$TREE, c('PLT_CN', 'CONDID', 'DIA', 'SPCD', 'TPA_UNADJ', 'SUBP', 'TREE', 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 = tpaHelper1_test, plts, db, grpBy, aGrpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = tpaHelper1_test, plts, db, grpBy, aGrpBy, 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, aGrpBy = aGrpBy, grpByOrig = grpByOrig)
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
a <- bind_rows(out[names(out) == 'a'])
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
aGrpBy <- c('YEAR', aGrpBy)
## 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 = tpaHelper2_test, popState, a, t, grpBy, aGrpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = tpaHelper2_test, popState, a, t, grpBy, aGrpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
aEst <- bind_rows(out[names(out) == 'aEst'])
tEst <- bind_rows(out[names(out) == 'tEst'])
out <- list(tEst = tEst, aEst = aEst, grpBy = grpBy, aGrpBy = aGrpBy, grpByOrig = grpByOrig)
}
return(out)
}
tpa_test <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
bySpecies = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
treeDomain = NULL,
areaDomain = NULL,
totals = FALSE,
byPlot = FALSE,
nCores = 1) {
## don't have to change original code
grpBy_quo <- rlang::enquo(grpBy)
areaDomain <- rlang::enquo(areaDomain)
treeDomain <- rlang::enquo(treeDomain)
### 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
}
### 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 = tpaStarter_test, db,
grpBy_quo = grpBy_quo, polys, returnSpatial,
bySpecies, bySizeClass,
landType, treeType, method,
lambda, treeDomain, areaDomain,
totals, byPlot, nCores, remote)
## Bring the results back
out <- unlist(out, recursive = FALSE)
aEst <- bind_rows(out[names(out) == 'aEst'])
tEst <- bind_rows(out[names(out) == 'tEst'])
grpBy <- out[names(out) == 'grpBy'][[1]]
aGrpBy <- out[names(out) == 'aGrpBy'][[1]]
grpByOrig <- out[names(out) == 'grpByOrig'][[1]]
if (byPlot){
tOut <- tEst
## Population estimates
} else {
##--------------------- TOTALS and RATIOS
# Area
aTotal <- aEst %>%
group_by(.dots = aGrpBy) %>%
summarize(AREA_TOTAL = sum(aEst, na.rm = TRUE),
aVar = sum(aVar, na.rm = TRUE),
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL * 100,
nPlots_AREA = sum(plotIn_AREA, na.rm = TRUE))
# # Tree
# tTotal <- tEst %>%
# group_by(.dots = grpBy) %>%
# #left_join(aTotal, by = c(aGrpBy)) %>%
# summarize(TREE_TOTAL = sum(tEst, na.rm = TRUE),
# BA_TOTAL = sum(bEst, na.rm = TRUE),
# ## Variances
# treeVar = sum(tVar, na.rm = TRUE),
# baVar = sum(bVar, na.rm = TRUE),
# #aVar = first(aVar),
# cvT = sum(cvEst_t, na.rm = TRUE),
# cvB = sum(cvEst_b, na.rm = TRUE),
# ## Sampling Errors
# TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
# BA_SE = sqrt(baVar) / BA_TOTAL * 100,
# nPlots_TREE = sum(plotIn_TREE, na.rm = TRUE)) #%>%
tTotal <- tEst %>%
group_by(.dots = grpBy) %>%
summarize_all(sum, na.rm = TRUE) %>%
mutate(TREE_TOTAL = tEst,
BA_TOTAL = bEst,
## Variances
treeVar = tVar,
baVar = bVar,
#aVar = first(aVar),
cvT = cvEst_t,
cvB = cvEst_b,
## Sampling Errors
TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
BA_SE = sqrt(baVar) / BA_TOTAL * 100,
nPlots_TREE = plotIn_TREE) %>%
select(grpBy, TREE_TOTAL, BA_TOTAL, treeVar, baVar, cvT, cvB, TREE_SE, BA_SE, nPlots_TREE)
## IF using polys, we treat each zone as a unique population
if (!is.null(polys)){
propGrp <- c('polyID', grpBy)
} else {
propGrp <- 'YEAR'
}
## Hand the proportions
tpTotal <- tEst %>%
group_by(.dots = unique(propGrp)) %>%
summarize(TREE_TOTAL_full = sum(tTEst, na.rm = TRUE), ## Need to sum this
BA_TOTAL_full = sum(bTEst, na.rm = TRUE), ## Need to sum this
tTVar = sum(tTVar, na.rm = TRUE),
bTVar = sum(bTVar, na.rm = TRUE),
cvTT = sum(cvEst_tT, na.rm = TRUE),
cvBT = sum(cvEst_bT, na.rm = TRUE))
suppressWarnings({
## Bring them together
tTotal <- tTotal %>%
left_join(aTotal, by = aGrpBy) %>%
left_join(tpTotal, by = unique(propGrp)) %>%
mutate(TPA = TREE_TOTAL / AREA_TOTAL,
BAA = BA_TOTAL / AREA_TOTAL,
tpaVar = (1/AREA_TOTAL^2) * (treeVar + (TPA^2 * aVar) - 2 * TPA * cvT),
baaVar = (1/AREA_TOTAL^2) * (baVar + (BAA^2 * aVar) - (2 * BAA * cvB)),
TPA_SE = sqrt(tpaVar) / TPA * 100,
BAA_SE = sqrt(baaVar) / BAA * 100,
TPA_PERC = TREE_TOTAL / (TREE_TOTAL_full) * 100,
BAA_PERC = BA_TOTAL / (BA_TOTAL_full) * 100,
tpVar = (1/TREE_TOTAL_full^2) * (treeVar + (TPA_PERC^2 * tTVar) - 2 * TPA_PERC * cvTT),
bpVar = (1/BA_TOTAL_full^2) * (baVar + (BAA_PERC^2 * bTVar) - (2 * BAA_PERC * cvBT)),
TPA_PERC_SE = sqrt(tpVar) / TPA_PERC * 100,
BAA_PERC_SE = sqrt(bpVar) / BAA_PERC * 100)
})
if (totals) {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TREE_TOTAL, BA_TOTAL, AREA_TOTAL, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, TREE_SE, BA_SE, AREA_TOTAL_SE, nPlots_TREE, nPlots_AREA)
} else {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, nPlots_TREE, nPlots_AREA)
}
# 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)
}
tpa_backup <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
bySpecies = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
treeDomain = NULL,
areaDomain = NULL,
totals = FALSE,
byPlot = 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 = '_'))
## 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).'))
}
# 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', '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') %>%
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')
## Add species to groups
if (bySpecies) {
db$TREE <- db$TREE %>%
left_join(select(intData$REF_SPECIES_2018, c('SPCD','COMMON_NAME', 'GENUS', 'SPECIES')), by = 'SPCD') %>%
mutate(SCIENTIFIC_NAME = paste(GENUS, SPECIES, sep = ' ')) %>%
mutate_if(is.factor,
as.character)
grpBy <- c(grpBy, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
grpByOrig <- c(grpByOrig, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
}
## 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),]
}
# Seperate area grouping names, (ex. TPA red oak in total land area of ingham county, rather than only area where red oak occurs)
if (!is.null(polys)){
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND) | grpBy %in% names(pltSF)])
} else {
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND)])
}
## 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 = tpaHelper1, plts, db, grpBy, aGrpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = tpaHelper1, plts, db, grpBy, aGrpBy, 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)
a <- bind_rows(out[names(out) == 'a'])
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
aGrpBy <- c('YEAR', aGrpBy)
## 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 = tpaHelper2, popState, a, t, grpBy, aGrpBy, method, lambda)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = tpaHelper2, popState, a, t, grpBy, aGrpBy, method, lambda, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
aEst <- bind_rows(out[names(out) == 'aEst'])
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 = 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 = 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
# Area
aEst <- aEst %>%
mutate_at(vars(aEst), ~(.*wgt)) %>%
mutate_at(vars(aVar), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = aGrpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
tEst <- tEst %>%
mutate_at(vars(tEst:bTEst), ~(.*wgt)) %>%
mutate_at(vars(tVar:cvEst_bT), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(tEst:cvEst_bT), sum, na.rm = TRUE)
}
##--------------------- TOTALS and RATIOS
# Area
aTotal <- aEst %>%
group_by(.dots = aGrpBy) %>%
summarize(AREA_TOTAL = sum(aEst, na.rm = TRUE),
aVar = sum(aVar, na.rm = TRUE),
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL * 100,
nPlots_AREA = sum(plotIn_AREA, na.rm = TRUE))
# # Tree
# tTotal <- tEst %>%
# group_by(.dots = grpBy) %>%
# #left_join(aTotal, by = c(aGrpBy)) %>%
# summarize(TREE_TOTAL = sum(tEst, na.rm = TRUE),
# BA_TOTAL = sum(bEst, na.rm = TRUE),
# ## Variances
# treeVar = sum(tVar, na.rm = TRUE),
# baVar = sum(bVar, na.rm = TRUE),
# #aVar = first(aVar),
# cvT = sum(cvEst_t, na.rm = TRUE),
# cvB = sum(cvEst_b, na.rm = TRUE),
# ## Sampling Errors
# TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
# BA_SE = sqrt(baVar) / BA_TOTAL * 100,
# nPlots_TREE = sum(plotIn_TREE, na.rm = TRUE)) #%>%
tTotal <- tEst %>%
group_by(.dots = grpBy) %>%
summarize_all(sum, na.rm = TRUE) %>%
mutate(TREE_TOTAL = tEst,
BA_TOTAL = bEst,
## Variances
treeVar = tVar,
baVar = bVar,
#aVar = first(aVar),
cvT = cvEst_t,
cvB = cvEst_b,
## Sampling Errors
TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
BA_SE = sqrt(baVar) / BA_TOTAL * 100,
nPlots_TREE = plotIn_TREE) %>%
select(grpBy, TREE_TOTAL, BA_TOTAL, treeVar, baVar, cvT, cvB, TREE_SE, BA_SE, nPlots_TREE)
## IF using polys, we treat each zone as a unique population
if (!is.null(polys)){
propGrp <- c('polyID', grpBy)
} else {
propGrp <- 'YEAR'
}
## Hand the proportions
tpTotal <- tEst %>%
group_by(.dots = unique(propGrp)) %>%
summarize(TREE_TOTAL_full = sum(tTEst, na.rm = TRUE), ## Need to sum this
BA_TOTAL_full = sum(bTEst, na.rm = TRUE), ## Need to sum this
tTVar = sum(tTVar, na.rm = TRUE),
bTVar = sum(bTVar, na.rm = TRUE),
cvTT = sum(cvEst_tT, na.rm = TRUE),
cvBT = sum(cvEst_bT, na.rm = TRUE))
suppressWarnings({
## Bring them together
tTotal <- tTotal %>%
left_join(aTotal, by = aGrpBy) %>%
left_join(tpTotal, by = unique(propGrp)) %>%
mutate(TPA = TREE_TOTAL / AREA_TOTAL,
BAA = BA_TOTAL / AREA_TOTAL,
tpaVar = (1/AREA_TOTAL^2) * (treeVar + (TPA^2 * aVar) - 2 * TPA * cvT),
baaVar = (1/AREA_TOTAL^2) * (baVar + (BAA^2 * aVar) - (2 * BAA * cvB)),
TPA_SE = sqrt(tpaVar) / TPA * 100,
BAA_SE = sqrt(baaVar) / BAA * 100,
TPA_PERC = TREE_TOTAL / (TREE_TOTAL_full) * 100,
BAA_PERC = BA_TOTAL / (BA_TOTAL_full) * 100,
tpVar = (1/TREE_TOTAL_full^2) * (treeVar + (TPA_PERC^2 * tTVar) - 2 * TPA_PERC * cvTT),
bpVar = (1/BA_TOTAL_full^2) * (baVar + (BAA_PERC^2 * bTVar) - (2 * BAA_PERC * cvBT)),
TPA_PERC_SE = sqrt(tpVar) / TPA_PERC * 100,
BAA_PERC_SE = sqrt(bpVar) / BAA_PERC * 100)
})
if (totals) {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TREE_TOTAL, BA_TOTAL, AREA_TOTAL, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, TREE_SE, BA_SE, AREA_TOTAL_SE, nPlots_TREE, nPlots_AREA)
} else {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, nPlots_TREE, nPlots_AREA)
}
# 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 tpa_backup tpa_test tpaStarter_test tpa tpaStarter
Documented in tpa
tpaStarter <- function(x,
db,
grpBy_quo = NULL,
polys = NULL,
returnSpatial = FALSE,
bySpecies = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
treeDomain = NULL,
areaDomain = NULL,
totals = FALSE,
byPlot = 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]
}
## 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)
library(sf)
})
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')
## 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') %>%
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')
## Add species to groups
if (bySpecies) {
db$TREE <- db$TREE %>%
left_join(select(intData$REF_SPECIES_2018, c('SPCD','COMMON_NAME', 'GENUS', 'SPECIES')), by = 'SPCD') %>%
mutate(SCIENTIFIC_NAME = paste(GENUS, SPECIES, sep = ' ')) %>%
mutate_if(is.factor,
as.character)
grpBy <- c(grpBy, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
grpByOrig <- c(grpByOrig, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
}
## 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),]
}
# Seperate area grouping names, (ex. TPA red oak in total land area of ingham county, rather than only area where red oak occurs)
if (!is.null(polys)){
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND) | grpBy %in% names(pltSF)])
} else {
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND)])
}
## Only the necessary plots for EVAL of interest
db$PLOT <- filter(db$PLOT, PLT_CN %in% pops$PLT_CN)
## Narrow up the tables to the necessary variables, reduces memory load
## sent out the cores
## 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]
### Only joining tables necessary to produce plot level estimates, adjusted for non-response
db$PLOT <- select(db$PLOT, c('PLT_CN', 'STATECD', 'MACRO_BREAKPOINT_DIA', 'INVYR', 'MEASYEAR', 'PLOT_STATUS_CD', grpP, 'aD_p', 'sp'))
db$COND <- select(db$COND, c('PLT_CN', 'CONDPROP_UNADJ', 'PROP_BASIS', 'COND_STATUS_CD', 'CONDID', grpC, 'aD_c', 'landD')) %>%
filter(PLT_CN %in% db$PLOT$PLT_CN)
db$TREE <- select(db$TREE, c('PLT_CN', 'CONDID', 'DIA', 'SPCD', 'TPA_UNADJ', 'SUBP', 'TREE', grpT, 'tD', 'typeD')) %>%
filter(PLT_CN %in% db$PLOT$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 = tpaHelper1, plts,
db[names(db) %in% c('COND', 'TREE')],
grpBy, aGrpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = tpaHelper1, plts,
db[names(db) %in% c('COND', 'TREE')],
grpBy, aGrpBy, 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, aGrpBy = aGrpBy, grpByOrig = grpByOrig)
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
a <- bind_rows(out[names(out) == 'a'])
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
aGrpBy <- c('YEAR', aGrpBy)
## 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 = tpaHelper2, popState, a, t, grpBy, aGrpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = tpaHelper2, popState, a, t, grpBy, aGrpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
aEst <- bind_rows(out[names(out) == 'aEst'])
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
# Area
aEst <- aEst %>%
mutate_at(vars(aEst), ~(.*wgt)) %>%
mutate_at(vars(aVar), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = aGrpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
tEst <- tEst %>%
mutate_at(vars(tEst:bTEst), ~(.*wgt)) %>%
mutate_at(vars(tVar:cvEst_bT), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(tEst:cvEst_bT), sum, na.rm = TRUE)
}
out <- list(tEst = tEst, aEst = aEst, grpBy = grpBy, aGrpBy = aGrpBy, grpByOrig = grpByOrig)
}
return(out)
}
#' @export
tpa <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
bySpecies = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
treeDomain = NULL,
areaDomain = NULL,
totals = FALSE,
variance = FALSE,
byPlot = FALSE,
nCores = 1) {
## don't have to change original code
grpBy_quo <- rlang::enquo(grpBy)
areaDomain <- rlang::enquo(areaDomain)
treeDomain <- rlang::enquo(treeDomain)
### 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 = tpaStarter, db,
grpBy_quo = grpBy_quo, polys, returnSpatial,
bySpecies, bySizeClass,
landType, treeType, method,
lambda, treeDomain, areaDomain,
totals, byPlot, nCores, remote, mr)
## Bring the results back
out <- unlist(out, recursive = FALSE)
aEst <- bind_rows(out[names(out) == 'aEst'])
tEst <- bind_rows(out[names(out) == 'tEst'])
grpBy <- out[names(out) == 'grpBy'][[1]]
aGrpBy <- out[names(out) == 'aGrpBy'][[1]]
grpByOrig <- out[names(out) == 'grpByOrig'][[1]]
if (byPlot){
tOut <- tEst
## Population estimates
} 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))
maxyearsA <- aEst %>%
select(aGrpBy) %>%
group_by(.dots = aGrpBy[!c(aGrpBy %in% 'YEAR')]) %>%
summarise(YEAR = max(YEAR, na.rm = TRUE))
# Combine estimates
aEst <- aEst %>%
ungroup() %>%
select(-c(YEAR)) %>%
left_join(maxyearsA, by = aGrpBy[!c(aGrpBy %in% 'YEAR')])
tEst <- tEst %>%
ungroup() %>%
select(-c(YEAR)) %>%
left_join(maxyearsT, by = grpBy[!c(grpBy %in% 'YEAR')])
}
})})
##--------------------- TOTALS and RATIOS
# Area
aTotal <- aEst %>%
group_by(.dots = aGrpBy) %>%
summarize(AREA_TOTAL = sum(aEst, na.rm = TRUE),
aVar = sum(aVar, na.rm = TRUE),
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL * 100,
AREA_TOTAL_VAR = aVar,
nPlots_AREA = sum(plotIn_AREA, na.rm = TRUE))
# # Tree
# tTotal <- tEst %>%
# group_by(.dots = grpBy) %>%
# #left_join(aTotal, by = c(aGrpBy)) %>%
# summarize(TREE_TOTAL = sum(tEst, na.rm = TRUE),
# BA_TOTAL = sum(bEst, na.rm = TRUE),
# ## Variances
# treeVar = sum(tVar, na.rm = TRUE),
# baVar = sum(bVar, na.rm = TRUE),
# #aVar = first(aVar),
# cvT = sum(cvEst_t, na.rm = TRUE),
# cvB = sum(cvEst_b, na.rm = TRUE),
# ## Sampling Errors
# TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
# BA_SE = sqrt(baVar) / BA_TOTAL * 100,
# nPlots_TREE = sum(plotIn_TREE, na.rm = TRUE)) #%>%
tTotal <- tEst %>%
group_by(.dots = grpBy) %>%
summarize_all(sum, na.rm = TRUE) %>%
mutate(TREE_TOTAL = tEst,
BA_TOTAL = bEst,
## Variances
treeVar = tVar,
baVar = bVar,
#aVar = first(aVar),
cvT = cvEst_t,
cvB = cvEst_b,
## Sampling Errors
TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
TREE_VAR = treeVar,
BA_SE = sqrt(baVar) / BA_TOTAL * 100,
BA_VAR = baVar,
N = sum(N, na.rm = TRUE),
nPlots_TREE = plotIn_TREE) %>%
select(grpBy, TREE_TOTAL, BA_TOTAL, treeVar, baVar, cvT, cvB, TREE_SE, BA_SE, TREE_VAR, BA_VAR,
nPlots_TREE, N)
## IF using polys, we treat each zone as a unique population
if (!is.null(polys)){
propGrp <- c('polyID', grpBy)
} else {
propGrp <- 'YEAR'
}
## Hand the proportions
tpTotal <- tEst %>%
group_by(.dots = unique(propGrp)) %>%
summarize(TREE_TOTAL_full = sum(tTEst, na.rm = TRUE), ## Need to sum this
BA_TOTAL_full = sum(bTEst, na.rm = TRUE), ## Need to sum this
tTVar = sum(tTVar, na.rm = TRUE),
bTVar = sum(bTVar, na.rm = TRUE),
cvTT = sum(cvEst_tT, na.rm = TRUE),
cvBT = sum(cvEst_bT, na.rm = TRUE))
suppressWarnings({
## Bring them together
tTotal <- tTotal %>%
left_join(aTotal, by = aGrpBy) %>%
left_join(tpTotal, by = unique(propGrp)) %>%
mutate(TPA = TREE_TOTAL / AREA_TOTAL,
BAA = BA_TOTAL / AREA_TOTAL,
tpaVar = (1/AREA_TOTAL^2) * (treeVar + (TPA^2 * aVar) - 2 * TPA * cvT),
baaVar = (1/AREA_TOTAL^2) * (baVar + (BAA^2 * aVar) - (2 * BAA * cvB)),
TPA_SE = sqrt(tpaVar) / TPA * 100,
BAA_SE = sqrt(baaVar) / BAA * 100,
TPA_VAR = tpaVar,
BAA_VAR = baaVar,
TPA_PERC = TREE_TOTAL / (TREE_TOTAL_full) * 100,
BAA_PERC = BA_TOTAL / (BA_TOTAL_full) * 100,
tpVar = (1/TREE_TOTAL_full^2) * (treeVar + (TPA_PERC^2 * tTVar) - 2 * TPA_PERC * cvTT),
bpVar = (1/BA_TOTAL_full^2) * (baVar + (BAA_PERC^2 * bTVar) - (2 * BAA_PERC * cvBT)),
TPA_PERC_SE = sqrt(tpVar) / TPA_PERC * 100,
BAA_PERC_SE = sqrt(bpVar) / BAA_PERC * 100,
TPA_PERC_VAR = tpVar,
BAA_PERC_VAR = bpVar)
})
if (totals) {
if (variance){
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TREE_TOTAL, BA_TOTAL, AREA_TOTAL, TPA_VAR, BAA_VAR,
TPA_PERC_VAR, BAA_PERC_VAR, TREE_VAR, BA_VAR, AREA_TOTAL_VAR, nPlots_TREE, nPlots_AREA, N)
} else {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TREE_TOTAL, BA_TOTAL, AREA_TOTAL, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, TREE_SE, BA_SE, AREA_TOTAL_SE, nPlots_TREE, nPlots_AREA)
}
} else {
if (variance) {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TPA_VAR, BAA_VAR,
TPA_PERC_VAR, BAA_PERC_VAR, nPlots_TREE, nPlots_AREA, N)
} else {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, nPlots_TREE, nPlots_AREA)
}
}
# 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_AREA))})})
## 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)
}
tpaStarter_test <- function(x,
db,
grpBy_quo = NULL,
polys = NULL,
returnSpatial = FALSE,
bySpecies = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
treeDomain = NULL,
areaDomain = NULL,
totals = FALSE,
byPlot = FALSE,
nCores = 1,
remote){
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]
}
## 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)
}
### 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(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') %>%
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', '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') %>%
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')
## Add species to groups
if (bySpecies) {
db$TREE <- db$TREE %>%
left_join(select(intData$REF_SPECIES_2018, c('SPCD','COMMON_NAME', 'GENUS', 'SPECIES')), by = 'SPCD') %>%
mutate(SCIENTIFIC_NAME = paste(GENUS, SPECIES, sep = ' ')) %>%
mutate_if(is.factor,
as.character)
grpBy <- c(grpBy, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
grpByOrig <- c(grpByOrig, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
}
## 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),]
}
# Seperate area grouping names, (ex. TPA red oak in total land area of ingham county, rather than only area where red oak occurs)
if (!is.null(polys)){
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND) | grpBy %in% names(pltSF)])
} else {
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND)])
}
## Only the necessary plots for EVAL of interest
db$PLOT <- filter(db$PLOT, PLT_CN %in% pops$PLT_CN)
## Narrow up the tables to the necessary variables, reduces memory load
## sent out the cores
## 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]
### Only joining tables necessary to produce plot level estimates, adjusted for non-response
db$PLOT <- select(db$PLOT, c('PLT_CN', 'STATECD', 'MACRO_BREAKPOINT_DIA', 'INVYR', 'MEASYEAR', 'PLOT_STATUS_CD', grpP, 'aD_p', 'sp'))
db$COND <- select(db$COND, c('PLT_CN', 'CONDPROP_UNADJ', 'PROP_BASIS', 'COND_STATUS_CD', 'CONDID', grpC, 'aD_c', 'landD'))
db$TREE <- select(db$TREE, c('PLT_CN', 'CONDID', 'DIA', 'SPCD', 'TPA_UNADJ', 'SUBP', 'TREE', 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 = tpaHelper1_test, plts, db, grpBy, aGrpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = tpaHelper1_test, plts, db, grpBy, aGrpBy, 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, aGrpBy = aGrpBy, grpByOrig = grpByOrig)
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
a <- bind_rows(out[names(out) == 'a'])
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
aGrpBy <- c('YEAR', aGrpBy)
## 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 = tpaHelper2_test, popState, a, t, grpBy, aGrpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = tpaHelper2_test, popState, a, t, grpBy, aGrpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
aEst <- bind_rows(out[names(out) == 'aEst'])
tEst <- bind_rows(out[names(out) == 'tEst'])
out <- list(tEst = tEst, aEst = aEst, grpBy = grpBy, aGrpBy = aGrpBy, grpByOrig = grpByOrig)
}
return(out)
}
tpa_test <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
bySpecies = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
treeDomain = NULL,
areaDomain = NULL,
totals = FALSE,
byPlot = FALSE,
nCores = 1) {
## don't have to change original code
grpBy_quo <- rlang::enquo(grpBy)
areaDomain <- rlang::enquo(areaDomain)
treeDomain <- rlang::enquo(treeDomain)
### 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
}
### 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 = tpaStarter_test, db,
grpBy_quo = grpBy_quo, polys, returnSpatial,
bySpecies, bySizeClass,
landType, treeType, method,
lambda, treeDomain, areaDomain,
totals, byPlot, nCores, remote)
## Bring the results back
out <- unlist(out, recursive = FALSE)
aEst <- bind_rows(out[names(out) == 'aEst'])
tEst <- bind_rows(out[names(out) == 'tEst'])
grpBy <- out[names(out) == 'grpBy'][[1]]
aGrpBy <- out[names(out) == 'aGrpBy'][[1]]
grpByOrig <- out[names(out) == 'grpByOrig'][[1]]
if (byPlot){
tOut <- tEst
## Population estimates
} else {
##--------------------- TOTALS and RATIOS
# Area
aTotal <- aEst %>%
group_by(.dots = aGrpBy) %>%
summarize(AREA_TOTAL = sum(aEst, na.rm = TRUE),
aVar = sum(aVar, na.rm = TRUE),
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL * 100,
nPlots_AREA = sum(plotIn_AREA, na.rm = TRUE))
# # Tree
# tTotal <- tEst %>%
# group_by(.dots = grpBy) %>%
# #left_join(aTotal, by = c(aGrpBy)) %>%
# summarize(TREE_TOTAL = sum(tEst, na.rm = TRUE),
# BA_TOTAL = sum(bEst, na.rm = TRUE),
# ## Variances
# treeVar = sum(tVar, na.rm = TRUE),
# baVar = sum(bVar, na.rm = TRUE),
# #aVar = first(aVar),
# cvT = sum(cvEst_t, na.rm = TRUE),
# cvB = sum(cvEst_b, na.rm = TRUE),
# ## Sampling Errors
# TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
# BA_SE = sqrt(baVar) / BA_TOTAL * 100,
# nPlots_TREE = sum(plotIn_TREE, na.rm = TRUE)) #%>%
tTotal <- tEst %>%
group_by(.dots = grpBy) %>%
summarize_all(sum, na.rm = TRUE) %>%
mutate(TREE_TOTAL = tEst,
BA_TOTAL = bEst,
## Variances
treeVar = tVar,
baVar = bVar,
#aVar = first(aVar),
cvT = cvEst_t,
cvB = cvEst_b,
## Sampling Errors
TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
BA_SE = sqrt(baVar) / BA_TOTAL * 100,
nPlots_TREE = plotIn_TREE) %>%
select(grpBy, TREE_TOTAL, BA_TOTAL, treeVar, baVar, cvT, cvB, TREE_SE, BA_SE, nPlots_TREE)
## IF using polys, we treat each zone as a unique population
if (!is.null(polys)){
propGrp <- c('polyID', grpBy)
} else {
propGrp <- 'YEAR'
}
## Hand the proportions
tpTotal <- tEst %>%
group_by(.dots = unique(propGrp)) %>%
summarize(TREE_TOTAL_full = sum(tTEst, na.rm = TRUE), ## Need to sum this
BA_TOTAL_full = sum(bTEst, na.rm = TRUE), ## Need to sum this
tTVar = sum(tTVar, na.rm = TRUE),
bTVar = sum(bTVar, na.rm = TRUE),
cvTT = sum(cvEst_tT, na.rm = TRUE),
cvBT = sum(cvEst_bT, na.rm = TRUE))
suppressWarnings({
## Bring them together
tTotal <- tTotal %>%
left_join(aTotal, by = aGrpBy) %>%
left_join(tpTotal, by = unique(propGrp)) %>%
mutate(TPA = TREE_TOTAL / AREA_TOTAL,
BAA = BA_TOTAL / AREA_TOTAL,
tpaVar = (1/AREA_TOTAL^2) * (treeVar + (TPA^2 * aVar) - 2 * TPA * cvT),
baaVar = (1/AREA_TOTAL^2) * (baVar + (BAA^2 * aVar) - (2 * BAA * cvB)),
TPA_SE = sqrt(tpaVar) / TPA * 100,
BAA_SE = sqrt(baaVar) / BAA * 100,
TPA_PERC = TREE_TOTAL / (TREE_TOTAL_full) * 100,
BAA_PERC = BA_TOTAL / (BA_TOTAL_full) * 100,
tpVar = (1/TREE_TOTAL_full^2) * (treeVar + (TPA_PERC^2 * tTVar) - 2 * TPA_PERC * cvTT),
bpVar = (1/BA_TOTAL_full^2) * (baVar + (BAA_PERC^2 * bTVar) - (2 * BAA_PERC * cvBT)),
TPA_PERC_SE = sqrt(tpVar) / TPA_PERC * 100,
BAA_PERC_SE = sqrt(bpVar) / BAA_PERC * 100)
})
if (totals) {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TREE_TOTAL, BA_TOTAL, AREA_TOTAL, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, TREE_SE, BA_SE, AREA_TOTAL_SE, nPlots_TREE, nPlots_AREA)
} else {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, nPlots_TREE, nPlots_AREA)
}
# 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)
}
tpa_backup <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
bySpecies = FALSE,
bySizeClass = FALSE,
landType = 'forest',
treeType = 'live',
method = 'TI',
lambda = .5,
treeDomain = NULL,
areaDomain = NULL,
totals = FALSE,
byPlot = 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 = '_'))
## 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).'))
}
# 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', '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') %>%
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')
## Add species to groups
if (bySpecies) {
db$TREE <- db$TREE %>%
left_join(select(intData$REF_SPECIES_2018, c('SPCD','COMMON_NAME', 'GENUS', 'SPECIES')), by = 'SPCD') %>%
mutate(SCIENTIFIC_NAME = paste(GENUS, SPECIES, sep = ' ')) %>%
mutate_if(is.factor,
as.character)
grpBy <- c(grpBy, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
grpByOrig <- c(grpByOrig, 'SPCD', 'COMMON_NAME', 'SCIENTIFIC_NAME')
}
## 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),]
}
# Seperate area grouping names, (ex. TPA red oak in total land area of ingham county, rather than only area where red oak occurs)
if (!is.null(polys)){
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND) | grpBy %in% names(pltSF)])
} else {
aGrpBy <- c(grpBy[grpBy %in% names(db$PLOT) | grpBy %in% names(db$COND)])
}
## 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 = tpaHelper1, plts, db, grpBy, aGrpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = tpaHelper1, plts, db, grpBy, aGrpBy, 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)
a <- bind_rows(out[names(out) == 'a'])
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
aGrpBy <- c('YEAR', aGrpBy)
## 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 = tpaHelper2, popState, a, t, grpBy, aGrpBy, method, lambda)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = tpaHelper2, popState, a, t, grpBy, aGrpBy, method, lambda, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
aEst <- bind_rows(out[names(out) == 'aEst'])
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 = 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 = 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
# Area
aEst <- aEst %>%
mutate_at(vars(aEst), ~(.*wgt)) %>%
mutate_at(vars(aVar), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = aGrpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
tEst <- tEst %>%
mutate_at(vars(tEst:bTEst), ~(.*wgt)) %>%
mutate_at(vars(tVar:cvEst_bT), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(tEst:cvEst_bT), sum, na.rm = TRUE)
}
##--------------------- TOTALS and RATIOS
# Area
aTotal <- aEst %>%
group_by(.dots = aGrpBy) %>%
summarize(AREA_TOTAL = sum(aEst, na.rm = TRUE),
aVar = sum(aVar, na.rm = TRUE),
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL * 100,
nPlots_AREA = sum(plotIn_AREA, na.rm = TRUE))
# # Tree
# tTotal <- tEst %>%
# group_by(.dots = grpBy) %>%
# #left_join(aTotal, by = c(aGrpBy)) %>%
# summarize(TREE_TOTAL = sum(tEst, na.rm = TRUE),
# BA_TOTAL = sum(bEst, na.rm = TRUE),
# ## Variances
# treeVar = sum(tVar, na.rm = TRUE),
# baVar = sum(bVar, na.rm = TRUE),
# #aVar = first(aVar),
# cvT = sum(cvEst_t, na.rm = TRUE),
# cvB = sum(cvEst_b, na.rm = TRUE),
# ## Sampling Errors
# TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
# BA_SE = sqrt(baVar) / BA_TOTAL * 100,
# nPlots_TREE = sum(plotIn_TREE, na.rm = TRUE)) #%>%
tTotal <- tEst %>%
group_by(.dots = grpBy) %>%
summarize_all(sum, na.rm = TRUE) %>%
mutate(TREE_TOTAL = tEst,
BA_TOTAL = bEst,
## Variances
treeVar = tVar,
baVar = bVar,
#aVar = first(aVar),
cvT = cvEst_t,
cvB = cvEst_b,
## Sampling Errors
TREE_SE = sqrt(treeVar) / TREE_TOTAL * 100,
BA_SE = sqrt(baVar) / BA_TOTAL * 100,
nPlots_TREE = plotIn_TREE) %>%
select(grpBy, TREE_TOTAL, BA_TOTAL, treeVar, baVar, cvT, cvB, TREE_SE, BA_SE, nPlots_TREE)
## IF using polys, we treat each zone as a unique population
if (!is.null(polys)){
propGrp <- c('polyID', grpBy)
} else {
propGrp <- 'YEAR'
}
## Hand the proportions
tpTotal <- tEst %>%
group_by(.dots = unique(propGrp)) %>%
summarize(TREE_TOTAL_full = sum(tTEst, na.rm = TRUE), ## Need to sum this
BA_TOTAL_full = sum(bTEst, na.rm = TRUE), ## Need to sum this
tTVar = sum(tTVar, na.rm = TRUE),
bTVar = sum(bTVar, na.rm = TRUE),
cvTT = sum(cvEst_tT, na.rm = TRUE),
cvBT = sum(cvEst_bT, na.rm = TRUE))
suppressWarnings({
## Bring them together
tTotal <- tTotal %>%
left_join(aTotal, by = aGrpBy) %>%
left_join(tpTotal, by = unique(propGrp)) %>%
mutate(TPA = TREE_TOTAL / AREA_TOTAL,
BAA = BA_TOTAL / AREA_TOTAL,
tpaVar = (1/AREA_TOTAL^2) * (treeVar + (TPA^2 * aVar) - 2 * TPA * cvT),
baaVar = (1/AREA_TOTAL^2) * (baVar + (BAA^2 * aVar) - (2 * BAA * cvB)),
TPA_SE = sqrt(tpaVar) / TPA * 100,
BAA_SE = sqrt(baaVar) / BAA * 100,
TPA_PERC = TREE_TOTAL / (TREE_TOTAL_full) * 100,
BAA_PERC = BA_TOTAL / (BA_TOTAL_full) * 100,
tpVar = (1/TREE_TOTAL_full^2) * (treeVar + (TPA_PERC^2 * tTVar) - 2 * TPA_PERC * cvTT),
bpVar = (1/BA_TOTAL_full^2) * (baVar + (BAA_PERC^2 * bTVar) - (2 * BAA_PERC * cvBT)),
TPA_PERC_SE = sqrt(tpVar) / TPA_PERC * 100,
BAA_PERC_SE = sqrt(bpVar) / BAA_PERC * 100)
})
if (totals) {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TREE_TOTAL, BA_TOTAL, AREA_TOTAL, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, TREE_SE, BA_SE, AREA_TOTAL_SE, nPlots_TREE, nPlots_AREA)
} else {
tOut <- tTotal %>%
select(grpBy, TPA, BAA, TPA_PERC, BAA_PERC, TPA_SE, BAA_SE,
TPA_PERC_SE, BAA_PERC_SE, nPlots_TREE, nPlots_AREA)
}
# 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.