georest.usfs - FIA and LCMS, made usable

Two public USDA Forest Service APIs that answer complementary halves of the same question and look nothing alike:

LCMS

FIA

Nature

Wall-to-wall 30 m map

Probability sample of plots

Span

1985-2025, annual

1984-present, multi-year panels

Answers

What changed, and where

What is there, and how much

Uncertainty

Map accuracy

Design-based standard error

Neither needs authentication. Earth Engine is optional - only the arbitrary-geometry path uses it (pip install "georest[ee]").

Two things to know before trusting a number from here:

  1. Every FIA estimate carries its sampling error, and cells too thin to report are flagged. An FIA estimate without its error is not a fact.

  2. LCMS and FIA are not directly comparable. This package makes the comparison easy and labels it; it will not hand you a blended number.

github github

1. Discovery works offline

FIA’s /fullreport takes an attribute, two groupings, and an evaluation: 752 x 96 x 96 x 1143. Nobody memorizes that, and the docs still say ‘under construction’.

All three catalogs ship in the wheel, so search is instant and needs no network. That is not theoretical - FIADB-API returned HTTP 200 with an HTML error page across every endpoint while this package was being written, and these cells kept working.

from georest import usfs
from georest.usfs import align, fia, lcms, lcms_ee, vocab

# 752 estimate attributes - search them
fia.find_attributes('carbon', limit=8)
snum description group land_basis units eval_typ basis
0 48 Aboveground carbon in live seedlings, shrubs, ... Carbon Forest land short tons EXPCURR COND
1 49 Belowground carbon in live seedlings, shrubs, ... Carbon Forest land short tons EXPCURR COND
2 50 Carbon in stumps, coarse roots, and coarse woo... Carbon Forest land short tons EXPCURR COND
3 51 Carbon in litter, in short tons, on forest land Carbon Forest land short tons EXPCURR COND
4 52 Carbon in organic soil, in short tons, on fore... Carbon Forest land short tons EXPCURR COND
5 62 Aboveground carbon in live seedlings, shrubs, ... Carbon Timberland short tons EXPCURR COND
6 63 Belowground carbon in live seedlings, shrubs, ... Carbon Timberland short tons EXPCURR COND
7 64 Carbon in stumps, coarse roots, and coarse woo... Carbon Timberland short tons EXPCURR COND 
# Narrow by land basis and the evaluation type an attribute requires
fia.find_attributes('net growth volume', land_basis='Timberland', limit=6)
snum description group land_basis units eval_typ basis
0 208 Average annual net growth of merchantable bole... Annual net growth volume Timberland cubic feet EXPGROW PTREE
1 209 Average annual net growth of sawlog wood volum... Annual net growth volume Timberland board feet (International 1/4-inch rule) EXPGROW PTREE
2 210 Average annual net growth of sawlog wood volum... Annual net growth volume Timberland cubic feet EXPGROW PTREE
3 211 Average annual net growth of merchantable bole... Annual net growth volume Timberland cubic feet EXPGROW PTREE
4 212 Average annual net growth of merchantable bole... Annual net growth volume Timberland cubic feet EXPGROW PTREE
5 956 Average annual net growth of sawlog wood volum... Annual net growth volume Timberland board feet (International 1/4-inch rule) EXPGROW PTREE 
# 1,143 evaluations. growth_only=True keeps those supporting growth,
# removals and mortality - a prerequisite the docs do not surface.
fia.find_evaluations('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 
# 96 grouping variables, with the real database column behind each
fia.find_groupings('species', limit=6)
label db_column has_metadata
0 Artificial regen species C.STDORGSP True
1 Species T.SPCD True
2 Species (seedling) S.SEEDLING True
3 Species group T.SPGRPCD True
4 Species group (seedling) S.SPGRPCD True
5 Species group - Major T.SPCD True

describe_grouping returns FIA’s own prose for a variable, including what each code value means - the difference between a column of integers and a column you can interpret.

print(fia.describe_grouping('Aspect')[:600])

Aspect. The direction of slope, (measured to the nearest degree) for most of the condition. North is from 316 to 360 degrees and 0 to 45 degrees East is from 46 to 135 degrees South is from 136 to 225 degrees West is from 226 to 315 degrees

2. FIA estimates, with their error attached

This is what the package is built around. Every cell comes back with se_pct and plots, and is flagged when SE% > 30 or n < 30.

n = 0 and n = 1 get their own reason codes rather than folding into ‘high error’, because a zero-plot cell reports a zero standard error - which reads as maximum precision when it means no information.

# Local validation first: this pairing is impossible and never leaves.
try:
    fia.validate(wc=12025, snum=999999)
except fia.FIAValidationError as e:
    print('caught locally:', e)

caught locally: no attribute with snum=999999. Use find_attributes('<term>') to search the 752 available.

A note on FIADB-API availability

The estimate endpoint is intermittently unavailable, and it fails in an unusual way: HTTP 200 with an HTML error page rather than a 5xx (‘Internal Server Error: list index out of range’). During development the same query succeeded and failed minutes apart.

usfs detects that, retries it like the 5xx it is, and raises UpstreamError with the real message. The next cell probes once so the rest of the notebook can degrade gracefully instead of erroring out - everything above this point is offline and unaffected.

FIA_UP = True
try:
    _probe = fia.estimate(wc=12025, snum=2)
    print('FIADB-API is up -', len(_probe), 'cells from the probe')
except usfs.UpstreamError as e:
    FIA_UP = False
    print('FIADB-API is down right now:')
    print(' ', e)
    print('\nEstimate cells below will skip. Discovery still works.')

FIADB-API is up - 2 cells from the probe

# Forest area by county x forest type group, Alabama
df = None
if FIA_UP:
    df = fia.estimate(wc=12025, snum=2,
                     rselected='County code and name',
                     cselected='Forest type group')
    print(f'{len(df)} cells, units = ' + str(df['units'].iloc[0]))
    cols = ['column', 'estimate', 'se_pct', 'plots', 'unreliable']
    display(df[df['row'] == 'Total'][cols].head(8))
else:
    print('skipped - FIADB-API unavailable')

499 cells, units = acres
column estimate se_pct plots unreliable
485 `0100 White / red / jack pine group 1.574836e+04 54.919458 4 True
486 `0140 Longleaf / slash pine group 1.184334e+06 6.070555 258 False
487 `0160 Loblolly / shortleaf pine group 9.542519e+06 1.837364 2012 False
488 `0170 Other eastern softwoods group 8.855775e+04 23.601668 21 True
489 `0400 Oak / pine group 2.571904e+06 4.135748 610 False
490 `0500 Oak / hickory group 6.688526e+06 2.195050 1556 False
491 `0600 Oak / gum / cypress group 2.089333e+06 4.551432 498 False
492 `0700 Elm / ash / cottonwood group 5.677498e+05 9.034413 151 False

Look at the flagged rows. On a real cross-tabulation like this one, a large fraction of cells are too thin to report - the honest shape of a sample split many ways, and worth seeing rather than smoothing over.

When this was last run against a live API: 380 of 499 cells flagged.

if df is not None:
    flagged = df[df['unreliable']]
    print(f'{len(flagged)} of {len(df)} cells flagged unreliable')
    display(flagged[['column', 'estimate', 'se_pct', 'plots',
                     'unreliable_reason']].head(5))
else:
    print('skipped - no estimate to inspect')

380 of 499 cells flagged unreliable
column estimate se_pct plots unreliable_reason
0 `0140 Longleaf / slash pine group 7286.993798 73.248072 2 only 2 plots (< 30); normal approximation for ...
1 `0160 Loblolly / shortleaf pine group 143287.569561 19.304803 29 only 29 plots (< 30); normal approximation for...
2 `0400 Oak / pine group 37934.790844 34.726526 12 only 12 plots (< 30); normal approximation for...
3 `0500 Oak / hickory group 69119.604724 27.358213 15 only 15 plots (< 30); normal approximation for...
4 `0600 Oak / gum / cypress group 41159.074801 36.131380 9 only 9 plots (< 30); normal approximation for ... 
# reliable() drops them - a caller decision, never done silently
if df is not None:
    clean = fia.reliable(df)
    print(f'{len(clean)} cells survive the default thresholds')
else:
    print('skipped')

119 cells survive the default thresholds

Several states in one estimate

Pass a list of evaluations and FIADB-API estimates across all of them in one request. That is the only way to get the combined standard error right: adding separate state estimates gives the right total, but their errors cannot be combined from the separate results. Delaware alone is 4.51% and Maryland 2.37%; together, 2.15%.

if FIA_UP:
    de = fia.find_evaluations('Delaware')['wc'].iloc[0]
    md = fia.find_evaluations('Maryland')['wc'].iloc[0]
    both = fia.estimate([de, md], snum=2, rselected='State code')
    display(both[['row', 'estimate', 'se_pct', 'plots', 'unreliable']])
    print('API warnings:', both.attrs['warnings'] or 'none')
else:
    print('skipped - FIADB-API unavailable')
row estimate se_pct plots unreliable
0 `0010 10 Delaware 3.440177e+05 4.513635 114 False
1 `0024 24 Maryland 2.390953e+06 2.367253 413 False
2 Total 2.734971e+06 2.145955 527 False 
API warnings: none

3. LCMS - land cover, land use, and change

Small, clean, release-versioned. 3,643 precomputed summary areas: 3,137 counties, 502 ranger districts, and CONUS / All-Lands rollups.

Releases are not interchangeable

Two traps worth seeing before you pin a release:

  1. Products differ. 2025-6 is a tree canopy release carrying only NLCD_Percent_Tree_Canopy_Cover - no Land_Cover at all. So ‘the latest release’ is ambiguous unless you say latest of what.

  2. The API’s summaries cover less than the imagery. The API’s pre-computed 2025-11 summaries list only CONUS and AK, while the 2025-11 Earth Engine asset covers CONUS, AK, HI and PRUSVI through 2025. For Hawaii or Puerto Rico, summarize the asset over a geometry instead of pinning an older release.

for r in lcms.lcms_releases():
    prods = [p['Name'] for p in r['Products']]
    sas = [s['StudyArea'] for s in r['StudyAreas']]
    print(f"{r['VersionNumber']:9s} {r['StartYear']}-{r['EndYear']} "
          f"areas={r['SummaryAreaCount']:5d}")
    print(f"          products={prods}")
    print(f"          study_areas={sas}")

2025-11   1985-2025 areas= 3643
          products=['Change', 'Land_Cover', 'Land_Use']
          study_areas=['AK', 'CONUS']
2025-6    1985-2025 areas= 3795
          products=['NLCD_Percent_Tree_Canopy_Cover']
          study_areas=['CONUS', 'AK']
2024-10   1985-2024 areas= 3730
          products=['Change', 'Land_Cover', 'Land_Use']
          study_areas=['PRUSVI', 'HAWAII', 'AK', 'CONUS']
2023-9    1985-2023 areas= 3716
          products=['Change', 'Land_Cover', 'Land_Use']
          study_areas=['HAWAII', 'PRUSVI', 'CONUS', 'SEAK']
2022-8    1985-2022 areas=    0
          products=['Change', 'Land_Cover', 'Land_Use']
          study_areas=['PRUSVI', 'CONUS', 'SEAK']

# 'latest' resolves per product
print('latest with Land_Cover :', lcms.latest_release('Land_Cover'))
print('latest with TCC        :',
      lcms.latest_release('NLCD_Percent_Tree_Canopy_Cover'))

latest with Land_Cover : 2025-11
latest with TCC        : 2025-6

# Asking a tree-canopy release for land cover is caught locally,
# and the error names a release that would work.
try:
    lcms.lcms_summary('Land_Cover', state='Oregon', county='Crook',
                    year=2025, release='2025-6')
except ValueError as e:
    print(e)

release '2025-6' does not publish 'Land_Cover' - it has ['NLCD_Percent_Tree_Canopy_Cover']. Releases with 'Land_Cover': ['2025-11', '2024-10', '2023-9', '2022-8'] (pass release='2025-11')

# Class areas for one county, one year. 2025 is the last year in the
# current Land_Cover release (2025-11); see latest_release above.
lc = lcms.lcms_summary(state='Oregon', county='Crook', year=2025)
lc[['year', 'class_name', 'acres', 'source']]
year class_name acres source
0 2025 Barren & Grass/Forb/Herb Mix 14.010875 lcms-api
1 2025 Barren or Impervious 25920.341384 lcms-api
2 2025 Barren & Shrubs Mix 1128.876225 lcms-api
3 2025 Barren & Trees Mix 32.914437 lcms-api
4 2025 Grass/Forb/Herb 434991.192239 lcms-api
5 2025 Grass/Forb/Herb & Shrubs Mix 632294.559757 lcms-api
6 2025 Grass/Forb/Herb & Trees Mix 31487.329102 lcms-api
7 2025 Non-Processing Area Mask 0.444790 lcms-api
8 2025 Shrubs 214845.650364 lcms-api
9 2025 Shrubs & Trees Mix 1810.071630 lcms-api
10 2025 Trees 563582.114595 lcms-api
11 2025 Water 5595.454258 lcms-api 
# Omit `year` for the full 1985-2025 series (~64 KB - fetch eagerly)
ts = lcms.lcms_summary('Change', state='Oregon', county='Crook')
print(f'{len(ts)} rows across ' + str(ts['year'].nunique()) + ' years')
ts[ts['class_name'] == 'Wildfire'].nlargest(5, 'acres')[
    ['year', 'class_name', 'acres']]

369 rows across 41 years
year class_name acres
368 2025 Wildfire 22036.215445
149 2001 Wildfire 7242.288073
358 2024 Wildfire 2187.030889
332 2021 Wildfire 1711.773109
302 2018 Wildfire 1371.508999

Palettes come from the API

Each class arrives with its ClassValue and ClassPalette, so a legend cannot drift out of sync with a release. Any hardcoded LCMS palette is a maintenance liability this removes.

lcms.lcms_classes('Change').head(8)
class_name class_value palette
0 Wind 1 ff09f3
1 Hurricane 2 541aff
2 Snow or Ice Transition 3 e4f5fd
3 Desiccation 4 cc982e
4 Inundation 5 0adaff
5 Prescribed Fire 6 a10018
6 Wildfire 7 d54309
7 Mechanical Land Transformation 8 fafa4b 
vp = lcms.lcms_vis_params('Change')
print('min=', vp['min'], 'max=', vp['max'],
      'palette entries=', len(vp['palette']))
list(vp['classLegendDict'].items())[:5]

min= 1 max= 16 palette entries= 16

[('Wind', 'ff09f3'),
 ('Hurricane', '541aff'),
 ('Snow or Ice Transition', 'e4f5fd'),
 ('Desiccation', 'cc982e'),
 ('Inundation', '0adaff')]

4. Arbitrary geometry (needs Earth Engine)

Install with pip install "georest[ee]", then authenticate and initialize Earth Engine as usual.

The API serves only its 3,643 precomputed areas - bbox, geojson and geometry are all rejected and POST returns 403. So lcms_summary dispatches a geometry to Earth Engine instead, computing the same class areas over any polygon.

Same call, same output shape, and the frame records which backend produced it. Skip this section if Earth Engine is not initialized - everything above works without it.

RUN_EE = False  # flip to True once ee is initialized

if RUN_EE:
    import ee
    aoi = ee.Geometry.Rectangle([-120.9, 44.1, -120.5, 44.4])
    ee_lc = lcms.lcms_summary('Land_Cover', geometry=aoi,
                            year=2025, scale=90)
    display(ee_lc[['year', 'class_name', 'acres', 'source']])
else:
    print('RUN_EE is False - skipping the Earth Engine path')

RUN_EE is False - skipping the Earth Engine path

For display of LCMS in the geeViz viewer, prefer the built-in thematic handling: the assets carry class names, values and palettes, so autoViz reads them directly.

# The map uses geeViz (pip install geeViz), which reads the class names
# and palette straight off the LCMS asset.
if RUN_EE:
    import geeViz.geeView as gv
    Map = gv.Map
    Map.clearMap()
    lc2025 = lcms_ee.lcms_ee_collection('Land_Cover').filter(
        ee.Filter.eq('year', 2025))
    Map.addLayer(lc2025.mosaic(), {'autoViz': True},
                 'LCMS Land Cover 2025')
    Map.centerObject(aoi)
    Map.view()
else:
    print('RUN_EE is False - skipping the map')

RUN_EE is False - skipping the map

5. Putting them side by side - carefully

align.compare_area reports both figures and their difference, and returns the caveats with the numbers.

The difference is an observation, not an error term. A 15% gap does not mean either source is 15% wrong: FIA’s ‘forest land’ is a land-use definition based on stocking and potential, while LCMS land cover describes present canopy - so a recently harvested stand stays forest land in FIA while LCMS may map it as grass the same year.

result = None
if FIA_UP:
    result = align.compare_area(wc=412022, state='Oregon',
                             county='Crook', year=2025)
    print(align.summarize_comparison(result))
else:
    print('skipped - the FIA half of the comparison is unavailable')

# The LCMS half stands on its own either way:
treed = align.lcms_tree_area(state='Oregon', county='Crook', year=2025)
treed[['year', 'acres', 'source', 'estimator']]

LCMS maps 596,912 acres of tree-bearing land cover. FIA estimates 1,034,367 acres of forest land (SE 5.9%, 303 plots), a difference of -42.3%. These are different estimators of related but distinct quantities - map-derived cover versus a design-based land-use estimate - so the gap is not an error in either. FIA counts recently harvested stands as forest land; LCMS maps present canopy.
year acres source estimator
0 2025 596912.429764 lcms wall-to-wall map (map accuracy)

The caveats travel with the numbers. A number this easy to quote is a number that travels without its footnotes unless you attach them.

# Available with or without a live FIA call - these are static
# facts about the two datasets, not properties of one comparison.
for c in align.COMPARISON_CAVEATS:
    print('-', c)
    print()

- LCMS area is map-derived; FIA area is a design-based estimate. The difference mixes map accuracy with sampling error and is not an error term for either.

- FIA 'forest land' is a land-use definition based on stocking and potential; LCMS land cover describes present canopy. A recently harvested stand stays forest land in FIA while LCMS may map it as grass or barren the same year.

- Reference periods differ: an FIA evaluation spans several years of panels, while an LCMS year is a single annual map.

- Minimum mapping unit and edge handling differ, which matters most in fragmented landscapes.

6. Cache management

Vocabularies are cached under ~/.georest/usfs with a 30-day stale-while-revalidate, falling back to the bundled snapshot and then to an empty result. A lookup never blocks on the network and never raises.

refresh_all() is the escape hatch for ‘a new release just dropped’ - not the normal path.

from pathlib import Path
print('cache dir: ~/' + vocab.cache_dir().relative_to(Path.home()).as_posix())
# vocab.refresh_all()   # uncomment to force a refresh

cache dir: ~/.georest/usfs