feat: implement JD proxy's with CROPGRIDS area-share allocation

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).

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<br/><i>sources</i>"] --> harmonize["harmonize<br/><i>common grid</i>"]
   harmonize --> emit["emit<br/><i>per-pixel LUC</i>"]
-  emit --> attribute["attribute<br/><i>per-crop rollup;<br/>JD or statistical</i>"]
+  emit --> attribute["attribute<br/><i>per-crop rollup;<br/>JD, JD proxy or statistical</i>"]
   attribute --> trace["trace<br/><i>emissions factors</i>"]
 ```
 
 - **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
 

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.