feat: Ingest the FDaP 2026a dataset for cocoa coffee

docs/architecture.md

@@ -89,7 +89,7 @@ Harmonization is **one GDAL VRT per source band per tile**. Each VRT declares th
 There is no reprojection step and no mosaic step:
 
 - **No reprojection.** Every ingested COG is already `EPSG:4326` — `geo.validate_geotiff` asserts it at ingest — so the VRT only ever rescales and offsets within a single CRS. An earlier design added a second VRT pass (`rasterio.vrt.WarpedVRT`) to reproject; it was removed once the 4326-at-ingest invariant was enforced. (A `WarpedVRT` does survive in `geo.py`, where the statistical leg downsamples to the MapSPAM grid; it too is a same-CRS rescale.)
-- **No mosaic.** The unit of work is one tile, so there is never more than one source per band to combine. Two source layouts are handled: ten-degree-tiled datasets (GLAD, GNW peatlands, Harris AGB, Huang BGB, SoilGrids, CDL, TCL, GACED30, GPW grassland, GPW livestock, Descals oil palm) map a whole source tile onto the whole destination tile, while whole-world datasets (IPCC climate zones, MapSPAM) use `SrcRect` to window the tile's footprint out of the global raster. A publisher whose own tiling or projection does not match the 10° graticule is mosaicked or warped onto it during ingest instead — GACED30 ships three-degree tiles, and its dataset module resamples the sixteen covering a tile onto that tile; Descals oil palm ships hundred-kilometre cells on a lattice that divides neither the graticule nor its own pixel grid evenly, and merges up to 66 of them per tile; GPW livestock ships one global raster per species and year on Interrupted Goode Homolosine, and its dataset module warps each tile's footprint onto a 30-arc-second grid — so harmonization still sees exactly one source per band.
+- **No mosaic.** The unit of work is one tile, so there is never more than one source per band to combine. Two source layouts are handled: ten-degree-tiled datasets (GLAD, GNW peatlands, Harris AGB, Huang BGB, SoilGrids, CDL, TCL, GACED30, GPW grassland, GPW livestock, Descals oil palm) map a whole source tile onto the whole destination tile, while whole-world datasets (IPCC climate zones, MapSPAM) use `SrcRect` to window the tile's footprint out of the global raster. A publisher whose own tiling or projection does not match the 10° graticule is mosaicked or warped onto it during ingest instead — GACED30 ships three-degree tiles, and its dataset module resamples the sixteen covering a tile onto that tile; Descals oil palm ships hundred-kilometre cells on a lattice that divides neither the graticule nor its own pixel grid evenly, and merges up to 66 of them per tile; the FDaP masks ship as one-degree files on their own 10 m grid, so their dataset module pieces up to 100 of them into a tile at that resolution, 111,319 or 111,320 pixels a side, without resampling; GPW livestock ships one global raster per species and year on Interrupted Goode Homolosine, and its dataset module warps each tile's footprint onto a 30-arc-second grid — so harmonization still sees exactly one source per band.
 
 A tile whose source was never published is tolerated when the caller opts in: the VRT is emitted with a band 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.
 

docs/further_research.md

@@ -29,7 +29,19 @@ Sampled over 1°×1° windows at 30 m on GLAD's land-cover map, not on GACED30,
 
 **Potential impact:** Potentially the largest single gap for the global extension. The tropical acacia and oil palm emission factors (78.6 and 45.0 t CO₂e ha⁻¹ yr⁻¹) motivate the front-loaded shape of the entire peatland model, yet the pipeline can attribute close to zero peat emissions to those crops. Note that the usual reassurance — that per-crop factors are protected because the numerator only counts pixels that are cropland in 2020 — does not hold for woody perennials, where that same restriction is the mechanism of the undercount. The question is not whether the forest layer *sees* the clearing, which tree-cover loss does; it is that the clearing’s destination is never treated as converted.
 
