diff --git a/README.md b/README.md
index f60f5fe..e5e0def 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
This repo contains an experimental methodology and data pipeline for estimating the land use change (LUC) related emissions associated with agricultural commodities. This methodology allocates LUC emissions to crops in proportion to their displacement of natural ecosystems, based on high-resolution satellite imagery. It's primarily intended for use in corporate GHG inventories, and follows the new GHGP Land Sector and Removals Standard. It supports two attribution methodologies over a shared per-pixel emissions core: the high-resolution "jurisdictional direct land use change" (jdLUC) calculation where detailed crop maps exist, and a coarser-resolution statistical approach for regions where they don't.
-As a proof of concept, the jurisdictional-direct leg covers eleven crops grown in the United States (barley, dry beans, cotton, maize, potato, rice, sorghum, soybean, sugarbeet, sugarcane and wheat). Although U.S. land use change emissions are relatively modest contributors to global totals, the U.S. agricultural sector is well studied and has strong data infrastructure, which makes it a good place to start testing methods. The statistical leg extends the same per-pixel emissions core to global crop coverage, using coarser sub-national statistics (IFPRI MapSPAM) where high-resolution crop maps aren't available.
+As a proof of concept, the jurisdictional-direct leg covers eleven crops grown in the United States (barley, dry beans, cotton, maize, potato, rice, sorghum, soybean, sugarbeet, sugarcane and wheat). Although U.S. land use change emissions are relatively modest contributors to global totals, the U.S. agricultural sector is well studied and has strong data infrastructure, which makes it a good place to start testing methods. The statistical leg extends the same per-pixel emissions core to global crop coverage, using coarser sub-national statistics (IFPRI MapSPAM) where high-resolution crop maps aren't available. A third leg attributes by CROPGRIDS crop area, and reports no production yet.
> **Assessing LSRS conformance?** The [executive summary](docs/executive_summary.md) gives a concise account of how this methodology maps to the GHGP Land Sector and Removals Standard — which requirements it meets, the key modeling choices behind them, and where it deviates or remains a work in progress.
@@ -40,7 +40,7 @@ The data is licensed [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/).
Once you're ready to look under the hood:
- `docs/executive_summary.md` — the short, non-technical overview: what jdLUC is, why LUC emissions matter, and the headline takeaways. Start here if you're new.
-- `docs/methodology.md` — the high-level overview of the methodology: the datasets behind it, how we quantify per-pixel emissions, and how the jurisdictional-direct and statistical attribution legs produce emissions factors. The *what* we compute and *why*, from a scientific standpoint.
+- `docs/methodology.md` — the high-level overview of the methodology: the datasets behind it, how we quantify per-pixel emissions, and how the attribution legs produce emissions factors. The *what* we compute and *why*, from a scientific standpoint.
- `docs/architecture.md` — the system architecture and the rationale behind it: why we build on the [Pangeo](https://pangeo.io/) stack, how the five-stage pipeline is structured, and the tooling, storage, and caching choices that make it reproducible on a single host. The *how* it's built and *why those choices*.
- `docs/data.md` — how to get the data products, grids, and full schemas for the harmonized inputs, per-pixel emissions, and emissions-factor table.
- `docs/coverage.md` — the full list of countries and crops the pipeline produces emissions factors for, by ISO 3166-1 alpha-3 code.
@@ -81,7 +81,7 @@ export GOOGLE_APPLICATION_CREDENTIALS="${HOME}/.config/gcloud/application_defaul
### Running the pipeline
-Every stage is scoped by one or more **ISO 3166 alpha-3 country codes**, or by `--backfill` for a global run. Stages 4–5 additionally take a `--methodology-name` (`STATISTICAL` or `JURISDICTIONAL_DIRECT`). The example below reproduces Honduras via the statistical leg. Honduras is the cheapest end-to-end run in the corpus: the unit of work is a 10° tile, and it falls inside a single one (`20N_090W`). Swap in `USA` for the jurisdictional-direct leg, which needs the US-only CDL and NASS inputs:
+Every stage is scoped by one or more **ISO 3166 alpha-3 country codes**, or by `--backfill` for a global run. Stage 4 additionally takes a `--methodology-name` (`STATISTICAL`, `JURISDICTIONAL_DIRECT` or `JURISDICTIONAL_DIRECT_PROXY`), and stage 5 takes the first two. The example below reproduces Honduras via the statistical leg. Honduras is the cheapest end-to-end run in the corpus: the unit of work is a 10° tile, and it falls inside a single one (`20N_090W`). Swap in `USA` for the jurisdictional-direct leg, which needs the US-only CDL and NASS inputs:
```bash
# 1. Ingest each source dataset for the countries (positional: dataset, then ISO codes).
diff --git a/docs/architecture.md b/docs/architecture.md
index 6ef56a6..6c6592b 100644
--- a/docs/architecture.md
+++ b/docs/architecture.md
@@ -36,19 +36,19 @@ The pipeline is a linear chain of five stages, each consuming the previous stage
flowchart LR
ingest["ingest
sources"] --> harmonize["harmonize
common grid"]
harmonize --> emit["emit
per-pixel LUC"]
- emit --> attribute["attribute
per-crop rollup;
JD or statistical"]
+ emit --> attribute["attribute
per-crop rollup;
JD, JD proxy or statistical"]
attribute --> trace["trace
emissions factors"]
```
- **Ingest** — download every external source dataset into `INGEST_ROOT` as tiled Cloud-Optimized GeoTIFFs (COGs, raster), FlatGeobuf (vector), or parquet (tabular), tagged with provenance metadata, and with minimal changes.
- **Harmonize** — align the ingested sources onto a common grid, one 10° tile at a time, yielding one lazy `xarray.Dataset` per tile backed by a zarr.
- **Emit** — compute per-pixel land conversion events and emissions as a lazy dask/xarray graph and persist to zarr.
-- **Attribute** — clip per-pixel emissions to jurisdiction polygons and crop masks, rolling up to per-(jurisdiction, crop) totals via the jurisdictional-direct or statistical leg.
-- **Trace** — join the attribution rollups to production (NASS yields for the direct leg, MapSPAM for the statistical leg) to produce the final emissions-factor table.
+- **Attribute** — clip per-pixel emissions to jurisdiction polygons and crop masks, rolling up to per-(jurisdiction, crop) totals via the jurisdictional-direct, statistical or jurisdictional-direct proxy leg.
+- **Trace** — join the attribution rollups to production (NASS yields for the direct leg, MapSPAM for the statistical leg) to produce the final emissions-factor table. The proxy leg has no production yet.
## Compute platform and libraries
-Raster computation is expressed as lazy `xarray.Dataset` / `dask.array` graphs built with rioxarray and rasterio. GDAL VRTs align and resample every source onto the common grid; the statistical leg's step down to the MapSPAM grid is an overlap-weighted mean computed in dask (`geo.downscale_darray`). Nothing materializes until a stage writes its result, at which point dask streams the computation tile-by-tile into the output store.
+Raster computation is expressed as lazy `xarray.Dataset` / `dask.array` graphs built with rioxarray and rasterio. GDAL VRTs align and resample every source onto the common grid; the statistical and proxy legs' steps down to the MapSPAM and CROPGRIDS grids are overlap-weighted means computed in dask (`geo.downscale_darray`). Nothing materializes until a stage writes its result, at which point dask streams the computation tile-by-tile into the output store.
Each library earns its place:
@@ -80,7 +80,7 @@ The unifying theme is **chunked formats** — every store can be read and writte
All rasters are harmonized onto the **GLAD tile grid**: global `EPSG:4326` at 0.00025° (~30 m), tiled on the 10° graticule. We chose it because most of the source datasets are already published on it — tree-cover loss and the grassland class are native to it exactly — so warping *to* it leaves the layers that decide a conversion unresampled. It is also an intuitive grid and compatible with a tiled processing approach: the unit of work is a single 10° tile, and every grid is anchored on the absolute 10° graticule rather than on the extent of whatever was requested.
-A second, coarser grid (the ~10 km MapSPAM resolution) exists for the statistical attribution leg, which downsamples per-pixel emissions to match the resolution of the MapSPAM crop statistics rather than implying a precision the crop data does not have.
+Two coarser grids exist for attribution, the ~10 km MapSPAM resolution for the statistical leg and CROPGRIDS' 0.05° for the proxy leg, each matching its crop data rather than implying a precision that data does not have.
### One VRT per dataset per tile — chunked end to end
@@ -95,13 +95,13 @@ There is no reprojection step and no mosaic step:
A tile whose source was never published is tolerated when the caller opts in: each of the dataset's bands is emitted with a header and no source, which GDAL reads as all-nodata. A minority of tiles need this for the core inputs, and the emissions core coerces those nodatas to zero, so that "no data here" does not propagate as NaN. For a dataset whose publisher maps one crop rather than the whole land surface it is the majority case — Descals oil palm covers 48 of the 280 tiles — and there the band header carries no `NoDataValue` at all, because 0 is a mapped class, so the sourceless band reads back as zero directly.
-The emissions core harmonizes **two dataset lists and merges the results**: the primary stack that every conversion depends on, and the per-crop masks that widen its destination evidence (`emit.CROP_SUPPLEMENT_DATASET_NAMES`). They are kept apart because `harmonize.workflow` hashes its `dataset_names` into the cache key, so a single list would re-materialize every primary band across every tile each time a mask is admitted; separated, admitting one costs only its own cache. Both are harmonized onto the same tile grid, so the merge is an exact join.
+The emissions core harmonizes **two dataset lists and merges the results**: the primary stack that every conversion depends on, and the per-crop masks that widen its destination evidence (`harmonize.Stack.CROP_SUPPLEMENT`). They are kept apart because `harmonize.workflow` hashes its `dataset_names` into the cache key, so a single list would re-materialize every primary band across every tile each time a mask is admitted; separated, admitting one costs only its own cache. Both are harmonized onto the same tile grid, so the merge is an exact join.
Using VRTs rather than materializing intermediates is central to the design: the pipeline must be **fully chunked from end to end**, so that at no point is an entire raster dumped into RAM — doing so would cause massive slowdowns and out-of-memory failures. The VRTs are opened with rioxarray and assembled into one `xarray.Dataset`, one variable per fully-qualified source band (`{source}:{product}:{band}`).
### Per-dataset nodata, dtype, and resampling
-Resampling is a per-dataset property: nearest-neighbor/mode for categorical layers (GLAD land classes, CDL, peatland mask, climate zones, GACED30 cropland, GPW grassland class, and Descals oil palm planting year) and bilinear/average for intensive layers (SoilGrids, Huang BGB, Harris AGB, and GPW livestock density). Extensive layers would be summed when downsampling, but this functionality isn't currently implemented. The nodata, dtype, and resampling conventions vary across the external datasets and must be **manually specified and validated** per dataset — these details matter, and getting them wrong silently corrupts downstream results, which is another reason ingestion is isolated and independently inspectable.
+Resampling is a per-dataset property: nearest-neighbor for categorical layers (GLAD land classes, CDL, peatland mask, climate zones, GACED30 cropland, GPW grassland class, Descals oil palm planting year, and the FDaP masks) and bilinear/average for intensive layers (SoilGrids, Huang BGB, Harris AGB, and GPW livestock density). Nearest applies when downsampling too, since a regular subsample keeps each class's area where mode would shrink minority classes. Extensive layers would be summed when downsampling, but this functionality isn't currently implemented. The nodata, dtype, and resampling conventions vary across the external datasets and must be **manually specified and validated** per dataset — these details matter, and getting them wrong silently corrupts downstream results, which is another reason ingestion is isolated and independently inspectable.
## Caching, versioning, and provenance
diff --git a/docs/further_research.md b/docs/further_research.md
index 92ae5c8..3ec6b50 100644
--- a/docs/further_research.md
+++ b/docs/further_research.md
@@ -8,9 +8,9 @@ The tooling described in [`validation.md`](validation.md) measures the entries b
## Woody perennial crops resolve to no destination
-**Status**: oil palm is resolved. Descals et al. (2024) is admitted as a per-crop destination layer, so forest cleared for oil palm and standing as oil palm at the assessment year now resolves to cropland, takes its conversion emissions, and carries the peat pulse and occupation that ride on it. That is 352 Mt of the 1,071 Mt in the gap below. Against WRI, Indonesian oil palm now reads 0.429 where the country's other crops read 0.32-0.33, so it sits above its siblings rather than below them: the expansion-share rule is not biased against woody perennials, and what is left for Indonesia is a country-level shortfall every crop shares. The entry stands for the rest.
+**Status**: oil palm is resolved. Descals et al. (2024) is admitted as a per-crop destination layer, so forest cleared for oil palm and standing as oil palm at the assessment year now resolves to cropland, takes its conversion emissions, and carries the peat pulse and occupation that ride on it. That is 352 Mt of the 1,071 Mt in the gap below. Against WRI, Indonesian oil palm now reads 0.429 where the country's other crops read 0.32-0.33, so it sits above its siblings rather than below them: the expansion-share rule is not biased against woody perennials, and what is left for Indonesia is a country-level shortfall every crop shares. Cocoa and coffee have the FDaP masks between 24°S and 24°N, not yet measured. The entry stands for the rest.
-**Issue**: apart from oil palm, the conversion matrix names cropland and pasture as destinations and nothing else, so forest cleared for cocoa, coconut, rubber, coffee, tea or plantation timber resolves to no destination at all. Tree-cover loss identifies the clearance and dates it; GACED30 and Global Pasture Watch are then asked what stands there at the assessment year, and neither recognises a young plantation. The pixel's source carbon is charged to nobody and reported in `dropped-emissions` — bucket 4 of the six in `methodology.md`, the one LSRS Requirement 11 names explicitly.
+**Issue**: apart from oil palm, cocoa and coffee, the conversion matrix names cropland and pasture as destinations and nothing else, so forest cleared for coconut, rubber, tea or plantation timber resolves to no destination at all. Tree-cover loss identifies the clearance and dates it; GACED30 and Global Pasture Watch are then asked what stands there at the assessment year, and neither recognises a young plantation. The pixel's source carbon is charged to nobody and reported in `dropped-emissions` — bucket 4 of the six in `methodology.md`, the one LSRS Requirement 11 names explicitly.
Two things follow, and the second is what makes this worth a commodity-level entry rather than a line in the conformance table:
@@ -33,13 +33,13 @@ Sampled over 1°×1° windows at 30 m on GLAD's land-cover map, not on GACED30,
## FDaP masks match FAOSTAT's national area, not where it lies
-**Status**: open. No stage reads the masks yet, so this bounds what they can carry once one does.
+**Status**: open. The masks are `emit`'s cocoa and coffee destinations, so this bounds what they carry.
The masks threshold FDaP's per-pixel probabilities per country and year (`jdluc/datasets/fdap.py`). Each commodity has a floor, the model's precision/recall crossover, and no threshold goes below it. Where a country's fitted threshold is above the floor, it is set so the masked area meets FAOSTAT's area harvested; that matches the national total and says nothing about which pixels carry it. Where the fit is below the floor, the floor applies and the mask maps less than FAOSTAT's area: 34 coffee entries across 11 countries, from Colombia to Vietnam, and Venezuela's cocoa in 2020. Other producers use the floor alone, so their area is not matched to FAOSTAT's at all, and countries with no FAOSTAT production mask nothing. The 2024 masks are fitted to FAOSTAT's 2023 area.
Where even the highest threshold maps more than FAOSTAT's area harvested, the mask stays larger than FAOSTAT's figure: cocoa in Ghana (2020–2024), Liberia (2018), Fiji (2020 and 2022), Guyana (2022 and 2024) and Suriname, and coffee in Burundi (2018, 2022 and 2024), Kenya (2022), Malawi, Paraguay and Zimbabwe. Mozambique's coffee is excluded nationally instead, since its mask stays over 100 times FAOSTAT's area even at that threshold.
-**Potential impact:** Once the masks are a destination, an over-mapped country attributes more of its clearance to cocoa or coffee than FAOSTAT's area supports, and a producer whose mask maps less than FAOSTAT's area leaves the rest in `dropped-emissions`.
+**Potential impact:** An over-mapped country claims more clearance as cropland than FAOSTAT's cocoa or coffee area supports, and a producer whose mask maps less than FAOSTAT's area leaves the rest in `dropped-emissions`.
**Potential improvement path:** `tools/fit-fdap-thresholds.py` refits any country against the current FAOSTAT release: the 2024 masks against FAOSTAT's 2024 area, and the thresholds fitted to earlier FAOSTAT figures. Where the pixels lie needs an independent check against field plots, which a fit to national totals cannot provide.
@@ -51,7 +51,7 @@ Where even the highest threshold maps more than FAOSTAT's area harvested, the ma
Two consequences follow, and the first is why this entry is separate from the one above:
-1. **No crop mask closes it.** Every candidate on the register — FDaP, Descals coconut, the forest-management map — maps a woody perennial. The fallow rule is a definitional exclusion rather than a detection failure, so the fix is a cropland definition that admits long-fallow rotation, not another crop.
+1. **No crop mask closes it.** Every per-crop mask, admitted or candidate — FDaP, Descals coconut, the forest-management map — maps a woody perennial. The fallow rule is a definitional exclusion rather than a detection failure, so the fix is a cropland definition that admits long-fallow rotation, not another crop.
2. **A low anchor ratio here is not a detection failure.** PNG's sLUC/WRI for oil palm is 0.074 against 0.429 for Indonesia and 0.664 for Malaysia. PNG's tree-cover loss is dominated by selective logging and subsistence rotation rather than plantation conversion, so most of it should not be charged to any commodity. Reading that figure as a palm-detection gap sends the next reader after the wrong layer.
**Potential impact:** Unquantified globally. PNG is one measured instance; the 43.0 Mha figure above is the 5-arcmin upper bound across all countries, and the crops carrying it are staples rather than traded commodities, so the exposure is concentrated in jurisdictional totals rather than in any export factor.
diff --git a/docs/methodology.md b/docs/methodology.md
index ed886bd..4901272 100644
--- a/docs/methodology.md
+++ b/docs/methodology.md
@@ -8,12 +8,13 @@ This document is a high-level overview of the datasets and emissions methodology
- `cdl_glad_comparison_supplement.md` — a historical record of the GLAD-cropland restriction, which the pipeline no longer applies.
- `further_research.md` — known limitations and areas for further research.
-The pipeline quantifies the land-use-change (LUC) emissions associated with agricultural commodities and allocates them to specific crops. It supports two attribution methodologies that share a single per-pixel emissions core:
+The pipeline quantifies the land-use-change (LUC) emissions associated with agricultural commodities and allocates them to specific crops. It supports three attribution methodologies that share a single per-pixel emissions core:
- **Jurisdictional direct** — high-resolution crop maps link emissions to crops by spatial intersection. Currently **United States only** (USDA Cropland Data Layer).
- **Statistical** — where only coarse crop statistics exist, emissions are attributed in proportion to each crop's share of local cropland expansion. This leg is **global** (IFPRI MapSPAM).
+- **Jurisdictional direct proxy** — emissions are attributed by each crop's share of CROPGRIDS' 2020 crop area. This leg is **global**, and reports no production yet.
-Both legs follow the GHGP Land Sector and Removals Standard's 20-year linearly-discounted lookback, and both roll up to World Bank administrative jurisdictions.
+All three legs follow the GHGP Land Sector and Removals Standard's 20-year linearly-discounted lookback, and all three roll up to World Bank administrative jurisdictions.
### At a glance
@@ -30,7 +31,7 @@ graph TD
bgb[(Huang BGB)]
soc[(SoilGrids OCS)]
clim[(IPCC climate zone)]
- masks[(crop maps
CDL · MapSPAM)]
+ masks[(crop maps
CDL · MapSPAM · CROPGRIDS)]
yields[(crop yields
NASS)]
jurr[(administrative boundaries
World Bank)]
@@ -78,45 +79,48 @@ ______________________________________________________________________
Every input is ingested from its upstream publisher into cloud storage as tiled Cloud-Optimized GeoTIFFs (COGs, raster), FlatGeobuf (vector), or parquet (tabular), tagged with provenance metadata. One module per dataset lives under `jdluc/datasets/`. The ingestion and harmonization machinery that produces these analysis-ready copies is described in `architecture.md`.
-| Dataset | Role | Source | Kind |
-| ---------------------------------- | -------------------------------------------------------------------------------------------------- | ------------------------------------------- | ------- |
-| **GNW Hansen tree-cover loss** | Annual year of gross tree-cover loss, 2001-2025 — the forest source, and the only dated layer | GNW/Hansen GFC-2025-v1.13 | raster |
-| **GACED30 cropland extent** | Annual cropland extent, 2000-2024 — the cropland destination | Zenodo (doi:10.5281/zenodo.18199675) | raster |
-| **GPW grassland class** | Annual grassland class, 2000-2024 — the rangeland and pasture sources, and the pasture destination | Zenodo (doi:10.5281/zenodo.13890400) | raster |
-| **Descals oil palm** | Per-pixel oil palm planting year, 1989-2022 — the oil palm destination | Zenodo (doi:10.5281/zenodo.13379129) | raster |
-| **GNW Harris AGB (2000)** | Forest above-ground woody biomass | GNW data-api (WHRC AGB v1.4) | raster |
-| **Huang BGB** | Forest below-ground (root) biomass | Figshare (doi:10.6084/m9.figshare.12199637) | raster |
-| **SoilGrids OCS (0–30 cm)** | Soil organic carbon stock | ISRIC SoilGrids (WCS) | raster |
-| **GNW Global Peatlands** | Binary peatland mask | GNW data-api (`gfw_peatlands` v20230315) | raster |
-| **IPCC climate zones** | Climate domain per pixel (drives IPCC factors) | Zenodo (doi:10.5281/zenodo.7303808) | raster |
-| **USDA NASS CDL** | US per-pixel crop identity (jurisdictional-direct leg) | USDA Cropland Data Layer | raster |
-| **IFPRI MapSPAM** | Global per-crop physical area + production, 2000/2005/2010/2020 (statistical leg) | Harvard Dataverse | raster |
-| **USDA NASS QuickStats** | State-level crop yields (jurisdictional-direct production) | NASS QuickStats API | tabular |
-| **World Bank Official Boundaries** | Admin-0/1/2 jurisdiction polygons | World Bank | vector |
-| **FAOSTAT Production, livestock** | National livestock stocks and carcass-weight meat production (statistical leg) | FAO bulk download | tabular |
-| **FAOSTAT Livestock Patterns** | National livestock units per head of each grazer (statistical leg) | FAO bulk download | tabular |
-| **FAOSTAT Production, crops** | National crop production and harvested area (validation yardstick) | FAO bulk download | tabular |
-| **GLAD GLCLUC v2** | Land cover / land-use time series (2000–2020) — ingestable, no longer read by the emissions core | GLAD/Hansen GeoTIFFs | raster |
-| **GPW livestock headcount** | Annual livestock heads per hectare, 2000/2005/2010/2020 (statistical leg) | Zenodo (doi:10.5281/zenodo.17491242) | raster |
-| **FDaP cocoa and coffee** | Annual cocoa and coffee masks, 2018/2020/2022/2024, made at ingest from per-pixel probabilities | Forest Data Partnership 2026a (GCS) | raster |
-| **CROPGRIDS crop area** | Crop area of 173 crops in 2020, at 0.05° | Figshare (doi:10.6084/m9.figshare.22491997) | raster |
-
-The first nine rasters feed the per-pixel emissions core. CDL and MapSPAM feed the two attribution legs respectively, and the statistical leg also reads GPW livestock and FAOSTAT's livestock table for its livestock commodities. NASS yields and World Bank boundaries are joined downstream when building emissions factors.
-
-FAOSTAT's crop table, GLAD GLCLUC, the FDaP masks and CROPGRIDS are ingested but feed neither leg, and are listed so the inventory matches `datasets.DatasetName`. FAOSTAT's crop production is the national yardstick the validation tooling measures against — the global analogue of NASS QuickStats, and deliberately outside the attribution path so that no leg can consume the number it is measured by (see `validation.md`). Its livestock table is the exception: it prices the livestock commodities' production (see `validation.md`). GLAD GLCLUC is the land-cover series the emissions core read before the three annual layers replaced it, kept ingestable because the comparisons in `further_research.md` measure against it. The FDaP masks are candidate destinations for cocoa and coffee, two of the woody crops that resolve to no destination (see `further_research.md`); `jdluc/datasets/fdap.py` describes how ingest turns the probabilities into masks. GPW livestock is Global Pasture Watch's FAOSTAT-adjusted annual headcount of buffalo, cattle, goats, horses and sheep. CROPGRIDS (Tang et al. 2024) combines 27 gridded crop maps, subnational statistics for 52 countries and FAOSTAT's national totals into each crop's 2020 crop area, the land the crop covers, in hectares per cell.
+| Dataset | Role | Source | Kind |
+| ---------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------- | ------------------------------------------- | ------- |
+| **GNW Hansen tree-cover loss** | Annual year of gross tree-cover loss, 2001-2025 — the forest source, and the only dated layer | GNW/Hansen GFC-2025-v1.13 | raster |
+| **GACED30 cropland extent** | Annual cropland extent, 2000-2024 — the cropland destination | Zenodo (doi:10.5281/zenodo.18199675) | raster |
+| **GPW grassland class** | Annual grassland class, 2000-2024 — the rangeland and pasture sources, and the pasture destination | Zenodo (doi:10.5281/zenodo.13890400) | raster |
+| **Descals oil palm** | Per-pixel oil palm planting year, 1989-2022 — the oil palm destination | Zenodo (doi:10.5281/zenodo.13379129) | raster |
+| **FDaP cocoa and coffee** | Annual cocoa and coffee masks, 2018/2020/2022/2024, made at ingest from per-pixel probabilities — the cocoa and coffee destinations | Forest Data Partnership 2026a (GCS) | raster |
+| **GNW Harris AGB (2000)** | Forest above-ground woody biomass | GNW data-api (WHRC AGB v1.4) | raster |
+| **Huang BGB** | Forest below-ground (root) biomass | Figshare (doi:10.6084/m9.figshare.12199637) | raster |
+| **SoilGrids OCS (0–30 cm)** | Soil organic carbon stock | ISRIC SoilGrids (WCS) | raster |
+| **GNW Global Peatlands** | Binary peatland mask | GNW data-api (`gfw_peatlands` v20230315) | raster |
+| **IPCC climate zones** | Climate domain per pixel (drives IPCC factors) | Zenodo (doi:10.5281/zenodo.7303808) | raster |
+| **USDA NASS CDL** | US per-pixel crop identity (jurisdictional-direct leg) | USDA Cropland Data Layer | raster |
+| **IFPRI MapSPAM** | Global per-crop physical area + production, 2000/2005/2010/2020 (statistical leg) | Harvard Dataverse | raster |
+| **USDA NASS QuickStats** | State-level crop yields (jurisdictional-direct production) | NASS QuickStats API | tabular |
+| **World Bank Official Boundaries** | Admin-0/1/2 jurisdiction polygons | World Bank | vector |
+| **FAOSTAT Production, livestock** | National livestock stocks and carcass-weight meat production (statistical and proxy legs) | FAO bulk download | tabular |
+| **FAOSTAT Livestock Patterns** | National livestock units per head of each grazer (statistical and proxy legs) | FAO bulk download | tabular |
+| **FAOSTAT Production, crops** | National crop production and harvested area (validation yardstick) | FAO bulk download | tabular |
+| **GLAD GLCLUC v2** | Land cover / land-use time series (2000–2020) — ingestable, not read by the emissions core | GLAD/Hansen GeoTIFFs | raster |
+| **GPW livestock headcount** | Annual livestock heads per hectare, 2000/2005/2010/2020 (statistical and proxy legs) | Zenodo (doi:10.5281/zenodo.17491242) | raster |
+| **CROPGRIDS crop area** | Physical area (not harvested area) of 173 crops in 2020, at 0.05° (jurisdictional-direct proxy leg) | Figshare (doi:10.6084/m9.figshare.22491997) | raster |
+
+The first ten rasters feed the per-pixel emissions core. CDL, MapSPAM and CROPGRIDS feed the jurisdictional-direct, statistical and proxy legs respectively, and the statistical and proxy legs also read GPW livestock and FAOSTAT's livestock tables for their livestock commodities. NASS yields and World Bank boundaries are joined downstream when building emissions factors.
+
+FAOSTAT's crop table and GLAD GLCLUC are ingested but feed no leg, and are listed so the inventory matches `datasets.DatasetName`. FAOSTAT's crop production is the national yardstick the validation tooling measures against — the global analogue of NASS QuickStats, and deliberately outside the attribution path so that no leg can consume the number it is measured by (see `validation.md`). Its livestock table is the exception: it prices the livestock commodities' production (see `validation.md`). GLAD GLCLUC is the land-cover series the emissions core read before the three annual layers replaced it, kept ingestable because the comparisons in `further_research.md` measure against it. GPW livestock is Global Pasture Watch's FAOSTAT-adjusted annual headcount of buffalo, cattle, goats, horses and sheep. CROPGRIDS (Tang et al. 2024) combines 27 gridded crop maps, subnational statistics for 52 countries and FAOSTAT's national totals into each crop's 2020 crop area, the land the crop covers, in hectares per cell.
GPW grassland is native to the GLAD tile grid: its published mosaics are 0.00025° with pixel edges on whole degrees, so a ten-degree tile is a windowed read at an integer offset with no resampling. Its dominant-class band separates cultivated grassland from natural/semi-natural grassland and, new in v2, open shrubland — a split the land-cover series it replaced could not express, since those codes encode cover fraction rather than vegetation type. Reported F1 from five-fold spatially blocked cross-validation is 0.64 for cultivated and 0.76 for natural/semi-natural; both are v1 figures, and open shrubland arrives in v2-beta unvalidated.
-Descals oil palm is the only per-crop extent layer in the inventory, and the only one that dates its own pixels: each 30 m pixel holds the year the plantation standing in 2021 was established, 0 where there is no oil palm. **It is read as extent, not as a date.** A pixel is an oil palm destination where its planting year is non-zero and no later than the assessment year; dating stays with tree-cover loss and the last grassland departure, as for every other conversion. The planting year does not date one, because a replant reads as the year it was replanted rather than as the year the land was first cleared — it bounds when the current crop was established, not when conversion happened. The upper bound is what tests the destination at the assessment year: a plantation established after 2020 is not yet evidence of what an earlier clearance was for, which excludes 2.2% of mapped palm globally. Extent accuracy is high on industrial plantations (producer's 91.0 ± 2.5%, user's 91.8 ± 1.2%) and much weaker on smallholders (71.4 ± 0.7% and 72.4 ± 1.8%); the planting year itself carries an RMSE of 2.65 years against field data. The publisher screens globally and then publishes only the 609 hundred-kilometre cells where oil palm was found, so 232 of the 280 ten-degree tiles touch no cell and are not ingested at all; a reader that opts into missing tiles sees them as zero, which is the right answer here, but it is the screen's recall that answers rather than a per-pixel classification. The companion 10 m industrial-versus-smallholder extent layer from the same record is not ingested: its palm footprint agrees with the planting-year layer to within 0.4%, so the only thing it adds is the two-way split, which would have to be mode-downsampled from 10 m onto the 27.8 m GLAD tile grid and would lose fragmented smallholder area in the process.
+Descals oil palm is one of three per-crop extent layers, beside FDaP's cocoa and coffee, and the only one that dates its own pixels: each 30 m pixel holds the year the plantation standing in 2021 was established, 0 where there is no oil palm. **It is read as extent, not as a date.** A pixel is an oil palm destination where its planting year is non-zero and no later than the assessment year; dating stays with tree-cover loss and the last grassland departure, as for every other conversion. The planting year does not date one, because a replant reads as the year it was replanted rather than as the year the land was first cleared — it bounds when the current crop was established, not when conversion happened. The upper bound is what tests the destination at the assessment year: a plantation established after 2020 is not yet evidence of what an earlier clearance was for, which excludes 2.2% of mapped palm globally. Extent accuracy is high on industrial plantations (producer's 91.0 ± 2.5%, user's 91.8 ± 1.2%) and much weaker on smallholders (71.4 ± 0.7% and 72.4 ± 1.8%); the planting year itself carries an RMSE of 2.65 years against field data. The publisher screens globally and then publishes only the 609 hundred-kilometre cells where oil palm was found, so 232 of the 280 ten-degree tiles touch no cell and are not ingested at all; a reader that opts into missing tiles sees them as zero, which is the right answer here, but it is the screen's recall that answers rather than a per-pixel classification. The companion 10 m industrial-versus-smallholder extent layer from the same record is not ingested: its palm footprint agrees with the planting-year layer to within 0.4%, so the only thing it adds is the two-way split.
+
+The FDaP masks are the cocoa and coffee destinations between 24°S and 24°N, read at 2020; `jdluc/datasets/fdap.py` describes how ingest thresholds them. Harmonize subsamples the 10 m masks onto the GLAD tile grid by nearest, which keeps their area.
MapSPAM is the newest addition and the one that makes the statistical leg possible: it downscales sub-national crop statistics to a ~10 km (5 arc-minute) grid via a cross-entropy allocation, publishing physical area and production per crop for 2000, 2005, 2010, and 2020 (note: **no 2015 snapshot**). Its crop taxonomy is coarser in earlier years (the 2000 snapshot reports only 21, partly grouped, crops), which the statistical leg reconciles to a common per-crop taxonomy before use (see §3).
-### Two grids
+### Three grids
-Every source is ingested in `EPSG:4326` and resampled onto one of two common grids, depending on which stage consumes them:
+Every source is ingested in `EPSG:4326` and resampled onto one of three common grids, depending on which stage consumes them:
- **GLAD tile grid** — `EPSG:4326` at 0.00025° (~30 m), tiled on the 10° graticule. Everything in the emissions core and the jurisdictional-direct leg lives here.
- **MapSPAM grid** — the coarser ~10 km MapSPAM resolution. The statistical leg downsamples per-pixel emissions to this grid to match the resolution of the MapSPAM crop statistics, rather than implying a 30 m precision the underlying crop data does not have.
+- **CROPGRIDS grid** — 0.05° (~5.5 km), which the proxy leg downsamples to likewise.
The choice of the GLAD tile grid as the common backbone, and the mechanics of resampling every source onto it one 10° tile at a time, are covered in `architecture.md`.
@@ -132,15 +136,15 @@ Three annual layers decide what happened on a pixel, each answering one question
A pixel's **source** is forest where tree-cover loss falls inside the lookback. Otherwise it is the grassland class the pixel most recently left: rangeland (natural or open shrubland) or pasture (cultivated grassland), whichever departure is the later. A tree-cover loss outranks a grassland departure, so a pixel that lost forest is never also charged grassland carbon.
-A pixel's **destination** is read at the assessment year alone, against three layers in a fixed order: cropland where Descals maps oil palm established by 2020, then cropland where GACED30 calls it cropland in 2020, then pasture where Global Pasture Watch calls it cultivated grassland. A per-crop layer outranks both general ones, being the more specific evidence, and cropland outranks pasture, so the resolved destinations are disjoint. The layers themselves are not: every layer that claims a pixel is recorded, so a contested pixel stays visible in the emitted `destination-dataset` band even though only one layer resolves it.
+A pixel's **destination** is read at the assessment year alone, against five layers in a fixed order: cropland where Descals maps oil palm established by 2020, then cropland where FDaP's 2020 cocoa or coffee mask holds it, then cropland where GACED30 calls it cropland in 2020, then pasture where Global Pasture Watch calls it cultivated grassland. A per-crop layer outranks both general ones, being the more specific evidence, and cropland outranks pasture, so the resolved destinations are disjoint. The layers themselves are not: every layer that claims a pixel is recorded, so a contested pixel stays visible in the emitted `destination-dataset` band even though only one layer resolves it.
Source and destination together name one of five conversions, each releasing the pools its row and column carry:
-| | **→ cropland** (Descals or GACED30 at 2020) | **→ pasture** (GPW cultivated at 2020) |
-| ----------------------------------------- | ------------------------------------------- | -------------------------------------- |
-| **forest →** (dated by tree-cover loss) | forest biomass + mineral soil | forest biomass |
-| **rangeland →** (dated by last departure) | grassland biomass + mineral soil | grassland biomass |
-| **pasture →** (dated by last departure) | grassland biomass + mineral soil | *churn, not a conversion* |
+| | **→ cropland** (Descals, FDaP or GACED30 at 2020) | **→ pasture** (GPW cultivated at 2020) |
+| ----------------------------------------- | ------------------------------------------------- | -------------------------------------- |
+| **forest →** (dated by tree-cover loss) | forest biomass + mineral soil | forest biomass |
+| **rangeland →** (dated by last departure) | grassland biomass + mineral soil | grassland biomass |
+| **pasture →** (dated by last departure) | grassland biomass + mineral soil | *churn, not a conversion* |
A pixel matching none of these emits nothing — there is no default row to fire on it. Pasture that stays pasture is churn rather than conversion, and rangeland that becomes natural rangeland again is not a conversion to any commodity. The soil column is where a pasture destination differs: it leaves the ground at its grassland reference state, so the mineral term is accounted and zero rather than absent, and on peat the drainage pulse replaces that term rather than adding to it (§ Carbon stocks and fluxes).
@@ -195,20 +199,20 @@ ______________________________________________________________________
## 3. Attribution and emissions factors
-Attribution turns the crop-agnostic per-pixel emissions layer into per-(jurisdiction, commodity) totals, and then into emissions factors. `attribute.py` dispatches each country to one of two methodologies, fanning out over the ten-degree tiles that country touches and summing the per-tile partials (see [`coverage.md`](coverage.md) for the full country list and the territorial extent rule), and `trace.py` converts the rollups into the final emissions-factor table. Both legs clip to provincial (World Bank admin-1) polygons.
+Attribution turns the crop-agnostic per-pixel emissions layer into per-(jurisdiction, commodity) totals, and then into emissions factors. `attribute.py` dispatches each country to one of three methodologies, fanning out over the ten-degree tiles that country touches and summing the per-tile partials (see [`coverage.md`](coverage.md) for the full country list and the territorial extent rule), and `trace.py` converts the rollups into the final emissions-factor table. All three legs clip to provincial (World Bank admin-1) polygons.
### Shared framing
-Both legs quantify and allocate LUC emissions to a crop from a specific region in the same way. The crop's sourcing year is specified and a 20-year lookback window is defined, then that window is split into spans for which land use is known at both the start and end. Within each span, conversions from high- to low-carbon-density states are identified and the associated carbon emissions are quantified (§2), attributing them to all production within the lookback window. To more strongly penalize production that follows land conversion more quickly, emissions are discounted using a linear temporal ramp. Because this step operates independently on every pixel, it applies globally and is shared between the direct and statistical methodologies.
+All three legs quantify and allocate LUC emissions to a crop from a specific region in the same way. The crop's sourcing year is specified and a 20-year lookback window is defined, then that window is split into spans for which land use is known at both the start and end. Within each span, conversions from high- to low-carbon-density states are identified and the associated carbon emissions are quantified (§2), attributing them to all production within the lookback window. To more strongly penalize production that follows land conversion more quickly, emissions are discounted using a linear temporal ramp. Because this step operates independently on every pixel, it applies globally and is shared by all three methodologies.
### Rows that name no commodity
-Each jurisdiction's table, in both legs, has two rows that are not commodities: `DROPPED` and `UNATTRIBUTED`. Neither has a production figure, so neither has an emissions factor, and neither changes any crop's factor. They are easy to mix up. The difference is whether their emissions count toward the jurisdiction's emissions total, which is the `emissions-per-hectare` layer added up over the jurisdiction:
+Each jurisdiction's table, in every leg, has two rows that are not commodities: `DROPPED` and `UNATTRIBUTED`. Neither has a production figure, so neither has an emissions factor, and neither changes any crop's factor. They are easy to mix up. The difference is whether their emissions count toward the jurisdiction's emissions total, which is the `emissions-per-hectare` layer added up over the jurisdiction:
- **`DROPPED` does not count toward the total.** It is carbon released by clearing forest or grassland where, as far as our maps show, the land did not become cropland or pasture: for example, it was built on, or it became a tree crop we have no map for (§2, buckets 3 to 5). Its row has only a total, `emissions_mt`.
- **`UNATTRIBUTED` counts toward the total.** It is emissions from land that did become cropland or pasture, but that none of the commodities we model takes, mostly land used for crops or grazing we don't model (§2, bucket 2). Its row breaks the emissions down the same way a commodity row does.
-So the commodity rows plus `UNATTRIBUTED` add up to the emissions total, and `DROPPED` is extra, on top of it. What ends up in `UNATTRIBUTED` differs between the two legs, and each leg's section below says what it is.
+So the commodity rows plus `UNATTRIBUTED` add up to the emissions total, and `DROPPED` is extra, on top of it. What ends up in `UNATTRIBUTED` differs between the legs, and each leg's section below says what it is.
### Jurisdictional direct (`jurisdictional_direct.py`)
@@ -255,13 +259,17 @@ Each quantity is split by its own distribution — production by the pooled prod
Where no reference year places the group in a pixel there is no basis for a split at all. Two of the six 2000 groups have a catch-all constituent (`GROUP_TO_RESIDUAL_NAME`: "other oil crops" and "other pulses"), which absorbs the whole unattributable remainder rather than handing a slice to every named sibling. The other four have none, so the total is split evenly across constituents as a last resort. That last resort is not rare: measured across the MapSPAM snapshots it divides 18.1 Mt of 2000 production and 4.84 Mha of area — about 7.7% of the four affected groups — among crops that may not grow in the pixel at all. Its effect on a finished factor is nonetheless bounded well under a percent, because the 2000 snapshot is the only thing the decomposition touches and it carries 6.25% of the emissions numerator and 3.13% of the production denominator. See "MapSPAM group decomposition falls back to an even split" in `further_research.md`.
-**Tropical woody-perennial factors other than oil palm are lower bounds.** GACED30 admits woody crops in principle, but measured against independent references it recovers only a small fraction of them — 1.1% and 0.3% of mapped oil palm on two Riau windows — so a woody perennial with no layer of its own is one of the destinations bucket 4 above is blind to, and the conversion emissions reaching coffee, coconut and rubber are far below what external references report. Those crops' factors should be read as lower bounds rather than estimates until a per-crop destination layer supplies what the general cropland mask cannot. Oil palm has one, so its factor is no longer bounded that way; `further_research.md` covers the mechanism for the rest, and separately a smaller allocation gap specific to oil palm. The asymmetry worth naming is that MapSPAM's expansion denominator *does* see these crops, so they dilute every other crop's share in the cells they grow in while collecting almost none of the emissions themselves.
+**Tropical woody-perennial factors other than oil palm are lower bounds.** GACED30 admits woody crops in principle, but measured against independent references it recovers only a small fraction of them — 1.1% and 0.3% of mapped oil palm on two Riau windows — so a woody perennial with no layer of its own is one of the destinations bucket 4 above is blind to, and the conversion emissions reaching coconut and rubber are far below what external references report. Those crops' factors should be read as lower bounds rather than estimates until a per-crop destination layer supplies what the general cropland mask cannot. Oil palm has one, and coffee has FDaP's between 24°S and 24°N; `further_research.md` covers the mechanism for the rest, and separately a smaller allocation gap specific to oil palm. The asymmetry worth naming is that MapSPAM's expansion denominator *does* see these crops, so they dilute every other crop's share in the cells they grow in while collecting almost none of the emissions themselves.
**Windowed production denominator (a departure from WRI).** The emissions numerator is time-resolved: it links conversion emissions to each crop's *historical* expansion, weighted by recency. Dividing that by a single *current-year* production snapshot — as in WRI's methodology — would place the numerator and denominator on inconsistent time bases. A crop that expanded early in the window and then contracted would carry real conversion emissions against a shrunken (or zero) present-day production, yielding an inflated or undefined emissions factor. Instead, production is reduced over the same 20-year lookback window using the same linear temporal ramp applied to emissions: each span's production is the mean of its two MapSPAM snapshot years, weighted by that span's discount weight and normalized by the sum of the weights. This ties emission allocation and production to the same years with the same recency weighting. When production is flat across the window the result reduces exactly to the current-production value; it diverges only when production actually changed. This consistent linkage is a material improvement over attributing historical-expansion emissions to present-day production alone.
+### Jurisdictional direct proxy (`jurisdictional_direct_proxy.py`)
+
+Each pixel's emissions are split by the land class it became, cropland or pasture, and averaged onto the CROPGRIDS grid. In each cell, each crop takes its share of the 2020 crop area of all 173 CROPGRIDS crops of what cropland is charged, and each livestock commodity takes its 2020 grazing share, as in the statistical leg, of what pasture is charged. CROPGRIDS has one year, so every span is divided by the same shares. A crop's hectares are its CROPGRIDS crop area, and `UNATTRIBUTED` is cropland's charge in cells with no CROPGRIDS crop.
+
### From rollups to emissions factors (`trace.py`)
-`trace.py` takes the attribution rollup and derives the emissions factor identically for both methodologies — the only methodology-specific step is where production comes from (NASS yield × area for direct; for statistical, MapSPAM production for crops and GPW heads × FAOSTAT carcass weight per head for livestock):
+`trace.py` takes the attribution rollup and derives the emissions factor identically for the jurisdictional-direct and statistical methodologies — the only methodology-specific step is where production comes from (NASS yield × area for direct; for statistical, MapSPAM production for crops and GPW heads × FAOSTAT carcass weight per head for livestock). The proxy leg has no production yet:
```
emissions_factor_kgco2e_per_kg = Σ emissions / Σ production
diff --git a/jdluc/__tests__/emit_test.py b/jdluc/__tests__/emit_test.py
index 87dc466..ac73a48 100644
--- a/jdluc/__tests__/emit_test.py
+++ b/jdluc/__tests__/emit_test.py
@@ -6,6 +6,9 @@
import xarray
from jdluc.datasets.descals_oil_palm import DATASET as OIL_PALM
+from jdluc.datasets.fdap import COCOA_DATASET as COCOA
+from jdluc.datasets.fdap import COFFEE_DATASET as COFFEE
+from jdluc.datasets.fdap import YEARS as FDAP_YEARS
from jdluc.datasets.gnw_global_peatlands import DATASET as PEATLANDS
from jdluc.datasets.gnw_tcl import DATASET as TREE_COVER_LOSS
from jdluc.datasets.gnw_tcl import LOSS_YEAR_OFFSET
@@ -639,6 +642,8 @@ def test_get_last_departure_year(is_rangeland_by_year: list[bool], year: int) ->
def get_dset_for_one_pixel(
grassland_by_year: list[float],
+ is_cocoa: bool,
+ is_coffee: bool,
is_cropland: bool,
loss_year: float | None,
planting_year: float,
@@ -656,6 +661,12 @@ def get_dset_for_one_pixel(
]: get_darray_for_data(
data=[[Cropland.CROPLAND if is_cropland else Cropland.NOT_CROPLAND]]
),
+ COCOA.fully_qualified_band_names[
+ FDAP_YEARS.index(ASSESSMENT_YEAR)
+ ]: get_darray_for_data(data=[[int(is_cocoa)]]),
+ COFFEE.fully_qualified_band_names[
+ FDAP_YEARS.index(ASSESSMENT_YEAR)
+ ]: get_darray_for_data(data=[[int(is_coffee)]]),
OIL_PALM.fully_qualified_band_name: get_darray_for_data(
data=[[planting_year]]
),
@@ -893,6 +904,8 @@ def test_get_conversion_record(
result = get_conversion_record(
dset=get_dset_for_one_pixel(
grassland_by_year=grassland_by_year,
+ is_cocoa=False,
+ is_coffee=False,
is_cropland=is_cropland,
loss_year=loss_year,
planting_year=0,
@@ -904,6 +917,8 @@ def test_get_conversion_record(
@pytest.mark.parametrize(
(
+ "is_cocoa",
+ "is_coffee",
"is_cropland",
"is_cultivated",
"planting_year",
@@ -913,6 +928,8 @@ def test_get_conversion_record(
),
(
pytest.param(
+ False,
+ False,
True,
False,
0,
@@ -922,6 +939,8 @@ def test_get_conversion_record(
id="GACED30 alone claims it",
),
pytest.param(
+ False,
+ False,
False,
True,
0,
@@ -931,6 +950,8 @@ def test_get_conversion_record(
id="GPW alone claims it",
),
pytest.param(
+ False,
+ False,
True,
True,
0,
@@ -939,8 +960,12 @@ def test_get_conversion_record(
False,
id="both claim it, and the earlier member resolves it",
),
- pytest.param(False, False, 0, 0, False, False, id="neither claims it"),
pytest.param(
+ False, False, False, False, 0, 0, False, False, id="neither claims it"
+ ),
+ pytest.param(
+ False,
+ False,
False,
False,
2012,
@@ -950,6 +975,8 @@ def test_get_conversion_record(
id="oil palm claims what no other layer sees",
),
pytest.param(
+ False,
+ False,
False,
True,
2012,
@@ -959,6 +986,8 @@ def test_get_conversion_record(
id="oil palm outranks pasture on a pixel both claim",
),
pytest.param(
+ False,
+ False,
False,
False,
ASSESSMENT_YEAR,
@@ -968,6 +997,8 @@ def test_get_conversion_record(
id="a plantation established in the assessment year still counts",
),
pytest.param(
+ False,
+ False,
False,
False,
ASSESSMENT_YEAR + 1,
@@ -976,9 +1007,46 @@ def test_get_conversion_record(
False,
id="a plantation established after it does not",
),
+ pytest.param(
+ True,
+ False,
+ False,
+ False,
+ 0,
+ DestinationDataset.FDAP_COCOA,
+ True,
+ False,
+ id="cocoa claims what no other layer sees",
+ ),
+ pytest.param(
+ False,
+ True,
+ False,
+ True,
+ 0,
+ DestinationDataset.FDAP_COFFEE | DestinationDataset.GPW_GRASSLAND,
+ True,
+ False,
+ id="coffee outranks pasture on a pixel both claim",
+ ),
+ pytest.param(
+ True,
+ True,
+ True,
+ False,
+ 0,
+ DestinationDataset.FDAP_COCOA
+ | DestinationDataset.FDAP_COFFEE
+ | DestinationDataset.LIAO_GACED30,
+ True,
+ False,
+ id="cocoa, coffee and GACED30 all keep their bits",
+ ),
),
)
def test_get_conversion_record_destination(
+ is_cocoa: bool,
+ is_coffee: bool,
is_cropland: bool,
is_cultivated: bool,
planting_year: float,
@@ -993,6 +1061,8 @@ def test_get_conversion_record_destination(
*NO_GRASSLAND[:-1],
float(Grassland.CULTIVATED if is_cultivated else Grassland.OTHER),
],
+ is_cocoa=is_cocoa,
+ is_coffee=is_coffee,
is_cropland=is_cropland,
loss_year=2012,
planting_year=planting_year,
@@ -1108,6 +1178,8 @@ def test_the_halves_a_conversion_record_resolves(
result = get_conversion_record(
dset=get_dset_for_one_pixel(
grassland_by_year=grassland_by_year,
+ is_cocoa=False,
+ is_coffee=False,
is_cropland=is_cropland,
loss_year=loss_year,
planting_year=0,
@@ -1152,6 +1224,8 @@ def test_a_conversion_record_tells_pasture_from_rangeland_sources(
result = get_conversion_record(
dset=get_dset_for_one_pixel(
grassland_by_year=grassland_by_year,
+ is_cocoa=False,
+ is_coffee=False,
is_cropland=True,
loss_year=loss_year,
planting_year=0,
diff --git a/jdluc/__tests__/harmonize_test.py b/jdluc/__tests__/harmonize_test.py
index 32e8725..4e0e0b5 100644
--- a/jdluc/__tests__/harmonize_test.py
+++ b/jdluc/__tests__/harmonize_test.py
@@ -92,7 +92,7 @@ def test_grid_get_resolution_for_world() -> None:
(
(BandType.CATEGORICAL, tiling.XY(1, 1), "nearest"),
(BandType.CATEGORICAL, tiling.XY(2, 2), "nearest"),
- (BandType.CATEGORICAL, tiling.XY(3, 3), "mode"),
+ (BandType.CATEGORICAL, tiling.XY(3, 3), "nearest"),
(BandType.INTENSIVE, tiling.XY(1, 1), "bilinear"),
(BandType.INTENSIVE, tiling.XY(2, 2), "nearest"),
(BandType.INTENSIVE, tiling.XY(3, 3), "average"),
diff --git a/jdluc/__tests__/jurisdictional_direct_proxy_test.py b/jdluc/__tests__/jurisdictional_direct_proxy_test.py
new file mode 100644
index 0000000..7c26424
--- /dev/null
+++ b/jdluc/__tests__/jurisdictional_direct_proxy_test.py
@@ -0,0 +1,263 @@
+import numpy
+import pytest
+import xarray
+
+from jdluc.datasets.gnw_global_peatlands import DATASET as PEATLANDS
+from jdluc.datasets.gpw_grassland import DATASET as GRASSLAND
+from jdluc.datasets.gpw_grassland import YEARS as GRASSLAND_YEARS
+from jdluc.datasets.gpw_grassland import Grassland
+from jdluc.datasets.tang_cropgrids import DATASET as CROPGRIDS
+from jdluc.datasets.tang_cropgrids import YEAR as CROPGRIDS_YEAR
+from jdluc.datasets.tang_cropgrids import Crop
+from jdluc.emit import (
+ CROPLAND,
+ PASTURELAND,
+ PEATLAND_EMISSIONS_ANNUAL_TCO2E_PER_HA,
+ SPAN_TO_LINEAR_DISCOUNT_WEIGHT,
+ Conversion,
+ EmissionComponent,
+ NonCommodity,
+)
+from jdluc.jurisdictional_direct_proxy import (
+ get_commodity_name_to_totals,
+ get_crop_to_area_share,
+ get_name_to_band,
+)
+from jdluc.statistical import Livestock
+
+# The discount weights add up to 0.2, so 50 and 500 per span discount to 10 and 100
+VEGETATION_PER_SPAN = 50.0
+SOIL_PER_SPAN = 500.0
+LAND_CLASS_COMPONENT_NAMES = [
+ f"{land_class:s}-{component!s}-emissions:tco2e-per-ha"
+ for land_class in (CROPLAND, PASTURELAND)
+ for component in EmissionComponent
+]
+
+
+def get_darray_for_cells(data: list[list[float]]) -> xarray.DataArray:
+ arr = numpy.array(data, dtype=numpy.float64)
+ y, x = arr.shape
+ return xarray.DataArray(
+ coords={"y": range(y), "x": range(x)}, data=arr, dims=("y", "x")
+ )
+
+
+def get_dset_for_pixel(
+ conversion: Conversion, grassland: float, is_peat: bool
+) -> xarray.Dataset:
+ """One pixel, charged the same emissions in every span."""
+ values = {
+ "conversion": float(conversion),
+ f"{CROPLAND:s}-peatland-occupation:tco2e-per-ha": 0.0,
+ "dropped-emissions:tco2e-per-ha": 0.0,
+ GRASSLAND.fully_qualified_band_names[
+ GRASSLAND_YEARS.index(CROPGRIDS_YEAR)
+ ]: grassland,
+ f"{PASTURELAND:s}-peatland-occupation:tco2e-per-ha": 0.0,
+ PEATLANDS.fully_qualified_band_name: float(is_peat),
+ } | {
+ name: value
+ for before, after in SPAN_TO_LINEAR_DISCOUNT_WEIGHT
+ for name, value in (
+ (f"soil-emissions:tco2e-per-ha:{before:d}-{after:d}", SOIL_PER_SPAN),
+ (
+ f"vegetation-emissions:tco2e-per-ha:{before:d}-{after:d}",
+ VEGETATION_PER_SPAN,
+ ),
+ )
+ }
+ return xarray.Dataset(
+ {name: get_darray_for_cells(data=[[value]]) for name, value in values.items()}
+ )
+
+
+def get_dset_for_crop_areas(
+ crop_to_hectares: dict[Crop, list[list[float]]],
+) -> xarray.Dataset:
+ """Every CROPGRIDS crop band, zero unless given."""
+ zeros = numpy.zeros_like(next(iter(crop_to_hectares.values()))).tolist()
+ crop_to_data = dict.fromkeys(Crop, zeros) | crop_to_hectares
+ return xarray.Dataset(
+ {
+ band_name: get_darray_for_cells(data=crop_to_data[crop])
+ for crop, band_name in zip(
+ Crop, CROPGRIDS.fully_qualified_band_names, strict=True
+ )
+ }
+ )
+
+
+@pytest.mark.parametrize(
+ ("conversion", "is_peat", "name_to_expected"),
+ (
+ pytest.param(
+ Conversion.FOREST_TO_CROPLAND,
+ False,
+ {f"{CROPLAND:s}-forest-emissions:tco2e-per-ha": 110.0},
+ id="forest-to-cropland-charges-cropland",
+ ),
+ pytest.param(
+ Conversion.FOREST_TO_PASTURE,
+ False,
+ {f"{PASTURELAND:s}-forest-emissions:tco2e-per-ha": 110.0},
+ id="forest-to-pasture-charges-pastureland",
+ ),
+ pytest.param(
+ Conversion.RANGELAND_TO_CROPLAND,
+ True,
+ {
+ f"{CROPLAND:s}-grassland-emissions:tco2e-per-ha": 10.0,
+ f"{CROPLAND:s}-peatland-conversion-emissions:tco2e-per-ha": 100.0,
+ },
+ id="peat-soil-stays-with-the-cropland-it-became",
+ ),
+ pytest.param(
+ Conversion.RANGELAND_TO_PASTURE,
+ True,
+ {
+ f"{PASTURELAND:s}-grassland-emissions:tco2e-per-ha": 10.0,
+ f"{PASTURELAND:s}-peatland-conversion-emissions:tco2e-per-ha": 100.0,
+ },
+ id="peat-soil-stays-with-the-pasture-it-became",
+ ),
+ ),
+)
+def test_get_name_to_band_splits_each_component_by_land_class(
+ conversion: Conversion, is_peat: bool, name_to_expected: dict[str, float]
+) -> None:
+ name_to_band = get_name_to_band(
+ dset=get_dset_for_pixel(
+ conversion=conversion, grassland=float(Grassland.OTHER), is_peat=is_peat
+ )
+ )
+ # Every other land class and component is charged nothing
+ assert {
+ name: name_to_band[name].item() for name in LAND_CLASS_COMPONENT_NAMES
+ } == pytest.approx(
+ dict.fromkeys(LAND_CLASS_COMPONENT_NAMES, 0.0) | name_to_expected
+ )
+
+
+@pytest.mark.parametrize("conversion", Conversion)
+@pytest.mark.parametrize("is_peat", (False, True))
+def test_get_name_to_band_charges_every_conversion_exactly_once(
+ conversion: Conversion, is_peat: bool
+) -> None:
+ name_to_band = get_name_to_band(
+ dset=get_dset_for_pixel(
+ conversion=conversion, grassland=float(Grassland.OTHER), is_peat=is_peat
+ )
+ )
+ charged = sum(name_to_band[name].item() for name in LAND_CLASS_COMPONENT_NAMES)
+ assert charged == pytest.approx(0.0 if conversion == Conversion.NONE else 110.0)
+
+
+@pytest.mark.parametrize(
+ ("grassland", "expected"),
+ (
+ pytest.param(
+ float(Grassland.CULTIVATED), 1.0, id="cultivated-grassland-is-pasture"
+ ),
+ pytest.param(float(Grassland.NATURAL), 0.0, id="natural-grassland-is-not"),
+ pytest.param(numpy.nan, 0.0, id="no-observation-reads-as-not-pasture"),
+ ),
+)
+def test_get_name_to_band_pasture_fraction(grassland: float, expected: float) -> None:
+ name_to_band = get_name_to_band(
+ dset=get_dset_for_pixel(
+ conversion=Conversion.NONE, grassland=grassland, is_peat=False
+ )
+ )
+ assert name_to_band[f"pasture:fraction:{CROPGRIDS_YEAR:d}"].item() == expected
+
+
+def test_get_crop_to_area_share_divides_by_every_crop_cropgrids_maps() -> None:
+ result = get_crop_to_area_share(
+ crops=(Crop.COCOA,),
+ dset=get_dset_for_crop_areas(
+ crop_to_hectares={
+ # Cocoa beside coffee, a cell with no crop at all, and the ocean
+ Crop.COCOA: [[1.0, 0.0, numpy.nan]],
+ Crop.COFFEE: [[3.0, 0.0, numpy.nan]],
+ }
+ ),
+ )
+ assert set(result) == {Crop.COCOA}
+ assert result[Crop.COCOA].values.tolist() == [[0.25, 0.0, 0.0]]
+
+
+def test_get_commodity_name_to_totals() -> None:
+ """Two cells of 2 and 3 hectares, one row as a province one cell tall clips to.
+
+ Cocoa grows on the first cell and no crop on the second; half of each is pasture. Every band
+ charges both cells alike, and each land class and component charges its own amount, so a total
+ reading the wrong band is visible. Cropland's peat occupation drains the whole cell and
+ pasture's half of it.
+ """
+ per_hectare = {
+ f"{CROPLAND:s}-forest-emissions:tco2e-per-ha": 100.0,
+ f"{CROPLAND:s}-grassland-emissions:tco2e-per-ha": 20.0,
+ f"{CROPLAND:s}-peatland-conversion-emissions:tco2e-per-ha": 3.0,
+ f"{CROPLAND:s}-peatland-occupation:tco2e-per-ha": PEATLAND_EMISSIONS_ANNUAL_TCO2E_PER_HA,
+ "dropped-emissions:tco2e-per-ha": 9.0,
+ f"pasture:fraction:{CROPGRIDS_YEAR:d}": 0.5,
+ f"{PASTURELAND:s}-forest-emissions:tco2e-per-ha": 4000.0,
+ f"{PASTURELAND:s}-grassland-emissions:tco2e-per-ha": 500.0,
+ f"{PASTURELAND:s}-peatland-conversion-emissions:tco2e-per-ha": 60.0,
+ f"{PASTURELAND:s}-peatland-occupation:tco2e-per-ha": PEATLAND_EMISSIONS_ANNUAL_TCO2E_PER_HA
+ / 2,
+ }
+ dset = get_dset_for_crop_areas(crop_to_hectares={Crop.COCOA: [[50.0, 0.0]]}).assign(
+ {
+ name: get_darray_for_cells(data=[[value, value]])
+ for name, value in per_hectare.items()
+ }
+ | {"hectares-per-cell:ha": get_darray_for_cells(data=[[2.0, 3.0]])}
+ )
+ result = get_commodity_name_to_totals(
+ commodities=(Crop.COCOA, Livestock.BEEF_CATTLE),
+ commodity_to_share={
+ Crop.COCOA: get_darray_for_cells(data=[[1.0, 0.0]]),
+ # PASTURE's quarter is not asked for
+ Livestock.BEEF_CATTLE: get_darray_for_cells(data=[[0.75, 0.75]]),
+ },
+ dset=dset,
+ )
+ assert result == {
+ # Cropland's charge on the first cell's 2 ha, and CROPGRIDS' area rather than a share of it
+ Crop.COCOA.name: pytest.approx(
+ {
+ "commodity_hectares": 50.0,
+ "emissions_mt": 320.6,
+ "forest_emissions_mt": 200.0,
+ "grassland_emissions_mt": 40.0,
+ "peatland_commodity_hectares": 2.0,
+ "peatland_conversion_emissions_mt": 6.0,
+ "peatland_occupation_emissions_mt": 74.6,
+ }
+ ),
+ # Pasture's charge on three quarters of both cells' 5 ha, and of their 2.5 ha of pasture
+ Livestock.BEEF_CATTLE.name: pytest.approx(
+ {
+ "commodity_hectares": 1.875,
+ "emissions_mt": 17169.9375,
+ "forest_emissions_mt": 15000.0,
+ "grassland_emissions_mt": 1875.0,
+ "peatland_commodity_hectares": 1.875,
+ "peatland_conversion_emissions_mt": 225.0,
+ "peatland_occupation_emissions_mt": 69.9375,
+ }
+ ),
+ NonCommodity.DROPPED.name: pytest.approx({"emissions_mt": 45.0}),
+ # Cropland's charge on the second cell's 3 ha, plus the quarter of pasture's no row takes
+ NonCommodity.UNATTRIBUTED.name: pytest.approx(
+ {
+ "emissions_mt": 6204.2125,
+ "forest_emissions_mt": 5300.0,
+ "grassland_emissions_mt": 685.0,
+ "peatland_conversion_emissions_mt": 84.0,
+ "peatland_occupation_emissions_mt": 135.2125,
+ }
+ ),
+ }
diff --git a/jdluc/attribute.py b/jdluc/attribute.py
index 448ecbb..ac36074 100644
--- a/jdluc/attribute.py
+++ b/jdluc/attribute.py
@@ -5,7 +5,10 @@
- STATISTICAL (default): downscales the per-pixel emissions to the MAPSPAM grid and
attributes them across commodities by their shares of each cell's expansion.
- JURISDICTIONAL_DIRECT: masks the per-pixel emissions to each crop's CDL codes.
-Both clip to provincial (World Bank admin-1) polygons and sum commodity area, peatland
+ - JURISDICTIONAL_DIRECT_PROXY: downscales the per-pixel emissions to the CROPGRIDS grid, split
+ by whether each pixel became cropland or pasture, and divides cropland's part among CROPGRIDS'
+ crops by their 2020 crop-area shares and pasture's part among grazers.
+All three clip to provincial (World Bank admin-1) polygons and sum commodity area, peatland
commodity area, peatland-occupation emissions, and total emissions. Returns a pandas.DataFrame indexed by
(admin_level, admin_id, commodity_name, methodology); the per-(country, tile) sub-workflows are
cached and `--concurrency` bounds how many tiles are in flight at once.
@@ -28,15 +31,17 @@
from jdluc import (
jurisdictional_direct,
+ jurisdictional_direct_proxy,
statistical,
)
-from jdluc.datasets import worldbank_jurisdictions
+from jdluc.datasets import tang_cropgrids, worldbank_jurisdictions
logger = logging.getLogger(__name__)
class Methodology(enum.IntEnum):
JURISDICTIONAL_DIRECT = enum.auto()
+ JURISDICTIONAL_DIRECT_PROXY = enum.auto()
STATISTICAL = enum.auto()
@@ -59,11 +64,15 @@ def workflow(
iso_3166s: collections.abc.Iterable[str],
methodology: Methodology,
) -> pandas.DataFrame:
- workflow_for_tile = (
- jurisdictional_direct.workflow
- if methodology == Methodology.JURISDICTIONAL_DIRECT
- else statistical.workflow
- )
+ match methodology:
+ case Methodology.JURISDICTIONAL_DIRECT:
+ workflow_for_tile = jurisdictional_direct.workflow
+ case Methodology.JURISDICTIONAL_DIRECT_PROXY:
+ workflow_for_tile = jurisdictional_direct_proxy.workflow
+ case Methodology.STATISTICAL:
+ workflow_for_tile = statistical.workflow
+ case _:
+ raise ValueError(methodology)
# NB: have each thread handle a tile, to avoid races in harmonize, emit, or downscale
tile_id_to_iso_3166s: dict[str, list[str]] = collections.defaultdict(list)
@@ -124,13 +133,19 @@ def workflow_for_tile_id(tile_id: str) -> list[pandas.DataFrame]:
def get_commodity_names(methodology: Methodology) -> tuple[str, ...]:
- if methodology == Methodology.STATISTICAL:
- return tuple(
- sorted(e.name for e in (*statistical.Crop, *statistical.Livestock))
- )
- else:
- assert methodology == Methodology.JURISDICTIONAL_DIRECT
- return tuple(sorted(e.name for e in jurisdictional_direct.Crop))
+ match methodology:
+ case Methodology.JURISDICTIONAL_DIRECT:
+ return tuple(sorted(e.name for e in jurisdictional_direct.Crop))
+ case Methodology.JURISDICTIONAL_DIRECT_PROXY:
+ return tuple(
+ sorted(e.name for e in (*tang_cropgrids.Crop, *statistical.Livestock))
+ )
+ case Methodology.STATISTICAL:
+ return tuple(
+ sorted(e.name for e in (*statistical.Crop, *statistical.Livestock))
+ )
+ case _:
+ raise ValueError(methodology)
DEFAULT_CONCURRENCY = 4
diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py
index 71154a0..ae9d68e 100644
--- a/jdluc/datasets/fdap.py
+++ b/jdluc/datasets/fdap.py
@@ -331,8 +331,7 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None:
return base.RasterDataset(
band_names=[f"{year=:d}" for year in YEARS],
- # Averaged rather than taken by mode, so each harmonized pixel keeps the share of it masked
- band_type=base.BandType.INTENSIVE,
+ band_type=base.BandType.CATEGORICAL,
dtype="uint8",
no_data=NO_DATA,
partitioning=tiling.Partitioning.TEN_DEGREE_TILE,
diff --git a/jdluc/datasets/tang_cropgrids.py b/jdluc/datasets/tang_cropgrids.py
index bedb58c..df29c01 100644
--- a/jdluc/datasets/tang_cropgrids.py
+++ b/jdluc/datasets/tang_cropgrids.py
@@ -213,6 +213,7 @@ class Crop(enum.StrEnum):
RESOLUTION = tiling.XY(x=7_200, y=3_600)
+YEAR = 2020
def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None:
diff --git a/jdluc/emit.py b/jdluc/emit.py
index f4af306..6a0c112 100644
--- a/jdluc/emit.py
+++ b/jdluc/emit.py
@@ -34,6 +34,7 @@
from jdluc import geo, harmonize, storage, tiling, utils
from jdluc.datasets import (
descals_oil_palm,
+ fdap,
gnw_global_peatlands,
gnw_harris_agb,
gnw_tcl,
@@ -289,6 +290,8 @@ class DestinationDataset(enum.IntFlag):
"""Which datasets (the can superpose) claim a pixel's land class at the assessment year."""
DESCALS_OIL_PALM = enum.auto()
+ FDAP_COCOA = enum.auto()
+ FDAP_COFFEE = enum.auto()
LIAO_GACED30 = enum.auto()
GPW_GRASSLAND = enum.auto()
@@ -304,6 +307,22 @@ class DestinationDataset(enum.IntFlag):
)
& (planting_year <= ASSESSMENT_YEAR)
),
+ DestinationDataset.FDAP_COCOA: lambda dset: (
+ dset[
+ fdap.COCOA_DATASET.fully_qualified_band_names[
+ fdap.YEARS.index(ASSESSMENT_YEAR)
+ ]
+ ]
+ == 1
+ ),
+ DestinationDataset.FDAP_COFFEE: lambda dset: (
+ dset[
+ fdap.COFFEE_DATASET.fully_qualified_band_names[
+ fdap.YEARS.index(ASSESSMENT_YEAR)
+ ]
+ ]
+ == 1
+ ),
DestinationDataset.LIAO_GACED30: lambda dset: (
dset[
liao_gaced30.DATASET.fully_qualified_band_names[
@@ -324,7 +343,10 @@ class DestinationDataset(enum.IntFlag):
assert set(DestinationDataset) == set(DESTINATION_DATASET_TO_PREDICATE)
CROPLAND_DATASETS = (
- DestinationDataset.DESCALS_OIL_PALM | DestinationDataset.LIAO_GACED30
+ DestinationDataset.DESCALS_OIL_PALM
+ | DestinationDataset.FDAP_COCOA
+ | DestinationDataset.FDAP_COFFEE
+ | DestinationDataset.LIAO_GACED30
)
PASTURE_DATASETS = DestinationDataset.GPW_GRASSLAND
# Each member belongs to one group, and the cropland members come first -- so a cropland bit
@@ -334,6 +356,7 @@ class DestinationDataset(enum.IntFlag):
# How band names spell each destination's land class
CROPLAND = "cropland"
PASTURELAND = "pastureland"
+PEATLAND_OCCUPATION = "peatland-occupation"
def get_last_departure_year(is_source: xarray.DataArray) -> xarray.DataArray:
@@ -675,7 +698,7 @@ def merge_name_units(name: str, units: typing.Hashable | None) -> str:
)
-@storage.cache_to_zarr(version=1)
+@storage.cache_to_zarr(version=2)
def workflow(tile_id: str) -> xarray.Dataset:
logger.info(f"Running the land conversion and emissions worflow for {tile_id=:s}")
dset = geo.exact_merge(
@@ -789,12 +812,12 @@ def workflow(tile_id: str) -> xarray.Dataset:
for (before, after), darray in span_to_vegetation_emissions.items()
}
| {
- f"{CROPLAND:s}-peatland-occupation": cropland_occupation_emissions,
+ f"{CROPLAND:s}-{PEATLAND_OCCUPATION:s}": cropland_occupation_emissions,
# NB: charged to nobody, so it is reported beside the total rather than inside it
"dropped-emissions": dropped_emissions,
"emissions-per-hectare": emissions_per_hectare,
"hectares-per-pixel": get_hectares_per_pixel(darray=emissions_per_hectare),
- f"{PASTURELAND:s}-peatland-occupation": pastureland_occupation_emissions,
+ f"{PASTURELAND:s}-{PEATLAND_OCCUPATION:s}": pastureland_occupation_emissions,
}
)
diff --git a/jdluc/harmonize.py b/jdluc/harmonize.py
index 923307f..6bc2d1a 100644
--- a/jdluc/harmonize.py
+++ b/jdluc/harmonize.py
@@ -168,7 +168,10 @@ def get_resampling_for_band_type(
):
match band_type:
case base.BandType.CATEGORICAL:
- return rasterio.enums.Resampling.mode
+ # Although mode is the natural choice, we choose nearest because
+ # it is an unbiased estimate over area -- mode systematically shrinks
+ # the minority class
+ return rasterio.enums.Resampling.nearest
case base.BandType.EXTENSIVE:
raise NotImplementedError("GDAL doesn't implement sum resampling")
case base.BandType.INTENSIVE:
@@ -386,7 +389,11 @@ class Stack(enum.Enum):
DatasetName.LIAO_GACED30,
DatasetName.SOILGRIDS_OCS,
)
- CROP_SUPPLEMENT = (DatasetName.DESCALS_OIL_PALM,)
+ CROP_SUPPLEMENT = (
+ DatasetName.DESCALS_OIL_PALM,
+ DatasetName.FDAP_COCOA,
+ DatasetName.FDAP_COFFEE,
+ )
def all_unique(it: collections.abc.Iterable[str]) -> bool:
diff --git a/jdluc/jurisdictional_direct.py b/jdluc/jurisdictional_direct.py
index c23aed6..f2d81a3 100644
--- a/jdluc/jurisdictional_direct.py
+++ b/jdluc/jurisdictional_direct.py
@@ -44,11 +44,11 @@ class Crop(enum.Enum):
# Harmonized
gnw_global_peatlands.DATASET.fully_qualified_band_name,
# Emissions
- f"{emit.CROPLAND:s}-peatland-occupation:tco2e-per-ha",
+ f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha",
"dropped-emissions:tco2e-per-ha",
"emissions-per-hectare:tco2e-per-ha",
"hectares-per-pixel:ha",
- f"{emit.PASTURELAND:s}-peatland-occupation:tco2e-per-ha",
+ f"{emit.PASTURELAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha",
*(f"{component!s}-emissions:tco2e-per-ha" for component in emit.EmissionComponent),
# Crop
usda_nass_cdl.DATASET.fully_qualified_band_name,
@@ -137,7 +137,7 @@ def get_commodity_name_to_totals(
}
-@storage.cache_to_parquet(version=0)
+@storage.cache_to_parquet(version=2)
def workflow(
commodity_names: tuple[str, ...],
iso_3166: str,
@@ -217,10 +217,10 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]:
]
== 1,
pastureland_occupation_per_hectare=clipped[
- f"{emit.PASTURELAND:s}-peatland-occupation:tco2e-per-ha"
+ f"{emit.PASTURELAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha"
],
peatland_occupation_per_hectare=clipped[
- f"{emit.CROPLAND:s}-peatland-occupation:tco2e-per-ha"
+ f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha"
],
)
for commodity_name, totals in commodity_name_to_totals.items():
diff --git a/jdluc/jurisdictional_direct_proxy.py b/jdluc/jurisdictional_direct_proxy.py
new file mode 100644
index 0000000..e4c01ff
--- /dev/null
+++ b/jdluc/jurisdictional_direct_proxy.py
@@ -0,0 +1,346 @@
+import collections
+import collections.abc
+import logging
+
+import numpy
+import pandas
+import rasterio.enums
+import xarray
+from rioxarray.exceptions import NoDataInBounds
+
+from jdluc import emit, geo, harmonize, statistical, storage, tiling, utils
+from jdluc.datasets import (
+ DatasetName,
+ faostat,
+ gpw_grassland,
+ gpw_livestock,
+ tang_cropgrids,
+ worldbank_jurisdictions,
+)
+
+logger = logging.getLogger(__name__)
+
+
+assert tang_cropgrids.YEAR in gpw_grassland.YEARS
+assert tang_cropgrids.YEAR in gpw_livestock.YEARS
+
+type Commodity = tang_cropgrids.Crop | statistical.Livestock
+assert not {e.name for e in tang_cropgrids.Crop} & {
+ e.name for e in statistical.Livestock
+}
+
+DATASET_NAMES = (DatasetName.GPW_LIVESTOCK, DatasetName.TANG_CROPGRIDS_CROP_AREA)
+
+
+def get_name_to_band(dset: xarray.Dataset) -> dict[str, xarray.DataArray]:
+ """The 30 m bands this leg divides, each per hectare or a fraction of the pixel.
+
+ Each emissions component is split by the land class its pixel became, and discounted over
+ the spans: the shares that divide it come from one year, so no span needs its own band.
+ """
+ span_to_component_to_emissions = emit.get_span_to_component_to_emissions(dset=dset)
+ return (
+ {
+ f"{land_class:s}-{component!s}-emissions:tco2e-per-ha": emit.get_linear_discounted_total(
+ span_to_value={
+ span: component_to_emissions[component].where(
+ dset["conversion"].isin(conversions), other=0
+ )
+ for span, component_to_emissions in span_to_component_to_emissions.items()
+ }
+ )
+ for land_class, conversions in (
+ (emit.CROPLAND, emit.TO_CROPLAND),
+ (emit.PASTURELAND, emit.TO_PASTURE),
+ )
+ for component in emit.EmissionComponent
+ }
+ | {
+ name: dset[name]
+ for name in (
+ f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha",
+ "dropped-emissions:tco2e-per-ha",
+ f"{emit.PASTURELAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha",
+ )
+ }
+ | {
+ # NB: a mask here and an averaging downscale after, so this arrives as a cell fraction
+ f"pasture:fraction:{tang_cropgrids.YEAR:d}": (
+ dset[
+ gpw_grassland.DATASET.fully_qualified_band_names[
+ gpw_grassland.YEARS.index(tang_cropgrids.YEAR)
+ ]
+ ]
+ == gpw_grassland.Grassland.CULTIVATED
+ ).astype(numpy.float32)
+ }
+ )
+
+
+@storage.cache_to_zarr(version=0)
+def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset:
+ logger.info("Computing emissions on the GLAD grid")
+ luc_and_emissions = geo.exact_merge(
+ harmonize.workflow(
+ dataset_names=harmonize.Stack.LUC_AND_EMISSIONS.value,
+ ignore_missing_tiles=True,
+ skip_ingest=False,
+ tile_id=tile_id,
+ tile_resolution=tiling.TileResolution.GLAD,
+ ),
+ # NB: this should call the harmonize workflow with identical args and hit
+ # the cache from the preceding call
+ emit.workflow(tile_id=tile_id),
+ )
+
+ logger.info("Downsampling emissions from the GLAD grid to the CROPGRIDS grid")
+ grid = harmonize.Grid.from_tile_id_resolution(
+ resolution=tiling.TileResolution.CROPGRIDS, tile_id=tile_id
+ )
+ return xarray.Dataset(
+ {
+ name: geo.downscale_darray(
+ darray=band.rename(name),
+ epsg=grid.epsg,
+ height=grid.resolution.y,
+ # NB: without an equal-area CRS, this presents a small error
+ resampling=rasterio.enums.Resampling.average,
+ transform=grid.transform,
+ width=grid.resolution.x,
+ )
+ for name, band in get_name_to_band(dset=luc_and_emissions).items()
+ }
+ )
+
+
+def get_crop_hectares(
+ crop: tang_cropgrids.Crop, dset: xarray.Dataset
+) -> xarray.DataArray:
+ """`crop`'s physical area in each cell, in hectares; zero over the ocean, which CROPGRIDS
+ leaves missing."""
+ return dset[
+ tang_cropgrids.DATASET.fully_qualified_band_names[
+ list(tang_cropgrids.Crop).index(crop)
+ ]
+ ].fillna(0)
+
+
+def get_crop_to_area_share(
+ crops: tuple[tang_cropgrids.Crop, ...], dset: xarray.Dataset
+) -> dict[tang_cropgrids.Crop, xarray.DataArray]:
+ """Each crop's share of a cell's crop area, over every crop CROPGRIDS maps."""
+ total_area = sum(
+ (get_crop_hectares(crop=crop, dset=dset) for crop in tang_cropgrids.Crop),
+ start=xarray.DataArray(numpy.float32(0)),
+ )
+ total_area = total_area.where(total_area > 0)
+ return {
+ # Share is zero when no crop occupies the cell at all
+ crop: (get_crop_hectares(crop=crop, dset=dset) / total_area).fillna(0)
+ for crop in crops
+ }
+
+
+def get_commodity_name_to_totals(
+ commodities: tuple[Commodity, ...],
+ commodity_to_share: dict[Commodity, xarray.DataArray],
+ dset: xarray.Dataset,
+) -> dict[str, dict[str, float]]:
+ """Each commodity's emissions and area, summed over a jurisdiction's cells.
+
+ A crop takes its crop-area share of everything charged to cropland, and a livestock commodity
+ its grazing share of everything charged to pasture, in every component and in peat occupation
+ alike. Beside them sit `DROPPED` and `UNATTRIBUTED` (see `emit.NonCommodity`); here
+ `UNATTRIBUTED` is what the shares leave of the charge: cropland's where CROPGRIDS maps no crop,
+ and the shares of commodities not asked for.
+ """
+ hectares = dset["hectares-per-cell:ha"]
+ land_class_to_component_to_emissions = {
+ land_class: {
+ component: dset[f"{land_class:s}-{component!s}-emissions:tco2e-per-ha"]
+ * hectares
+ for component in emit.EmissionComponent
+ }
+ for land_class in (emit.CROPLAND, emit.PASTURELAND)
+ }
+ land_class_to_occupation = {
+ land_class: dset[f"{land_class:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha"]
+ * hectares
+ for land_class in (emit.CROPLAND, emit.PASTURELAND)
+ }
+
+ commodity_to_totals: dict[
+ Commodity | emit.NonCommodity, dict[str, xarray.DataArray]
+ ] = collections.defaultdict(dict)
+ for commodity in commodities:
+ totals = commodity_to_totals[commodity]
+ share = commodity_to_share[commodity]
+ if isinstance(commodity, statistical.Livestock):
+ land_class = emit.PASTURELAND
+ totals["commodity_hectares"] = (
+ dset[f"pasture:fraction:{tang_cropgrids.YEAR:d}"] * hectares * share
+ )
+ else:
+ land_class = emit.CROPLAND
+ totals["commodity_hectares"] = get_crop_hectares(crop=commodity, dset=dset)
+ component_to_emissions = {
+ component: emissions * share
+ for component, emissions in land_class_to_component_to_emissions[
+ land_class
+ ].items()
+ }
+ occupation = land_class_to_occupation[land_class] * share
+ totals |= {
+ component.column: component_to_emissions[component]
+ for component in emit.EmissionComponent
+ }
+ # NB: occupation is a flat annual rate over the area drained, so dividing it back out
+ # recovers that area on `emit`'s 30 m intersection
+ totals["peatland_commodity_hectares"] = (
+ occupation / emit.PEATLAND_EMISSIONS_ANNUAL_TCO2E_PER_HA
+ )
+ totals["peatland_occupation_emissions_mt"] = occupation
+ totals["emissions_mt"] = sum(component_to_emissions.values(), start=occupation)
+
+ commodity_to_totals[emit.NonCommodity.DROPPED]["emissions_mt"] = (
+ dset["dropped-emissions:tco2e-per-ha"] * hectares
+ )
+
+ # The whole charge each component carries, as a commodity holding every share would take it
+ component_to_charged = {
+ component: land_class_to_component_to_emissions[emit.CROPLAND][component]
+ + land_class_to_component_to_emissions[emit.PASTURELAND][component]
+ for component in emit.EmissionComponent
+ }
+ unattributed = commodity_to_totals[emit.NonCommodity.UNATTRIBUTED]
+ unattributed |= {
+ component.column: component_to_charged[component]
+ - sum(
+ commodity_to_totals[commodity][component.column]
+ for commodity in commodities
+ )
+ for component in emit.EmissionComponent
+ }
+ unattributed["peatland_occupation_emissions_mt"] = (
+ land_class_to_occupation[emit.CROPLAND]
+ + land_class_to_occupation[emit.PASTURELAND]
+ - sum(
+ commodity_to_totals[commodity]["peatland_occupation_emissions_mt"]
+ for commodity in commodities
+ )
+ )
+ unattributed["emissions_mt"] = sum(
+ (unattributed[component.column] for component in emit.EmissionComponent),
+ start=unattributed["peatland_occupation_emissions_mt"],
+ )
+
+ return utils.get_sum_totals(enum_to_name_to_darray=commodity_to_totals)
+
+
+SCHEMA = {
+ "admin_id": str,
+ "admin_level": str,
+ "commodity_hectares": float,
+ "commodity_name": str,
+ "emissions_mt": float,
+ "forest_emissions_mt": float,
+ "grassland_emissions_mt": float,
+ "jurisdiction_name": str,
+ "peatland_commodity_hectares": float,
+ "peatland_conversion_emissions_mt": float,
+ "peatland_occupation_emissions_mt": float,
+}
+
+
+@storage.cache_to_parquet(version=0)
+def workflow(
+ commodity_names: tuple[str, ...],
+ iso_3166: str,
+ tile_id: str,
+) -> pandas.DataFrame:
+ commodities: tuple[Commodity, ...] = tuple(
+ tang_cropgrids.Crop[name]
+ if name in tang_cropgrids.Crop.__members__
+ else statistical.Livestock[name]
+ for name in commodity_names
+ )
+ crops = tuple(
+ commodity
+ for commodity in commodities
+ if isinstance(commodity, tang_cropgrids.Crop)
+ )
+ logger.info(f"Computing emissions for {commodities=:} and {tile_id=:s}")
+ merged = geo.exact_merge(
+ get_downscaled_luc_emissions(tile_id=tile_id),
+ harmonize.workflow(
+ dataset_names=DATASET_NAMES,
+ ignore_missing_tiles=True,
+ skip_ingest=False,
+ tile_id=tile_id,
+ tile_resolution=tiling.TileResolution.CROPGRIDS,
+ ),
+ )
+ # Each grazer's livestock units per head, as FAOSTAT counts them for this country
+ merged = merged.assign(
+ statistical.get_band_name_to_livestock_units(
+ dset=merged,
+ species_to_units_per_head={
+ gpw_livestock.Species[species.name]: units_per_head
+ for species, units_per_head in faostat.get_species_to_livestock_units_per_head(
+ iso_3166=iso_3166,
+ livestock_units=faostat.load(
+ dataset=faostat.LIVESTOCK_PATTERNS_DATASET
+ ),
+ stocks=faostat.load(dataset=faostat.LIVESTOCK_DATASET),
+ ).items()
+ },
+ )
+ )
+ # NB: on the whole tile, since a clip one cell wide has no spacing to measure a cell by
+ merged["hectares-per-cell:ha"] = emit.get_hectares_per_pixel(
+ darray=merged["dropped-emissions:tco2e-per-ha"]
+ )
+ # NB: small enough to hold in memory, which speeds up the many reads for each jurisdiction
+ merged = merged.load()
+
+ def it() -> collections.abc.Iterator[dict[str, float | str]]:
+ for (
+ jurisdiction
+ ) in worldbank_jurisdictions.iter_jurisdiction_for_iso_3166_tile_id(
+ admin_level=worldbank_jurisdictions.AdminLevel.PROVINCIAL,
+ iso_3166=iso_3166,
+ tile_id=tile_id,
+ ):
+ logger.info(f"Clipping to provincial geometry from {jurisdiction.id=:s}")
+ try:
+ clipped = geo.clip_dset(dset=merged, geometry=jurisdiction.geometry)
+ except NoDataInBounds as exc:
+ logger.warning(repr(exc))
+ else:
+ logger.info(
+ f"Populating emissions for {jurisdiction.id=!s} over "
+ f"{len(commodities):d} commodities"
+ )
+ commodity_name_to_totals = get_commodity_name_to_totals(
+ commodities=commodities,
+ commodity_to_share=get_crop_to_area_share(crops=crops, dset=clipped)
+ | statistical.get_livestock_to_grazing_share(
+ dset=clipped, year=tang_cropgrids.YEAR
+ ),
+ dset=clipped,
+ )
+ for commodity_name, totals in commodity_name_to_totals.items():
+ yield totals | {
+ "admin_id": jurisdiction.id,
+ "admin_level": jurisdiction.level,
+ "commodity_name": commodity_name,
+ "jurisdiction_name": jurisdiction.name,
+ }
+
+ if data := list(it()):
+ assert set(data[0]) == set(SCHEMA)
+ return (
+ pandas.DataFrame.from_records(columns=list(SCHEMA), data=data)
+ .astype(SCHEMA)
+ .set_index(["admin_level", "admin_id", "commodity_name"])
+ )
diff --git a/jdluc/statistical.py b/jdluc/statistical.py
index 140bde0..514333e 100644
--- a/jdluc/statistical.py
+++ b/jdluc/statistical.py
@@ -100,9 +100,9 @@ class Livestock(enum.StrEnum):
EMIT_VARIABLE_NAMES = [
- f"{emit.CROPLAND:s}-peatland-occupation:tco2e-per-ha",
+ f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha",
"dropped-emissions:tco2e-per-ha",
- f"{emit.PASTURELAND:s}-peatland-occupation:tco2e-per-ha",
+ f"{emit.PASTURELAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha",
]
# Components divided by the pooled expansion share. The peat pulse is not one of them: it is stored
@@ -113,7 +113,7 @@ class Livestock(enum.StrEnum):
}
-@storage.cache_to_zarr(version=0)
+@storage.cache_to_zarr(version=2)
def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset:
logger.info("Computing emissions on the GLAD grid")
luc_and_emissions = geo.exact_merge(
@@ -144,9 +144,7 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset:
(emit.CROPLAND, emit.TO_CROPLAND),
(emit.PASTURELAND, emit.TO_PASTURE),
):
- name = (
- f"{land_class:s}-peatland-conversion:tco2e-per-ha:{before:d}-{after:d}"
- )
+ name = f"{land_class:s}-{emit.EmissionComponent.PEATLAND_CONVERSION!s}:tco2e-per-ha:{before:d}-{after:d}"
luc_and_emissions[name] = component_to_darray[
emit.EmissionComponent.PEATLAND_CONVERSION
].where(luc_and_emissions["conversion"].isin(conversions), other=0)
@@ -405,12 +403,13 @@ def get_commodity_name_to_totals(
pasture expanded, and the cropland peat pulse where no crop expanded.
"""
cropland_occupation_per_hectare = dset[
- f"{emit.CROPLAND:s}-peatland-occupation:tco2e-per-ha"
+ f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha"
]
hectares = emit.get_hectares_per_pixel(darray=cropland_occupation_per_hectare)
cropland_occupation = cropland_occupation_per_hectare * hectares
pastureland_occupation = (
- dset[f"{emit.PASTURELAND:s}-peatland-occupation:tco2e-per-ha"] * hectares
+ dset[f"{emit.PASTURELAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha"]
+ * hectares
)
commodity_to_totals: dict[
@@ -439,7 +438,7 @@ def get_commodity_name_to_totals(
emit.EmissionComponent.PEATLAND_CONVERSION: emit.get_linear_discounted_total(
span_to_value={
(before, after): dset[
- f"{land_class:s}-peatland-conversion:tco2e-per-ha:{before:d}-{after:d}"
+ f"{land_class:s}-{emit.EmissionComponent.PEATLAND_CONVERSION!s}:tco2e-per-ha:{before:d}-{after:d}"
]
* hectares
* span_to_commodity_to_peatland_conversion_share[
@@ -532,7 +531,7 @@ def get_commodity_name_to_totals(
span_to_value={
(before, after): sum(
dset[
- f"{land_class:s}-peatland-conversion:tco2e-per-ha:{before:d}-{after:d}"
+ f"{land_class:s}-{emit.EmissionComponent.PEATLAND_CONVERSION!s}:tco2e-per-ha:{before:d}-{after:d}"
]
for land_class in (emit.CROPLAND, emit.PASTURELAND)
)
@@ -621,7 +620,7 @@ def get_crop_to_area_share(
}
-@storage.cache_to_parquet(version=0)
+@storage.cache_to_parquet(version=2)
def workflow(
commodity_names: tuple[str, ...],
iso_3166: str,
diff --git a/jdluc/tiling.py b/jdluc/tiling.py
index 8902a80..43edbbd 100644
--- a/jdluc/tiling.py
+++ b/jdluc/tiling.py
@@ -27,6 +27,7 @@ class TileResolution(XY, enum.Enum):
# The dataset is provided in ESRI:54052, so no EPSG:4326 grid fits exactly; this is
# the coarsest one which is finer than the source grid and tiles into MAPSPAM
GPW_LIVESTOCK = 1_200, 1_200
+ CROPGRIDS = 200, 200
MAPSPAM = 120, 120
@property
diff --git a/jdluc/trace.py b/jdluc/trace.py
index 335fcab..56a8747 100644
--- a/jdluc/trace.py
+++ b/jdluc/trace.py
@@ -116,7 +116,7 @@ def safe_div(numer: float, denom: float) -> float:
} # type: ignore
-@storage.cache_to_parquet(ignored_args=["concurrency"], version=0)
+@storage.cache_to_parquet(ignored_args=["concurrency"], version=2)
def workflow(
concurrency: int,
commodity_names: tuple[str, ...],
diff --git a/validation/__tests__/schema_test.py b/validation/__tests__/schema_test.py
index 22b5786..b972424 100644
--- a/validation/__tests__/schema_test.py
+++ b/validation/__tests__/schema_test.py
@@ -38,9 +38,12 @@ def test_methodology_names_match_the_enum() -> None:
def test_every_methodology_is_named_here() -> None:
- """A third leg would otherwise be silently absent from every comparison."""
+ """A new leg would otherwise be silently absent from every comparison."""
assert {schema.STATISTICAL, schema.JURISDICTIONAL_DIRECT} == {
- methodology.name for methodology in attribute.Methodology
+ methodology.name
+ for methodology in attribute.Methodology
+ # Reports no production, so it has no emissions factor to compare
+ if methodology != attribute.Methodology.JURISDICTIONAL_DIRECT_PROXY
}