R/standStructNew.R
In hunter-stanke/rFIA_TX: Duplicate rFIA for TX
Defines functions
standStruct_backup
standStruct
standStructStarter
Documented in
standStruct
standStructStarter <- function(x,
db,
grpBy_quo = NULL,
polys = NULL,
returnSpatial = FALSE,
landType = 'forest',
method = 'TI',
lambda = .5,
areaDomain = NULL,
totals = FALSE,
byPlot = FALSE,
tidy = TRUE,
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') %>%
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 or COND tables. Did you accidentally quote the variables names? e.g. use grpBy = ECOSUBCD (correct) instead of grpBy = "ECOSUBCD". ', collapse = ', '))
} else {
# Convert to character
grpBy <- names(d_quo)
}
} else {
grpBy <- NULL
}
reqTables <- c('PLOT', 'TREE', 'COND', 'POP_PLOT_STRATUM_ASSGN', 'POP_ESTN_UNIT', 'POP_EVAL',
'POP_STRATUM', 'POP_EVAL_TYP', 'POP_EVAL_GRP')
if (!is.null(polys) & first(class(polys)) %in% c('sf', 'SpatialPolygons', 'SpatialPolygonsDataFrame') == FALSE){
stop('polys must be spatial polygons object of class sp or sf. ')
}
if (landType %in% c('timber', 'forest') == FALSE){
stop('landType must be one of: "forest" or "timber".')
}
if (any(reqTables %in% names(db) == FALSE)){
missT <- reqTables[reqTables %in% names(db) == FALSE]
stop(paste('Tables', paste (as.character(missT), collapse = ', '), 'not found in object db.'))
}
if (str_to_upper(method) %in% c('TI', 'SMA', 'LMA', 'EMA', 'ANNUAL') == FALSE) {
warning(paste('Method', method, 'unknown. Defaulting to Temporally Indifferent (TI).'))
}
# I like a unique ID for a plot through time
if (byPlot) {grpBy <- c('pltID', grpBy)}
# Save original grpBy for pretty return with spatial objects
grpByOrig <- grpBy
## IF the object was clipped
if ('prev' %in% names(db$PLOT)){
## Only want the current plots, no grm
db$PLOT <- filter(db$PLOT, prev == 0)
}
### AREAL SUMMARY PREP
if(!is.null(polys)) {
# # Convert polygons to an sf object
# polys <- polys %>%
# as('sf')%>%
# mutate_if(is.factor,
# as.character)
# ## A unique ID
# polys$polyID <- 1:nrow(polys)
#
# # Add shapefile names to grpBy
grpBy = c(grpBy, 'polyID')
## Make plot data spatial, projected same as polygon layer
pltSF <- select(db$PLOT, c('LON', 'LAT', pltID)) %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
distinct(pltID, .keep_all = TRUE)
coordinates(pltSF) <- ~LON+LAT
proj4string(pltSF) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
pltSF <- as(pltSF, 'sf') %>%
st_transform(crs = st_crs(polys))
## Split up polys
polyList <- split(polys, as.factor(polys$polyID))
suppressWarnings({suppressMessages({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
cl <- makeCluster(nCores)
clusterEvalQ(cl, {
library(dplyr)
library(stringr)
library(rFIA)
})
out <- parLapply(cl, X = names(polyList), fun = areal_par, pltSF, polyList)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(polyList), FUN = areal_par, pltSF, polyList, mc.cores = nCores)
}
})})
pltSF <- bind_rows(out)
# A warning
if (length(unique(pltSF$pltID)) < 1){
stop('No plots in db overlap with polys.')
}
## Add polygon names to PLOT
db$PLOT <- db$PLOT %>%
left_join(select(pltSF, polyID, pltID), by = 'pltID')
# Test if any polygons cross state boundaries w/ different recent inventory years (continued w/in loop)
if ('mostRecent' %in% names(db) & length(unique(db$POP_EVAL$STATECD)) > 1){
mergeYears <- pltSF %>%
right_join(select(db$PLOT, PLT_CN, pltID), by = 'pltID') %>%
inner_join(select(db$POP_PLOT_STRATUM_ASSGN, c('PLT_CN', 'EVALID', 'STATECD')), by = 'PLT_CN') %>%
inner_join(select(db$POP_EVAL, c('EVALID', 'END_INVYR')), by = 'EVALID') %>%
group_by(polyID) %>%
summarize(maxYear = max(END_INVYR, na.rm = TRUE))
}
## TO RETURN SPATIAL PLOTS
}
if (byPlot & returnSpatial){
grpBy <- c(grpBy, 'LON', 'LAT')
} # END AREAL
## Build domain indicator function which is 1 if observation meets criteria, and 0 otherwise
# Land type domain indicator
if (tolower(landType) == 'forest'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1, 1, 0)
} else if (tolower(landType) == 'timber'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1 & db$COND$SITECLCD %in% c(1, 2, 3, 4, 5, 6) & db$COND$RESERVCD == 0, 1, 0)
}
# 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)
### 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 == 'EXPCURR') %>%
filter(!is.na(END_INVYR) & !is.na(EVALID) & END_INVYR >= 2003) %>%
distinct(END_INVYR, EVALID, .keep_all = TRUE)# %>%
#group_by(END_INVYR) %>%
#summarise(id = list(EVALID)
## If a most-recent subset, make sure that we don't get two reporting years in
## western states
if (mr) {
db$POP_EVAL <- db$POP_EVAL %>%
group_by(EVAL_TYP, STATECD) %>%
filter(END_INVYR == max(END_INVYR, na.rm = TRUE)) %>%
ungroup()
}
## Cut STATECD
db$POP_EVAL <- select(db$POP_EVAL, -c(STATECD))
### The population tables
pops <- select(db$POP_EVAL, c('EVALID', 'ESTN_METHOD', 'CN', 'END_INVYR', 'EVAL_TYP')) %>%
rename(EVAL_CN = CN) %>%
left_join(select(db$POP_ESTN_UNIT, c('CN', 'EVAL_CN', 'AREA_USED', 'P1PNTCNT_EU')), by = c('EVAL_CN')) %>%
rename(ESTN_UNIT_CN = CN) %>%
left_join(select(db$POP_STRATUM, c('ESTN_UNIT_CN', 'EXPNS', 'P2POINTCNT', 'CN', 'P1POINTCNT', 'ADJ_FACTOR_SUBP', 'ADJ_FACTOR_MICR', "ADJ_FACTOR_MACR")), by = c('ESTN_UNIT_CN')) %>%
rename(STRATUM_CN = CN) %>%
left_join(select(db$POP_PLOT_STRATUM_ASSGN, c('STRATUM_CN', 'PLT_CN', 'INVYR', 'STATECD')), by = 'STRATUM_CN') %>%
mutate_if(is.factor,
as.character)
### Which estimator to use?
if (str_to_upper(method) %in% c('ANNUAL')){
## Want to use the year where plots are measured, no repeats
## Breaking this up into pre and post reporting becuase
## Estimation units get weird on us otherwise
popOrig <- pops
pops <- pops %>%
group_by(STATECD) %>%
filter(END_INVYR == INVYR) %>%
ungroup()
prePops <- popOrig %>%
group_by(STATECD) %>%
filter(INVYR < min(END_INVYR, na.rm = TRUE)) %>%
distinct(PLT_CN, .keep_all = TRUE) %>%
ungroup()
pops <- bind_rows(pops, prePops) %>%
mutate(YEAR = INVYR)
} else { # Otherwise temporally indifferent
pops <- rename(pops, YEAR = END_INVYR)
}
## P2POINTCNT column is NOT consistent for annnual estimates, plots
## within individual strata and est units are related to different INVYRs
p2_INVYR <- pops %>%
group_by(ESTN_UNIT_CN, STRATUM_CN, INVYR) %>%
summarize(P2POINTCNT_INVYR = length(unique(PLT_CN)))
## Want a count of p2 points / eu, gets screwed up with grouping below
p2eu_INVYR <- p2_INVYR %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, INVYR, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN, INVYR) %>%
summarize(p2eu_INVYR = sum(P2POINTCNT_INVYR, na.rm = TRUE))
p2eu <- pops %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(p2eu = sum(P2POINTCNT, na.rm = TRUE))
## Rejoin
pops <- pops %>%
left_join(p2_INVYR, by = c('ESTN_UNIT_CN', 'STRATUM_CN', 'INVYR')) %>%
left_join(p2eu_INVYR, by = c('ESTN_UNIT_CN', 'INVYR')) %>%
left_join(p2eu, by = 'ESTN_UNIT_CN')
## Recode a few of the estimation methods to make things easier below
pops$ESTN_METHOD = recode(.x = pops$ESTN_METHOD,
`Post-Stratification` = 'strat',
`Stratified random sampling` = 'strat',
`Double sampling for stratification` = 'double',
`Simple random sampling` = 'simple',
`Subsampling units of unequal size` = 'simple')
## Only the necessary plots for EVAL of interest
db$PLOT <- filter(db$PLOT, PLT_CN %in% pops$PLT_CN)
## Which grpByNames are in which table? Helps us subset below
grpP <- names(db$PLOT)[names(db$PLOT) %in% grpBy]
grpC <- names(db$COND)[names(db$COND) %in% grpBy & names(db$COND) %in% grpP == FALSE]
### 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', 'TPA_UNADJ', 'CCLCD', 'SUBP', 'TREE'))
#filter(DIA >= 5)
## 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 = ssHelper1, plts, db, grpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = ssHelper1, plts, db, grpBy, byPlot, mc.cores = nCores)
}
})
if (byPlot){
## back to dataframes
out <- unlist(out, recursive = FALSE)
tOut <- bind_rows(out[names(out) == 't'])
## Make it spatial
if (returnSpatial){
tOut <- tOut %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
st_as_sf(coords = c('LON', 'LAT'),
crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
}
out <- list(tEst = tOut, grpBy = grpBy, grpByOrig = grpByOrig)
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
## Splitting up by STATECD and groups of 25 ESTN_UNIT_CNs
#estunit <- distinct(pops, ESTN_UNIT_CN) #%>%
#mutate(estID)
#estID <- seq(1, nrow(estunit), 50)
#estunit$estID <- rep_len(estID, length.out = nrow(estunit))
#pops <- pops %>%
# left_join(estunit, by = 'ESTN_UNIT_CN') #%>%
#mutate(estBreaks = )
#popState <- split(pops, as.factor(pops$estID))
popState <- split(pops, as.factor(pops$STATECD))
#
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
## Use the same cluster as above
# cl <- makeCluster(nCores)
# clusterEvalQ(cl, {
# library(dplyr)
# library(stringr)
# library(rFIA)
# })
out <- parLapply(cl, X = names(popState), fun = ssHelper2, popState, t, grpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = ssHelper2, popState, t, grpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
##### ----------------- MOVING AVERAGES
if (str_to_upper(method) %in% c("SMA", 'EMA', 'LMA')){
### ---- SIMPLE MOVING AVERAGE
if (str_to_upper(method) == 'SMA'){
## Assuming a uniform weighting scheme
wgts <- pops %>%
group_by(ESTN_UNIT_CN) %>%
summarize(wgt = 1 / length(unique(INVYR)))
#aEst <- left_join(aEst, wgts, by = 'ESTN_UNIT_CN')
tEst <- left_join(tEst, wgts, by = 'ESTN_UNIT_CN')
#### ----- Linear MOVING AVERAGE
} else if (str_to_upper(method) == 'LMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
## Want a number of INVYRs per EU
neu <- wgts %>%
group_by(ESTN_UNIT_CN) %>%
summarize(n = sum(rank, na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = rank / n) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
#### ----- EXPONENTIAL MOVING AVERAGE
} else if (str_to_upper(method) == 'EMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
if (length(lambda) < 2){
## Want sum of weighitng functions
neu <- wgts %>%
mutate(l = lambda) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
} else {
grpBy <- c('lambda', grpBy)
#aGrpBy <- c('lambda', aGrpBy)
## Duplicate weights for each level of lambda
yrWgts <- list()
for (i in 1:length(unique(lambda))) {
yrWgts[[i]] <- mutate(wgts, lambda = lambda[i])
}
wgts <- bind_rows(yrWgts)
## Want sum of weighitng functions
neu <- wgts %>%
group_by(lambda, ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = c('lambda', 'ESTN_UNIT_CN')) %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(lambda, ESTN_UNIT_CN, INVYR, wgt)
}
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
}
### Applying the weights
# Area
tEst <- tEst %>%
mutate_at(vars(aEst:moEst), ~(.*wgt)) %>%
mutate_at(vars(aVar:cvEst_mo), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
}
out <- list(tEst = tEst, grpBy = grpBy, grpByOrig = grpByOrig)
}
return(out)
}
#' @export
standStruct <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
landType = 'forest',
method = 'TI',
lambda = .5,
areaDomain = NULL,
totals = FALSE,
variance = FALSE,
byPlot = FALSE,
tidy = TRUE,
nCores = 1) {
## don't have to change original code
grpBy_quo <- rlang::enquo(grpBy)
areaDomain <- rlang::enquo(areaDomain)
### 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 = standStructStarter, db,
grpBy_quo = grpBy_quo, polys, returnSpatial,
landType, method,
lambda, areaDomain,
totals, byPlot, tidy, nCores, remote, mr)
## Bring the results back
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
grpBy <- out[names(out) == 'grpBy'][[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))
# Combine estimates
tEst <- tEst %>%
ungroup() %>%
select(-c(YEAR)) %>%
left_join(maxyearsT, by = grpBy[!c(grpBy %in% 'YEAR')])
}
})})
##--------------------- TOTALS and RATIOS
suppressWarnings({
## Bring them together
tOut <- ungroup(tEst) %>%
group_by(.dots = grpBy) %>%
summarize_all(sum,na.rm = TRUE) %>%
# Renaming, computing ratios, and SE
mutate(POLE_AREA = pEst,
MATURE_AREA = maEst,
LATE_AREA = lEst,
MOSAIC_AREA = moEst,
AREA_TOTAL = aEst,
## Ratios
POLE_PERC = POLE_AREA / AREA_TOTAL * 100,
MATURE_PERC = MATURE_AREA / AREA_TOTAL * 100,
LATE_PERC = LATE_AREA / AREA_TOTAL * 100,
MOSAIC_PERC = MOSAIC_AREA / AREA_TOTAL * 100,
## Ratio Var
ppVar = (1/AREA_TOTAL^2) * (pVar + (POLE_PERC^2 * aVar) - 2 * POLE_PERC * cvEst_p),
mapVar = (1/AREA_TOTAL^2) * (maVar + (MATURE_PERC^2 * aVar) - 2 * MATURE_PERC * cvEst_ma),
lpVar = (1/AREA_TOTAL^2) * (lVar + (LATE_PERC^2 * aVar) - 2 * LATE_PERC * cvEst_l),
mopVar = (1/AREA_TOTAL^2) * (moVar + (MOSAIC_PERC^2 * aVar) - 2 * MOSAIC_PERC * cvEst_mo),
## SE RATIO
POLE_PERC_SE = sqrt(ppVar) / POLE_PERC * 100,
MATURE_PERC_SE = sqrt(mapVar) / MATURE_PERC * 100,
LATE_PERC_SE = sqrt(lpVar) / LATE_PERC * 100,
MOSAIC_PERC_SE = sqrt(mopVar) / MOSAIC_PERC * 100,
## VAR RATIO
POLE_PERC_VAR = ppVar,
MATURE_PERC_VAR = mapVar,
LATE_PERC_VAR = lpVar,
MOSAIC_PERC_VAR = mopVar,
## SE TOTAL
POLE_AREA_SE = sqrt(pVar) / POLE_AREA * 100,
MATURE_AREA_SE = sqrt(maVar) / MATURE_AREA * 100,
LATE_AREA_SE = sqrt(lVar) / LATE_AREA * 100,
MOSAIC_AREA_SE = sqrt(moVar) / MOSAIC_AREA * 100,
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL *100,
## VAR TOTAL
POLE_AREA_VAR = pVar,
MATURE_AREA_VAR = maVar,
LATE_AREA_VAR = lVar,
MOSAIC_AREA_VAR = moVar,
AREA_TOTAL_VAR = aVar,
## nPlots
nPlots_AREA = plotIn_AREA)
})
if (totals) {
if (variance) {
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_AREA","MATURE_AREA","LATE_AREA","MOSAIC_AREA","AREA_TOTAL",
"POLE_PERC_VAR","MATURE_PERC_VAR","LATE_PERC_VAR", "MOSAIC_PERC_VAR",
"POLE_AREA_VAR", "MATURE_AREA_VAR", "LATE_AREA_VAR", "MOSAIC_AREA_VAR",
"AREA_TOTAL_VAR","nPlots_AREA", 'N')
} else {
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_AREA","MATURE_AREA","LATE_AREA","MOSAIC_AREA","AREA_TOTAL",
"POLE_PERC_SE","MATURE_PERC_SE","LATE_PERC_SE", "MOSAIC_PERC_SE",
"POLE_AREA_SE", "MATURE_AREA_SE", "LATE_AREA_SE", "MOSAIC_AREA_SE",
"AREA_TOTAL_SE","nPlots_AREA")
}
} else {
if (variance){
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_PERC_VAR","MATURE_PERC_VAR","LATE_PERC_VAR", "MOSAIC_PERC_VAR",
"nPlots_AREA", 'N')
} else {
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_PERC_SE","MATURE_PERC_SE","LATE_PERC_SE", "MOSAIC_PERC_SE",
"nPlots_AREA")
}
}
# Snag the names
tNames <- names(tOut)[names(tOut) %in% grpBy == FALSE]
## Tidy things up if they didn't specify polys, returnSpatial
if (tidy & is.null(polys) & returnSpatial == FALSE){
## Writing the variance options after the below, and don't want to change below
## So instead, do a temporary name swap
if (variance) names(tOut) <- str_replace(names(tOut), '_VAR', '_SE')
## pivot longer
stagePERC <- pivot_longer(select(tOut, grpBy, POLE_PERC:MOSAIC_PERC, nPlots_AREA), names_to = 'STAGE', values_to = 'PERC', cols = POLE_PERC:MOSAIC_PERC) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
stagePERC_SE <- pivot_longer(select(tOut, grpBy, POLE_PERC_SE:MOSAIC_PERC_SE), names_to = 'STAGE', values_to = 'PERC_SE', cols = POLE_PERC_SE:MOSAIC_PERC_SE) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
## rejoin
stage <- left_join(stagePERC, stagePERC_SE, by = c(grpBy, 'STAGE')) %>%
select(grpBy, STAGE, PERC, PERC_SE, nPlots_AREA)
if (totals){
stageAREA <- pivot_longer(select(tOut, grpBy, POLE_AREA:MOSAIC_AREA), names_to = 'STAGE', values_to = 'AREA', cols = POLE_AREA:MOSAIC_AREA) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
stageAREA_SE <- pivot_longer(select(tOut, grpBy, POLE_AREA_SE:MOSAIC_AREA_SE), names_to = 'STAGE', values_to = 'AREA_SE', cols = POLE_AREA_SE:MOSAIC_AREA_SE) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
## Rejoin
stage <- stage %>%
left_join(stageAREA, by = c(grpBy, 'STAGE')) %>%
left_join(stageAREA_SE, by = c(grpBy, 'STAGE'))%>%
select(grpBy, STAGE, PERC, AREA, PERC_SE, AREA_SE, nPlots_AREA)
}
tOut <- stage
## Writing the variance options after the below, and don't want to change below
## So instead, do a temporary name swap
if (variance) names(tOut) <- str_replace(names(tOut), '_SE', '_VAR')
}
} # End byPlot
## 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)
}
standStruct_backup <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
landType = 'forest',
method = 'TI',
lambda = .5,
areaDomain = NULL,
byPlot = FALSE,
totals = FALSE,
tidy = TRUE,
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') %>%
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 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(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', 'all') == FALSE){
stop('landType must be one of: "forest", "timber", or "all".')
}
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)
} else if (tolower(landType) == 'all') {
db$COND$landD <- 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)
### 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 == 'EXPCURR') %>%
filter(!is.na(END_INVYR) & !is.na(EVALID) & END_INVYR >= 2003) %>%
distinct(END_INVYR, EVALID, .keep_all = TRUE)# %>%
#group_by(END_INVYR) %>%
#summarise(id = list(EVALID)
## Make an annual panel ID, associated with an INVYR
### The population tables
pops <- select(db$POP_EVAL, c('EVALID', 'ESTN_METHOD', 'CN', 'END_INVYR')) %>%
rename(EVAL_CN = CN) %>%
left_join(select(db$POP_ESTN_UNIT, c('CN', 'EVAL_CN', 'AREA_USED', 'P1PNTCNT_EU')), by = c('EVAL_CN')) %>%
rename(ESTN_UNIT_CN = CN) %>%
left_join(select(db$POP_STRATUM, c('ESTN_UNIT_CN', 'EXPNS', 'P2POINTCNT', 'CN', 'P1POINTCNT', 'ADJ_FACTOR_SUBP', 'ADJ_FACTOR_MICR', "ADJ_FACTOR_MACR")), by = c('ESTN_UNIT_CN')) %>%
rename(STRATUM_CN = CN) %>%
left_join(select(db$POP_PLOT_STRATUM_ASSGN, c('STRATUM_CN', 'PLT_CN', 'INVYR', 'STATECD')), by = 'STRATUM_CN') %>%
ungroup() %>%
mutate_if(is.factor,
as.character)
### Which estimator to use?
if (str_to_upper(method) %in% c('ANNUAL')){
## Want to use the year where plots are measured, no repeats
## Breaking this up into pre and post reporting becuase
## Estimation units get weird on us otherwise
popOrig <- pops
pops <- pops %>%
group_by(STATECD) %>%
filter(END_INVYR == INVYR) %>%
ungroup()
prePops <- popOrig %>%
group_by(STATECD) %>%
filter(INVYR < min(END_INVYR, na.rm = TRUE)) %>%
distinct(PLT_CN, .keep_all = TRUE) %>%
ungroup()
pops <- bind_rows(pops, prePops) %>%
mutate(YEAR = INVYR)
} else { # Otherwise temporally indifferent
pops <- rename(pops, YEAR = END_INVYR)
}
## P2POINTCNT column is NOT consistent for annnual estimates, plots
## within individual strata and est units are related to different INVYRs
p2_INVYR <- pops %>%
group_by(ESTN_UNIT_CN, STRATUM_CN, INVYR) %>%
summarize(P2POINTCNT_INVYR = length(unique(PLT_CN)))
## Want a count of p2 points / eu, gets screwed up with grouping below
p2eu_INVYR <- p2_INVYR %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, INVYR, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN, INVYR) %>%
summarize(p2eu_INVYR = sum(P2POINTCNT_INVYR, na.rm = TRUE))
p2eu <- pops %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(p2eu = sum(P2POINTCNT, na.rm = TRUE))
## Rejoin
pops <- pops %>%
left_join(p2_INVYR, by = c('ESTN_UNIT_CN', 'STRATUM_CN', 'INVYR')) %>%
left_join(p2eu_INVYR, by = c('ESTN_UNIT_CN', 'INVYR')) %>%
left_join(p2eu, by = 'ESTN_UNIT_CN')
## Recode a few of the estimation methods to make things easier below
pops$ESTN_METHOD = recode(.x = pops$ESTN_METHOD,
`Post-Stratification` = 'strat',
`Stratified random sampling` = 'strat',
`Double sampling for stratification` = 'double',
`Simple random sampling` = 'simple',
`Subsampling units of unequal size` = 'simple')
## 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 = ssHelper1, plts, db, grpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = ssHelper1, plts, db, grpBy, byPlot, mc.cores = nCores)
}
})
if (byPlot){
## back to dataframes
out <- unlist(out, recursive = FALSE)
tOut <- bind_rows(out[names(out) == 't'])
## Make it spatial
if (returnSpatial){
tOut <- tOut %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
st_as_sf(coords = c('LON', 'LAT'),
crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
}
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
## Splitting up by ESTN_UNIT
popState <- split(pops, as.factor(pops$STATECD))
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
## Use the same cluster as above
# cl <- makeCluster(nCores)
# clusterEvalQ(cl, {
# library(dplyr)
# library(stringr)
# library(rFIA)
# })
out <- parLapply(cl, X = names(popState), fun = ssHelper2, popState, t, grpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = ssHelper2, popState, t, grpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
##### ----------------- MOVING AVERAGES
if (str_to_upper(method) %in% c("SMA", 'EMA', 'LMA')){
### ---- SIMPLE MOVING AVERAGE
if (str_to_upper(method) == 'SMA'){
## Assuming a uniform weighting scheme
wgts <- pops %>%
group_by(ESTN_UNIT_CN) %>%
summarize(wgt = 1 / length(unique(INVYR)))
#aEst <- left_join(aEst, wgts, by = 'ESTN_UNIT_CN')
tEst <- left_join(tEst, wgts, by = 'ESTN_UNIT_CN')
#### ----- Linear MOVING AVERAGE
} else if (str_to_upper(method) == 'LMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
## Want a number of INVYRs per EU
neu <- wgts %>%
group_by(ESTN_UNIT_CN) %>%
summarize(n = sum(rank, na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = rank / n) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
#### ----- EXPONENTIAL MOVING AVERAGE
} else if (str_to_upper(method) == 'EMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
if (length(lambda) < 2){
## Want sum of weighitng functions
neu <- wgts %>%
mutate(l = lambda) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
} else {
grpBy <- c('lambda', grpBy)
#aGrpBy <- c('lambda', aGrpBy)
## Duplicate weights for each level of lambda
yrWgts <- list()
for (i in 1:length(unique(lambda))) {
yrWgts[[i]] <- mutate(wgts, lambda = lambda[i])
}
wgts <- bind_rows(yrWgts)
## Want sum of weighitng functions
neu <- wgts %>%
group_by(lambda, ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = c('lambda', 'ESTN_UNIT_CN')) %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(lambda, ESTN_UNIT_CN, INVYR, wgt)
}
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
}
### Applying the weights
# Area
tEst <- tEst %>%
mutate_at(vars(aEst:moEst), ~(.*wgt)) %>%
mutate_at(vars(aVar:cvEst_mo), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
}
##--------------------- TOTALS and RATIOS
suppressWarnings({
## Bring them together
tOut <- ungroup(tEst) %>%
group_by(.dots = grpBy) %>%
summarize_all(sum,na.rm = TRUE) %>%
# Renaming, computing ratios, and SE
mutate(POLE_AREA = pEst,
MATURE_AREA = maEst,
LATE_AREA = lEst,
MOSAIC_AREA = moEst,
AREA_TOTAL = aEst,
## Ratios
POLE_PERC = POLE_AREA / AREA_TOTAL * 100,
MATURE_PERC = MATURE_AREA / AREA_TOTAL * 100,
LATE_PERC = LATE_AREA / AREA_TOTAL * 100,
MOSAIC_PERC = MOSAIC_AREA / AREA_TOTAL * 100,
## Ratio Var
ppVar = (1/AREA_TOTAL^2) * (pVar + (POLE_PERC^2 * aVar) - 2 * POLE_PERC * cvEst_p),
mapVar = (1/AREA_TOTAL^2) * (maVar + (MATURE_PERC^2 * aVar) - 2 * MATURE_PERC * cvEst_ma),
lpVar = (1/AREA_TOTAL^2) * (lVar + (LATE_PERC^2 * aVar) - 2 * LATE_PERC * cvEst_l),
mopVar = (1/AREA_TOTAL^2) * (moVar + (MOSAIC_PERC^2 * aVar) - 2 * MOSAIC_PERC * cvEst_mo),
## SE RATIO
POLE_PERC_SE = sqrt(ppVar) / POLE_PERC * 100,
MATURE_PERC_SE = sqrt(mapVar) / MATURE_PERC * 100,
LATE_PERC_SE = sqrt(lpVar) / LATE_PERC * 100,
MOSAIC_PERC_SE = sqrt(mopVar) / MOSAIC_PERC * 100,
## SE TOTAL
POLE_AREA_SE = sqrt(pVar) / POLE_AREA * 100,
MATURE_AREA_SE = sqrt(maVar) / MATURE_AREA * 100,
LATE_AREA_SE = sqrt(lVar) / LATE_AREA * 100,
MOSAIC_AREA_SE = sqrt(moVar) / MOSAIC_AREA * 100,
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL *100,
## nPlots
nPlots_AREA = plotIn_AREA)
})
if (totals) {
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_AREA","MATURE_AREA","LATE_AREA","MOSAIC_AREA","AREA_TOTAL",
"POLE_PERC_SE","MATURE_PERC_SE","LATE_PERC_SE", "MOSAIC_PERC_SE",
"POLE_AREA_SE", "MATURE_AREA_SE", "LATE_AREA_SE", "MOSAIC_AREA_SE",
"AREA_TOTAL_SE","nPlots_AREA")
} else {
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_PERC_SE","MATURE_PERC_SE","LATE_PERC_SE", "MOSAIC_PERC_SE",
"nPlots_AREA")
}
# Snag the names
tNames <- names(tOut)[names(tOut) %in% grpBy == FALSE]
## Tidy things up if they didn't specify polys, returnSpatial
if (tidy & is.null(polys) & returnSpatial == FALSE){
## pivot longer
stagePERC <- pivot_longer(select(tOut, grpBy, POLE_PERC:MOSAIC_PERC, nPlots_AREA), names_to = 'STAGE', values_to = 'PERC', cols = POLE_PERC:MOSAIC_PERC) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
stagePERC_SE <- pivot_longer(select(tOut, grpBy, POLE_PERC_SE:MOSAIC_PERC_SE), names_to = 'STAGE', values_to = 'PERC_SE', cols = POLE_PERC_SE:MOSAIC_PERC_SE) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
## rejoin
stage <- left_join(stagePERC, stagePERC_SE, by = c(grpBy, 'STAGE')) %>%
select(grpBy, STAGE, PERC, PERC_SE, nPlots_AREA)
if (totals){
stageAREA <- pivot_longer(select(tOut, grpBy, POLE_AREA:MOSAIC_AREA), names_to = 'STAGE', values_to = 'AREA', cols = POLE_AREA:MOSAIC_AREA) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
stageAREA_SE <- pivot_longer(select(tOut, grpBy, POLE_AREA_SE:MOSAIC_AREA_SE), names_to = 'STAGE', values_to = 'AREA_SE', cols = POLE_AREA_SE:MOSAIC_AREA_SE) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
## Rejoin
stage <- stage %>%
left_join(stageAREA, by = c(grpBy, 'STAGE')) %>%
left_join(stageAREA_SE, by = c(grpBy, 'STAGE'))%>%
select(grpBy, STAGE, PERC, AREA, PERC_SE, AREA_SE, nPlots_AREA)
}
tOut <- stage
}
} # End byPlot
## 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 standStruct_backup standStruct standStructStarter
Documented in standStruct
standStructStarter <- function(x,
db,
grpBy_quo = NULL,
polys = NULL,
returnSpatial = FALSE,
landType = 'forest',
method = 'TI',
lambda = .5,
areaDomain = NULL,
totals = FALSE,
byPlot = FALSE,
tidy = TRUE,
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') %>%
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 or COND tables. Did you accidentally quote the variables names? e.g. use grpBy = ECOSUBCD (correct) instead of grpBy = "ECOSUBCD". ', collapse = ', '))
} else {
# Convert to character
grpBy <- names(d_quo)
}
} else {
grpBy <- NULL
}
reqTables <- c('PLOT', 'TREE', 'COND', 'POP_PLOT_STRATUM_ASSGN', 'POP_ESTN_UNIT', 'POP_EVAL',
'POP_STRATUM', 'POP_EVAL_TYP', 'POP_EVAL_GRP')
if (!is.null(polys) & first(class(polys)) %in% c('sf', 'SpatialPolygons', 'SpatialPolygonsDataFrame') == FALSE){
stop('polys must be spatial polygons object of class sp or sf. ')
}
if (landType %in% c('timber', 'forest') == FALSE){
stop('landType must be one of: "forest" or "timber".')
}
if (any(reqTables %in% names(db) == FALSE)){
missT <- reqTables[reqTables %in% names(db) == FALSE]
stop(paste('Tables', paste (as.character(missT), collapse = ', '), 'not found in object db.'))
}
if (str_to_upper(method) %in% c('TI', 'SMA', 'LMA', 'EMA', 'ANNUAL') == FALSE) {
warning(paste('Method', method, 'unknown. Defaulting to Temporally Indifferent (TI).'))
}
# I like a unique ID for a plot through time
if (byPlot) {grpBy <- c('pltID', grpBy)}
# Save original grpBy for pretty return with spatial objects
grpByOrig <- grpBy
## IF the object was clipped
if ('prev' %in% names(db$PLOT)){
## Only want the current plots, no grm
db$PLOT <- filter(db$PLOT, prev == 0)
}
### AREAL SUMMARY PREP
if(!is.null(polys)) {
# # Convert polygons to an sf object
# polys <- polys %>%
# as('sf')%>%
# mutate_if(is.factor,
# as.character)
# ## A unique ID
# polys$polyID <- 1:nrow(polys)
#
# # Add shapefile names to grpBy
grpBy = c(grpBy, 'polyID')
## Make plot data spatial, projected same as polygon layer
pltSF <- select(db$PLOT, c('LON', 'LAT', pltID)) %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
distinct(pltID, .keep_all = TRUE)
coordinates(pltSF) <- ~LON+LAT
proj4string(pltSF) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
pltSF <- as(pltSF, 'sf') %>%
st_transform(crs = st_crs(polys))
## Split up polys
polyList <- split(polys, as.factor(polys$polyID))
suppressWarnings({suppressMessages({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
cl <- makeCluster(nCores)
clusterEvalQ(cl, {
library(dplyr)
library(stringr)
library(rFIA)
})
out <- parLapply(cl, X = names(polyList), fun = areal_par, pltSF, polyList)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(polyList), FUN = areal_par, pltSF, polyList, mc.cores = nCores)
}
})})
pltSF <- bind_rows(out)
# A warning
if (length(unique(pltSF$pltID)) < 1){
stop('No plots in db overlap with polys.')
}
## Add polygon names to PLOT
db$PLOT <- db$PLOT %>%
left_join(select(pltSF, polyID, pltID), by = 'pltID')
# Test if any polygons cross state boundaries w/ different recent inventory years (continued w/in loop)
if ('mostRecent' %in% names(db) & length(unique(db$POP_EVAL$STATECD)) > 1){
mergeYears <- pltSF %>%
right_join(select(db$PLOT, PLT_CN, pltID), by = 'pltID') %>%
inner_join(select(db$POP_PLOT_STRATUM_ASSGN, c('PLT_CN', 'EVALID', 'STATECD')), by = 'PLT_CN') %>%
inner_join(select(db$POP_EVAL, c('EVALID', 'END_INVYR')), by = 'EVALID') %>%
group_by(polyID) %>%
summarize(maxYear = max(END_INVYR, na.rm = TRUE))
}
## TO RETURN SPATIAL PLOTS
}
if (byPlot & returnSpatial){
grpBy <- c(grpBy, 'LON', 'LAT')
} # END AREAL
## Build domain indicator function which is 1 if observation meets criteria, and 0 otherwise
# Land type domain indicator
if (tolower(landType) == 'forest'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1, 1, 0)
} else if (tolower(landType) == 'timber'){
db$COND$landD <- ifelse(db$COND$COND_STATUS_CD == 1 & db$COND$SITECLCD %in% c(1, 2, 3, 4, 5, 6) & db$COND$RESERVCD == 0, 1, 0)
}
# 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)
### 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 == 'EXPCURR') %>%
filter(!is.na(END_INVYR) & !is.na(EVALID) & END_INVYR >= 2003) %>%
distinct(END_INVYR, EVALID, .keep_all = TRUE)# %>%
#group_by(END_INVYR) %>%
#summarise(id = list(EVALID)
## If a most-recent subset, make sure that we don't get two reporting years in
## western states
if (mr) {
db$POP_EVAL <- db$POP_EVAL %>%
group_by(EVAL_TYP, STATECD) %>%
filter(END_INVYR == max(END_INVYR, na.rm = TRUE)) %>%
ungroup()
}
## Cut STATECD
db$POP_EVAL <- select(db$POP_EVAL, -c(STATECD))
### The population tables
pops <- select(db$POP_EVAL, c('EVALID', 'ESTN_METHOD', 'CN', 'END_INVYR', 'EVAL_TYP')) %>%
rename(EVAL_CN = CN) %>%
left_join(select(db$POP_ESTN_UNIT, c('CN', 'EVAL_CN', 'AREA_USED', 'P1PNTCNT_EU')), by = c('EVAL_CN')) %>%
rename(ESTN_UNIT_CN = CN) %>%
left_join(select(db$POP_STRATUM, c('ESTN_UNIT_CN', 'EXPNS', 'P2POINTCNT', 'CN', 'P1POINTCNT', 'ADJ_FACTOR_SUBP', 'ADJ_FACTOR_MICR', "ADJ_FACTOR_MACR")), by = c('ESTN_UNIT_CN')) %>%
rename(STRATUM_CN = CN) %>%
left_join(select(db$POP_PLOT_STRATUM_ASSGN, c('STRATUM_CN', 'PLT_CN', 'INVYR', 'STATECD')), by = 'STRATUM_CN') %>%
mutate_if(is.factor,
as.character)
### Which estimator to use?
if (str_to_upper(method) %in% c('ANNUAL')){
## Want to use the year where plots are measured, no repeats
## Breaking this up into pre and post reporting becuase
## Estimation units get weird on us otherwise
popOrig <- pops
pops <- pops %>%
group_by(STATECD) %>%
filter(END_INVYR == INVYR) %>%
ungroup()
prePops <- popOrig %>%
group_by(STATECD) %>%
filter(INVYR < min(END_INVYR, na.rm = TRUE)) %>%
distinct(PLT_CN, .keep_all = TRUE) %>%
ungroup()
pops <- bind_rows(pops, prePops) %>%
mutate(YEAR = INVYR)
} else { # Otherwise temporally indifferent
pops <- rename(pops, YEAR = END_INVYR)
}
## P2POINTCNT column is NOT consistent for annnual estimates, plots
## within individual strata and est units are related to different INVYRs
p2_INVYR <- pops %>%
group_by(ESTN_UNIT_CN, STRATUM_CN, INVYR) %>%
summarize(P2POINTCNT_INVYR = length(unique(PLT_CN)))
## Want a count of p2 points / eu, gets screwed up with grouping below
p2eu_INVYR <- p2_INVYR %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, INVYR, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN, INVYR) %>%
summarize(p2eu_INVYR = sum(P2POINTCNT_INVYR, na.rm = TRUE))
p2eu <- pops %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(p2eu = sum(P2POINTCNT, na.rm = TRUE))
## Rejoin
pops <- pops %>%
left_join(p2_INVYR, by = c('ESTN_UNIT_CN', 'STRATUM_CN', 'INVYR')) %>%
left_join(p2eu_INVYR, by = c('ESTN_UNIT_CN', 'INVYR')) %>%
left_join(p2eu, by = 'ESTN_UNIT_CN')
## Recode a few of the estimation methods to make things easier below
pops$ESTN_METHOD = recode(.x = pops$ESTN_METHOD,
`Post-Stratification` = 'strat',
`Stratified random sampling` = 'strat',
`Double sampling for stratification` = 'double',
`Simple random sampling` = 'simple',
`Subsampling units of unequal size` = 'simple')
## Only the necessary plots for EVAL of interest
db$PLOT <- filter(db$PLOT, PLT_CN %in% pops$PLT_CN)
## Which grpByNames are in which table? Helps us subset below
grpP <- names(db$PLOT)[names(db$PLOT) %in% grpBy]
grpC <- names(db$COND)[names(db$COND) %in% grpBy & names(db$COND) %in% grpP == FALSE]
### 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', 'TPA_UNADJ', 'CCLCD', 'SUBP', 'TREE'))
#filter(DIA >= 5)
## 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 = ssHelper1, plts, db, grpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = ssHelper1, plts, db, grpBy, byPlot, mc.cores = nCores)
}
})
if (byPlot){
## back to dataframes
out <- unlist(out, recursive = FALSE)
tOut <- bind_rows(out[names(out) == 't'])
## Make it spatial
if (returnSpatial){
tOut <- tOut %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
st_as_sf(coords = c('LON', 'LAT'),
crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
}
out <- list(tEst = tOut, grpBy = grpBy, grpByOrig = grpByOrig)
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
## Splitting up by STATECD and groups of 25 ESTN_UNIT_CNs
#estunit <- distinct(pops, ESTN_UNIT_CN) #%>%
#mutate(estID)
#estID <- seq(1, nrow(estunit), 50)
#estunit$estID <- rep_len(estID, length.out = nrow(estunit))
#pops <- pops %>%
# left_join(estunit, by = 'ESTN_UNIT_CN') #%>%
#mutate(estBreaks = )
#popState <- split(pops, as.factor(pops$estID))
popState <- split(pops, as.factor(pops$STATECD))
#
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
## Use the same cluster as above
# cl <- makeCluster(nCores)
# clusterEvalQ(cl, {
# library(dplyr)
# library(stringr)
# library(rFIA)
# })
out <- parLapply(cl, X = names(popState), fun = ssHelper2, popState, t, grpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = ssHelper2, popState, t, grpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
##### ----------------- MOVING AVERAGES
if (str_to_upper(method) %in% c("SMA", 'EMA', 'LMA')){
### ---- SIMPLE MOVING AVERAGE
if (str_to_upper(method) == 'SMA'){
## Assuming a uniform weighting scheme
wgts <- pops %>%
group_by(ESTN_UNIT_CN) %>%
summarize(wgt = 1 / length(unique(INVYR)))
#aEst <- left_join(aEst, wgts, by = 'ESTN_UNIT_CN')
tEst <- left_join(tEst, wgts, by = 'ESTN_UNIT_CN')
#### ----- Linear MOVING AVERAGE
} else if (str_to_upper(method) == 'LMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
## Want a number of INVYRs per EU
neu <- wgts %>%
group_by(ESTN_UNIT_CN) %>%
summarize(n = sum(rank, na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = rank / n) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
#### ----- EXPONENTIAL MOVING AVERAGE
} else if (str_to_upper(method) == 'EMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
if (length(lambda) < 2){
## Want sum of weighitng functions
neu <- wgts %>%
mutate(l = lambda) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
} else {
grpBy <- c('lambda', grpBy)
#aGrpBy <- c('lambda', aGrpBy)
## Duplicate weights for each level of lambda
yrWgts <- list()
for (i in 1:length(unique(lambda))) {
yrWgts[[i]] <- mutate(wgts, lambda = lambda[i])
}
wgts <- bind_rows(yrWgts)
## Want sum of weighitng functions
neu <- wgts %>%
group_by(lambda, ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = c('lambda', 'ESTN_UNIT_CN')) %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(lambda, ESTN_UNIT_CN, INVYR, wgt)
}
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
}
### Applying the weights
# Area
tEst <- tEst %>%
mutate_at(vars(aEst:moEst), ~(.*wgt)) %>%
mutate_at(vars(aVar:cvEst_mo), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
}
out <- list(tEst = tEst, grpBy = grpBy, grpByOrig = grpByOrig)
}
return(out)
}
#' @export
standStruct <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
landType = 'forest',
method = 'TI',
lambda = .5,
areaDomain = NULL,
totals = FALSE,
variance = FALSE,
byPlot = FALSE,
tidy = TRUE,
nCores = 1) {
## don't have to change original code
grpBy_quo <- rlang::enquo(grpBy)
areaDomain <- rlang::enquo(areaDomain)
### 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 = standStructStarter, db,
grpBy_quo = grpBy_quo, polys, returnSpatial,
landType, method,
lambda, areaDomain,
totals, byPlot, tidy, nCores, remote, mr)
## Bring the results back
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
grpBy <- out[names(out) == 'grpBy'][[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))
# Combine estimates
tEst <- tEst %>%
ungroup() %>%
select(-c(YEAR)) %>%
left_join(maxyearsT, by = grpBy[!c(grpBy %in% 'YEAR')])
}
})})
##--------------------- TOTALS and RATIOS
suppressWarnings({
## Bring them together
tOut <- ungroup(tEst) %>%
group_by(.dots = grpBy) %>%
summarize_all(sum,na.rm = TRUE) %>%
# Renaming, computing ratios, and SE
mutate(POLE_AREA = pEst,
MATURE_AREA = maEst,
LATE_AREA = lEst,
MOSAIC_AREA = moEst,
AREA_TOTAL = aEst,
## Ratios
POLE_PERC = POLE_AREA / AREA_TOTAL * 100,
MATURE_PERC = MATURE_AREA / AREA_TOTAL * 100,
LATE_PERC = LATE_AREA / AREA_TOTAL * 100,
MOSAIC_PERC = MOSAIC_AREA / AREA_TOTAL * 100,
## Ratio Var
ppVar = (1/AREA_TOTAL^2) * (pVar + (POLE_PERC^2 * aVar) - 2 * POLE_PERC * cvEst_p),
mapVar = (1/AREA_TOTAL^2) * (maVar + (MATURE_PERC^2 * aVar) - 2 * MATURE_PERC * cvEst_ma),
lpVar = (1/AREA_TOTAL^2) * (lVar + (LATE_PERC^2 * aVar) - 2 * LATE_PERC * cvEst_l),
mopVar = (1/AREA_TOTAL^2) * (moVar + (MOSAIC_PERC^2 * aVar) - 2 * MOSAIC_PERC * cvEst_mo),
## SE RATIO
POLE_PERC_SE = sqrt(ppVar) / POLE_PERC * 100,
MATURE_PERC_SE = sqrt(mapVar) / MATURE_PERC * 100,
LATE_PERC_SE = sqrt(lpVar) / LATE_PERC * 100,
MOSAIC_PERC_SE = sqrt(mopVar) / MOSAIC_PERC * 100,
## VAR RATIO
POLE_PERC_VAR = ppVar,
MATURE_PERC_VAR = mapVar,
LATE_PERC_VAR = lpVar,
MOSAIC_PERC_VAR = mopVar,
## SE TOTAL
POLE_AREA_SE = sqrt(pVar) / POLE_AREA * 100,
MATURE_AREA_SE = sqrt(maVar) / MATURE_AREA * 100,
LATE_AREA_SE = sqrt(lVar) / LATE_AREA * 100,
MOSAIC_AREA_SE = sqrt(moVar) / MOSAIC_AREA * 100,
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL *100,
## VAR TOTAL
POLE_AREA_VAR = pVar,
MATURE_AREA_VAR = maVar,
LATE_AREA_VAR = lVar,
MOSAIC_AREA_VAR = moVar,
AREA_TOTAL_VAR = aVar,
## nPlots
nPlots_AREA = plotIn_AREA)
})
if (totals) {
if (variance) {
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_AREA","MATURE_AREA","LATE_AREA","MOSAIC_AREA","AREA_TOTAL",
"POLE_PERC_VAR","MATURE_PERC_VAR","LATE_PERC_VAR", "MOSAIC_PERC_VAR",
"POLE_AREA_VAR", "MATURE_AREA_VAR", "LATE_AREA_VAR", "MOSAIC_AREA_VAR",
"AREA_TOTAL_VAR","nPlots_AREA", 'N')
} else {
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_AREA","MATURE_AREA","LATE_AREA","MOSAIC_AREA","AREA_TOTAL",
"POLE_PERC_SE","MATURE_PERC_SE","LATE_PERC_SE", "MOSAIC_PERC_SE",
"POLE_AREA_SE", "MATURE_AREA_SE", "LATE_AREA_SE", "MOSAIC_AREA_SE",
"AREA_TOTAL_SE","nPlots_AREA")
}
} else {
if (variance){
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_PERC_VAR","MATURE_PERC_VAR","LATE_PERC_VAR", "MOSAIC_PERC_VAR",
"nPlots_AREA", 'N')
} else {
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_PERC_SE","MATURE_PERC_SE","LATE_PERC_SE", "MOSAIC_PERC_SE",
"nPlots_AREA")
}
}
# Snag the names
tNames <- names(tOut)[names(tOut) %in% grpBy == FALSE]
## Tidy things up if they didn't specify polys, returnSpatial
if (tidy & is.null(polys) & returnSpatial == FALSE){
## Writing the variance options after the below, and don't want to change below
## So instead, do a temporary name swap
if (variance) names(tOut) <- str_replace(names(tOut), '_VAR', '_SE')
## pivot longer
stagePERC <- pivot_longer(select(tOut, grpBy, POLE_PERC:MOSAIC_PERC, nPlots_AREA), names_to = 'STAGE', values_to = 'PERC', cols = POLE_PERC:MOSAIC_PERC) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
stagePERC_SE <- pivot_longer(select(tOut, grpBy, POLE_PERC_SE:MOSAIC_PERC_SE), names_to = 'STAGE', values_to = 'PERC_SE', cols = POLE_PERC_SE:MOSAIC_PERC_SE) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
## rejoin
stage <- left_join(stagePERC, stagePERC_SE, by = c(grpBy, 'STAGE')) %>%
select(grpBy, STAGE, PERC, PERC_SE, nPlots_AREA)
if (totals){
stageAREA <- pivot_longer(select(tOut, grpBy, POLE_AREA:MOSAIC_AREA), names_to = 'STAGE', values_to = 'AREA', cols = POLE_AREA:MOSAIC_AREA) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
stageAREA_SE <- pivot_longer(select(tOut, grpBy, POLE_AREA_SE:MOSAIC_AREA_SE), names_to = 'STAGE', values_to = 'AREA_SE', cols = POLE_AREA_SE:MOSAIC_AREA_SE) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
## Rejoin
stage <- stage %>%
left_join(stageAREA, by = c(grpBy, 'STAGE')) %>%
left_join(stageAREA_SE, by = c(grpBy, 'STAGE'))%>%
select(grpBy, STAGE, PERC, AREA, PERC_SE, AREA_SE, nPlots_AREA)
}
tOut <- stage
## Writing the variance options after the below, and don't want to change below
## So instead, do a temporary name swap
if (variance) names(tOut) <- str_replace(names(tOut), '_SE', '_VAR')
}
} # End byPlot
## 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)
}
standStruct_backup <- function(db,
grpBy = NULL,
polys = NULL,
returnSpatial = FALSE,
landType = 'forest',
method = 'TI',
lambda = .5,
areaDomain = NULL,
byPlot = FALSE,
totals = FALSE,
tidy = TRUE,
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') %>%
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 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(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', 'all') == FALSE){
stop('landType must be one of: "forest", "timber", or "all".')
}
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)
} else if (tolower(landType) == 'all') {
db$COND$landD <- 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)
### 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 == 'EXPCURR') %>%
filter(!is.na(END_INVYR) & !is.na(EVALID) & END_INVYR >= 2003) %>%
distinct(END_INVYR, EVALID, .keep_all = TRUE)# %>%
#group_by(END_INVYR) %>%
#summarise(id = list(EVALID)
## Make an annual panel ID, associated with an INVYR
### The population tables
pops <- select(db$POP_EVAL, c('EVALID', 'ESTN_METHOD', 'CN', 'END_INVYR')) %>%
rename(EVAL_CN = CN) %>%
left_join(select(db$POP_ESTN_UNIT, c('CN', 'EVAL_CN', 'AREA_USED', 'P1PNTCNT_EU')), by = c('EVAL_CN')) %>%
rename(ESTN_UNIT_CN = CN) %>%
left_join(select(db$POP_STRATUM, c('ESTN_UNIT_CN', 'EXPNS', 'P2POINTCNT', 'CN', 'P1POINTCNT', 'ADJ_FACTOR_SUBP', 'ADJ_FACTOR_MICR', "ADJ_FACTOR_MACR")), by = c('ESTN_UNIT_CN')) %>%
rename(STRATUM_CN = CN) %>%
left_join(select(db$POP_PLOT_STRATUM_ASSGN, c('STRATUM_CN', 'PLT_CN', 'INVYR', 'STATECD')), by = 'STRATUM_CN') %>%
ungroup() %>%
mutate_if(is.factor,
as.character)
### Which estimator to use?
if (str_to_upper(method) %in% c('ANNUAL')){
## Want to use the year where plots are measured, no repeats
## Breaking this up into pre and post reporting becuase
## Estimation units get weird on us otherwise
popOrig <- pops
pops <- pops %>%
group_by(STATECD) %>%
filter(END_INVYR == INVYR) %>%
ungroup()
prePops <- popOrig %>%
group_by(STATECD) %>%
filter(INVYR < min(END_INVYR, na.rm = TRUE)) %>%
distinct(PLT_CN, .keep_all = TRUE) %>%
ungroup()
pops <- bind_rows(pops, prePops) %>%
mutate(YEAR = INVYR)
} else { # Otherwise temporally indifferent
pops <- rename(pops, YEAR = END_INVYR)
}
## P2POINTCNT column is NOT consistent for annnual estimates, plots
## within individual strata and est units are related to different INVYRs
p2_INVYR <- pops %>%
group_by(ESTN_UNIT_CN, STRATUM_CN, INVYR) %>%
summarize(P2POINTCNT_INVYR = length(unique(PLT_CN)))
## Want a count of p2 points / eu, gets screwed up with grouping below
p2eu_INVYR <- p2_INVYR %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, INVYR, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN, INVYR) %>%
summarize(p2eu_INVYR = sum(P2POINTCNT_INVYR, na.rm = TRUE))
p2eu <- pops %>%
distinct(ESTN_UNIT_CN, STRATUM_CN, .keep_all = TRUE) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(p2eu = sum(P2POINTCNT, na.rm = TRUE))
## Rejoin
pops <- pops %>%
left_join(p2_INVYR, by = c('ESTN_UNIT_CN', 'STRATUM_CN', 'INVYR')) %>%
left_join(p2eu_INVYR, by = c('ESTN_UNIT_CN', 'INVYR')) %>%
left_join(p2eu, by = 'ESTN_UNIT_CN')
## Recode a few of the estimation methods to make things easier below
pops$ESTN_METHOD = recode(.x = pops$ESTN_METHOD,
`Post-Stratification` = 'strat',
`Stratified random sampling` = 'strat',
`Double sampling for stratification` = 'double',
`Simple random sampling` = 'simple',
`Subsampling units of unequal size` = 'simple')
## 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 = ssHelper1, plts, db, grpBy, byPlot)
#stopCluster(cl) # Keep the cluster active for the next run
} else { # Unix systems
out <- mclapply(names(plts), FUN = ssHelper1, plts, db, grpBy, byPlot, mc.cores = nCores)
}
})
if (byPlot){
## back to dataframes
out <- unlist(out, recursive = FALSE)
tOut <- bind_rows(out[names(out) == 't'])
## Make it spatial
if (returnSpatial){
tOut <- tOut %>%
filter(!is.na(LAT) & !is.na(LON)) %>%
st_as_sf(coords = c('LON', 'LAT'),
crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
grpBy <- grpBy[grpBy %in% c('LAT', 'LON') == FALSE]
}
## Population estimation
} else {
## back to dataframes
out <- unlist(out, recursive = FALSE)
t <- bind_rows(out[names(out) == 't'])
## Adding YEAR to groups
grpBy <- c('YEAR', grpBy)
## Splitting up by ESTN_UNIT
popState <- split(pops, as.factor(pops$STATECD))
suppressWarnings({
## Compute estimates in parallel -- Clusters in windows, forking otherwise
if (Sys.info()['sysname'] == 'Windows'){
## Use the same cluster as above
# cl <- makeCluster(nCores)
# clusterEvalQ(cl, {
# library(dplyr)
# library(stringr)
# library(rFIA)
# })
out <- parLapply(cl, X = names(popState), fun = ssHelper2, popState, t, grpBy, method)
stopCluster(cl)
} else { # Unix systems
out <- mclapply(names(popState), FUN = ssHelper2, popState, t, grpBy, method, mc.cores = nCores)
}
})
## back to dataframes
out <- unlist(out, recursive = FALSE)
tEst <- bind_rows(out[names(out) == 'tEst'])
##### ----------------- MOVING AVERAGES
if (str_to_upper(method) %in% c("SMA", 'EMA', 'LMA')){
### ---- SIMPLE MOVING AVERAGE
if (str_to_upper(method) == 'SMA'){
## Assuming a uniform weighting scheme
wgts <- pops %>%
group_by(ESTN_UNIT_CN) %>%
summarize(wgt = 1 / length(unique(INVYR)))
#aEst <- left_join(aEst, wgts, by = 'ESTN_UNIT_CN')
tEst <- left_join(tEst, wgts, by = 'ESTN_UNIT_CN')
#### ----- Linear MOVING AVERAGE
} else if (str_to_upper(method) == 'LMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
## Want a number of INVYRs per EU
neu <- wgts %>%
group_by(ESTN_UNIT_CN) %>%
summarize(n = sum(rank, na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = rank / n) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
#### ----- EXPONENTIAL MOVING AVERAGE
} else if (str_to_upper(method) == 'EMA'){
wgts <- pops %>%
distinct(YEAR, ESTN_UNIT_CN, INVYR, .keep_all = TRUE) %>%
arrange(YEAR, ESTN_UNIT_CN, INVYR) %>%
group_by(as.factor(YEAR), as.factor(ESTN_UNIT_CN)) %>%
mutate(rank = min_rank(INVYR))
if (length(lambda) < 2){
## Want sum of weighitng functions
neu <- wgts %>%
mutate(l = lambda) %>%
group_by(ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = 'ESTN_UNIT_CN') %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(ESTN_UNIT_CN, INVYR, wgt)
} else {
grpBy <- c('lambda', grpBy)
#aGrpBy <- c('lambda', aGrpBy)
## Duplicate weights for each level of lambda
yrWgts <- list()
for (i in 1:length(unique(lambda))) {
yrWgts[[i]] <- mutate(wgts, lambda = lambda[i])
}
wgts <- bind_rows(yrWgts)
## Want sum of weighitng functions
neu <- wgts %>%
group_by(lambda, ESTN_UNIT_CN) %>%
summarize(l = 1-first(lambda),
sumwgt = sum(l*(1-l)^(1-rank), na.rm = TRUE))
## Rejoining and computing wgts
wgts <- wgts %>%
left_join(neu, by = c('lambda', 'ESTN_UNIT_CN')) %>%
mutate(wgt = l*(1-l)^(1-rank) / sumwgt) %>%
ungroup() %>%
select(lambda, ESTN_UNIT_CN, INVYR, wgt)
}
#aEst <- left_join(aEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
tEst <- left_join(tEst, wgts, by = c('ESTN_UNIT_CN', 'INVYR'))
}
### Applying the weights
# Area
tEst <- tEst %>%
mutate_at(vars(aEst:moEst), ~(.*wgt)) %>%
mutate_at(vars(aVar:cvEst_mo), ~(.*(wgt^2))) %>%
group_by(ESTN_UNIT_CN, .dots = grpBy) %>%
summarize_at(vars(aEst:plotIn_AREA), sum, na.rm = TRUE)
}
##--------------------- TOTALS and RATIOS
suppressWarnings({
## Bring them together
tOut <- ungroup(tEst) %>%
group_by(.dots = grpBy) %>%
summarize_all(sum,na.rm = TRUE) %>%
# Renaming, computing ratios, and SE
mutate(POLE_AREA = pEst,
MATURE_AREA = maEst,
LATE_AREA = lEst,
MOSAIC_AREA = moEst,
AREA_TOTAL = aEst,
## Ratios
POLE_PERC = POLE_AREA / AREA_TOTAL * 100,
MATURE_PERC = MATURE_AREA / AREA_TOTAL * 100,
LATE_PERC = LATE_AREA / AREA_TOTAL * 100,
MOSAIC_PERC = MOSAIC_AREA / AREA_TOTAL * 100,
## Ratio Var
ppVar = (1/AREA_TOTAL^2) * (pVar + (POLE_PERC^2 * aVar) - 2 * POLE_PERC * cvEst_p),
mapVar = (1/AREA_TOTAL^2) * (maVar + (MATURE_PERC^2 * aVar) - 2 * MATURE_PERC * cvEst_ma),
lpVar = (1/AREA_TOTAL^2) * (lVar + (LATE_PERC^2 * aVar) - 2 * LATE_PERC * cvEst_l),
mopVar = (1/AREA_TOTAL^2) * (moVar + (MOSAIC_PERC^2 * aVar) - 2 * MOSAIC_PERC * cvEst_mo),
## SE RATIO
POLE_PERC_SE = sqrt(ppVar) / POLE_PERC * 100,
MATURE_PERC_SE = sqrt(mapVar) / MATURE_PERC * 100,
LATE_PERC_SE = sqrt(lpVar) / LATE_PERC * 100,
MOSAIC_PERC_SE = sqrt(mopVar) / MOSAIC_PERC * 100,
## SE TOTAL
POLE_AREA_SE = sqrt(pVar) / POLE_AREA * 100,
MATURE_AREA_SE = sqrt(maVar) / MATURE_AREA * 100,
LATE_AREA_SE = sqrt(lVar) / LATE_AREA * 100,
MOSAIC_AREA_SE = sqrt(moVar) / MOSAIC_AREA * 100,
AREA_TOTAL_SE = sqrt(aVar) / AREA_TOTAL *100,
## nPlots
nPlots_AREA = plotIn_AREA)
})
if (totals) {
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_AREA","MATURE_AREA","LATE_AREA","MOSAIC_AREA","AREA_TOTAL",
"POLE_PERC_SE","MATURE_PERC_SE","LATE_PERC_SE", "MOSAIC_PERC_SE",
"POLE_AREA_SE", "MATURE_AREA_SE", "LATE_AREA_SE", "MOSAIC_AREA_SE",
"AREA_TOTAL_SE","nPlots_AREA")
} else {
tOut <- tOut %>%
select(grpBy,"POLE_PERC","MATURE_PERC", "LATE_PERC", "MOSAIC_PERC",
"POLE_PERC_SE","MATURE_PERC_SE","LATE_PERC_SE", "MOSAIC_PERC_SE",
"nPlots_AREA")
}
# Snag the names
tNames <- names(tOut)[names(tOut) %in% grpBy == FALSE]
## Tidy things up if they didn't specify polys, returnSpatial
if (tidy & is.null(polys) & returnSpatial == FALSE){
## pivot longer
stagePERC <- pivot_longer(select(tOut, grpBy, POLE_PERC:MOSAIC_PERC, nPlots_AREA), names_to = 'STAGE', values_to = 'PERC', cols = POLE_PERC:MOSAIC_PERC) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
stagePERC_SE <- pivot_longer(select(tOut, grpBy, POLE_PERC_SE:MOSAIC_PERC_SE), names_to = 'STAGE', values_to = 'PERC_SE', cols = POLE_PERC_SE:MOSAIC_PERC_SE) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
## rejoin
stage <- left_join(stagePERC, stagePERC_SE, by = c(grpBy, 'STAGE')) %>%
select(grpBy, STAGE, PERC, PERC_SE, nPlots_AREA)
if (totals){
stageAREA <- pivot_longer(select(tOut, grpBy, POLE_AREA:MOSAIC_AREA), names_to = 'STAGE', values_to = 'AREA', cols = POLE_AREA:MOSAIC_AREA) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
stageAREA_SE <- pivot_longer(select(tOut, grpBy, POLE_AREA_SE:MOSAIC_AREA_SE), names_to = 'STAGE', values_to = 'AREA_SE', cols = POLE_AREA_SE:MOSAIC_AREA_SE) %>%
mutate(STAGE = str_split(STAGE,pattern= '_', simplify = TRUE,)[,1])
## Rejoin
stage <- stage %>%
left_join(stageAREA, by = c(grpBy, 'STAGE')) %>%
left_join(stageAREA_SE, by = c(grpBy, 'STAGE'))%>%
select(grpBy, STAGE, PERC, AREA, PERC_SE, AREA_SE, nPlots_AREA)
}
tOut <- stage
}
} # End byPlot
## 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.