-**Potential improvement path:** A dedicated perennial-plantation layer is the route in, and the repo ingests one: Descals et al. (2024) global oil palm, whose 30 m pixels carry the planting year as well as the extent, so the peat pulse can be dated directly rather than inferred from a five-year GLAD span. Two properties bound what it can settle. Its planting year dates the plantation standing in 2021, so a replant reads as the replant year rather than as the year the land was first cleared — it says when the current crop was established, not when conversion happened, and pairing it with TCL or a GLAD span is what separates the two. And its accuracy is asymmetric: producer's and user's accuracy are 91.0 ± 2.5% and 91.8 ± 1.2% on industrial plantations against 71.4 ± 0.7% and 72.4 ± 1.8% on smallholders, so a smallholder frontier is the weak end of the layer as well as of GLAD. That has now been done for oil palm: the layer is admitted to the emissions numerator alongside the general cropland mask, and mapped plantation on peat is drained for the occupation term. It is read as extent at the assessment year and dates nothing, precisely because of the replant property above. One consequence is accepted rather than avoided: where the planting year precedes the tree-cover loss, the pixel still resolves to oil palm, so plantation-cycle canopy loss inside an existing estate is charged as conversion — 5.3-8.6% of the conversions the layer admits on Riau. Vetoing those would be dating with the layer, which it is not trusted to do. It also covers only oil palm, which is 352 Mt of the 1,071 Mt in the gap above; the remaining 719 Mt sits on crops with no dedicated extent layer here, cocoa, coconut and rubber among them, and coffee and tea with no open global map at all. That list is woody crops only, which is narrower than "perennial": sugarcane, banana and plantain are perennial but herbaceous, so GLAD's definition admits them. The classifier does not honor that admission where the canopy is dense — the screen above puts sugarcane at 64.2% coverage, banana at 33.9% and plantain at 28.6% — which brings those 71 (country, crop) pairs into scope alongside the woody crops. A general woody-cropland upgrade is also available: GLC_FCS30D (Zhang et al., 2024) carries a separable tree/shrub orchard cropland class, annual for 2000–2022, which would admit woody cropland to the numerator under a different carbon treatment than annual cropland — something GACED30's binary cropland class cannot support. Establish its per-class accuracy for that class first, since the published global overall accuracy of 80.88% spans 35 classes and says nothing about it, along with whether the 25.7% palm recall holds outside Riau and whether the 24% of that window it calls shrubland is immature palm. Adopting it also forces a regrid: its 30 m lattice does not nest into GLAD's 0.00025°, on which TCL, SDPT v2 and Global Pasture Watch are all native.
+**Potential improvement path:** A dedicated perennial-plantation layer is the route in, and the repo ingests one: Descals et al. (2024) global oil palm, whose 30 m pixels carry the planting year as well as the extent, so the peat pulse can be dated directly rather than inferred from a five-year GLAD span. Two properties bound what it can settle. Its planting year dates the plantation standing in 2021, so a replant reads as the replant year rather than as the year the land was first cleared — it says when the current crop was established, not when conversion happened, and pairing it with TCL or a GLAD span is what separates the two. And its accuracy is asymmetric: producer's and user's accuracy are 91.0 ± 2.5% and 91.8 ± 1.2% on industrial plantations against 71.4 ± 0.7% and 72.4 ± 1.8% on smallholders, so a smallholder frontier is the weak end of the layer as well as of GLAD. That has now been done for oil palm: the layer is admitted to the emissions numerator alongside the general cropland mask, and mapped plantation on peat is drained for the occupation term. It is read as extent at the assessment year and dates nothing, precisely because of the replant property above. One consequence is accepted rather than avoided: where the planting year precedes the tree-cover loss, the pixel still resolves to oil palm, so plantation-cycle canopy loss inside an existing estate is charged as conversion — 5.3-8.6% of the conversions the layer admits on Riau. Vetoing those would be dating with the layer, which it is not trusted to do. It also covers only oil palm, which is 352 Mt of the 1,071 Mt in the gap above; the remaining 719 Mt sits on crops with no dedicated extent layer the emissions core reads, cocoa, coffee, coconut and rubber among them. For cocoa and coffee the repo ingests one, the Forest Data Partnership's 10 m masks for 2018–2024, which no stage reads yet; tea has no open global map at all. That list is woody crops only, which is narrower than "perennial": sugarcane, banana and plantain are perennial but herbaceous, so GLAD's definition admits them. The classifier does not honor that admission where the canopy is dense — the screen above puts sugarcane at 64.2% coverage, banana at 33.9% and plantain at 28.6% — which brings those 71 (country, crop) pairs into scope alongside the woody crops. A general woody-cropland upgrade is also available: GLC_FCS30D (Zhang et al., 2024) carries a separable tree/shrub orchard cropland class, annual for 2000–2022, which would admit woody cropland to the numerator under a different carbon treatment than annual cropland — something GACED30's binary cropland class cannot support. Establish its per-class accuracy for that class first, since the published global overall accuracy of 80.88% spans 35 classes and says nothing about it, along with whether the 25.7% palm recall holds outside Riau and whether the 24% of that window it calls shrubland is immature palm. Adopting it also forces a regrid: its 30 m lattice does not nest into GLAD's 0.00025°, on which TCL, SDPT v2 and Global Pasture Watch are all native.
+
+## 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.
+
+The masks threshold FDaP's per-pixel probabilities per country and year (`jdluc/datasets/fdap.py`). Where a country has a fitted threshold, 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. Other producers use the model's precision/recall crossover, 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 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.
 
 ## Shifting cultivation resolves to no destination
 

docs/methodology.md

