How much forest does Oregon have, and how sure are we?

Forest Inventory and Analysis estimates with georest.usfs

The Forest Service’s Forest Inventory and Analysis (FIA) program measures trees on tens of thousands of field plots across the country and turns them into statewide numbers: forest area, timber volume, carbon, growth, mortality. Because those numbers come from a sample, every one of them has a standard error. This notebook shows how to get the numbers and how to read their uncertainty:

  1. find what can be estimated and for which inventory,

  2. get Oregon’s forest area, statewide and by county, with error bars,

  3. see why splitting a sample many ways makes some numbers unusable,

  4. combine two states correctly, and

  5. see how much the definition of “forest” matters.

The data come from FIADB-API, the public service behind the Forest Service’s EVALIDator tool.

pip install georest plotly

import pandas as pd
import plotly.express as px
import plotly.io as pio

from georest.usfs import fia

pd.set_option("display.max_colwidth", None)   # show text in full, never cut with "..."
pio.renderers.default = "notebook_connected"   # interactive here and on georest.io


def label(code_and_name):
    """'`41001 41001 OR Baker' -> 'Baker'; '`0200 Douglas-fir group' -> 'Douglas-fir group'."""
    text = str(code_and_name).lstrip("`")
    parts = text.split(" ")
    if len(parts) > 3 and len(parts[2]) == 2 and parts[2].isupper():   # county labels
        return " ".join(parts[3:])
    return " ".join(parts[1:]) if parts[0].isdigit() else text

1. What can be estimated, and for which inventory?

FIA offers 752 estimate attributes, each with a number (snum). Searching them runs offline against a catalog that ships with georest.

fia.find_attributes("area of forest land", land_basis="Forest land").head(5)
snum description group land_basis units eval_typ basis
0 2 Area of forest land, in acres Area Forest land acres EXPCURR COND
1 127 Area change of forest land, in acres, on remeasured conditions where both measurements are forest land Area change total Forest land acres EXPCHNG PCOND
2 128 Area change of forest land, in acres, on remeasured conditions where either measurement is forest land Area change total Forest land acres EXPCHNG PCOND
3 136 Average annual area change of forest land, in acres, on remeasured conditions where both measurements are forest land Annual area change Forest land acres EXPCHNG PCOND
4 137 Average annual area change of forest land, in acres, on remeasured conditions where either measurement is forest land Annual area change Forest land acres EXPCHNG PCOND

Attribute 2 is plain forest area in acres. Next, the evaluation (wc): a state’s inventory over a span of years. Oregon’s current one pools plots measured from 2015 to 2024.

oregon = fia.find_evaluations("Oregon")
OREGON = int(oregon["wc"].iloc[0])
oregon
wc state statecd years growth_acct most_recent
0 412024 Oregon 41 2015;2016;2017;2018;2019;2020;2021;2022;2023;2024 Y Y

2. Oregon’s forest area, statewide and by county

estimate returns one row per cell with its estimate, its standard error as a percentage (se_pct), and the number of plots behind it. Totals arrive as the row labelled Total.

by_county = fia.estimate(OREGON, 2, rselected="County code and name")
total = by_county[by_county["row"] == "Total"].iloc[0]

half_width = 1.96 * total["estimate"] * total["se_pct"] / 100
print(f"Oregon has {total['estimate'] / 1e6:.2f} million acres of forest land, "
      f"with a standard error of {total['se_pct']:.2f}% from {total['plots']:,} plots.")
print(f"95% confidence interval: {(total['estimate'] - half_width) / 1e6:.2f} to "
      f"{(total['estimate'] + half_width) / 1e6:.2f} million acres.")

Oregon has 29.78 million acres of forest land, with a standard error of 0.43% from 10,592 plots.
95% confidence interval: 29.53 to 30.03 million acres.

Statewide, the estimate is tight. County by county it is not: each county has far fewer plots. The bars below show each county’s estimate with its 95% interval. Hover for the plot count and standard error; a county in orange is flagged as too uncertain to report.

counties = by_county[by_county["row"] != "Total"].copy()
counties["county"] = counties["row"].map(label)
counties["ci95"] = 1.96 * counties["estimate"] * counties["se_pct"] / 100
counties["status"] = counties["unreliable"].map({True: "too uncertain", False: "reportable"})
counties = counties.sort_values("estimate")

fig = px.bar(counties, x="estimate", y="county", orientation="h", error_x="ci95",
             color="status", color_discrete_map={"reportable": "#2f6f4f", "too uncertain": "#d9480f"},
             custom_data=["plots", "se_pct", "unreliable_reason"],
             labels={"estimate": "Forest land (acres)", "county": ""},
             title="Oregon forest land by county, with 95% intervals")
fig.update_traces(hovertemplate=("<b>%{y}</b><br>%{x:,.0f} acres<br>"
                                 "%{customdata[0]} plots, SE %{customdata[1]:.1f}%"
                                 "<br>%{customdata[2]}<extra></extra>"))
fig.update_layout(height=820, legend_title_text="")
fig.show()

3. Split a sample many ways and some numbers stop meaning anything

Grouping by forest type divides the same plots among 16 groups. The common types are still precise; the rare ones rest on a handful of plots. estimate flags any cell with fewer than 30 plots or a standard error above 30%, and says why.

