Using R and Python to get forest resource data through the EVALIDator API

New ways to access FIA estimates through R and Python.
FIA
forest analytics
Author

Matt Russell

Published

October 25, 2023

If you work with forestry data across the United States, you’ve probably had to turn to Forest Inventory and Analysis (FIA) data at one time or another to look at the status and trends in our forest resources. The FIA database has information contained from over 130,000 forest inventory plots across the country.

One of the tools that analysts often use to access FIA data is the the EVALIDator web interface. EVALIDator is a web-based Application Programming Interface (API) that generates population estimates for core forestry metrics. It provides the ability to query FIA information at the state level or a circular retrieval of information around a fixed geographic point.

There was an excellent webinar offered in August hosted by the FIA National User Group that highlighted some new features available through EVALIDator. For a long while, analysts have used it’s web interface to point and click their way to estimating forest resources after selecting the forestry metrics and geographic scale they’re interested in. As a long-time EVALIDator user, this was a great way to access FIA data. But dealing with the output in an HTML format and getting it into different software to analyze was always an additional step.

I was happy to learn in the webinar about new ways to directly access FIA estimates through R and Python using the EVALIDator API. For me, this can save time and effort given most of my analyses, graphing, and reporting are done through R.

All of the documentation for these steps can be found in the FIADB-API & EVALIDator user documentation. I’ll go through a few examples where I access data from the state of Maine.

Using R to access FIA data at the state level

This first example will access the total amount of carbon stored in the live aboveground portions of trees growing on forestland in Maine, measured in metric tonnes. I’ll ask for estimates to be provided by forest type group and ownership.

First, I’ll use R and will load a few packages to access the data:

library(tidyverse)
library(httr)
library(jsonlite)
library(rlist)
library(knitr)

If you look at the bottom of the documentation webpage, you’ll see example use cases for extracting FIA estimates using R and Python. There are two examples for each that use GET and POST scripts. The GET scripts allow you to enter a complete URL if you know the attributes you’re interested in. The POST script, which I copy here with the fiadb_api_POST() function, allows you to directly input the variables you’re interested in in R:

# fiadb_api_POST() will accept a FIADB-API full report URL and return data frames
# See descriptor: https://apps.fs.usda.gov/fiadb-api/

fiadb_api_POST <- function(argList){
  # make request
  resp <- POST(url = "https://apps.fs.usda.gov/fiadb-api/fullreport", 
               body = argList, 
               encode = "form")
  # parse response from JSON to R list
  respObj <- content(resp, "parsed", encoding = "ISO-8859-1")
  # create empty output list
  outputList <- list()
  # add estimates data frame to output list
  outputList[['estimates']] <- as.data.frame(do.call(rbind,respObj$estimates))

  # if estimate includes totals and subtotals, add those data frames to output list
  if ('subtotals' %in% names(respObj)){
    subtotals <- list()
    # one subtotal data frame for each grouping variable
    for (i in names(respObj$subtotals)){
      subtotals[[i]] <- as.data.frame(do.call(rbind,respObj$subtotals[[i]]))
    }
    outputList[['subtotals']] <- subtotals

    # totals data frame
    outputList[['totals']] <- as.data.frame(do.call(rbind,respObj$totals))
  }

  # add estimate metadata
  outputList[['metadata']] <- respObj$metadata

  return(outputList)
}

The first item to know before accessing FIA data is the numeric code corresponding to the variable you’re interested in. In my case, this is snum = 98, corresponding to the variable that represents “Forest carbon pool 1: live aboveground, in metric tonnes, on forest land.” A note of caution: there are thousands of variables to choose from, but I suspect the most popular ones are listed toward the top of the page.

The second item to know is which wc code you’re interested in. I have no idea why it’s abbreviated “wc”, but it connects the state FIPS code with the 4-digit FIA inventory year. For example, I’m interested in wc = 232021 to obtain data from Maine’s (FIPS code 23) most recent inventory, collected in 2021. You could go back to inventories collected decades ago if you were interested in looking at trends, and those codes are described here.

Finally, you can use the rselected and cselected statements to identify the variables you’d like to group by in rows and columns. In our case this is forest type group and ownership group. I’ll obtain the data in an NJSON format, but you can also obtain the data in HTML, XML, or other formats. You can specify these parameters in arg_list and then store the data in a data frame called all_rows. I love how the data are presented in a long and “tidy” format:

arg_list_maine <- list(snum = 98,
               wc = 232021,
               rselected = 'Forest type group', 
               cselected = 'Ownership group - Major',
               outputFormat = 'NJSON')

# submit list to POST request function
post_data_maine <- fiadb_api_POST(arg_list_maine)

# estimate data frame
all_rows_maine <- post_data_maine[['estimates']]

