library(knitr)
knitr::opts_chunk$set(message = F, warning = F)
# 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")
)

FIESTA Overview

The R package, FIESTA (Forest Inventory ESTimation and Analysis) is a research estimation tool for analysts that work with sample-based inventory data like that from the U.S. Department of Agriculture, Forest Service, Forest Inventory and Analysis (FIA) Program to accommodate: unique population boundaries, different evaluation time periods, customized stratification schemes, non-standard variance equations, integration of multi-scale remotely-sensed data and other auxiliary information, and interaction with other modeling and estimation tools from CRAN R's library of packages. FIESTA contains a collection of functions that can access FIA databases, summarize and compile plot and spatial data, and generate estimates with associated sampling errors.

Functions are organized by type or objective and are named with a corresponding prefix:

Core Functions

Estimation Modules

Analysis Tools

Overview of FIESTA Database (DB) tools

FIESTA's DB tools extract data from FIA's online publicly-available, comma-delimited files (*.csv or *.zip). FIA's CSV files are available by state from the FIA DataMart at the following link: https://apps.fs.usda.gov/fia/datamart/datamart.html. Because of FIA's confidentiality agreement to protect the privacy of landowners, as well as protecting the scientific integrity of FIA's sample design, the exact coordinates of the sample plot locations are not included in the public data. If the exact coordinates are necessary for your analysis, contact FIA's Spatial Data Services (https://www.fia.fs.fed.us/tools-data/spatial/index.php).

Objective of tutorialThe objective of this tutorial is to demonstrate the use of FIESTA's DB tools for accessing FIA data. These tools extract data from FIA Datamart using FIA's standard evaluations as well as customized evaluations.

An FIA Evaluation is a group of plots within the FIA database that is used for population estimates. An FIA Evaluation represents different inventory spans of data with different stratification and area adjustments for nonreponse. Each Evaluation is determined by the type of estimation (evalType) including: area and tree estimates, growth and mortality estimates, and area change estimates (evalType). These plots are identified by an evalid, which is a unique identifier in the format of a 2-digit State code, a 2-digit year code, and a 2-digit evaluation type code. For example, EVALID '491601' represents the Utah 2016 evaluation for current area estimates.

FUNCTION | DESCRIPTION -------------- | --------------------------------------------------------------- DBgetCSV() | Downloads comma-delimited file (.csv) or downloads and extracts a compressed csv file (.zip) from FIA's online DataMart. DBqryCSV() | Extracts and queries data from FIA's online DataMart, either CSV or ZIP files. DBgetEvalid() | Gets evalid for identifying an estimation group of plots for state or checks evalid. DBgetXY() | Extracts XY data from FIA database. DBgetPlots() | Extracts inventory plot data from FIA database. DBgetStrata() | Extracts strata information and total acres by estimation unit from FIA database, including plot-level assignment and a data frame with strata weights by estimation unit.

Set up

First, you'll need to load the FIESTA library:

library(FIESTA)

Next, you'll need to set up an "outfolder". This is just a file path to a folder where you'd like FIESTA to send your data output. For this vignette, we have saved our outfolder file path as the outfolder object in a temporary directory. We also set a few default options preferred for this vignette.

outfolder <- tempdir()

DB Examples

The following examples show how to extract data from FIA's publicly-available, online DataMart. Data can be returned as R objects or exported to CSV (.csv) files or a SQLite (.sqlite) database. The zip files are extracted on-the-fly from the online website. Web server connections will affect download speeds. We show examples for the following functions:

The following examples extract data from FIA's online DataMart (https://apps.fs.usda.gov/fia/datamart/datamart.html).

DBgetCSV()

The DBgetCSV function extracts data from FIA's publicly-available, online DataMart CSV/ZIP files. The zip files are extracted on-the-fly from the online website. Web server connections will affect download speeds.

Example 1: Extract PLOT data for Wyoming and Utah

View Example

WYUTplot <- DBgetCSV(DBtable = "PLOT", 
                     states = c("Wyoming", "Utah"))
table(WYUTplot$STATECD)

DBgetCSV()

## Get plot table for Wyoming
WYplots <- DBgetCSV("PLOT", "Wyoming")
dim(WYplots)

## Get plot table for Wyoming and Utah
WYUTplots <- DBgetCSV(DBtable = "PLOT", 
                      states = c("Wyoming", "Utah"))
table(WYUTplots$STATECD)

## Get survey table for Wyoming
WYsurvey <- DBgetCSV("SURVEY", "Wyoming")
WYsurvey

DBqryCSV()

The DBqryCSV function queries table from FIA' online publicly-available DataMart. The tables in the query must be specified in the sqltables parameter.

Example: Multiple Uses

View Example

DBqryCSV()

# Get number of plots by inventory year for the state of Wyoming
sql1 <- "select INVYR, count(*) AS NBRPLOTS 
         from PLOT 
         where statecd=56 group by INVYR"
DBqryCSV(sql = sql1, 
         states = "Wyoming", 
         sqltables = "PLOT")


# Get number of plots by inventory year for Utah and Wyoming
sql3 <- "select STATECD, INVYR, count(*) NBRPLOTS 
         from PLOT 
         where statecd in(50,33) 
         group by STATECD, INVYR"
DBqryCSV(sql = sql3, 
         states = c("Vermont", "New Hampshire"), 
         sqltables = "PLOT")


# Get number of plots by inventory year for Iowa (stcd=19) that have silver maple (SPCD=317)
sql4 <- "select p.STATECD, p.INVYR, count(*) NBRPLOTS 
         from PLOT p 
         join TREE t ON p.CN = t.PLT_CN 
         where p.statecd = 19 and t.SPCD = 317
         group by p.STATECD, p.INVYR"
DBqryCSV(sql = sql4, 
         states = "IOWA", 
         sqltables = c("PLOT", "TREE"))

DBgetEvalid()

The DBgetEvalid function gets information for FIA Evaluations.

Example 1: Get most current evalid and inventory years for Wyoming

View Example

WYeval <- DBgetEvalid(states = "Wyoming",
                      evalCur = TRUE)

names(WYeval)
WYeval$evalidlist
WYeval$invyrs
WYeval$invyrtab
WYeval$invtype

Example 2: Get most current evaluations for New York and Texas for VOL and GRM evalTypes

View Example

NYTXeval <- DBgetEvalid(states = c("New York", "Texas"), 
                        evalType = c("VOL", "GRM"))

names(NYTXeval)
NYTXeval$evalidlist
NYTXeval$evalTypelist

DBgetXY()

The DBgetXY function queries XY data from FIA' online publicly-available DataMart or SQLite database.

Example1: Get xy data for the state of Wyoming

View Example

DBgetXY()

xydat1 <- DBgetXY(states = "Wyoming", 
                  eval = "FIA",
                  eval_opts = eval_options(Cur = TRUE)
                  )
names(xydat1)
head(xydat1$xyCur_PUBLIC)

Example 2: Add variable in plot table, PLOT_STATUS_CD, and output spatial

View Example

xydat2 <- DBgetXY(states = "Wyoming", 
                  eval = "FIA",
                  eval_opts = eval_options(Cur = TRUE),
                  pvars2keep = c("PLOT_STATUS_CD"),
                  issp = TRUE
                  )
spxy2 <- xydat2$spxy

## Display points with by PLOT_STATUS_CD (1-light blue; 2-brown; 3-blue)
spxy2$color <- with(spxy2, 
          ifelse(PLOT_STATUS_CD == 2, "brown", 
                ifelse(PLOT_STATUS_CD == 3, "blue", "light blue")))
plot(sf::st_geometry(spxy2['PLOT_STATUS_CD']), pch = 16, cex = .5,
                     col = spxy2$color)

Example 3: Get XY data for Wyoming, inventory years 2015 to 2019

View Example

xydat3 <- DBgetXY(states = "Vermont", 
                  eval = "custom",
                  eval_opts = eval_options(invyrs = 2017:2019),
                  issp = TRUE
                  )
spxy3 <- xydat3$spxy

## Display points 
plot(sf::st_geometry(spxy3), pch = 16, cex = .5, col="grey")

## Now only include P2 plots only (intensity1 = TRUE)
xydat3b <- DBgetXY(states = "Vermont", 
                  eval = "custom",
                  eval_opts = eval_options(invyrs = 2017:2019),
                  intensity1 = TRUE,
                  issp = TRUE,
                  )
spxy3b <- xydat3b$spxy

## Display points 
plot(sf::st_geometry(spxy3b), pch = 16, cex = .5, add = TRUE)

DBgetPlots()

The DBgetPlots function extracts data from FIA's online DataMart or SQLite database.

Example 1: Get data for Idaho, most current FIA Evaluation, all plots, no trees

View Example

DBgetPlots()

dat1 <- DBgetPlots(states = "Rhode Island", 
                  eval = "FIA", 
                  eval_opts = eval_options(Cur = TRUE, 
                                           Type = "ALL"),
                  issp = TRUE
                  )
names(dat1)
plt1 <- dat1$tabs$plot
spxy1 <- dat1$xy_PUBLIC
table(plt1$INVYR)

# Display spatial output
plot(sf::st_geometry(spxy1), pch = 16, cex = .5)


# Add a filter to include only plots with Northern red oak forest type (FORTYPCD == 505)
# Note: *allFilter* filters for plots and/or conditions for all states specified.
dat1b <- DBgetPlots(states = "Rhode Island", 
                   eval = "FIA",
                   eval_opts = eval_options(Cur = TRUE, 
                                            Type = "ALL"),
                   issp = TRUE, 
                   allFilter = "FORTYPCD == 505")

names(dat1b)
spxy1b <- dat1b$xy_PUBLIC
dim(spxy1b)

# Display spatial output
plot(sf::st_geometry(spxy1b), pch = 16, cex = .5, add = TRUE, col="green")

Example 3: Get data for Delaware, most current FIA Evaluation, include plotgeom data and subplot tables

View Example

dat2 <- DBgetPlots(states = "Delaware", 
                   eval = "FIA",
                   eval_opts = eval_options(Cur = TRUE, 
                                            Type = "ALL"),
                   issubp = TRUE,
                   plotgeom = TRUE)
names(dat2)
tabs2 <- dat2$tabs
plt2 <- tabs2$plot

## subplot and subp_cond tables are added to tabs list
names(tabs2)

## PLOTGEOM data are appended to plt table (e.g., ALP_ADFORCD, FVS_VARIANT)
head(plt2)

Example 3: Get data for Delaware, most current FIA Evaluation, include pop tables

View Example

dat3 <- DBgetPlots(states = "Delaware", 
                   eval = "FIA",
                   eval_opts = eval_options(Cur = TRUE, 
                                            Type = "ALL"),
                   savePOP = TRUE,
                   othertables = c("POP_STRATUM", "POP_ESTN_UNIT"))

## savePOP = TRUE, saves the POP_PLOT_STRATUM_ASSGN table used to select plots 
names(dat3)

## pop_stratum and pop_estn_unit tables are added to tabs list
tabs3 <- dat3$tabs
names(tabs3)

Example 4: Export data to CSV files

View Example

DBgetPlots(states = "Rhode Island", 
           eval = "FIA",
           eval_opts = eval_options(Cur = TRUE, 
                                    Type = "ALL"),
           returndata = FALSE,
           savedata = TRUE,
           savedata_opts = savedata_options(outfolder = outfolder, 
                                            out_fmt = "csv",
                                            overwrite_layer = TRUE)
           )

## Read in data from outfolder
plt <- read.csv(file.path(outfolder, "plot.csv"), stringsAsFactors=FALSE)
head(plt)

Example 5: Export data to a SQLite database

View Example

DBgetPlots(states = "Rhode Island", 
           eval = "FIA",
           eval_opts = eval_options(Cur = TRUE, 
                                    Type = "ALL"),
           returndata = FALSE,
           savedata = TRUE,
           savedata_opts = savedata_options(outfolder = outfolder, 
                                            out_fmt = "sqlite",
                                            out_dsn = "data.db",
                                            overwrite_dsn = TRUE,
                                            overwrite_layer = TRUE)
           )

## Connect to database and list tables
sqlitefn <- file.path(outfolder, "data.db")
conn <- DBI::dbConnect(RSQLite::SQLite(), sqlitefn)
DBI::dbListTables(conn)

## Read in plot table
plt <- DBI::dbReadTable(conn, "plot")
dim(plt)

## List fields in plot
DBI::dbListFields(conn, "plot")

## Query plot data
DBI::dbGetQuery(conn, "select PLOT_STATUS_CD, count(*) from plot group by PLOT_STATUS_CD")

## Disconnect database connection
DBI::dbDisconnect(conn)

Example 6: Export data to a spatiaLite database

View Example

## Note: SpatiaLite is an extension to SQLite, providing vector geodatabase functionality.

DBgetPlots(states = "Delaware", 
           eval = "FIA",
           eval_opts = eval_options(Cur = TRUE, evalType = "ALL"),
           issp = TRUE,
           returndata = FALSE,
           savedata = TRUE,
           savedata_opts = savedata_options(outfolder = outfolder, 
                                            out_fmt = "sqlite",
                                            out_dsn = "datasp.db",
                                            overwrite_dsn = TRUE,
                                            overwrite_layer = TRUE)
           )

## Connect to database and list tables
splitefn <- file.path(outfolder, "datasp.db")
conn <- DBI::dbConnect(RSQLite::SQLite(), splitefn)
DBI::dbListTables(conn)

## Disconnect database connection
DBI::dbDisconnect(conn)

## Now, let's use sf functions to list and read data
sf::st_layers(splitefn)

## Import spatial xy coordinates (using sf package)
spxy6 <- sf::st_read(splitefn, layer="xy_public")
dim(spxy6)

## Display plots with public coordinates
plot(sf::st_geometry(spxy6), pch = 16, cex = .5)

Example 7: Most current evaluation for multiple evalTypes ('ALL', 'VOL', 'GRM')

View Example

dat7 <- DBgetPlots(states = "Rhode Island", 
                   eval = "FIA",
                   eval_opts = eval_options(Cur = TRUE, 
                                            Type = c("VOL", "CHNG", "P2VEG"))
                   )

names(dat7)
tabs7 <- dat7$tabs
names(tabs7)

ppsa7 <- dat7$pop_plot_stratum_assgn
table(ppsa7$EVALID)

Example 8: Get data for multiple evalids

View Example

dat8 <- DBgetPlots(eval = "FIA",
                   eval_opts = eval_options(Cur = TRUE, 
                                            evalid = c(101800, 101801, 101803))
                   )

names(dat8)
tabs8 <- dat8$tabs
names(tabs8)

ppsa8 <- dat8$pop_plot_stratum_assgn
table(ppsa8$EVALID)

Example 9: Get data for multiple states by Endyr

View Example

dat9 <- DBgetPlots(states = c("Connecticut", "South Carolina"), 
                   eval = "FIA",
                   eval_opts = eval_options(Cur = TRUE, 
                                            evalType = "ALL",
                                            Endyr = 2017),
                   )
names(dat9)
tabs9 <- dat9$tabs
names(tabs9)

ppsa9 <- dat9$pop_plot_stratum_assgn
table(ppsa9$EVALID)

Example 10: Get data for multiple eval Endyrs

View Example

dat10 <- DBgetPlots(states = "Vermont", 
                    eval = "FIA",
                    eval_opts = eval_options(Cur = TRUE, 
                                             evalType = "VOL",
                                             Endyr = 2015:2019),
                    )

names(dat10)
tabs10 <- dat10$tabs
names(tabs10)

ppsa10 <- dat10$pop_plot_stratum_assgn
table(ppsa10$EVALID)

Example 11: Get data for multiple inventory years

View Example

dat11 <- DBgetPlots(states = "Vermont", 
                    eval = "custom",
                    eval_opts = eval_options(invyrs = 2012:2014, 
                                             evalType = "ALL")
                    )

names(dat11)
tabs11 <- dat11$tabs
names(tabs11)
plt11 <- tabs11$plot

table(plt11$INVYR)

Example 12: Get data for Periodic inventory

View Example

dat12 <- DBgetPlots(states = "Wyoming", 
                    invtype = "PERIODIC",
                    eval = "FIA",
                    eval_opts = list(Cur = TRUE, 
                                     evalType = "VOL")
                    )

names(dat12)
tabs12 <- dat12$tabs
names(tabs12)
plt12 <- tabs12$plot

table(plt12$STATECD, plt12$INVYR)

Example 13: Intensity

View Example

The objective of this section is to understand the differences when using INTENSITY=1.

## With intensified plots
dat13 <- DBgetPlots(states = "Vermont", 
                    eval = "FIA",
                    eval_opts = list(Cur = TRUE, 
                                     Type = "ALL"),
                    issp = TRUE
                    )
tabs13 <- dat13$tabs
plt13 <- tabs13$plot
spxy13 <- dat13$xy_PUBLIC

## With only P2 plots (intensity1 = TRUE)
dat13b <- DBgetPlots(states = "Vermont", 
                     eval = "FIA",
                     eval_opts = list(Cur = TRUE, 
                                     Type = "ALL"),
                     intensity1 = TRUE,
                     issp = TRUE
                     )

tabs13b <- dat13b$tabs
plt13b <- tabs13b$plt

table(plt13$INVYR)
table(plt13b$INVYR)
spxy13b <- dat13b$xy_PUBLIC


# Display spatial output of public coordinates
plot(sf::st_geometry(spxy13), pch = 16, cex = .5)
plot(sf::st_geometry(spxy13b), pch = 16, cex = .5, add = TRUE, col = "grey")

DBgetStrata()

The DBgetStrata function queries the FIA database for post-stratification information.

Example1: Get strata for the most current evaluation for Wyoming

View Example

DBgetStrata()

DBgetStrata()

strat1 <- DBgetStrata(states = "Wyoming", 
                      eval_opts = eval_options(Cur = TRUE)
                      )
names(strat1)

## Look at plot assign data
pltassgn1 <- strat1$pltassgn
head(pltassgn1)

unique(pltassgn1$EVALID)
strat1$evalid  


## Look at area data for estimation unit
strat1$unitarea
strat1$unitvar
strat1$unitvar2
strat1$areavar

## Look at stratification data for estimation unit
strat1$stratalut
strat1$strvar
strat1$getwtvar

Example 2: Get strata information for a specific evaluation for Wyoming

View Example

strat2 <- DBgetStrata(eval_opts = eval_options(evalid = 561200))

unique(strat2$pltassgn$EVALID)
strat2$evalid  

Example 3: Get strata information for Wyoming, evaluation ending in 2014

View Example

strat3 <- DBgetStrata(states = "Wyoming",
                      eval_opts = eval_options(Endyr = 2014))

unique(strat3$pltassgn$EVALID)
strat3$evalid  

Example 4: Get strata information for a specific set of Wyoming plots

View Example

strat4 <- DBgetStrata(dat = WYplt)

head(strat4$unitarea)
head(WYunitarea)

Example 5: Get strata information for a specific set of Wyoming plots

View Example

strat5 <- DBgetStrata(states = c("Utah", "Wyoming"),
                      eval_opts = eval_options(Cur = TRUE))

table(strat5$pltassgn$EVALID)


tfrescino/FIESTA documentation built on Feb. 7, 2024, 7:09 a.m.