by_type = fia.estimate(OREGON, 2, rselected="Forest type group")
types = by_type[by_type["row"] != "Total"].copy()
types["forest type"] = types["row"].map(label)
types["ci95"] = 1.96 * types["estimate"] * types["se_pct"] / 100
types["status"] = types["unreliable"].map({True: "too uncertain", False: "reportable"})
types = types.sort_values("estimate")

fig = px.bar(types, x="estimate", y="forest type", orientation="h", error_x="ci95",
             color="status", color_discrete_map={"reportable": "#2f6f4f", "too uncertain": "#d9480f"},
             custom_data=["plots", "se_pct", "unreliable_reason"],
             labels={"estimate": "Forest land (acres)", "forest type": ""},
             title="Oregon forest land by forest type group")
fig.update_traces(hovertemplate=("<b>%{y}</b><br>%{x:,.0f} acres<br>"
                                 "%{customdata[0]} plots, SE %{customdata[1]:.1f}%"
                                 "<br>%{customdata[2]}<extra></extra>"))
fig.update_layout(height=560, legend_title_text="")
fig.show()

types.loc[types["unreliable"], ["forest type", "estimate", "se_pct", "plots", "unreliable_reason"]]
forest type estimate se_pct plots unreliable_reason
8 Oak / hickory group 1571.262810 94.134912 1 single plot - no variance estimate possible
2 Western white pine group 9614.857898 42.626107 7 only 7 plots (< 30); normal approximation for a confidence interval does not hold
10 Aspen / birch group 59825.644771 26.629777 24 only 24 plots (< 30); normal approximation for a confidence interval does not hold
9 Elm / ash / cottonwood group 77650.500373 22.950610 26 only 26 plots (< 30); normal approximation for a confidence interval does not hold

Those flagged rows are still returned, not silently dropped: how much of a table is too thin to use is worth knowing. When you want them gone, fia.reliable(df) removes them.

print(f"{len(fia.reliable(types))} of {len(types)} forest type groups are reportable.")

12 of 16 forest type groups are reportable.

4. Two states together

To estimate across states, pass a list of evaluations. The totals add up either way, but the standard error of a combined estimate cannot be worked out by adding the states’ errors. Here is Oregon plus Washington done three ways.

WASHINGTON = int(fia.find_evaluations("Washington")["wc"].iloc[0])

def state_total(df):
    return df[df["row"] == "Total"].iloc[0]

combined = state_total(fia.estimate([OREGON, WASHINGTON], 2))
separate = [state_total(fia.estimate(wc, 2)) for wc in (OREGON, WASHINGTON)]

se_acres = [r["estimate"] * r["se_pct"] / 100 for r in separate]
added = sum(r["estimate"] for r in separate)
pd.DataFrame([
    {"method": "one request, both evaluations (correct)",
     "acres": combined["estimate"], "se_pct": combined["se_pct"]},
    {"method": "add the errors (wrong: overstates)",
     "acres": added, "se_pct": 100 * sum(se_acres) / added},
    {"method": "combine errors in quadrature",
     "acres": added, "se_pct": 100 * sum(s ** 2 for s in se_acres) ** 0.5 / added},
]).style.format({"acres": "{:,.0f}", "se_pct": "{:.3f}%"}).hide(axis="index")
method acres se_pct
one request, both evaluations (correct) 51,804,091 0.331%
add the errors (wrong: overstates) 51,804,091 0.467%
combine errors in quadrature 51,804,091 0.331%

For two separate states, whose samples are independent, combining the errors in quadrature matches FIA’s own answer. That shortcut breaks inside a state, where the rows of an estimate share the same sample design: there, only a single request (with str_filter to choose the area) gives the right error.

5. What counts as forest?

FIA and the Resources Planning Act (RPA) assessment define forest land differently, and in a dry state with a lot of woodland the difference is large. forest_definition picks one; it is always sent explicitly so two people pulling “the same” number get the same answer.

NEVADA = int(fia.find_evaluations("Nevada")["wc"].iloc[0])
rows = []
for definition in ("FIADEF", "RPADEF"):
    t = state_total(fia.estimate(NEVADA, 2, forest_definition=definition))
    rows.append({"definition": definition, "acres": t["estimate"],
                 "ci95": 1.96 * t["estimate"] * t["se_pct"] / 100, "plots": t["plots"]})
defs = pd.DataFrame(rows)

fig = px.bar(defs, x="definition", y="acres", error_y="ci95", color="definition",
             color_discrete_sequence=["#2f6f4f", "#8c6d1f"], custom_data=["plots"],
             labels={"acres": "Forest land (acres)", "definition": ""},
             title="Nevada forest land under the FIA and RPA definitions")
fig.update_traces(hovertemplate="%{x}: %{y:,.0f} acres<br>%{customdata[0]:,} plots<extra></extra>")
fig.update_layout(showlegend=False)
fig.show()

Where to go next

  • Other attributes. fia.find_attributes("carbon"), "net growth", "mortality". Growth, removals and mortality need an evaluation with growth accounting: fia.find_evaluations("Oregon", growth_only=True).

  • Other breakdowns. fia.find_groupings("owner") lists grouping labels; pass one as rselected= and another as cselected= for a cross-table.

  • A sub-area. str_filter narrows the plots, for example str_filter="plot.unitcd in (0, 1, 2)" for western Oregon’s survey units.

  • Warnings. Anything FIADB-API flags about an estimate is in df.attrs["warnings"].

Full reference: georest.usfs.fia.