print(all_rows_maine)
    ESTIMATE                                  GRP1          GRP2 PLOT_COUNT
1    4181422   `0100 White / red / jack pine group  `0001 Public         24
2   34501610   `0100 White / red / jack pine group `0002 Private        246
3   14310434              `0120 Spruce / fir group  `0001 Public        110
4  105362576              `0120 Spruce / fir group `0002 Private       1021
5   62478.28 `0160 Loblolly / shortleaf pine group `0002 Private          1
6   200517.1          `0380 Exotic softwoods group  `0001 Public          1
7     531124          `0380 Exotic softwoods group `0002 Private          5
8   823381.8                `0400 Oak / pine group  `0001 Public          5
9   10607119                `0400 Oak / pine group `0002 Private         73
10  992844.6             `0500 Oak / hickory group  `0001 Public          8
11  11742364             `0500 Oak / hickory group `0002 Private         85
12  409060.3    `0700 Elm / ash / cottonwood group  `0001 Public          7
13   6009873    `0700 Elm / ash / cottonwood group `0002 Private         85
14  14766166     `0800 Maple / beech / birch group  `0001 Public        111
15 136494541     `0800 Maple / beech / birch group `0002 Private       1245
16   3373188             `0900 Aspen / birch group  `0001 Public         29
17  28979204             `0900 Aspen / birch group `0002 Private        329
18  5415.465           `0960 Other hardwoods group  `0001 Public          2
19  734850.1           `0960 Other hardwoods group `0002 Private         26
20  2668.444                      `0999 Nonstocked  `0001 Public          1
21  56632.09                      `0999 Nonstocked `0002 Private         10
         SE SE_PERCENT     VARIANCE
1   1013846   24.24643 1.027883e+12
2   2428370   7.038426 5.896982e+12
3   1346913   9.412105 1.814175e+12
4   3093652   2.936196 9.570681e+12
5  69535.69   111.2958   4835212844
6  197339.8   98.41544   3.8943e+10
7  322121.4     60.649 103762218840
8  383501.6    46.5764 147073471553
9   1405489   13.25043 1.975398e+12
10 376671.9   37.93866 141881751049
11  1392492    11.8587 1.939034e+12
12 210167.8    51.3782  44170517675
13 792626.1   13.18873 628256055243
14  1317475   8.922253  1.73574e+12
15  3597214   2.635427 1.293995e+13
16 724233.3   21.47029 524513911067
17  1928787   6.655763  3.72022e+12
18 4321.885   79.80635     18678690
19 191250.2   26.02575  36576657832
20 2380.922    89.2251      5668789
21 31032.09   54.79594    962990504

This makes the output easy to plot immediately. Here’s a graph of the output, where you can see the vast amount of carbon stored on private land in Maine, mostly in the spruce/fir and maple/beech/birch forest type groups:

ggplot(all_rows_maine, aes(x = as.character(GRP2), 
                           y = as.numeric(ESTIMATE))) +
  geom_bar(stat = "identity") +
  facet_wrap(~as.character(GRP1)) +
  labs(x = "Ownership",
       y = "Forest carbon (metric tonnes)")

Note you’ll need to play with the names of variables a bit to tidy them up, e.g., turning “`0001 Public” to simply “Public”. But the R functions allow you to quickly obtain the data of interest. You can also grab the subtotals of the output to sum all values within each forest type or ownership category:

subtotals_maine <- post_data_maine[['subtotals']]

print(subtotals_maine)
$GRP1
    ESTIMATE                                  GRP1 PLOT_COUNT       SE
1   38683032   `0100 White / red / jack pine group        269  2617501
2  119673009              `0120 Spruce / fir group       1126  3345294
3   62478.28 `0160 Loblolly / shortleaf pine group          1 69535.69
4   731641.2          `0380 Exotic softwoods group          6 377174.6
5   11430501                `0400 Oak / pine group         78  1452201
6   12735208             `0500 Oak / hickory group         92  1432097
7    6418933    `0700 Elm / ash / cottonwood group         92   815680
8  151260707     `0800 Maple / beech / birch group       1349  3791376
9   32352393             `0900 Aspen / birch group        357  2057692
10  740265.6           `0960 Other hardwoods group         28 191299.1
11  59300.54                      `0999 Nonstocked         11 31123.29
   SE_PERCENT     VARIANCE
1    6.766535 6.851312e+12
2    2.795362 1.119099e+13
3    111.2958   4835212844
4    51.55185 142260655272
5    12.70461 2.108887e+12
6    11.24518 2.050902e+12
7    12.70741 665333890414
8    2.506517 1.437453e+13
9    6.360247 4.234097e+12
10   25.84195  36595336521
11     52.484    968659293