@@ -98,10 +98,11 @@ Every input is ingested from its upstream publisher into cloud storage as tiled
 | **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  |
 
 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 and GLAD GLCLUC 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. GPW livestock is Global Pasture Watch's FAOSTAT-adjusted annual headcount of buffalo, cattle, goats, horses and sheep.
+FAOSTAT's crop table, GLAD GLCLUC and the FDaP masks 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.
 
 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.
 

jdluc/__tests__/datasets/base_test.py

@@ -0,0 +1,15 @@
+import pytest
+
+from jdluc.datasets import NAME_TO_CLS, base
+
+
+@pytest.mark.parametrize(
+    "dataset",
+    [
+        pytest.param(dataset, id=name)
+        for name, dataset in NAME_TO_CLS.items()
+        if isinstance(dataset, base.TabularDataset)
+    ],
+)
+def test_a_tabular_datasets_are_hashable(dataset: base.TabularDataset) -> None:
+    assert hash(dataset)

jdluc/__tests__/datasets/faostat_test.py

@@ -0,0 +1,85 @@
+import pandas
+import pytest
+
+from jdluc.datasets import faostat
+from jdluc.datasets.worldbank_jurisdictions import AdminLevel
+
+
+def get_faostat_frame(
+    column_name: str, key_to_value: dict[tuple[str, str, int], float]
+) -> pandas.DataFrame:
+    """A table shaped like `faostat.load`'s, one column over (country, species, year)."""
+    return pandas.DataFrame.from_records(
+        data=[
+            {
+                "admin_level": AdminLevel.NATIONAL.name,
+                "admin_id": iso_3166,
+                "jurisdiction_name": iso_3166,
+                "commodity_name": commodity_name,
+                "year": year,
+                column_name: value,
+            }
+            for (iso_3166, commodity_name, year), value in key_to_value.items()
+        ]
+    ).set_index(faostat.IDX_COLUMN_NAMES)
+
+
+NO_LIVESTOCK_UNITS = dict.fromkeys(faostat.Species, 0.0)
+CATTLE = faostat.Species.CATTLE
+
+
+@pytest.mark.parametrize(
+    ("livestock_units", "stocks", "expected"),
+    (
+        pytest.param(
+            {
+                ("BRA", "CATTLE", 2000): 70.0,
+                ("BRA", "CATTLE", 2005): 70.0,
+                ("BRA", "CATTLE", 2010): 95.0,
+            },
+            {
+                ("BRA", "CATTLE", 2000): 100.0,
+                ("BRA", "CATTLE", 2005): 100.0,
+                ("BRA", "CATTLE", 2010): 100.0,
+            },
+            NO_LIVESTOCK_UNITS | {CATTLE: 0.70},
+            id="a-year-whose-heads-were-since-revised-does-not-move-the-median",
+        ),
+        pytest.param(
+            # Buffalo as FAOSTAT gives South America's, at zero, and sheep with no figure at all:
+            # neither falls back to another coefficient
+            {("BRA", "BUFFALO", 2020): 0.0, ("BRA", "CATTLE", 2020): 70.0},
+            {
+                ("BRA", "BUFFALO", 2020): 100.0,
+                ("BRA", "CATTLE", 2020): 100.0,
+                ("BRA", "SHEEP", 2020): 100.0,
+            },
+            NO_LIVESTOCK_UNITS | {CATTLE: 0.70},
+            id="a-grazer-given-zero-or-no-livestock-units-counts-none",
+        ),
+        pytest.param(
+            {("BRA", "CATTLE", 2000): 10.0, ("BRA", "CATTLE", 2020): 70.0},
+            {("BRA", "CATTLE", 2000): 0.0, ("BRA", "CATTLE", 2020): 100.0},
+            NO_LIVESTOCK_UNITS | {CATTLE: 0.70},
+            id="a-year-without-heads-is-left-out-rather-than-divided-by",
+        ),
+        pytest.param(
+            {("ARG", "CATTLE", 2020): 50.0, ("BRA", "CATTLE", 2020): 70.0},
+            {("ARG", "CATTLE", 2020): 100.0, ("BRA", "CATTLE", 2020): 100.0},
+            NO_LIVESTOCK_UNITS | {CATTLE: 0.70},
+            id="only-the-countrys-own-rows-count",
+        ),
+    ),
+)
+def test_get_species_to_livestock_units_per_head(
+    livestock_units: dict[tuple[str, str, int], float],
+    stocks: dict[tuple[str, str, int], float],
+    expected: dict[faostat.Species, float],
+) -> None:
+    assert faostat.get_species_to_livestock_units_per_head(
+        iso_3166="BRA",
+        livestock_units=get_faostat_frame(
+            column_name="stocks_livestock_units", key_to_value=livestock_units
+        ),
+        stocks=get_faostat_frame(column_name="stocks_head", key_to_value=stocks),
+    ) == pytest.approx(expected)