Add pyfia.carbon: NSVB-based forest carbon estimation (all 6 IPCC pools + stock-change)

.gitignore

@@ -165,6 +165,12 @@ data/
 # Documentation build
 site/
 data/
+# Exception: package-bundled NSVB coefficient tables (vendored from GTR-WO-104 Supp1)
+!src/pyfia/carbon/nsvb/data/
+!src/pyfia/carbon/nsvb/data/**
+# Exception: NGHGI report-target CSVs (frozen EPA Chapter 6 published numbers)
+!src/pyfia/carbon/data/
+!src/pyfia/carbon/data/nghgi_*.csv
 
 # JOSS paper (kept separate from repo)
 paper.md

CHANGELOG.md

@@ -5,6 +5,36 @@ All notable changes to pyFIA will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+
+### Added
+- **`live_tree()` function** — NSVB live tree carbon estimation:
+  - Recomputes above-ground biomass from scratch using the NSVB framework (Westfall et al. 2023, GTR-WO-104)
+  - Species-specific S10a carbon fractions (0.40–0.55) replace the flat 0.47 multiplier used by `biomass()`
+  - 3-level coefficient lookup precedence: Bailey DIVISION, species-level, Jenkins fallback
+  - Cull adjustment using Harmon et al. (2011) DECAYCD=3 density proportions
+  - `pool='ag'|'bg'|'total'` — AG via NSVB, BG bridges to FIADB `TREE.CARBON_BG`
+  - Validated against FIADB `TREE.CARBON_AG` on Georgia EVALID 132401 (130,952 trees): median per-tree relative error 0.085%
+- **`standing_dead()` function** — NSVB standing dead carbon estimation:
+  - Same NSVB biomass pipeline as `live_tree()`, plus decay-class reductions from `REF_TREE_DECAY_PROP`
+  - `DENSITY_PROP` x wood, `BARK_LOSS_PROP` x bark, `BRANCH_LOSS_PROP` x branch by hardwood/softwood x DECAYCD
+  - S10b dead-tree carbon fractions by hardwood/softwood x DECAYCD
+  - No `TREE.CULL` adjustment for dead trees (per FIADB Appendix K)
+  - `pool='ag'|'bg'|'total'` — same pool semantics as `live_tree()`
+  - Broken-top corrections: crown-proportion adjustment (Appendix K) + paraboloid volume-ratio for trees with `ACTUALHT < HT`
+  - Vendored Table S11 (`REF_TREE_STND_DEAD_CR_PROP`) for mean intact crown ratios by ecoregion province
+  - Validated against FIADB on Georgia EVALID 132401 (6,870 trees): median 10.9% per-tree relative error
+- **`pyfia.carbon` subpackage** — NSVB equation library, coefficient loaders, carbon fractions:
+  - `pyfia.carbon.nsvb.equations` — Models 1, 2, 4, 5, harmonization, vectorized pipelines
+  - `pyfia.carbon.nsvb.coefficients` — S1a–S8b coefficient tables, Bailey DIVISION lookup
+  - `pyfia.carbon.nsvb.carbon_fractions` — S10a (live), S10b (dead), `REF_TREE_DECAY_PROP` loaders
+  - Vendored coefficient CSVs from GTR-WO-104 supplementary archive
+- **NSVB validation gate** — `tests/validation/test_live_tree_nsvb.py` and `test_standing_dead_nsvb.py`:
+  - Per-tree parity tests against FIADB `CARBON_AG` on real Georgia inventory data
+  - Layered diagnostics: carbon rel-error, biomass ratio, FIADB-implied carbon fraction
+  - Ratchet thresholds that detect regressions and reward improvements
+  - EVALID-scoped to the current annual evaluation (avoids legacy CRM data contamination)
+
 ## [1.2.0] - 2025-01-18
 
 ### Added

CLAUDE.md

@@ -70,6 +70,7 @@ pyfia/
 ├── src/pyfia/           # Library source code
 │   ├── core/            # Database, backends, and settings
 │   ├── estimation/      # Statistical estimation
+│   ├── carbon/          # Carbon estimation (all 6 IPCC/NGHGI pools + total_ecosystem)
 │   ├── filtering/       # Domain filtering
 │   ├── downloader/      # FIA data download from DataMart
 │   ├── evalidator/      # EVALIDator API client for validation

README.md

@@ -33,17 +33,19 @@ pip install pyfia
 ```
 
 ```python
-from pyfia import FIA, biomass, tpa, volume, area
+from pyfia import FIA, tpa, volume, area, total_ecosystem
 
 with FIA("path/to/FIA_database.duckdb") as db:
     db.clip_by_state(37)  # North Carolina
     db.clip_most_recent(eval_type="VOL")
 
     # Core estimates
     trees = tpa(db, tree_domain="STATUSCD == 1")
-    carbon = biomass(db, by_species=True)
     timber = volume(db, land_type="timber")
     forest = area(db, land_type="forest")
+
+    # Total ecosystem carbon (all 6 IPCC/NGHGI pools)
+    carbon = total_ecosystem(db)
 ```
 
 ## Core Functions
@@ -54,6 +56,14 @@ with FIA("path/to/FIA_database.duckdb") as db:
 | `biomass()` | Above/belowground biomass | `biomass(db, by_species=True)` |
 | `volume()` | Merchantable volume (ft³) | `volume(db, land_type="timber")` |
 | `area()` | Forest land area | `area(db, grp_by="FORTYPCD")` |
+| `live_tree()` | NSVB live tree carbon | `live_tree(db, pool="ag")` |
+| `standing_dead()` | NSVB standing dead carbon | `standing_dead(db, pool="ag")` |
+| `understory()` | Understory vegetation carbon | `understory(db, pool="total")` |
+| `downed_dead()` | Downed dead wood carbon | `downed_dead(db)` |
+| `litter()` | Litter carbon | `litter(db)` |
+| `soil_organic()` | Soil organic carbon | `soil_organic(db)` |
+| `total_ecosystem()` | Total ecosystem carbon (all 6 pools) | `total_ecosystem(db)` |
+| `stock_change()` | Carbon stock change between periods | `stock_change(db, pool="all")` |
 | `site_index()` | Site productivity index | `site_index(db, grp_by="COUNTYCD")` |
 | `mortality()` | Annual mortality rates | `mortality(db)` |
 | `growth()` | Net growth estimation | `growth(db)` |

docs/DEVELOPMENT.md

@@ -47,6 +47,17 @@ pyfia/
 │   ├── fia.py         # Main FIA database class
 │   ├── data_reader.py # Efficient data loading
 │   └── backends/      # DuckDB, SQLite, MotherDuck backends
+├── carbon/             # NSVB carbon estimation
+│   ├── live_tree.py   # Live tree carbon (NSVB AG + FIADB BG bridge)
+│   ├── standing_dead.py # Standing dead carbon (NSVB + decay reductions)
+│   ├── understory.py  # Understory vegetation carbon (condition-level)
+│   ├── downed_dead.py # Downed dead wood carbon (condition-level)
+│   ├── litter.py      # Litter carbon (condition-level)
+│   ├── soil_organic.py # Soil organic carbon (condition-level)
+│   ├── total_ecosystem.py # Sum of all 6 pools
+│   ├── stock_change.py # Carbon stock change between inventory periods
+│   ├── data/          # Non-NSVB coefficient tables (Birdsey/Smith & Heath)
+│   └── nsvb/          # NSVB equation library, coefficients, carbon fractions
 ├── estimation/         # Statistical estimation
 │   ├── base.py        # BaseEstimator with Template Method pattern
 │   ├── grm.py         # GRM data loading and adjustment
@@ -85,7 +96,7 @@ pyfia/
 - Key methods: `clip_by_evalid()`, `clip_by_state()`, `clip_most_recent()`
 
 **Estimation Functions**
-- Simple API: `area()`, `biomass()`, `volume()`, `tpa()`, `mortality()`, `growth()`, `removals()`, `area_change()`, `site_index()`, `tree_metrics()`, `carbon_pools()`
+- Simple API: `area()`, `biomass()`, `volume()`, `tpa()`, `live_tree()`, `standing_dead()`, `understory()`, `downed_dead()`, `litter()`, `soil_organic()`, `total_ecosystem()`, `stock_change()`, `mortality()`, `growth()`, `removals()`, `area_change()`, `site_index()`, `tree_metrics()`
 - All support domain filtering, grouping, variance calculations
 - BaseEstimator uses Template Method for consistent workflow
 

docs/api/index.md

@@ -25,6 +25,8 @@ result = estimator(db, **options)  # Returns pl.DataFrame
 | [`volume()`](volume.md) | Estimate tree volume |
 | [`tpa()`](tpa.md) | Trees per acre and basal area |
 | [`biomass()`](biomass.md) | Tree biomass and carbon |
+| [`live_tree()`](live_tree.md) | NSVB live tree carbon |
+| [`standing_dead()`](standing_dead.md) | NSVB standing dead carbon |
 | [`site_index()`](site_index.md) | Area-weighted mean site index |
 | [`mortality()`](mortality.md) | Annual tree mortality |
 | [`growth()`](growth.md) | Annual tree growth |

docs/api/live_tree.md

@@ -0,0 +1,119 @@
+# Live Tree Carbon Estimation
+
+Estimate live tree carbon using the NSVB biomass framework with species-specific carbon fractions.
+
+## Overview
+
+The `live_tree()` function recomputes above-ground live tree biomass from scratch using the National Scale Volume and Biomass (NSVB) framework of Westfall et al. (2023, GTR-WO-104) and converts to carbon via species-specific S10a carbon fractions. This produces carbon estimates that align with the EPA NGHGI LULUCF live tree pool and match FIADB's pre-computed `CARBON_AG` column for NSVB-era inventories (September 2023 onward).
+
+```python
+import pyfia
+
+db = pyfia.FIA("georgia.duckdb")
+db.clip_by_state("GA")
+db.clip_most_recent(eval_type="VOL")
+
+# Above-ground live tree carbon
+result = pyfia.live_tree(db, pool="ag")
+
+# Total carbon (AG + BG bridge)
+total = pyfia.live_tree(db, pool="total")
+```
+
+## Function Reference
+
+::: pyfia.live_tree
+    options:
+      show_root_heading: true
+      show_source: true
+
+## Carbon Pools
+
+| Pool | Description | Method |
+|------|-------------|--------|
+| `"ag"` | Above-ground (default) | NSVB pipeline: stem wood + bark + branches, harmonized to total AGB, then multiplied by species-specific S10a carbon fractions |
+| `"bg"` | Below-ground (coarse roots) | Bridge to FIADB `TREE.CARBON_BG` (Phase 1 shortcut; native NSVB root model planned) |
+| `"total"` | AG + BG | NSVB above-ground + FIADB below-ground bridge |
+
+## How It Differs from `biomass()`
+
+| | `live_tree()` | `biomass()` |
+|---|---|---|
+| **Biomass source** | Recomputed from scratch via NSVB equations | Reads FIADB pre-computed `DRYBIO_*` columns |
+| **Carbon fraction** | Species-specific S10a (0.40-0.55) | Flat 0.47 multiplier |
+| **Coefficient lookup** | 3-level precedence (DIVISION, species, Jenkins) | N/A (pre-computed) |
+| **Cull adjustment** | NSVB cull formula with DECAYCD=3 density prop | Built into FIADB values |
+| **Transparency** | Full recompute, auditable | Black-box FIADB values |
+
+For NSVB-era inventories (2024+), both should agree closely. `live_tree()` is the preferred path for carbon accounting work that needs methodological transparency.
+
+## Technical Notes
+
+The NSVB pipeline predicts five biomass components per tree:
+
+1. Stem inside-bark wood volume (S1a) x wood density x 62.4 = gross wood weight
+2. Stem bark biomass (S6a)
+3. Branch biomass (S7a)
+4. Total above-ground biomass (S8a) - directly predicted
+
+The component sum is harmonized proportionally to the directly-predicted total AGB. Cull-reduced wood uses the Harmon et al. (2011) DECAYCD=3 density proportion (0.54 hardwood, 0.92 softwood). Carbon = harmonized AGB x species-specific S10a fraction.
+
+The optional `PLOTGEOM.ECOSUBCD` join activates Level 2 of the NSVB coefficient precedence (SPCD + Bailey DIVISION), closing a ~3% growing-stock biomass bias present in the species-level-only fallback. When `PLOTGEOM` is missing from older databases, the estimator falls back gracefully with a one-shot log warning.
+
+## Examples
+
+### Above-Ground Carbon Per Acre
+
+```python
+result = pyfia.live_tree(db, pool="ag")
+print(f"Carbon: {result['CARBON_ACRE'][0]:.2f} tons/acre")
+```
+
+### Carbon by Species
+
+```python
+result = pyfia.live_tree(db, pool="ag", by_species=True)
+result = pyfia.join_species_names(result, db)
+print(result.sort("CARBON_ACRE", descending=True).head(10))
+```
+
+### Carbon by Ownership Group
+
+```python
+result = pyfia.live_tree(
+    db,
+    pool="total",
+    grp_by="OWNGRPCD",
+    totals=True,
+    variance=True,
+)
+# OWNGRPCD: 10=National Forest, 20=Other Federal,
+#           30=State/Local, 40=Private
+print(result)
+```
+
+### Large Tree Carbon by Forest Type
+
+```python
+result = pyfia.live_tree(
+    db,
+    pool="ag",
+    grp_by="FORTYPCD",
+    tree_domain="DIA >= 20.0",
+)
+result = pyfia.join_forest_type_names(result, db)
+print(result)
+```
+
+### Carbon on Timberland with Standard Errors
+
+```python
+result = pyfia.live_tree(
+    db,
+    pool="ag",
+    land_type="timber",
+    variance=True,
+)
+print(f"Carbon: {result['CARBON_ACRE'][0]:.2f} +/- "
+      f"{result['CARBON_ACRE_SE'][0]:.2f} tons/acre")
+```