$GRP2
   ESTIMATE          GRP2 PLOT_COUNT      SE SE_PERCENT     VARIANCE
1  39065098  `0001 Public        277 1671841   4.279628 2.795052e+12
2 335082370 `0002 Private       2888 3374207   1.006978 1.138527e+13

Finally, you can grab the total number that sums all values. In this case, we learn that there’s about 374 million metric tonnes of carbon stored in the live aboveground portions of trees in Maine forests:

totals_maine <- post_data_maine[['totals']]
print(totals_maine)
   ESTIMATE PLOT_COUNT      SE SE_PERCENT    VARIANCE
1 374147468       3143 3430130  0.9167857 1.17658e+13

Using R to access FIA data around a fixed geographic point

A favorite use of EVALIDator by many is the ability to query FIA data around a specific location. For example, a user can generate population estimates for a 50-mile radius around a proposed mill that uses wood.

Here’s an example that queries FIA data in a 25-mile radius surrounding Bangor, Maine. It estimates the total forestland area (snum = 2) by stand age class in 20-year increments and stand size class (large-, medium-, or small-diameter trees). The latitude and longitude coordinates are specified in the function:

arg_list_bangor <- list(lat = 44.809,
                        lon = -68.771, 
                        radius = 25,
                        wc = 232021,
                        snum = 2,
                        rselected = 'Stand-size class', 
                        cselected = 'Stand age 20 yr classes (0 to 100 plus)',
                        outputFormat = 'NJSON')

# submit list to POST request function
post_data_bangor <- fiadb_api_POST(arg_list_bangor)

# estimate data frame
all_rows_bangor <- post_data_bangor[['estimates']]

print(all_rows_bangor)
   ESTIMATE                  GRP1               GRP2 PLOT_COUNT       SE
1  10445.11  `0001 Large diameter   `0001 0-20 years          2 7478.507
2  226522.4  `0001 Large diameter  `0002 21-40 years         50 34381.51
3  762825.5  `0001 Large diameter  `0003 41-60 years        162 62661.09
4   1906638  `0001 Large diameter  `0004 61-80 years        394 93680.07
5   1883745  `0001 Large diameter `0005 81-100 years        379 96030.66
6   1238973  `0001 Large diameter   `0006 100+ years        253 79016.52
7  32307.91 `0002 Medium diameter   `0001 0-20 years          7 13032.76
8   1351126 `0002 Medium diameter  `0002 21-40 years        276 82210.22
9   2250294 `0002 Medium diameter  `0003 41-60 years        456 104589.1
10  1895798 `0002 Medium diameter  `0004 61-80 years        375 97564.57
11 983302.6 `0002 Medium diameter `0005 81-100 years        207 71444.77
12 327157.7 `0002 Medium diameter   `0006 100+ years         65 42400.25
13   508578  `0003 Small diameter   `0001 0-20 years        124 48671.25
14  2459473  `0003 Small diameter  `0002 21-40 years        495 106273.7
15 995886.2  `0003 Small diameter  `0003 41-60 years        197 72323.27
16 350615.3  `0003 Small diameter  `0004 61-80 years         75 43021.92
17   200729  `0003 Small diameter `0005 81-100 years         41 32470.83
18 109848.2  `0003 Small diameter   `0006 100+ years         22 24218.61
19 24582.06      `0004 Nonstocked   `0001 0-20 years         11 8230.246
   SE_PERCENT    VARIANCE
1    71.59817    55928074
2    15.17797  1182088006
3    8.214342  3926412697
4    4.913365  8775955497
5    5.097859  9221888613
6     6.37758  6243610686
7    40.33923   169852837
8    6.084571  6758519770
9    4.647795 10938874239
10   5.146358  9518845131
11   7.265797  5104355424
12   12.96018  1797780913
13   9.570066  2368890284
14   4.320994 11294089519
15   7.262203  5230655949
16   12.27041  1850885767
17   16.17645  1054354747
18   22.04735   586541186
19   33.48071    67736949

You can plot these data directly to visualize the trends within the 25-mile radius:

ggplot(all_rows_bangor, aes(x = as.character(GRP2), 
                            y = as.numeric(ESTIMATE),
                            fill = as.character(GRP1))) +
  geom_bar(stat = "identity") +
  labs(x = "Stand age",
       y = "Forestland area (acres)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Conclusion

The new features available in the EVALIDator API that allow you to access data through R or Python can help to speed up data analysis using FIA data. You can access the entire history of FIA data and analyze forest resources data by state or a circular retrieval from a fixed geographic point. Special thanks to the USDA Forest Service staff that have made this available!

–

By Matt Russell. Email Matt with any questions or comments.