## ----setup, include = F-------------------------------------------------------
library(knitr)
knitr::opts_chunk$set(message = F, warning = F)
## ----include=FALSE------------------------------------------------------------
# Sets up output folding
hooks = knitr::knit_hooks$get()
hook_foldable = function(type) {
force(type)
function(x, options) {
res = hooks[[type]](x, options)
if (isFALSE(options[[paste0("fold.", type)]])) return(res)
paste0(
"<details><summary>", type, "</summary>\n\n",
res,
"\n\n</details>"
)
}
}
knitr::knit_hooks$set(
output = hook_foldable("output"),
plot = hook_foldable("plot")
)
## ----echo=-1------------------------------------------------------------------
data.table::setDTthreads(2)
## ----warning = F, message = F-------------------------------------------------
library(FIESTA)
## -----------------------------------------------------------------------------
outfolder <- tempdir()
## -----------------------------------------------------------------------------
# File names for spatial layers, stored as external data objects in FIESTA.
WYbhfn <- system.file("extdata", "sp_data/WYbighorn_adminbnd.shp", package = "FIESTA")
fornffn <- system.file("extdata", "sp_data/WYbighorn_forest_nonforest_250m.tif", package = "FIESTA")
# Get estimation unit and strata information for Bighorn National Forest.
stratdat.bh <- spGetStrata(xyplt = WYplt,
uniqueid = "CN",
unit_layer = WYbhfn,
strat_layer = fornffn,
spMakeSpatial_opts = list(xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269))
# Get names of output components
names(stratdat.bh)
# Plot assignment of strata and estimation unit (ONEUNIT, STRATUMCD)
head(stratdat.bh$pltassgn)
# Area by estimation unit
stratdat.bh$unitarea
# Pixel counts and strata weights (strwt) by strata and estimation unit
stratdat.bh$stratalut
# Variable names
stratdat.bh$unitvar # Estimation unit variable
stratdat.bh$strvar # Strata variable
stratdat.bh$areavar # Area variable
## -----------------------------------------------------------------------------
# File names for external spatial data
WYbhdistfn <- system.file("extdata", "sp_data/WYbighorn_districtbnd.shp", package = "FIESTA")
fornffn <- system.file("extdata", "sp_data/WYbighorn_forest_nonforest_250m.tif", package = "FIESTA")
# Get estimation unit and strata information for Bighorn National Forest Districts
stratdat.bhdist <- spGetStrata(xyplt = WYplt,
uniqueid = "CN",
unit_layer = WYbhdistfn,
unitvar = "DISTRICTNA",
strat_layer = fornffn,
spMakeSpatial_opts = list(xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269))
# Get names of output list components
names(stratdat.bhdist)
# Plot assignment of strata and estimation unit (DISTRICTNA, STRATUMCD)
head(stratdat.bhdist$pltassgn)
# Area by estimation units (Districts)
stratdat.bhdist$unitarea
# Pixel counts and strata weights (strwt) by strata and estimation unit
stratdat.bhdist$stratalut
# Variable names
stratdat.bhdist$unitvar # Estimation unit variable
stratdat.bhdist$strvar # Strata variable
stratdat.bhdist$areavar # Area variable
## ----results = FALSE, message = FALSE-----------------------------------------
GBpopdat <- modGBpop(popTabs = list(cond = FIESTA::WYcond, # FIA plot/condition data
tree = FIESTA::WYtree, # FIA tree data
seed = FIESTA::WYseed), # FIA seedling data
popTabIDs = list(cond = "PLT_CN"), # unique ID of plot in cond
pltassgn = FIESTA::WYpltassgn, # plot assignments
pltassgnid = "CN", # unique ID of plot in pltassgn
pjoinid = "PLT_CN", # plot id to join to pltassgn
unitarea = FIESTA::WYunitarea, # area by estimation units
unitvar = "ESTN_UNIT", # name of estimation unit variable
strata = TRUE, # if using post-stratification
stratalut = FIESTA::WYstratalut, # strata classes and pixels counts
strata_opts = strata_options(getwt = TRUE)) # strata options
## ----results = T--------------------------------------------------------------
names(GBpopdat)
## -----------------------------------------------------------------------------
# Number of plots by plot status
GBpopdat$plotsampcnt
# Number of conditions by condition status
GBpopdat$condsampcnt
# Number of plots and adjustment factors by strata
GBpopdat$stratalut
# Adjustment factors added to plot-level data
head(GBpopdat$pltidsadj)
# Adjustment factors added to tree data
head(GBpopdat$treex)
# Adjustment factors added to seedling data
head(GBpopdat$seedx)
## ----message = FALSE----------------------------------------------------------
GBpopdat.bh <- modGBpop(popTabs = popTables(plt = WYplt,
cond = WYcond,
tree = WYtree,
seed = WYseed),
stratdat = stratdat.bh)
# Get names of output list components
names(GBpopdat.bh)
## ----message=FALSE------------------------------------------------------------
## Using output as individual parameter inputs
GBpopdat.bh <- modGBpop(popTabs = popTables(plt = WYplt,
cond = WYcond,
tree = WYtree,
seed = WYseed),
popTabIDs = popTableIDs(plt = "CN"),
pltassgn = stratdat.bh$pltassgn,
pltassgnid = "CN",
unitvar = stratdat.bh$unitvar,
unitarea = stratdat.bh$unitarea,
areavar = stratdat.bh$areavar,
strata = TRUE,
stratalut =stratdat.bh$stratalut,
strvar = stratdat.bh$strvar)
## Get names of output list components
names(GBpopdat.bh)
## ----message = FALSE----------------------------------------------------------
# Bighorn National Forest District
# Using output list from spGetStrata()
GBpopdat.bhdist <- modGBpop(popTabs = popTables(plt = WYplt,
cond = WYcond,
tree = WYtree,
seed = WYseed),
stratdat = stratdat.bhdist)
## Get names of output list components
names(GBpopdat.bhdist)
## ----results = FALSE, message = F---------------------------------------------
SQLitefn <- system.file("extdata", "FIA_data/RIdat_eval2019.db", package="FIESTA")
conn <- DBI::dbConnect(RSQLite::SQLite(), SQLitefn)
DBI::dbListTables(conn)
DBI::dbDisconnect(conn)
## ----results = FALSE, message = F---------------------------------------------
GBpopdat.RI <- modGBpop(popTabs = popTables(plt = "plot",
cond = "cond",
tree = "tree",
seed = "seed"),
dsn = SQLitefn,
pltassgn = "pop_plot_stratum_assgn",
stratalut = "pop_stratum",
unitarea = "pop_estn_unit",
unitvar = "ESTN_UNIT",
areavar = "AREA_USED",
strata_opts = strata_options(getwt = TRUE,
getwtvar = "P1POINTCNT"))
names(GBpopdat.RI)
# Strata-level population data, including number of plots and adjustment factors
GBpopdat.RI$stratalut
## -----------------------------------------------------------------------------
area1.1 <- modGBarea(GBpopdat = GBpopdat,
landarea = "FOREST",
sumunits = TRUE)
## ----results = T--------------------------------------------------------------
names(area1.1)
## ----results = T--------------------------------------------------------------
area1.1$est
## ----results = TRUE-----------------------------------------------------------
raw1.1 <- area1.1$raw # extract raw data list object from output
names(raw1.1)
head(raw1.1$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
raw1.1$totest # estimates for population (i.e., WY)
## -----------------------------------------------------------------------------
## Area of forest land by forest type, Wyoming, 2011-2013
area1.2 <- modGBarea(GBpopdat = GBpopdat,
landarea = "FOREST", # est - forest land filter
rowvar = "FORTYPCD", # est - row domain
sumunits = TRUE, # est - sum estimation units to population
returntitle = TRUE) # out - return title information
## -----------------------------------------------------------------------------
names(area1.2)
## -----------------------------------------------------------------------------
## Estimate and percent sampling error of estimate
area1.2$est
## -----------------------------------------------------------------------------
raw1.2 <- area1.2$raw # extract raw data list object from output
names(raw1.2)
head(raw1.2$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
raw1.2$totest # estimates for population (i.e., WY)
head(raw1.2$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.2$rowest) # estimates by row for population (i.e., WY)
## Titles (list object) for estimate
titlelst1.2 <- area1.2$titlelst
names(titlelst1.2)
titlelst1.2
## ----message = FALSE----------------------------------------------------------
## Area of forest land by forest type and stand-size class, Wyoming, 2011-2013
area1.3 <- modGBarea(GBpopdat = GBpopdat, # pop - population calculations for WY, post-stratification
landarea = "FOREST", # est - forest land filter
rowvar = "FORTYPCD", # est - row domain
colvar = "STDSZCD", # est - column domain
sumunits = TRUE, # est - sum estimation units to population
savedata = TRUE, # out - save data to outfolder
returntitle = TRUE, # out - return title information
table_opts = table_options(row.FIAname = TRUE, # table - row domain names
col.FIAname = TRUE, # table - column domain names
allin1 = TRUE), # table - return output with est(pse)
savedata_opts = savedata_options(outfolder = outfolder, # save - outfolder for saving data
outfn.pre = "WY")) # save - prefix for output files
## -----------------------------------------------------------------------------
# Look at output list
names(area1.3)
# Estimate and percent sampling error of estimate
head(area1.3$est)
# Raw data (list object) for estimate
raw1.3 <- area1.3$raw # extract raw data list object from output
names(raw1.3)
head(raw1.3$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$totest) # estimates for population (i.e., WY)
head(raw1.3$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$rowest) # estimates by row for population (i.e., WY)
head(raw1.3$unit_colest) # estimates by column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$colest) # estimates by column for population (i.e., WY)
head(raw1.3$unit_grpest) # estimates by row and column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$grpest) # estimates by row and column for population (i.e., WY)
## Titles (list object) for estimate
titlelst1.3 <- area1.3$titlelst
names(titlelst1.3)
titlelst1.3
# List output files in outfolder
list.files(outfolder, pattern = "WY_area")
list.files(paste0(outfolder, "/rawdata"), pattern = "WY_area")
## ----messages = FALSE---------------------------------------------------------
area2.1 <- modGBarea(GBpopdat = GBpopdat.bh, # pop - population calculations for Bighorn NF, post-stratification
landarea = "FOREST", # est - forest land filter
rowvar = "FORTYPCD", # est - row domain
colvar = "STDSZCD", # est - column domain
returntitle = TRUE, # out - return title information
title_opts = title_options(title.ref = "Bighorn National Forest, 2011-2013"), # title - customize title reference
table_opts = table_options(row.FIAname = TRUE, # table - return FIA row names
col.FIAname = TRUE)) # table - return FIA column names
## ----results = T--------------------------------------------------------------
names(area2.1)
## ----results = T--------------------------------------------------------------
area2.1$est
## ----results = TRUE-----------------------------------------------------------
raw2.1 <- area2.1$raw # extract raw data list object from output
names(raw2.1)
head(raw2.1$unit_grpest) # estimates by row and group domains
# Titles (list object) for estimate
titlelst2.1 <- area2.1$titlelst
names(titlelst2.1)
titlelst2.1
## -----------------------------------------------------------------------------
area2.2 <- modGBarea(GBpopdat = GBpopdat.bh, # pop - population calculations for Bighorn NF, post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
rowvar = "FORTYPGRPCD", # est - row domain
colvar = "DSTRBCD1", # est - column domain
returntitle = TRUE, # out - return title information
title_opts = title_options(
title.ref = "Bighorn National Forest, 2011-2013" # title - customize title reference
),
table_opts = table_options(
row.FIAname = TRUE, # table - return FIA row names
col.FIAname = TRUE, # table - return FIA column names
estnull = 0,
rowlut = c(180, 200, 220, 260, 280, 900, 999),
raw.keep0 = TRUE
))
## ----results = T--------------------------------------------------------------
names(area2.2)
## ----results = T--------------------------------------------------------------
area2.2$est
## ----results = TRUE-----------------------------------------------------------
## Raw data (list object) for estimate
raw2.2 <- area2.2$raw # extract raw data list object from output
names(raw2.2)
head(raw2.2$unit_grpest) # estimates by row and group domains
## Titles (list object) for estimate
titlelst2.2 <- area2.2$titlelst
names(titlelst2.2)
titlelst2.2
## -----------------------------------------------------------------------------
area3.1 <- modGBarea(GBpopdat = GBpopdat.bhdist, # pop - population calculations for Bighorn NF, post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
pcfilter = "DSTRBCD1 > 0", # est - condition filter for table output
rowvar = "FORTYPGRPCD", # est - row domain
colvar = "DSTRBCD1", # est - column domain
returntitle = TRUE, # out - return title information
title_opts = title_options(
title.ref = "Bighorn National Forest, 2011-2013" # title - customize title reference
),
table_opts = table_options(
row.FIAname = TRUE, # table - return FIA row names
col.FIAname = TRUE # table - return FIA column names
))
## ----results = T--------------------------------------------------------------
names(area3.1)
## ----results = T--------------------------------------------------------------
area3.1$est
## ----results = TRUE-----------------------------------------------------------
raw3.1 <- area3.1$raw
names(raw3.1)
head(raw3.1$unit_rowest) # estimates by estimation unit for row domains
raw3.1$rowest # estimates for population for row domains
head(raw3.1$unit_colest) # estimates by estimation unit for column domains
raw3.1$colest # estimates for population for column domains
## Titles (list object) for estimate
titlelst3.1 <- area3.1$titlelst
names(titlelst3.1)
titlelst3.1
## -----------------------------------------------------------------------------
area4.1 <- modGBarea(GBpopdat = GBpopdat.RI, # pop - population calculations for Bighorn NF, post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
rowvar = "FORTYPCD", # est - row domain
colvar = "STDSZCD", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # table - return FIA row names
col.FIAname = TRUE # table - return FIA column names
))
## ----results = T--------------------------------------------------------------
names(area4.1)
## ----results = T--------------------------------------------------------------
area4.1$est
## ----results = TRUE-----------------------------------------------------------
raw4.1 <- area4.1$raw
names(raw4.1)
head(raw4.1$unit_grpest) # estimates by row and group domains
## Titles (list object) for estimate
titlelst4.1 <- area4.1$titlelst
names(titlelst4.1)
titlelst4.1
## -----------------------------------------------------------------------------
FIESTAutils::ref_estvar[, c("ESTTITLE", "ESTVAR", "ESTFILTER", "ESTUNITS")]
## -----------------------------------------------------------------------------
tree1.1 <- modGBtree(GBpopdat = GBpopdat, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "VOLCFNET", # est - net cubic-foot volume
estvar.filter = "STATUSCD == 1", # est - live trees only
returntitle = TRUE) # out - return title information
## -----------------------------------------------------------------------------
# Look at output list
names(tree1.1)
# Estimate and percent sampling error of estimate
tree1.1$est
# Raw data (list object) for estimate
raw1.1 <- tree1.1$raw # extract raw data list object from output
names(raw1.1)
head(raw1.1$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.1$totest) # estimates for population (i.e., WY)
# Titles (list object) for estimate
titlelst1.1 <- tree1.1$titlelst
names(titlelst1.1)
titlelst1.1
## -----------------------------------------------------------------------------
tree1.2 <- modGBtree(GBpopdat = GBpopdat, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "VOLCFNET", # est - net cubic-foot volume
estvar.filter = "STATUSCD == 1", # est - live trees only
rowvar = "FORTYPCD", # est - row domain
returntitle = TRUE) # out - return title information
## -----------------------------------------------------------------------------
# Look at output list
names(tree1.2)
# Estimate and percent sampling error of estimate
tree1.2$est
# Raw data (list object) for estimate
raw1.2 <- tree1.2$raw # extract raw data list object from output
names(raw1.2)
head(raw1.2$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.2$totest) # estimates for population (i.e., WY)
head(raw1.2$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.2$rowest) # estimates by row for population (i.e., WY)
# Titles (list object) for estimate
titlelst1.2 <- tree1.2$titlelst
names(titlelst1.2)
titlelst1.2
## -----------------------------------------------------------------------------
datBarplot(raw1.2$unit_rowest,
xvar = titlelst1.2$title.rowvar,
yvar = "est")
## -----------------------------------------------------------------------------
tree1.3 <- modGBtree(GBpopdat = GBpopdat, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "VOLCFNET", # est - net cubic-foot volume
estvar.filter = "STATUSCD == 1", # est - live trees only
rowvar = "FORTYPCD", # est - row domain
colvar = "STDSZCD", # est - column domain
returntitle = TRUE, # out - return title information
savedata = TRUE, # out - save data to outfolder
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
col.FIAname = TRUE, # est - column domain names
allin1 = TRUE # out - return output with est(pse)
),
savedata_opts = savedata_options(
outfolder = outfolder, # out - outfolder for saving data
outfn.pre = "WY" # out - prefix for output files
))
## -----------------------------------------------------------------------------
# Look at output list from modGBarea()
names(tree1.3)
# Estimate and percent sampling error of estimate
tree1.3$est
# Raw data (list object) for estimate
raw1.3 <- tree1.3$raw # extract raw data list object from output
names(raw1.3)
head(raw1.3$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$totest) # estimates for population (i.e., WY)
head(raw1.3$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$rowest) # estimates by row for population (i.e., WY)
head(raw1.3$unit_colest) # estimates by column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$colest) # estimates by column for population (i.e., WY)
head(raw1.3$unit_grpest) # estimates by row and column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$grpest) # estimates by row and column for population (i.e., WY)
# Titles (list object) for estimate
titlelst1.3 <- tree1.3$titlelst
names(titlelst1.3)
titlelst1.3
# List output files in outfolder
list.files(outfolder, pattern = "WY_tree")
list.files(paste0(outfolder, "/rawdata"), pattern = "WY_tree")
## -----------------------------------------------------------------------------
tree1.4 <- modGBtree(GBpopdat = GBpopdat, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "TPA_UNADJ", # est - number of trees per acre
estvar.filter = "STATUSCD == 1", # est - live trees only
rowvar = "SPCD", # est - row domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE # out - return output with est and pse
))
## -----------------------------------------------------------------------------
# Look at output list
names(tree1.4)
# Estimate and percent sampling error of estimate
tree1.4$est
## -----------------------------------------------------------------------------
tree1.5 <- modGBtree(GBpopdat = GBpopdat, # pop - population calculations
estseed = "add", # est - add seedling data
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "TPA_UNADJ", # est - number of trees per acre
estvar.filter = "STATUSCD == 1", # est - live trees only
rowvar = "SPCD", # est - row domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE, # out - return output with est and pse
))
## -----------------------------------------------------------------------------
# Look at output list
names(tree1.5)
# Estimate and percent sampling error of estimate
tree1.5$est
# Compare estimates with and without seedlings
head(tree1.4$est)
head(tree1.5$est)
## -----------------------------------------------------------------------------
tree1.6 <- modGBtree(GBpopdat = GBpopdat, # pop - population calculations
estseed = "only", # est - add seedling data
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "TPA_UNADJ", # est - number of trees per acre
rowvar = "SPCD", # est - row domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE # out - return output with est and pse
))
## -----------------------------------------------------------------------------
# Look at output list
names(tree1.6)
# Estimate and percent sampling error of estimate
tree1.6$est
# Compare estimates with, without, and only seedlings
head(tree1.4$est)
head(tree1.5$est)
head(tree1.6$est)
## -----------------------------------------------------------------------------
tree2.1 <- modGBtree(GBpopdat = GBpopdat.bh, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "TPA_UNADJ", # est - number of trees per acre
estvar.filter = "STATUSCD == 1", # est - live trees only
rowvar = "FORTYPCD", # est - row domain
colvar = "SPCD", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
col.FIAname = TRUE, # est - column domain names
allin1 = TRUE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Look at output list
names(tree2.1)
# Estimate and percent sampling error of estimate
head(tree2.1$est)
## -----------------------------------------------------------------------------
deadtree.filter <- "STATUSCD == 2 & STANDING_DEAD_CD == 1"
tree2.2 <- modGBtree(GBpopdat = GBpopdat.bh, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "VOLCFNET", # est - number of trees per acre
estvar.filter = deadtree.filter, # est - standing dead trees only
rowvar = "SPCD", # est - row domain
colvar = "AGENTCD", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
col.FIAname = TRUE, # est - column domain names
allin1 = TRUE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Look at output list
names(tree2.2)
# Estimate and percent sampling error of estimate
head(tree2.2$est)
## -----------------------------------------------------------------------------
# Look at tree data in GBpopdat.bh
head(GBpopdat.bh$treex)
# Use reference data frame stored as an R object in FIESTA
head(FIESTAutils::ref_diacl2in)
# Appends a new column to GBpopdat$treex classifying the DIA variable based on MIN and MAX columns in ref_diacl2in
dat <- datLUTclass(x = GBpopdat.bh$treex,
xvar = "DIA",
LUT = FIESTAutils::ref_diacl2in,
LUTclassnm = "DIACL2IN")
GBpopdat.bh$treex <- dat$xLUT
# Look at tree data, with new column (DIACL2IN)
head(GBpopdat.bh$treex)
# Look at table of new diameter classes (DIACL2IN)
table(GBpopdat.bh$treex$DIACL2IN)
# Another way to append diameter classes
# First, create a new variable using cut function to define 4 diameter classes
dat <- datLUTclass(x = GBpopdat.bh$treex,
xvar = "DIA",
cutbreaks = c(0, 5, 10, 20, 100))
GBpopdat.bh$treex <- dat$xLUT
# Look at tree data, with new column (DIACL2IN)
head(GBpopdat.bh$treex)
# Look at table of new diameter classes (DIACL)
table(GBpopdat.bh$treex$DIACL)
## -----------------------------------------------------------------------------
tree2.3 <- modGBtree(GBpopdat = GBpopdat.bh, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "TPA_UNADJ", # est - number of trees per acre
estvar.filter = "STATUSCD == 1", # est - live trees only
rowvar = "SPGRPCD", # est - row domain
colvar = "DIACL2IN", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = TRUE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Look at output list
names(tree2.3)
# Estimate and percent sampling error of estimate
head(tree2.3$est)
## -----------------------------------------------------------------------------
tree2.4 <- modGBtree(GBpopdat = GBpopdat.bh, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "TPA_UNADJ", # est - number of trees per acre
estvar.filter = "STATUSCD == 1", # est - live trees only
rowvar = "SPGRPCD", # est - row domain
colvar = "DIACL", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = TRUE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Look at output list
names(tree2.4)
# Estimate and percent sampling error of estimate
head(tree2.4$est)
## -----------------------------------------------------------------------------
tree2.5 <- modGBtree(GBpopdat = GBpopdat.bh, # pop - population calculations
estseed = "add", # est - add seedling data
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "TPA_UNADJ", # est - number of trees per acre
estvar.filter = "STATUSCD == 1", # est - live trees only
rowvar = "SPGRPCD", # est - row domain
colvar = "DIACL", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = TRUE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Look at output list
names(tree2.5)
# Estimate and percent sampling error of estimate
head(tree2.5$est)
## -----------------------------------------------------------------------------
deadtree.filter <- "STATUSCD == 2 & STANDING_DEAD_CD == 1"
tree3.1 <- modGBtree(GBpopdat = GBpopdat.bhdist, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = FALSE, # est - sum estimation units to population
estvar = "VOLCFNET", # est - number of trees per acre
estvar.filter = deadtree.filter, # est - only dead trees
rowvar = "FORTYPGRPCD", # est - row domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = TRUE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Look at output list
names(tree3.1)
# Estimate and percent sampling error of estimate
tree3.1$est
# Estimate and percent sampling error by district
tree3.1$raw$unit_rowest
## -----------------------------------------------------------------------------
tree4.1 <- modGBtree(GBpopdat = GBpopdat.RI, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "VOLCFNET", # est - net cubic-foot volume estimate
estvar.filter = "STATUSCD == 1", # est - live trees only
rowvar = "FORTYPCD", # est - row domain
colvar = "STDSZCD", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
col.FIAname = TRUE, # est - column domain names
allin1 = TRUE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Look at output list
names(tree4.1)
# Estimate and percent sampling error of estimate
head(tree4.1$est)
## -----------------------------------------------------------------------------
ratio1.1 <- modGBratio(GBpopdat = GBpopdat, # pop - population calculations
landarea = "TIMBERLAND", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "VOLCFNET", # est - net cubic-foot volume, numerator
estvarn.filter = "STATUSCD == 1", # est - live trees only, numerator
returntitle = TRUE) # out - return title information
## -----------------------------------------------------------------------------
# Look at output list
names(ratio1.1)
# Estimate and percent sampling error of estimate
head(ratio1.1$est)
# Raw data (list object) for estimate
raw1.1 <- ratio1.1$raw # extract raw data list object from output
names(raw1.1)
head(raw1.1$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.1$totest) # estimates for population (i.e., WY)
# Titles (list object) for estimate
titlelst <- ratio1.1$titlelst
names(titlelst)
titlelst
## -----------------------------------------------------------------------------
ratio1.2 <- modGBratio(GBpopdat = GBpopdat, # pop - population calculations
landarea = "TIMBERLAND", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "VOLCFNET", # est - net cubic-foot volume
estvarn.filter = "STATUSCD == 1", # est - live trees only
rowvar = "FORTYPCD", # est - row domain
returntitle = TRUE) # out - return title information
## -----------------------------------------------------------------------------
# Look at output list
names(ratio1.2)
# Estimate and percent sampling error of estimate
head(ratio1.2$est)
# Raw data (list object) for estimate
raw1.2 <- ratio1.2$raw # extract raw data list object from output
names(raw1.2)
head(raw1.2$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.2$totest) # estimates for population (i.e., WY)
head(raw1.2$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.2$rowest) # estimates by row for population (i.e., WY)
# Titles (list object) for estimate
titlelst <- ratio1.2$titlelst
names(titlelst)
titlelst
## -----------------------------------------------------------------------------
ratio1.3 <- modGBratio(GBpopdat = GBpopdat, # pop - population calculations
landarea = "TIMBERLAND", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "VOLCFNET", # est - net cubic-foot volume, numerator
estvarn.filter = "STATUSCD == 1", # est - live trees only, numerator
rowvar = "FORTYPCD", # est - row domain
colvar = "STDSZCD", # est - column domain
returntitle = TRUE, # out - return title information
savedata = TRUE, # out - save data to outfolder
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
col.FIAname = TRUE, # est - column domain names
allin1 = TRUE # out - return output with est(pse)
),
savedata_opts = savedata_options(
outfolder = outfolder, # out - outfolder for saving data
outfn.pre = "WY" # out - prefix for output files
))
## -----------------------------------------------------------------------------
# Look at output list from modGBarea()
names(ratio1.3)
# Estimate and percent sampling error of estimate
head(ratio1.3$est)
# Raw data (list object) for estimate
raw1.3 <- ratio1.3$raw # extract raw data list object from output
names(raw1.3)
head(raw1.3$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$totest) # estimates for population (i.e., WY)
head(raw1.3$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$rowest) # estimates by row for population (i.e., WY)
head(raw1.3$unit_colest) # estimates by column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$colest) # estimates by column for population (i.e., WY)
head(raw1.3$unit_grpest) # estimates by row and column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$grpest) # estimates by row and column for population (i.e., WY)
# Titles (list object) for estimate
titlelst <- ratio1.3$titlelst
names(titlelst)
titlelst
## -----------------------------------------------------------------------------
ratio2.1 <- modGBratio(GBpopdat = GBpopdat.bh, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "TPA_UNADJ", # est - number of trees per acre, numerator
estvarn.filter = "STATUSCD == 1", # est - live trees only, numerator
rowvar = "SPCD", # est - row domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE, # out - return output with est and pse
),
title_opts = title_options(
title.ref = "Bighorn National Forest"
))
## -----------------------------------------------------------------------------
# Look at output list
names(ratio2.1)
# Estimate and percent sampling error of estimate
head(ratio2.1$est)
## -----------------------------------------------------------------------------
ratio2.2 <- modGBratio(GBpopdat = GBpopdat.bh, # pop - population calculations
estseed = "add", # est - add seedling data
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "TPA_UNADJ", # est - number of trees per acre, numerator
estvarn.filter = "STATUSCD == 1", # est - live trees only, numerator
rowvar = "SPCD", # est - row domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE # out - return output with est and pse
),
title_opts = title_options(
title.ref = "Bighorn National Forest"
))
## -----------------------------------------------------------------------------
# Look at output list
names(ratio2.2)
# Estimate and percent sampling error of estimate
head(ratio2.2$est)
# Compare estimates with and without seedlings
head(ratio2.1$est)
head(ratio2.2$est)
## -----------------------------------------------------------------------------
ratio2.3 <- modGBratio(GBpopdat = GBpopdat.bh, # pop - population calculations
estseed = "only", # est - add seedling data
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "TPA_UNADJ", # est - number of trees per acre, numerator
rowvar = "SPCD", # est - row domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE # out - return output with est and pse
))
## -----------------------------------------------------------------------------
# Look at output list
names(ratio2.3)
# Estimate and percent sampling error of estimate
head(ratio2.3$est)
# Compare estimates with, without, and only seedlings
head(ratio2.1$est)
head(ratio2.2$est)
head(ratio2.3$est)
## -----------------------------------------------------------------------------
ratio2.4 <- modGBratio(GBpopdat = GBpopdat.bh, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "TPA_UNADJ", # est - number of trees per acre, numerator
estvarn.filter = "STATUSCD == 1", # est - live trees only, numerator
rowvar = "FORTYPCD", # est - row domain
colvar = "SPCD", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
col.FIAname = TRUE, # est - column domain names
allin1 = TRUE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Look at output list
names(ratio2.4)
# Estimate and percent sampling error of estimate
head(ratio2.4$est)
## -----------------------------------------------------------------------------
deadtree.filter <- "STATUSCD == 2 & STANDING_DEAD_CD == 1"
ratio2.5 <- modGBratio(GBpopdat = GBpopdat.bh, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "VOLCFNET", # est - number of trees per acre, numerator
estvarn.filter = deadtree.filter, # est - standing dead trees only, numerator
rowvar = "SPCD", # est - row domain
colvar = "AGENTCD", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
col.FIAname = TRUE, # est - column domain names
allin1 = TRUE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Look at output list
names(ratio2.5)
# Estimate and percent sampling error of estimate
head(ratio2.5$est)
## -----------------------------------------------------------------------------
ratio2.6 <- modGBratio(GBpopdat = GBpopdat.bh, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "TPA_UNADJ", # est - number of trees per acre, numerator
estvarn.filter = "STATUSCD == 1", # est - live trees only, numerator
rowvar = "SPGRPCD", # est - row domain
colvar = "DIACL2IN", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = list(
row.FIAname = TRUE, # est - row domain names
allin1 = TRUE # out - return output with est(pse)
),
title_opts = list(
title.ref = "Bighorn National Forest"
))
## -----------------------------------------------------------------------------
# Look at output list
names(ratio2.6)
# Estimate and percent sampling error of estimate
head(ratio2.6$est)
## -----------------------------------------------------------------------------
ratio3.2 <- modGBratio(GBpopdat = GBpopdat.bhdist, # pop - population calculations
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "TPA_UNADJ", # est - number of trees per acre, numerator
estvarn.filter = "STATUSCD == 1", # est - live trees only, numerator
rowvar = "SPGRPCD", # est - row domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = TRUE, # out - return output with est(pse)
),
title_opts = title_options(
title.ref="Bighorn National Forest Districts"
))
## -----------------------------------------------------------------------------
# Look at output list
names(ratio3.2)
# Estimate and percent sampling error of estimate
head(ratio3.2$est)
## -----------------------------------------------------------------------------
ratio1.4 <- modGBratio(GBpopdat = GBpopdat, # pop - population calculations for WY, post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "VOLCFNET", # est - net cubic-foot volume, numerator
estvarn.filter = "STATUSCD == 1", # est - live trees only
estvard = "VOLCFNET", # est - net cubic-foot volume, numerator
rowvar = "FORTYPCD", # est - row domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = TRUE, # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Look at output list
names(ratio1.4)
# Estimate and percent sampling error of estimate
head(ratio1.4$est)
## -----------------------------------------------------------------------------
## Get population data for Wyoming estimates, with post-stratification (for comparison)
GBpopdat.strat <- modGBpop(popTabs = popTables(
cond = WYcond, # FIA plot/condition data
tree = WYtree, # FIA tree data
seed = WYseed # FIA seedling data
),
pltassgn = WYpltassgn, # plot assignments
pltassgnid = "CN", # uniqueid of plots
unitarea = WYunitarea, # area by estimation units
unitvar = "ESTN_UNIT", # name of estimation unit
strata = TRUE, # if using post-stratification
stratalut = WYstratalut, # strata classes and pixels counts
strata_opts = strata_options(
getwt=TRUE, # calculate strata weights
getwtvar="P1POINTCNT" # use P1POINTCNT in stratalut to calculate weights
))
## Get population data for Wyoming estimates, with no post-stratification
GBpopdat.nostrat <- modGBpop(popTabs = popTables(
cond = WYcond, # FIA plot/condition data
tree = WYtree, # FIA tree data
seed = WYseed # FIA seedling data
),
pltassgn = WYpltassgn, # plot assignments
pltassgnid = "CN", # uniqueid of plots
unitarea = WYunitarea, # area by estimation units
unitvar = "ESTN_UNIT", # name of estimation unit
strata = FALSE) # if using post-stratification
## -----------------------------------------------------------------------------
## Area of forest land by forest type and stand-size class, Wyoming, 2011-2013, with post-stratification
area.strat <- modGBarea(GBpopdat = GBpopdat.strat, # pop - population calculations for WY, post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
rowvar = "FORTYPCD", # est - row domain
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE # out - return output with est(pse)
))
## Area of forest land by forest type and stand-size class, Wyoming, 2011-2013, no post-stratification
area.nostrat <- modGBarea(GBpopdat = GBpopdat.nostrat, # pop - population calculations for WY, no post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
rowvar = "FORTYPCD", # est - row domain
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Compare estimates and percent standard errors with and without post-stratification
head(area.strat$est)
head(area.nostrat$est)
## -----------------------------------------------------------------------------
## Number of live trees by species, Wyoming, 2011-2013, with post-stratification
tree.strat <- modGBtree(GBpopdat = GBpopdat.strat, # pop - population calculations for WY, post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "TPA_UNADJ", # est - number of trees per acre, numerator
estvar.filter = "STATUSCD == 1", # est - live trees only, numerator
rowvar = "FORTYPCD", # est - row domain
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE # out - return output with est(pse)
))
## Number of live trees by species, Wyoming, 2011-2013, no post-stratification
tree.nostrat <- modGBtree(GBpopdat = GBpopdat.nostrat, # pop - population calculations for WY, no post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvar = "TPA_UNADJ", # est - number of trees per acre, numerator
estvar.filter = "STATUSCD == 1", # est - live trees only, numerator
rowvar = "FORTYPCD", # est - row domain
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Compare estimates and percent standard errors with and without post-stratification
head(tree.strat$est)
head(tree.nostrat$est)
## -----------------------------------------------------------------------------
## Number of live trees per acre by species, Wyoming, 2011-2013, with post-stratification
ratio.strat <- modGBratio(GBpopdat = GBpopdat.strat, # pop - population calculations for WY, post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "TPA_UNADJ", # est - number of trees per acre, numerator
estvarn.filter = "STATUSCD == 1", # est - live trees only, numerator
rowvar = "FORTYPCD", # est - row domain
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE # out - return output with est(pse)
))
## Number of live trees per acre by species, Wyoming, 2011-2013, no post-stratification
ratio.nostrat <- modGBratio(GBpopdat = GBpopdat.nostrat, # pop - population calculations for WY, no post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = TRUE, # est - sum estimation units to population
estvarn = "TPA_UNADJ", # est - number of trees per acre, numerator
estvarn.filter = "STATUSCD == 1", # est - live trees only, numerator
rowvar = "FORTYPCD", # est - row domain
table_opts = table_options(
row.FIAname = TRUE, # est - row domain names
allin1 = FALSE # out - return output with est(pse)
))
## -----------------------------------------------------------------------------
# Compare estimates and percent standard errors with and without post-stratification
head(ratio.strat$est)
head(ratio.nostrat$est)
## -----------------------------------------------------------------------------
area.unit <- modGBarea(GBpopdat = GBpopdat, # pop - population calculations for WY, post-stratification
landarea = "FOREST", # est - forest land filter
sumunits = FALSE, # est - sum estimation units to population
rowvar = "FORTYPCD", # est - row domain
colvar = "STDSZCD", # est - column domain
returntitle = TRUE, # out - return title information
table_opts = table_options(
allin1 = TRUE # out - return output with est(pse)
))
# Estimate and percent sampling error of estimate (first 6 rows)
head(area.unit$est)
# Unique estimation units
unique(area.unit$est$unit)