## ----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( "

", type, "

\n\n", res, "\n\n

" ) } } 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)