From 2a0bf5c9fc6df6cc9f34241705bcf9b0f806e206 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Wed, 28 Oct 2026 14:09:06 -0700 Subject: [PATCH 01/21] Expand FAOSTAT to include coffee Since FAOSTAT groups both species of coffee into one item, track that --- jdluc/datasets/faostat.py | 27 ++++++++++++++++++++------- validation/prepare.py | 4 ++++ 2 files changed, 24 insertions(+), 7 deletions(-) diff --git a/jdluc/datasets/faostat.py b/jdluc/datasets/faostat.py index 7e2390b..2f6b44b 100644 --- a/jdluc/datasets/faostat.py +++ b/jdluc/datasets/faostat.py @@ -22,9 +22,10 @@ Two bulk archives, read into three tables: -- `CROP_DATASET`: area harvested and production of each MapSPAM crop FAOSTAT names one-for-one. - Yield is not carried, since a MapSPAM group's yield is not its members'; area and production are - additive, so a consumer sums the items it wants and divides +- `CROP_DATASET`: area harvested and production of each MapSPAM crop FAOSTAT names one-for-one, + plus coffee as one item, COFF, which MapSPAM splits into arabica and robusta. Yield is not + carried, since a MapSPAM group's yield is not its members'; area and production are additive, so + a consumer sums the items it wants and divides - `LIVESTOCK_DATASET`: stocks (live animals) and meat production (carcass weight, with the bone) of five grazers - `LIVESTOCK_PATTERNS_DATASET`: stocks in livestock units of the same five grazers, from Livestock @@ -152,6 +153,9 @@ class Crop(enum.IntEnum): ifpri_mapspam.CANONICAL_CROP_CLS.PMIL.name, # both are 79, Millet ifpri_mapspam.CANONICAL_CROP_CLS.SMIL.name, } +COMBINED_CROP_NAME_TO_ITEM_CODE = { + ifpri_mapspam.Crop2000.COFF.name: 656, +} # MapSPAM groups that aggregate many FAOSTAT items, by a member list that is MapSPAM's. SPAM_GROUP_CROP_NAMES = { ifpri_mapspam.CANONICAL_CROP_CLS.OCER.name, @@ -170,7 +174,13 @@ class Crop(enum.IntEnum): } assert not {crop.name for crop in Crop} & (SPLIT_CROP_NAMES | SPAM_GROUP_CROP_NAMES) assert set(CROP_TO_ITEM_CODE) == set(Crop) -assert len(set(CROP_TO_ITEM_CODE.values())) == len(CROP_TO_ITEM_CODE) +assert all( + ifpri_mapspam.GROUP_TO_CONSTITUENT_NAMES[combined_crop_name] <= SPLIT_CROP_NAMES + for combined_crop_name in COMBINED_CROP_NAME_TO_ITEM_CODE +) +assert len( + set(CROP_TO_ITEM_CODE.values()) | set(COMBINED_CROP_NAME_TO_ITEM_CODE.values()) +) == len(CROP_TO_ITEM_CODE) + len(COMBINED_CROP_NAME_TO_ITEM_CODE) class Species(enum.StrEnum): @@ -225,7 +235,8 @@ class ElementColumn: column_name="area_hectares", commodity_name_to_item_code={ crop.name: item_code for crop, item_code in CROP_TO_ITEM_CODE.items() - }, + } + | COMBINED_CROP_NAME_TO_ITEM_CODE, required=True, scale=1, ), @@ -233,7 +244,8 @@ class ElementColumn: column_name="production_kg", commodity_name_to_item_code={ crop.name: item_code for crop, item_code in CROP_TO_ITEM_CODE.items() - }, + } + | COMBINED_CROP_NAME_TO_ITEM_CODE, required=True, scale=KG_PER_TONNE, ), @@ -409,7 +421,8 @@ def _get_records_for_tile( product_name="production-crops", source_name="faostat", # v1: crop_name -> commodity_name - version="v1", + # v2: coffee as COFF + version="v2", ) LIVESTOCK_DATASET = base.TabularDataset( get_records_for_tile_id=functools.partial( diff --git a/validation/prepare.py b/validation/prepare.py index 4ce9a0a..708b5db 100644 --- a/validation/prepare.py +++ b/validation/prepare.py @@ -1509,6 +1509,10 @@ def read_faostat_production() -> pandas.DataFrame: frame = pandas.read_parquet(FAOSTAT_CACHE) if frame.index.names != [None]: frame = frame.reset_index() + # WRI's COFF is arabica only + frame = frame[ + ~frame["commodity_name"].isin(faostat.COMBINED_CROP_NAME_TO_ITEM_CODE) + ] return frame.rename(columns={"admin_id": "iso_3166", "commodity_name": "crop_name"}) From a6b01d8e59b967d7e236e49bee288e5c93021975 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Wed, 28 Oct 2026 15:07:13 -0700 Subject: [PATCH 02/21] Split dataset tests into their own modules --- jdluc/__tests__/datasets/base_test.py | 15 + jdluc/__tests__/datasets/faostat_test.py | 85 ++++ jdluc/__tests__/datasets/glad_glcluc_test.py | 5 + .../__tests__/datasets/ifpri_mapspam_test.py | 138 ++++++ .../datasets/worldbank_jurisdictions_test.py | 206 ++++++++ jdluc/__tests__/datasets_test.py | 440 ------------------ jdluc/datasets/ifpri_mapspam.py | 6 +- 7 files changed, 452 insertions(+), 443 deletions(-) create mode 100644 jdluc/__tests__/datasets/base_test.py create mode 100644 jdluc/__tests__/datasets/faostat_test.py create mode 100644 jdluc/__tests__/datasets/glad_glcluc_test.py create mode 100644 jdluc/__tests__/datasets/ifpri_mapspam_test.py create mode 100644 jdluc/__tests__/datasets/worldbank_jurisdictions_test.py delete mode 100644 jdluc/__tests__/datasets_test.py diff --git a/jdluc/__tests__/datasets/base_test.py b/jdluc/__tests__/datasets/base_test.py new file mode 100644 index 0000000..54a26b1 --- /dev/null +++ b/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) diff --git a/jdluc/__tests__/datasets/faostat_test.py b/jdluc/__tests__/datasets/faostat_test.py new file mode 100644 index 0000000..df7689b --- /dev/null +++ b/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) diff --git a/jdluc/__tests__/datasets/glad_glcluc_test.py b/jdluc/__tests__/datasets/glad_glcluc_test.py new file mode 100644 index 0000000..e3bb928 --- /dev/null +++ b/jdluc/__tests__/datasets/glad_glcluc_test.py @@ -0,0 +1,5 @@ +from jdluc.datasets.glad_glcluc import flatten_ranges + + +def test_flatten_ranges() -> None: + assert flatten_ranges(range(5), range(5, 10), range(10, 15)) == list(range(15)) diff --git a/jdluc/__tests__/datasets/ifpri_mapspam_test.py b/jdluc/__tests__/datasets/ifpri_mapspam_test.py new file mode 100644 index 0000000..4a3cddf --- /dev/null +++ b/jdluc/__tests__/datasets/ifpri_mapspam_test.py @@ -0,0 +1,138 @@ +import itertools + +import pytest + +from jdluc.datasets import ifpri_mapspam + + +@pytest.mark.parametrize("year", sorted(ifpri_mapspam.YEARS)) +def test_get_reported_crop_name_always_names_a_band_that_year_has(year: int) -> None: + # The load-bearing property: whatever it returns must be readable from that snapshot. A + # rename or a new group that broke this would otherwise surface as a KeyError deep in a + # pipeline run. + names = {e.name for e in ifpri_mapspam.YEAR_TO_CROP_CLS[year]} + for crop_name in sorted(ifpri_mapspam.RECOVERABLE_CROP_NAMES): + assert ( + ifpri_mapspam.get_reported_crop_name( + canonical_crop_name=crop_name, year=year + ) + in names + ) + + +@pytest.mark.parametrize("year", sorted(ifpri_mapspam.YEARS)) +def test_reported_names_collapse_only_in_the_decompose_year(year: int) -> None: + # 19 constituents share 6 groups in 2000, so 13 names disappear; every other year names all + # 32 separately. This is the whole reason the decomposition exists, asserted directly. + reported = [ + ifpri_mapspam.get_reported_crop_name(canonical_crop_name=name, year=year) + for name in sorted(ifpri_mapspam.RECOVERABLE_CROP_NAMES) + ] + collapsed = len(reported) - len(set(reported)) + assert collapsed == (13 if year == ifpri_mapspam.DECOMPOSITION_YEAR else 0) + + +@pytest.mark.parametrize( + "year", sorted(set(ifpri_mapspam.YEARS) - {ifpri_mapspam.DECOMPOSITION_YEAR}) +) +def test_get_reported_crop_name_defers_to_the_rename_outside_2000( + year: int, +) -> None: + for crop_name in sorted(ifpri_mapspam.RECOVERABLE_CROP_NAMES): + assert ifpri_mapspam.get_reported_crop_name( + canonical_crop_name=crop_name, year=year + ) == ifpri_mapspam.get_renamed_crop_name( + canonical_crop_name=crop_name, year=year + ) + + +def test_is_reported_as_group_is_exactly_the_constituents_in_2000() -> None: + grouped = { + (crop_name, year) + for year in ifpri_mapspam.YEARS + for crop_name in ifpri_mapspam.RECOVERABLE_CROP_NAMES + if ifpri_mapspam.is_reported_as_group(canonical_crop_name=crop_name, year=year) + } + assert grouped == { + (crop_name, ifpri_mapspam.DECOMPOSITION_YEAR) + for crop_name in ifpri_mapspam.CONSTITUENT_TO_GROUP_NAME + } + + +def test_comparable_spans_agree_with_the_vocabularies() -> None: + # Not a restatement of the constant: it checks the constant against how much of the + # unrecoverable vocabulary the two snapshots actually share. Listing a pair as comparable + # when their names barely overlap is what would make name-matching measure the taxonomy + # instead of the land. + for before, after in itertools.combinations(sorted(ifpri_mapspam.YEARS), 2): + one, two = ( + ifpri_mapspam.YEAR_TO_UNRECOVERABLE_CROP_NAMES[before], + ifpri_mapspam.YEAR_TO_UNRECOVERABLE_CROP_NAMES[after], + ) + shared = len(one & two) / len(one | two) + if (before, after) in ifpri_mapspam.SPANS_WITH_COMPARABLE_CROP_NAMES: + assert shared > 0.5, ( + f"{before}->{after} called comparable but shares {shared:.0%}" + ) + else: + assert shared < 0.2, ( + f"{before}->{after} called incomparable but shares {shared:.0%}" + ) + + +@pytest.mark.parametrize("quantity", sorted(ifpri_mapspam.Quantity)) +@pytest.mark.parametrize("year", sorted(ifpri_mapspam.YEARS)) +def test_get_band_name_is_positionally_consistent( + year: int, quantity: ifpri_mapspam.Quantity +) -> None: + # `get_band_name` indexes the band list by the crop's position in the taxonomy + # enum. Nothing about a band name reveals a reordering, so a crop silently reading another + # crop's raster is the failure this exists to catch. Checked by matching each band against + # the slugified enum *value*, which is what the band name is built from. + for crop in ifpri_mapspam.YEAR_TO_CROP_CLS[year]: + band = ifpri_mapspam.get_band_name( + quantity=quantity, reported_crop_name=crop.name, year=year + ) + assert crop.value.lower().replace(" ", "-") in band, ( + f"{year} {quantity} {crop.name} resolved to {band}" + ) + + +@pytest.mark.parametrize("quantity", sorted(ifpri_mapspam.Quantity)) +@pytest.mark.parametrize("year", sorted(ifpri_mapspam.YEARS)) +def test_get_band_name_is_injective( + year: int, quantity: ifpri_mapspam.Quantity +) -> None: + # Two crops resolving to one band would silently double-count it + bands = [ + ifpri_mapspam.get_band_name( + quantity=quantity, reported_crop_name=crop.name, year=year + ) + for crop in ifpri_mapspam.YEAR_TO_CROP_CLS[year] + ] + assert len(set(bands)) == len(bands) + + +@pytest.mark.parametrize("year", sorted(ifpri_mapspam.YEARS)) +def test_unrecoverable_crop_names_partition_the_taxonomy(year: int) -> None: + # Every crop a snapshot reports is either reachable from a recoverable crop or unrecoverable, + # never both and never neither -- otherwise the share denominator would drop or double-count + # a band. + reported = {e.name for e in ifpri_mapspam.YEAR_TO_CROP_CLS[year]} + unrecoverable = ifpri_mapspam.YEAR_TO_UNRECOVERABLE_CROP_NAMES[year] + reachable = { + ifpri_mapspam.get_reported_crop_name(canonical_crop_name=name, year=year) + for name in ifpri_mapspam.RECOVERABLE_CROP_NAMES + } + assert unrecoverable | reachable == reported + assert not (unrecoverable & reachable) + + +def test_unrecoverable_crop_names_grow_with_the_taxonomy() -> None: + # 2000 leaves only its two catch-alls unreachable; later releases name crops the 2000 + # taxonomy has no route to, and 2020 adds four more again. + counts = { + year: len(ifpri_mapspam.YEAR_TO_UNRECOVERABLE_CROP_NAMES[year]) + for year in ifpri_mapspam.YEARS + } + assert counts == {2000: 2, 2005: 10, 2010: 10, 2020: 14} diff --git a/jdluc/__tests__/datasets/worldbank_jurisdictions_test.py b/jdluc/__tests__/datasets/worldbank_jurisdictions_test.py new file mode 100644 index 0000000..88c6f79 --- /dev/null +++ b/jdluc/__tests__/datasets/worldbank_jurisdictions_test.py @@ -0,0 +1,206 @@ +import geopandas +import pytest +import shapely + +from jdluc.datasets.worldbank_jurisdictions import ( + AdminLevel, + get_jurisdiction_for_admin_level, + get_ten_degree_tile_ids_for_admin_id, + iter_jurisdiction_for_iso_3166_tile_id, + with_taiwan_carved_out_of_china, + with_taiwan_province_separated, +) +from jdluc.tiling import get_box_for_tile_id + + +@pytest.mark.integration +def test_get_ten_degree_tile_ids_for_admin_id() -> None: + iso_3166 = "BRA" # Brazil + expected = [ + "00N_040W", + "00N_050W", + "00N_060W", + "00N_070W", + "00N_080W", + "10N_050W", + "10N_060W", + "10N_070W", + "10N_080W", + "10S_040W", + "10S_050W", + "10S_060W", + "10S_070W", + "10S_080W", + # This tile is not in the GNW tileset + # "20S_030W", + "20S_050W", + "20S_060W", + "30S_060W", + ] + assert ( + sorted( + get_ten_degree_tile_ids_for_admin_id( + admin_id=iso_3166, admin_level=AdminLevel.NATIONAL + ) + ) + == expected + ) + + +@pytest.mark.integration +def test_iter_jurisdiction_for_iso_3166_tile_id_loses_no_province() -> None: + # The tile prefilter is only safe because every province overlaps at least one of the tiles + # its country covers: that is what keeps attribute.merge_dfs summing over an unchanged + # index, and it is why the prefilter needed no cache version bump. A province is allowed to + # be absent only when it covers no published tile at all. + iso_3166 = "BRA" + seen = { + jurisdiction.id + for tile_id in get_ten_degree_tile_ids_for_admin_id( + admin_id=iso_3166, admin_level=AdminLevel.NATIONAL + ) + for jurisdiction in iter_jurisdiction_for_iso_3166_tile_id( + admin_level=AdminLevel.PROVINCIAL, iso_3166=iso_3166, tile_id=tile_id + ) + } + expected = { + str(admin_id) + for admin_id in get_jurisdiction_for_admin_level( + admin_level=AdminLevel.PROVINCIAL + ).index + if str(admin_id).startswith(iso_3166) + and get_ten_degree_tile_ids_for_admin_id( + admin_id=str(admin_id), admin_level=AdminLevel.PROVINCIAL + ) + } + assert seen == expected + + +@pytest.mark.integration +def test_iter_jurisdiction_for_iso_3166_tile_id_only_yields_overlapping() -> None: + # And the filter must actually filter: a province yielded for a tile has to intersect it, + # otherwise the prefilter is a no-op and the speedup is imaginary. + tile_id = "10S_050W" + box = get_box_for_tile_id(tile_id=tile_id) + yielded = list( + iter_jurisdiction_for_iso_3166_tile_id( + admin_level=AdminLevel.PROVINCIAL, iso_3166="BRA", tile_id=tile_id + ) + ) + assert yielded, f"{tile_id} is a Brazilian tile, so something should overlap it" + assert all(box.intersects(j.geometry) for j in yielded) + + +@pytest.mark.integration +def test_every_admin_1_prefix_has_an_admin_0_row() -> None: + # Provinces are selected by their country's `admin_id` prefix, so a province whose prefix + # names no national row is missing from every output. Upstream, Taiwan was one. + national = set( + get_jurisdiction_for_admin_level(admin_level=AdminLevel.NATIONAL).index + ) + prefixes = { + str(admin_id)[:3] + for admin_id in get_jurisdiction_for_admin_level( + admin_level=AdminLevel.PROVINCIAL + ).index + } + assert prefixes - national == set() + + +MAINLAND_CHINA = shapely.Polygon( + shell=[ + (122.7, 37.4), + (121.6, 28.3), + (109.9, 20.2), + (100.1, 21.5), + (89.0, 27.3), + (78.4, 32.5), + (73.5, 39.4), + (86.9, 49.1), + (105.0, 41.6), + (120.9, 53.3), + (134.8, 48.4), + (130.6, 42.4), + ] +) +TAIWAN = shapely.Polygon( + shell=[ + (120.9, 21.9), + (120.0, 23.1), + (120.2, 23.8), + (121.0, 25.0), + (121.5, 25.3), + (121.9, 25.1), + (121.9, 24.5), + (121.4, 23.1), + ] +) +JAPAN = shapely.Polygon( + shell=[ + (130.7, 31.0), + (129.7, 32.8), + (130.0, 33.8), + (132.6, 35.5), + (136.8, 37.5), + (139.0, 37.9), + (140.0, 40.6), + (140.0, 41.6), + (140.4, 43.3), + (141.7, 45.5), + (145.3, 44.3), + (145.8, 43.3), + (143.2, 41.9), + (141.5, 41.4), + (142.1, 39.5), + (141.0, 38.2), + (140.9, 35.7), + (138.8, 34.6), + (135.8, 33.4), + (133.0, 32.7), + (131.9, 31.6), + ] +) + + +def test_with_taiwan_carved_out_of_china() -> None: + from_world_bank = geopandas.GeoDataFrame( + crs=4326, + data={"ISO_A3": ["CHN", "JPN"], "NAM_0": ["China", "Japan"]}, + geometry=[shapely.union(MAINLAND_CHINA, TAIWAN), JAPAN], + ) + result = with_taiwan_carved_out_of_china(admin_0=from_world_bank).set_index( + "ISO_A3" + ) + assert result.crs == from_world_bank.crs + assert len(result) == 3 + assert list(result.columns) == ["NAM_0", "geometry"] + assert set(result.index) == {"CHN", "JPN", "TWN"} + assert result.loc["TWN"]["NAM_0"] == "Taiwan" + assert shapely.equals(result.loc["CHN"].geometry, MAINLAND_CHINA) + assert shapely.equals(result.loc["JPN"].geometry, JAPAN) + assert shapely.equals(result.loc["TWN"].geometry, TAIWAN) + + +def test_with_taiwan_province_separated() -> None: + from_world_bank = geopandas.GeoDataFrame( + crs=4326, + data={ + "ADM1CD_c": ["CHNXXX", "JPNXXX", "TWN001"], + "ISO_A3": ["CHN", "JPN", "CHN"], + "NAM_0": ["China", "Japan", "China"], + "NAM_1": ["merged", "merged", "Taiwan Sheng"], + }, + geometry=[MAINLAND_CHINA, JAPAN, TAIWAN], + ) + result = with_taiwan_province_separated(admin_1=from_world_bank).set_index( + "ADM1CD_c" + ) + assert result.crs == from_world_bank.crs + assert len(result) == 3 + assert list(result.columns) == ["ISO_A3", "NAM_0", "NAM_1", "geometry"] + assert set(result.index) == {"CHNXXX", "JPNXXX", "TWN001"} + assert result.loc["TWN001"]["ISO_A3"] == "TWN" + assert result.loc["TWN001"]["NAM_0"] == "Taiwan" + assert shapely.equals(result.loc["CHNXXX"].geometry, MAINLAND_CHINA) + assert shapely.equals(result.loc["JPNXXX"].geometry, JAPAN) + assert shapely.equals(result.loc["TWN001"].geometry, TAIWAN) diff --git a/jdluc/__tests__/datasets_test.py b/jdluc/__tests__/datasets_test.py deleted file mode 100644 index bbd770e..0000000 --- a/jdluc/__tests__/datasets_test.py +++ /dev/null @@ -1,440 +0,0 @@ -import itertools - -import geopandas -import pandas -import pytest -import shapely - -from jdluc.datasets import NAME_TO_CLS, base, faostat, ifpri_mapspam -from jdluc.datasets.glad_glcluc import flatten_ranges -from jdluc.datasets.worldbank_jurisdictions import ( - AdminLevel, - get_jurisdiction_for_admin_level, - get_ten_degree_tile_ids_for_admin_id, - iter_jurisdiction_for_iso_3166_tile_id, - with_taiwan_carved_out_of_china, - with_taiwan_province_separated, -) -from jdluc.tiling import get_box_for_tile_id - - -def test_flatten_ranges() -> None: - assert flatten_ranges(range(5), range(5, 10), range(10, 15)) == list(range(15)) - - -@pytest.mark.integration -def test_get_ten_degree_tile_ids_for_admin_id() -> None: - iso_3166 = "BRA" # Brazil - expected = [ - "00N_040W", - "00N_050W", - "00N_060W", - "00N_070W", - "00N_080W", - "10N_050W", - "10N_060W", - "10N_070W", - "10N_080W", - "10S_040W", - "10S_050W", - "10S_060W", - "10S_070W", - "10S_080W", - # This tile is not in the GNW tileset - # "20S_030W", - "20S_050W", - "20S_060W", - "30S_060W", - ] - assert ( - sorted( - get_ten_degree_tile_ids_for_admin_id( - admin_id=iso_3166, admin_level=AdminLevel.NATIONAL - ) - ) - == expected - ) - - -@pytest.mark.integration -def test_iter_jurisdiction_for_iso_3166_tile_id_loses_no_province() -> None: - # The tile prefilter is only safe because every province overlaps at least one of the tiles - # its country covers: that is what keeps attribute.merge_dfs summing over an unchanged - # index, and it is why the prefilter needed no cache version bump. A province is allowed to - # be absent only when it covers no published tile at all. - iso_3166 = "BRA" - seen = { - jurisdiction.id - for tile_id in get_ten_degree_tile_ids_for_admin_id( - admin_id=iso_3166, admin_level=AdminLevel.NATIONAL - ) - for jurisdiction in iter_jurisdiction_for_iso_3166_tile_id( - admin_level=AdminLevel.PROVINCIAL, iso_3166=iso_3166, tile_id=tile_id - ) - } - expected = { - str(admin_id) - for admin_id in get_jurisdiction_for_admin_level( - admin_level=AdminLevel.PROVINCIAL - ).index - if str(admin_id).startswith(iso_3166) - and get_ten_degree_tile_ids_for_admin_id( - admin_id=str(admin_id), admin_level=AdminLevel.PROVINCIAL - ) - } - assert seen == expected - - -@pytest.mark.integration -def test_iter_jurisdiction_for_iso_3166_tile_id_only_yields_overlapping() -> None: - # And the filter must actually filter: a province yielded for a tile has to intersect it, - # otherwise the prefilter is a no-op and the speedup is imaginary. - tile_id = "10S_050W" - box = get_box_for_tile_id(tile_id=tile_id) - yielded = list( - iter_jurisdiction_for_iso_3166_tile_id( - admin_level=AdminLevel.PROVINCIAL, iso_3166="BRA", tile_id=tile_id - ) - ) - assert yielded, f"{tile_id} is a Brazilian tile, so something should overlap it" - assert all(box.intersects(j.geometry) for j in yielded) - - -@pytest.mark.integration -def test_every_admin_1_prefix_has_an_admin_0_row() -> None: - # Provinces are selected by their country's `admin_id` prefix, so a province whose prefix - # names no national row is missing from every output. Upstream, Taiwan was one. - national = set( - get_jurisdiction_for_admin_level(admin_level=AdminLevel.NATIONAL).index - ) - prefixes = { - str(admin_id)[:3] - for admin_id in get_jurisdiction_for_admin_level( - admin_level=AdminLevel.PROVINCIAL - ).index - } - assert prefixes - national == set() - - -@pytest.mark.parametrize("year", sorted(ifpri_mapspam.YEARS)) -def test_get_reported_crop_name_always_names_a_band_that_year_has(year: int) -> None: - # The load-bearing property: whatever it returns must be readable from that snapshot. A - # rename or a new group that broke this would otherwise surface as a KeyError deep in a - # pipeline run. - names = {e.name for e in ifpri_mapspam.YEAR_TO_CROP_CLS[year]} - for crop_name in sorted(ifpri_mapspam.RECOVERABLE_CROP_NAMES): - assert ( - ifpri_mapspam.get_reported_crop_name( - canonical_crop_name=crop_name, year=year - ) - in names - ) - - -@pytest.mark.parametrize("year", sorted(ifpri_mapspam.YEARS)) -def test_reported_names_collapse_only_in_the_decompose_year(year: int) -> None: - # 19 constituents share 6 groups in 2000, so 13 names disappear; every other year names all - # 32 separately. This is the whole reason the decomposition exists, asserted directly. - reported = [ - ifpri_mapspam.get_reported_crop_name(canonical_crop_name=name, year=year) - for name in sorted(ifpri_mapspam.RECOVERABLE_CROP_NAMES) - ] - collapsed = len(reported) - len(set(reported)) - assert collapsed == (13 if year == ifpri_mapspam.DECOMPOSITION_YEAR else 0) - - -@pytest.mark.parametrize( - "year", sorted(set(ifpri_mapspam.YEARS) - {ifpri_mapspam.DECOMPOSITION_YEAR}) -) -def test_get_reported_crop_name_defers_to_the_rename_outside_2000( - year: int, -) -> None: - for crop_name in sorted(ifpri_mapspam.RECOVERABLE_CROP_NAMES): - assert ifpri_mapspam.get_reported_crop_name( - canonical_crop_name=crop_name, year=year - ) == ifpri_mapspam.get_renamed_crop_name( - canonical_crop_name=crop_name, year=year - ) - - -def test_is_reported_as_group_is_exactly_the_constituents_in_2000() -> None: - grouped = { - (crop_name, year) - for year in ifpri_mapspam.YEARS - for crop_name in ifpri_mapspam.RECOVERABLE_CROP_NAMES - if ifpri_mapspam.is_reported_as_group(canonical_crop_name=crop_name, year=year) - } - assert grouped == { - (crop_name, ifpri_mapspam.DECOMPOSITION_YEAR) - for crop_name in ifpri_mapspam.CONSTITUENT_TO_GROUP_NAME - } - - -def test_comparable_spans_agree_with_the_vocabularies() -> None: - # Not a restatement of the constant: it checks the constant against how much of the - # unrecoverable vocabulary the two snapshots actually share. Listing a pair as comparable - # when their names barely overlap is what would make name-matching measure the taxonomy - # instead of the land. - for before, after in itertools.combinations(sorted(ifpri_mapspam.YEARS), 2): - one, two = ( - ifpri_mapspam.YEAR_TO_UNRECOVERABLE_CROP_NAMES[before], - ifpri_mapspam.YEAR_TO_UNRECOVERABLE_CROP_NAMES[after], - ) - shared = len(one & two) / len(one | two) - if (before, after) in ifpri_mapspam.SPANS_WITH_COMPARABLE_CROP_NAMES: - assert shared > 0.5, ( - f"{before}->{after} called comparable but shares {shared:.0%}" - ) - else: - assert shared < 0.2, ( - f"{before}->{after} called incomparable but shares {shared:.0%}" - ) - - -@pytest.mark.parametrize("quantity", sorted(ifpri_mapspam.Quantity)) -@pytest.mark.parametrize("year", sorted(ifpri_mapspam.YEARS)) -def test_get_band_name_is_positionally_consistent( - year: int, quantity: ifpri_mapspam.Quantity -) -> None: - # `get_band_name` indexes the band list by the crop's position in the taxonomy - # enum. Nothing about a band name reveals a reordering, so a crop silently reading another - # crop's raster is the failure this exists to catch. Checked by matching each band against - # the slugified enum *value*, which is what the band name is built from. - for crop in ifpri_mapspam.YEAR_TO_CROP_CLS[year]: - band = ifpri_mapspam.get_band_name( - quantity=quantity, reported_crop_name=crop.name, year=year - ) - assert crop.value.lower().replace(" ", "-") in band, ( - f"{year} {quantity} {crop.name} resolved to {band}" - ) - - -@pytest.mark.parametrize("quantity", sorted(ifpri_mapspam.Quantity)) -@pytest.mark.parametrize("year", sorted(ifpri_mapspam.YEARS)) -def test_get_band_name_is_injective( - year: int, quantity: ifpri_mapspam.Quantity -) -> None: - # Two crops resolving to one band would silently double-count it - bands = [ - ifpri_mapspam.get_band_name( - quantity=quantity, reported_crop_name=crop.name, year=year - ) - for crop in ifpri_mapspam.YEAR_TO_CROP_CLS[year] - ] - assert len(set(bands)) == len(bands) - - -@pytest.mark.parametrize("year", sorted(ifpri_mapspam.YEARS)) -def test_unrecoverable_crop_names_partition_the_taxonomy(year: int) -> None: - # Every crop a snapshot reports is either reachable from a recoverable crop or unrecoverable, - # never both and never neither -- otherwise the share denominator would drop or double-count - # a band. - reported = {e.name for e in ifpri_mapspam.YEAR_TO_CROP_CLS[year]} - unrecoverable = ifpri_mapspam.YEAR_TO_UNRECOVERABLE_CROP_NAMES[year] - reachable = { - ifpri_mapspam.get_reported_crop_name(canonical_crop_name=name, year=year) - for name in ifpri_mapspam.RECOVERABLE_CROP_NAMES - } - assert unrecoverable | reachable == reported - assert not (unrecoverable & reachable) - - -def test_unrecoverable_crop_names_grow_with_the_taxonomy() -> None: - # 2000 leaves only its two catch-alls unreachable; later releases name crops the 2000 - # taxonomy has no route to, and 2020 adds four more again. - counts = { - year: len(ifpri_mapspam.YEAR_TO_UNRECOVERABLE_CROP_NAMES[year]) - for year in ifpri_mapspam.YEARS - } - assert counts == {2000: 2, 2005: 10, 2010: 10, 2020: 14} - - -@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) - - -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) - - -MAINLAND_CHINA = shapely.Polygon( - shell=[ - (122.7, 37.4), - (121.6, 28.3), - (109.9, 20.2), - (100.1, 21.5), - (89.0, 27.3), - (78.4, 32.5), - (73.5, 39.4), - (86.9, 49.1), - (105.0, 41.6), - (120.9, 53.3), - (134.8, 48.4), - (130.6, 42.4), - ] -) -TAIWAN = shapely.Polygon( - shell=[ - (120.9, 21.9), - (120.0, 23.1), - (120.2, 23.8), - (121.0, 25.0), - (121.5, 25.3), - (121.9, 25.1), - (121.9, 24.5), - (121.4, 23.1), - ] -) -JAPAN = shapely.Polygon( - shell=[ - (130.7, 31.0), - (129.7, 32.8), - (130.0, 33.8), - (132.6, 35.5), - (136.8, 37.5), - (139.0, 37.9), - (140.0, 40.6), - (140.0, 41.6), - (140.4, 43.3), - (141.7, 45.5), - (145.3, 44.3), - (145.8, 43.3), - (143.2, 41.9), - (141.5, 41.4), - (142.1, 39.5), - (141.0, 38.2), - (140.9, 35.7), - (138.8, 34.6), - (135.8, 33.4), - (133.0, 32.7), - (131.9, 31.6), - ] -) - - -def test_with_taiwan_carved_out_of_china() -> None: - from_world_bank = geopandas.GeoDataFrame( - crs=4326, - data={"ISO_A3": ["CHN", "JPN"], "NAM_0": ["China", "Japan"]}, - geometry=[shapely.union(MAINLAND_CHINA, TAIWAN), JAPAN], - ) - result = with_taiwan_carved_out_of_china(admin_0=from_world_bank).set_index( - "ISO_A3" - ) - assert result.crs == from_world_bank.crs - assert len(result) == 3 - assert list(result.columns) == ["NAM_0", "geometry"] - assert set(result.index) == {"CHN", "JPN", "TWN"} - assert result.loc["TWN"]["NAM_0"] == "Taiwan" - assert shapely.equals(result.loc["CHN"].geometry, MAINLAND_CHINA) - assert shapely.equals(result.loc["JPN"].geometry, JAPAN) - assert shapely.equals(result.loc["TWN"].geometry, TAIWAN) - - -def test_with_taiwan_province_separated() -> None: - from_world_bank = geopandas.GeoDataFrame( - crs=4326, - data={ - "ADM1CD_c": ["CHNXXX", "JPNXXX", "TWN001"], - "ISO_A3": ["CHN", "JPN", "CHN"], - "NAM_0": ["China", "Japan", "China"], - "NAM_1": ["merged", "merged", "Taiwan Sheng"], - }, - geometry=[MAINLAND_CHINA, JAPAN, TAIWAN], - ) - result = with_taiwan_province_separated(admin_1=from_world_bank).set_index( - "ADM1CD_c" - ) - assert result.crs == from_world_bank.crs - assert len(result) == 3 - assert list(result.columns) == ["ISO_A3", "NAM_0", "NAM_1", "geometry"] - assert set(result.index) == {"CHNXXX", "JPNXXX", "TWN001"} - assert result.loc["TWN001"]["ISO_A3"] == "TWN" - assert result.loc["TWN001"]["NAM_0"] == "Taiwan" - assert shapely.equals(result.loc["CHNXXX"].geometry, MAINLAND_CHINA) - assert shapely.equals(result.loc["JPNXXX"].geometry, JAPAN) - assert shapely.equals(result.loc["TWN001"].geometry, TAIWAN) diff --git a/jdluc/datasets/ifpri_mapspam.py b/jdluc/datasets/ifpri_mapspam.py index e99662a..92346ba 100644 --- a/jdluc/datasets/ifpri_mapspam.py +++ b/jdluc/datasets/ifpri_mapspam.py @@ -315,8 +315,8 @@ def get_renamed_crop_name(canonical_crop_name: str, year: int) -> str: """`canonical_crop_name` as `year`'s taxonomy spells it, matched on the shared crop label. The rename alone -- it cannot see groups, so it is wrong for a constituent in the decompose - year. `get_reported_crop_name` is the one to call; this stays separate because `datasets_test` - checks that function against it. + year. `get_reported_crop_name` is the one to call; this stays separate because + `ifpri_mapspam_test` checks that function against it. """ return YEAR_TO_CROP_CLS[year](CANONICAL_CROP_CLS[canonical_crop_name].value).name @@ -605,7 +605,7 @@ def get_band_name(quantity: Quantity, reported_crop_name: str, year: int) -> str renames or groups is a `KeyError`. Resolve with `get_reported_crop_name` first. Positional: the band list is built from the same enum in the same order, so a crop's index in - the taxonomy is its index among the bands. `datasets_test` pins that for every year and + the taxonomy is its index among the bands. `ifpri_mapspam_test` pins that for every year and quantity, since nothing in the band names themselves would reveal a reordering. """ crop_cls = YEAR_TO_CROP_CLS[year] From d205a9930450ebb7d2818c2ad757be417eb5d319 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Wed, 28 Oct 2026 15:31:19 -0700 Subject: [PATCH 03/21] Add ingestion logic for FDaP 2026a for cocoa, including pre-processing which binarizes using per-country/year thresholds, applies jurisdictional exclusion lists, and smooths temporally --- jdluc/__tests__/datasets/fdap_test.py | 140 ++++++++ jdluc/datasets/__init__.py | 3 + jdluc/datasets/fdap.py | 451 ++++++++++++++++++++++++++ jdluc/ingest.py | 1 + 4 files changed, 595 insertions(+) create mode 100644 jdluc/__tests__/datasets/fdap_test.py create mode 100644 jdluc/datasets/fdap.py diff --git a/jdluc/__tests__/datasets/fdap_test.py b/jdluc/__tests__/datasets/fdap_test.py new file mode 100644 index 0000000..c1b4c34 --- /dev/null +++ b/jdluc/__tests__/datasets/fdap_test.py @@ -0,0 +1,140 @@ +import math + +import geopandas +import numpy +import pytest +import rasterio.transform +import shapely + +from jdluc.datasets import fdap +from jdluc.datasets.worldbank_jurisdictions import AdminLevel + + +def test_get_pixel_range_gives_each_degree_the_pixels_centred_in_it() -> None: + pixel_ranges = [ + fdap.get_pixel_range(start_degrees=degrees, stop_degrees=degrees + 1) + for degrees in range(-180, 180) + ] + # Read along the pixels, the degrees their centres fall in are every degree once, in order + assert [ + math.floor((pixel + 0.5) * fdap.PIXEL_DEGREES) + for pixels in pixel_ranges + for pixel in (pixels.start - 1, pixels.start, pixels.stop - 1, pixels.stop) + ] == [ + centre_degrees + for degrees in range(-180, 180) + for centre_degrees in (degrees - 1, degrees, degrees, degrees + 1) + ] + + +@pytest.mark.parametrize("start_degrees", range(-180, 180)) +def test_get_pixel_range_reads_each_source_inside_its_file(start_degrees: int) -> None: + # A source file is 11133 pixels a side, from the pixel holding its west or south edge + pixels = fdap.get_pixel_range( + start_degrees=start_degrees, stop_degrees=start_degrees + 1 + ) + first_pixel = math.floor(start_degrees / fdap.PIXEL_DEGREES) + assert first_pixel <= pixels.start + assert pixels.stop <= first_pixel + 11_133 + + +# Three countries side by side, one degree each: AAA has fitted thresholds, BBB is a producer +# without them, and CCC is not a producer. AAA's western half is an excluded province. +FDAP_NATIONAL = geopandas.GeoSeries( + crs=4326, + data=[shapely.box(0, 0, 1, 1), shapely.box(1, 0, 2, 1), shapely.box(2, 0, 3, 1)], + index=["AAA", "BBB", "CCC"], +) +FDAP_PROVINCIAL = geopandas.GeoSeries( + crs=4326, data=[shapely.box(0, 0, 0.5, 1)], index=["AAA001"] +) +FDAP_FLOOR = 0.5 +# Four pixels per country, of which AAA's first two are in the excluded province +FDAP_TRANSFORM = rasterio.transform.from_origin(0, 1, 0.25, 1) +FDAP_SHAPE = (1, 12) +FDAP_ISO_3166_TO_UNEXCLUDED_PIXELS = {"AAA": slice(2, 4), "BBB": slice(4, 8)} + + +def get_fdap_mask(fill_probability: float, year: int) -> numpy.ndarray: + return fdap.get_mask( + admin_level_to_geometries={ + AdminLevel.NATIONAL: FDAP_NATIONAL, + AdminLevel.PROVINCIAL: FDAP_PROVINCIAL, + AdminLevel.DISTRICT: geopandas.GeoSeries(crs=4326), + }, + excluded_admin_ids={ + AdminLevel.NATIONAL: set(), + AdminLevel.PROVINCIAL: {"AAA001"}, + AdminLevel.DISTRICT: set(), + }, + floor=FDAP_FLOOR, + iso_3166_year_to_threshold={ + ("AAA", 2018): 0.2, + ("AAA", 2020): 0.8, + ("AAA", 2022): fdap.HIGHEST_THRESHOLD, + ("AAA", 2024): fdap.LOWEST_THRESHOLD, + }, + probability=numpy.full( + dtype="float32", fill_value=fill_probability, shape=FDAP_SHAPE + ), + producer_iso_3166s={"AAA", "BBB"}, + transform=FDAP_TRANSFORM, + year=year, + ) + + +@pytest.mark.parametrize("year", fdap.YEARS) +def test_get_mask_leaves_out_exclusions_and_non_producers_even_at_1( + year: int, +) -> None: + assert get_fdap_mask(fill_probability=1.0, year=year).tolist() == [ + [0] * 2 + [1] * 6 + [0] * 4 + ] + + +@pytest.mark.parametrize( + ("iso_3166", "year", "threshold"), + [ + pytest.param("AAA", 2018, FDAP_FLOOR, id="fitted-below-the-floor"), + pytest.param("AAA", 2020, 0.8, id="fitted-above-the-floor"), + pytest.param( + "AAA", 2022, fdap.HIGHEST_THRESHOLD, id="the-upper-end-of-the-sweep" + ), + pytest.param( + "AAA", 2024, FDAP_FLOOR, id="the-lower-end-of-the-sweep-below-the-floor" + ), + pytest.param( + "BBB", 2018, FDAP_FLOOR, id="a-producer-without-fitted-thresholds" + ), + ], +) +def test_get_mask_compares_strictly_with_the_threshold( + iso_3166: str, year: int, threshold: float +) -> None: + pixels = FDAP_ISO_3166_TO_UNEXCLUDED_PIXELS[iso_3166] + at = numpy.float32(threshold) + above = numpy.nextafter(at, numpy.float32(1)) + assert not get_fdap_mask(fill_probability=float(at), year=year)[:, pixels].any() + assert get_fdap_mask(fill_probability=float(above), year=year)[:, pixels].all() + + +@pytest.mark.parametrize( + ("masks", "expected"), + [ + pytest.param((1, 0, 1, 1), (1, 1, 1, 1), id="a-one-year-gap-in-2020-fills"), + pytest.param((0, 1, 0, 1), (0, 1, 1, 1), id="a-one-year-gap-in-2022-fills"), + pytest.param((1, 0, 0, 1), (1, 0, 0, 1), id="a-two-year-gap-stays"), + pytest.param((0, 1, 1, 1), (0, 1, 1, 1), id="2018-never-fills"), + pytest.param((1, 1, 1, 0), (1, 1, 1, 0), id="2024-never-fills"), + ], +) +def test_get_smoothed_masks_fills_one_year_gaps_in_the_inner_years( + masks: tuple[int, ...], expected: tuple[int, ...] +) -> None: + smoothed = fdap.get_smoothed_masks( + year_to_mask={ + year: numpy.full(dtype="uint8", fill_value=mask, shape=(1, 1)) + for year, mask in zip(fdap.YEARS, masks, strict=True) + } + ) + assert tuple(int(smoothed[year][0, 0]) for year in fdap.YEARS) == expected diff --git a/jdluc/datasets/__init__.py b/jdluc/datasets/__init__.py index cae9ac4..634f4dc 100644 --- a/jdluc/datasets/__init__.py +++ b/jdluc/datasets/__init__.py @@ -4,6 +4,7 @@ base, descals_oil_palm, faostat, + fdap, glad_glcluc, gnw_global_peatlands, gnw_harris_agb, @@ -33,6 +34,7 @@ def _generate_next_value_( FAOSTAT_LIVESTOCK_PATTERNS = enum.auto() FAOSTAT_PRODUCTION_CROPS = enum.auto() FAOSTAT_PRODUCTION_LIVESTOCK = enum.auto() + FDAP_COCOA = enum.auto() GLAD_GLCLUC = enum.auto() GNW_GLOBAL_PEATLANDS = enum.auto() GNW_HARRIS_AGB = enum.auto() @@ -65,6 +67,7 @@ def _generate_next_value_( DatasetName.FAOSTAT_LIVESTOCK_PATTERNS: faostat.LIVESTOCK_PATTERNS_DATASET, DatasetName.FAOSTAT_PRODUCTION_CROPS: faostat.CROP_DATASET, DatasetName.FAOSTAT_PRODUCTION_LIVESTOCK: faostat.LIVESTOCK_DATASET, + DatasetName.FDAP_COCOA: fdap.COCOA_DATASET, DatasetName.GLAD_GLCLUC: glad_glcluc.DATASET, DatasetName.GNW_GLOBAL_PEATLANDS: gnw_global_peatlands.DATASET, DatasetName.GNW_HARRIS_AGB: gnw_harris_agb.DATASET, diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py new file mode 100644 index 0000000..63f0bbf --- /dev/null +++ b/jdluc/datasets/fdap.py @@ -0,0 +1,451 @@ +"""Forest Data Partnership | Commodity probabilities, 2018-2024 + +license: CC BY 4.0 + +year: 2018, 2020, 2022, 2024 + +Google, for the Forest Data Partnership (2026). Commodity probabilities, release 2026a [Data set]. gs://earth-engine-public-requester-pays/forestdatapartnership/2026a -- the record ingested here, whose licence asks for the attribution "Produced by Google for the Forest Data Partnership" +TODO: add reference for the thresholds + +https://github.com/google/forest-data-partnership + +# Methodology + +- A multilayer perceptron probe, one per commodity, reads each 10 m pixel's Google AlphaEarth + Foundations satellite embedding, a yearly summary of satellite imagery, and gives the + probability that the pixel grows the commodity. It covers 24°S to 24°N +- The authors say the model maps too much of each tree crop on average, and that its results + change a lot from year to year + +We don't keep the probabilities. Ingest turns them into a binary mask for each year, in four steps: + +1. Merge: the one-degree source files are pieced together into ten-degree tiles at their own 10 m + resolution, without resampling. Neighbouring files overlap by one pixel, which is taken from + the file it falls in +2. Threshold: a pixel is in the mask when its probability is above its country's threshold for + that year, or above the commodity's floor if that is higher. A country without thresholds uses + the floor if FAOSTAT says it grows the commodity, and is left out otherwise +3. Exclude: pixels in the excluded areas are left out +4. Smooth: a pixel missing from the 2020 or 2022 mask is added back when it is in the masks of both + the year before and the year after. 2018 and 2024 are not smoothed +""" + +import enum +import logging +import math +import typing + +import geopandas +import google.auth +import numpy +import pandas +import rasterio +import rasterio.features +import rasterio.transform +import rasterio.windows +import shapely + +from jdluc import storage, tiling, utils +from jdluc.datasets import base, worldbank_jurisdictions + +logger = logging.getLogger(__name__) + +RELEASE = "2026a" +YEARS: typing.Final[tuple[int, ...]] = (2018, 2020, 2022, 2024) +# 10 m on the equator of the WGS 84 ellipsoid, whose semi-major axis is 6,378,137 m +PIXEL_DEGREES = 10 / (math.pi * 6_378_137 / 180) +LOWEST_THRESHOLD = 0.0025 +HIGHEST_THRESHOLD = 0.99 +NO_DATA = None + + +class Commodity(enum.StrEnum): + COCOA = enum.auto() + + +def get_pixel_range(start_degrees: int, stop_degrees: int) -> range: + """Global indices of the pixels whose centres fall in [start, stop) degrees, on either axis.""" + return range( + math.ceil(start_degrees / PIXEL_DEGREES - 0.5), + math.ceil(stop_degrees / PIXEL_DEGREES - 0.5), + ) + + +def get_mask( + admin_level_to_geometries: dict[ + worldbank_jurisdictions.AdminLevel, geopandas.GeoSeries + ], + excluded_admin_ids: dict[worldbank_jurisdictions.AdminLevel, set[str]], + floor: float, + iso_3166_year_to_threshold: dict[tuple[str, int], float], + probability: numpy.ndarray, + producer_iso_3166s: set[str], + transform: rasterio.Affine, + year: int, +) -> numpy.ndarray: + """1 where `probability` is strictly greater than its pixel's threshold, else 0.""" + box = shapely.box( + *rasterio.transform.array_bounds(*probability.shape, transform=transform) + ) + national = admin_level_to_geometries[worldbank_jurisdictions.AdminLevel.NATIONAL] + defaults = [ + (national[iso_3166], floor) + for iso_3166 in sorted(producer_iso_3166s) + if national[iso_3166].intersects(box) + ] + fitted = [ + (national[iso_3166], max(threshold, floor)) + for (iso_3166, the_year), threshold in sorted( + iso_3166_year_to_threshold.items() + ) + if the_year == year and national[iso_3166].intersects(box) + ] + exclusions = [ + (geometry, 1.0) + for admin_level, admin_ids in sorted(excluded_admin_ids.items()) + for geometry in admin_level_to_geometries[admin_level][sorted(admin_ids)] + if geometry.intersects(box) + ] + # Later shapes overwrite earlier ones, and nothing is greater than 1.0 + threshold = rasterio.features.rasterize( + dtype="float32", + fill=1, + out_shape=probability.shape, + shapes=defaults + fitted + exclusions, + transform=transform, + ) + return (probability > threshold).astype("uint8") + + +def get_smoothed_masks( + year_to_mask: dict[int, numpy.ndarray], +) -> dict[int, numpy.ndarray]: + """Each year's mask, with each inner year also 1 where the years either side are both 1.""" + years = sorted(year_to_mask) + first_year, *_, last_year = years + return { + year: ( + year_to_mask[year] + if year in {first_year, last_year} + else year_to_mask[year] + | (year_to_mask[years[idx - 1]] & year_to_mask[years[idx + 1]]) + ) + for idx, year in enumerate(years) + } + + +@utils.threadsafe_cache +def get_billing_project() -> str: + _, project = google.auth.default() + assert project is not None + return project + + +@utils.threadsafe_cache +def get_index() -> pandas.DataFrame: + return pandas.read_parquet( + path=f"gs://earth-engine-public-requester-pays/forestdatapartnership/{RELEASE:s}/fdap_index.parquet", + storage_options={"requester_pays": get_billing_project()}, + ) + + +def get_raster_dataset( + commodity: Commodity, + excluded_admin_ids: dict[worldbank_jurisdictions.AdminLevel, set[str]], + floor: float, + iso_3166_year_to_threshold: dict[tuple[str, int], float], + producer_iso_3166s: set[str], + version: str, +) -> base.RasterDataset: + iso_3166s = {iso_3166 for iso_3166, _ in iso_3166_year_to_threshold} + assert set(iso_3166_year_to_threshold) == { + (iso_3166, year) for iso_3166 in iso_3166s for year in YEARS + } + assert all( + LOWEST_THRESHOLD <= threshold <= HIGHEST_THRESHOLD + for threshold in iso_3166_year_to_threshold.values() + ) + assert iso_3166s <= producer_iso_3166s + assert set(excluded_admin_ids) == set(worldbank_jurisdictions.AdminLevel) + assert ( + not iso_3166s & excluded_admin_ids[worldbank_jurisdictions.AdminLevel.NATIONAL] + ) + + def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: + admin_level_to_geometries = { + admin_level: worldbank_jurisdictions.get_jurisdiction_for_admin_level( + admin_level=admin_level + ).geometry + for admin_level in worldbank_jurisdictions.AdminLevel + } + max_latitude, min_longitude = tiling.get_lat_lon_for_tile_id(tile_id=tile_id) + columns = get_pixel_range( + start_degrees=min_longitude, stop_degrees=min_longitude + 10 + ) + rows = get_pixel_range( + start_degrees=max_latitude - 10, stop_degrees=max_latitude + ) + index = get_index() + sources = index[ + (index["commodity"] == commodity) + & (index["lng"] >= min_longitude) + & (index["lng"] < min_longitude + 10) + & (index["lat"] >= max_latitude - 10) + & (index["lat"] < max_latitude) + ] + lat_lons = sorted( + set(zip(sources["lat"].tolist(), sources["lng"].tolist(), strict=True)) + ) + logger.info(f"Thresholding {len(lat_lons):d} source cells onto {tile_id=:s}") + with ( + rasterio.Env( + CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif", + GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", + GS_USER_PROJECT=get_billing_project(), + ), + # Blocks no source writes are zeros when the file closes + rasterio.open( + BIGTIFF="YES", + NUM_THREADS="ALL_CPUS", + compress="deflate", + count=len(YEARS), + crs=4326, + driver="GTiff", + dtype="uint8", + fp=local_path, + height=len(rows), + interleave="band", + mode="w", + nodata=NO_DATA, + tiled=True, + transform=rasterio.transform.from_origin( + columns.start * PIXEL_DEGREES, + rows.stop * PIXEL_DEGREES, + PIXEL_DEGREES, + PIXEL_DEGREES, + ), + width=len(columns), + ) as dataset, + ): + for latitude, longitude in lat_lons: + source_columns = get_pixel_range( + start_degrees=longitude, stop_degrees=longitude + 1 + ) + source_rows = get_pixel_range( + start_degrees=latitude, stop_degrees=latitude + 1 + ) + window = rasterio.windows.Window( + col_off=source_columns.start - columns.start, + height=len(source_rows), + row_off=rows.stop - source_rows.stop, + width=len(source_columns), + ) + year_to_mask: dict[int, numpy.ndarray] = {} + for year in YEARS: + (uri,) = sources[ + (sources["lat"] == latitude) + & (sources["lng"] == longitude) + & (sources["year"] == year) + ]["path"] + with rasterio.open(fp=storage.to_gdal_path(uri=uri)) as source: + left = round(source.transform.c / PIXEL_DEGREES) + top = round(source.transform.f / PIXEL_DEGREES) + # On the global grid, a whole number of pixels from (0, 0) + assert source.transform.almost_equals( + rasterio.transform.from_origin( + left * PIXEL_DEGREES, + top * PIXEL_DEGREES, + PIXEL_DEGREES, + PIXEL_DEGREES, + ), + precision=1e-12, + ) + assert ( + left <= source_columns.start + and source_columns.stop <= left + source.width + ) + assert ( + top - source.height <= source_rows.start + and source_rows.stop <= top + ) + # Neighbouring sources share a column or row, which only one of them reads + probability = source.read( + 1, + window=rasterio.windows.Window( + col_off=source_columns.start - left, + height=len(source_rows), + row_off=top - source_rows.stop, + width=len(source_columns), + ), + ) + year_to_mask[year] = get_mask( + admin_level_to_geometries=admin_level_to_geometries, + excluded_admin_ids=excluded_admin_ids, + floor=floor, + iso_3166_year_to_threshold=iso_3166_year_to_threshold, + probability=probability, + producer_iso_3166s=producer_iso_3166s, + transform=rasterio.windows.transform( + transform=dataset.transform, window=window + ), + year=year, + ) + year_to_smoothed_mask = get_smoothed_masks(year_to_mask=year_to_mask) + for band_idx, year in enumerate(YEARS, start=1): + dataset.write(year_to_smoothed_mask[year], band_idx, window=window) + + 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, + dtype="uint8", + no_data=NO_DATA, + partitioning=tiling.Partitioning.TEN_DEGREE_TILE, + product_name=f"{commodity:s}-{RELEASE:s}", + save_tile_id_to_local_path=save_tile_id_to_local_path, + source_name="fdap", + version=version, + ) + + +COCOA_DATASET = get_raster_dataset( + commodity=Commodity.COCOA, + excluded_admin_ids={ + worldbank_jurisdictions.AdminLevel.NATIONAL: set(), + worldbank_jurisdictions.AdminLevel.PROVINCIAL: set(), + worldbank_jurisdictions.AdminLevel.DISTRICT: set(), + }, + floor=0.54, + iso_3166_year_to_threshold={ + ("BRA", 2018): 0.923, + ("BRA", 2020): 0.946, + ("BRA", 2022): 0.976, + ("BRA", 2024): 0.958, + ("CIV", 2018): 0.972, + ("CIV", 2020): 0.982, + ("CIV", 2022): 0.969, + ("CIV", 2024): 0.961, + ("CMR", 2018): 0.372, + ("CMR", 2020): 0.434, + ("CMR", 2022): 0.573, + ("CMR", 2024): 0.616, + ("COD", 2018): 0.628, + ("COD", 2020): 0.81, + ("COD", 2022): 0.91, + ("COD", 2024): 0.919, + ("COL", 2018): 0.919, + ("COL", 2020): 0.851, + ("COL", 2022): 0.986, + ("COL", 2024): 0.971, + ("DOM", 2018): 0.016, + ("DOM", 2020): LOWEST_THRESHOLD, + ("DOM", 2022): 0.004, + ("DOM", 2024): 0.017, + ("ECU", 2018): 0.16, + ("ECU", 2020): 0.292, + ("ECU", 2022): 0.862, + ("ECU", 2024): 0.812, + ("GHA", 2018): 0.972, + ("GHA", 2020): HIGHEST_THRESHOLD, + ("GHA", 2022): HIGHEST_THRESHOLD, + ("GHA", 2024): HIGHEST_THRESHOLD, + ("IDN", 2018): 0.257, + ("IDN", 2020): 0.357, + ("IDN", 2022): 0.355, + ("IDN", 2024): 0.418, + ("IND", 2018): 0.066, + ("IND", 2020): 0.043, + ("IND", 2022): 0.552, + ("IND", 2024): 0.078, + ("LBR", 2018): HIGHEST_THRESHOLD, + ("LBR", 2020): 0.988, + ("LBR", 2022): 0.986, + ("LBR", 2024): 0.984, + ("NGA", 2018): 0.415, + ("NGA", 2020): 0.319, + ("NGA", 2022): 0.544, + ("NGA", 2024): 0.533, + ("PER", 2018): 0.558, + ("PER", 2020): 0.866, + ("PER", 2022): 0.91, + ("PER", 2024): 0.812, + ("PNG", 2018): 0.579, + ("PNG", 2020): 0.589, + ("PNG", 2022): 0.779, + ("PNG", 2024): 0.568, + ("TZA", 2018): 0.249, + ("TZA", 2020): 0.329, + ("TZA", 2022): 0.58, + ("TZA", 2024): 0.431, + ("UGA", 2018): LOWEST_THRESHOLD, + ("UGA", 2020): LOWEST_THRESHOLD, + ("UGA", 2022): 0.014, + ("UGA", 2024): 0.028, + ("VEN", 2018): 0.849, + ("VEN", 2020): 0.466, + ("VEN", 2022): 0.958, + ("VEN", 2024): 0.819, + }, + # Countries with FAOSTAT area harvested of cocoa (item 661) in any year from 2018 to 2024 + producer_iso_3166s={ + "AGO", + "BLZ", + "BOL", + "BRA", + "CAF", + "CIV", + "CMR", + "COD", + "COG", + "COL", + "COM", + "CRI", + "CUB", + "DMA", + "DOM", + "ECU", + "FJI", + "FSM", + "GAB", + "GHA", + "GIN", + "GNQ", + "GRD", + "GTM", + "GUY", + "HND", + "HTI", + "IDN", + "IND", + "JAM", + "LBR", + "LCA", + "LKA", + "MDG", + "MEX", + "MYS", + "NGA", + "NIC", + "PAN", + "PER", + "PHL", + "PNG", + "SLB", + "SLE", + "SLV", + "STP", + "SUR", + "TGO", + "THA", + "TLS", + "TTO", + "TZA", + "UGA", + "VCT", + "VEN", + "VNM", + "VUT", + "WSM", + }, + version="v0", +) diff --git a/jdluc/ingest.py b/jdluc/ingest.py index 8a35db1..c1a0dd9 100644 --- a/jdluc/ingest.py +++ b/jdluc/ingest.py @@ -19,6 +19,7 @@ uv run python -m jdluc.ingest LIAO_GACED30 IDN uv run python -m jdluc.ingest GPW_GRASSLAND BRA uv run python -m jdluc.ingest DESCALS_OIL_PALM IDN MYS + uv run python -m jdluc.ingest FDAP_COCOA CIV GHA """ import argparse From 7b1a833acc80d64256d91237fb2435a2f47cf175 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Wed, 28 Oct 2026 15:54:59 -0700 Subject: [PATCH 04/21] Add ingestion logic for FDaP 2026a for coffee --- jdluc/datasets/__init__.py | 2 + jdluc/datasets/fdap.py | 160 +++++++++++++++++++++++++++++++++++++ 2 files changed, 162 insertions(+) diff --git a/jdluc/datasets/__init__.py b/jdluc/datasets/__init__.py index 634f4dc..cb0e48c 100644 --- a/jdluc/datasets/__init__.py +++ b/jdluc/datasets/__init__.py @@ -35,6 +35,7 @@ def _generate_next_value_( FAOSTAT_PRODUCTION_CROPS = enum.auto() FAOSTAT_PRODUCTION_LIVESTOCK = enum.auto() FDAP_COCOA = enum.auto() + FDAP_COFFEE = enum.auto() GLAD_GLCLUC = enum.auto() GNW_GLOBAL_PEATLANDS = enum.auto() GNW_HARRIS_AGB = enum.auto() @@ -68,6 +69,7 @@ def _generate_next_value_( DatasetName.FAOSTAT_PRODUCTION_CROPS: faostat.CROP_DATASET, DatasetName.FAOSTAT_PRODUCTION_LIVESTOCK: faostat.LIVESTOCK_DATASET, DatasetName.FDAP_COCOA: fdap.COCOA_DATASET, + DatasetName.FDAP_COFFEE: fdap.COFFEE_DATASET, DatasetName.GLAD_GLCLUC: glad_glcluc.DATASET, DatasetName.GNW_GLOBAL_PEATLANDS: gnw_global_peatlands.DATASET, DatasetName.GNW_HARRIS_AGB: gnw_harris_agb.DATASET, diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index 63f0bbf..83c9cbb 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -61,6 +61,7 @@ class Commodity(enum.StrEnum): COCOA = enum.auto() + COFFEE = enum.auto() def get_pixel_range(start_degrees: int, stop_degrees: int) -> range: @@ -449,3 +450,162 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: }, version="v0", ) + +COFFEE_DATASET = get_raster_dataset( + commodity=Commodity.COFFEE, + excluded_admin_ids={ + worldbank_jurisdictions.AdminLevel.NATIONAL: set(), + worldbank_jurisdictions.AdminLevel.PROVINCIAL: set(), + worldbank_jurisdictions.AdminLevel.DISTRICT: set(), + }, + floor=0.33, + iso_3166_year_to_threshold={ + ("BRA", 2018): 0.554, + ("BRA", 2020): 0.658, + ("BRA", 2022): 0.763, + ("BRA", 2024): 0.593, + ("CIV", 2018): LOWEST_THRESHOLD, + ("CIV", 2020): 0.005, + ("CIV", 2022): 0.007, + ("CIV", 2024): 0.008, + ("COL", 2018): 0.285, + ("COL", 2020): 0.196, + ("COL", 2022): 0.845, + ("COL", 2024): 0.674, + ("ETH", 2018): 0.019, + ("ETH", 2020): 0.187, + ("ETH", 2022): 0.965, + ("ETH", 2024): 0.705, + ("GTM", 2018): 0.134, + ("GTM", 2020): 0.083, + ("GTM", 2022): 0.128, + ("GTM", 2024): 0.073, + ("HND", 2018): 0.189, + ("HND", 2020): 0.205, + ("HND", 2022): 0.524, + ("HND", 2024): 0.202, + ("IDN", 2018): 0.679, + ("IDN", 2020): 0.337, + ("IDN", 2022): 0.309, + ("IDN", 2024): 0.752, + ("IND", 2018): 0.055, + ("IND", 2020): 0.191, + ("IND", 2022): 0.316, + ("IND", 2024): 0.131, + ("KEN", 2018): 0.472, + ("KEN", 2020): 0.669, + ("KEN", 2022): HIGHEST_THRESHOLD, + ("KEN", 2024): 0.989, + ("MEX", 2018): 0.322, + ("MEX", 2020): 0.196, + ("MEX", 2022): 0.341, + ("MEX", 2024): 0.136, + ("NIC", 2018): 0.317, + ("NIC", 2020): 0.172, + ("NIC", 2022): 0.223, + ("NIC", 2024): 0.124, + ("PER", 2018): 0.024, + ("PER", 2020): 0.081, + ("PER", 2022): 0.337, + ("PER", 2024): 0.063, + ("PNG", 2018): 0.03, + ("PNG", 2020): 0.092, + ("PNG", 2022): 0.783, + ("PNG", 2024): 0.309, + ("SLV", 2018): 0.005, + ("SLV", 2020): LOWEST_THRESHOLD, + ("SLV", 2022): 0.024, + ("SLV", 2024): 0.006, + ("UGA", 2018): 0.751, + ("UGA", 2020): 0.741, + ("UGA", 2022): 0.979, + ("UGA", 2024): 0.956, + ("VNM", 2018): 0.208, + ("VNM", 2020): 0.049, + ("VNM", 2022): 0.52, + ("VNM", 2024): 0.296, + }, + # Countries with FAOSTAT area harvested of coffee (item 656) in any year from 2018 to 2024 + producer_iso_3166s={ + "AGO", + "BDI", + "BEN", + "BLZ", + "BOL", + "BRA", + "CAF", + "CHN", + "CIV", + "CMR", + "COD", + "COG", + "COK", + "COL", + "COM", + "CPV", + "CRI", + "CUB", + "DMA", + "DOM", + "ECU", + "ETH", + "FJI", + "GAB", + "GHA", + "GIN", + "GNQ", + "GTM", + "GUY", + "HND", + "HTI", + "IDN", + "IND", + "JAM", + "KEN", + "KHM", + "LAO", + "LBR", + "LKA", + "MDG", + "MEX", + "MMR", + "MOZ", + "MWI", + "MYS", + "NCL", + "NGA", + "NIC", + "NPL", + "PAN", + "PER", + "PHL", + "PNG", + "PRI", + "PRY", + "PYF", + "RWA", + "SAU", + "SLE", + "SLV", + "STP", + "SUR", + "TGO", + "THA", + "TLS", + "TON", + "TTO", + "TWN", + "TZA", + "UGA", + "USA", + "VCT", + "VEN", + "VNM", + "VUT", + "WSM", + "YEM", + "ZMB", + "ZWE", + }, + version="v0", +) From 260cd91a1961b378227899ce9a3980e1ad2874a7 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Wed, 28 Oct 2026 16:03:16 -0700 Subject: [PATCH 05/21] update docs --- docs/architecture.md | 2 +- docs/further_research.md | 2 +- docs/methodology.md | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/architecture.md b/docs/architecture.md index 1933fdc..8c882bb 100644 --- a/docs/architecture.md +++ b/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. diff --git a/docs/further_research.md b/docs/further_research.md index d8c7689..38d3ecd 100644 --- a/docs/further_research.md +++ b/docs/further_research.md @@ -29,7 +29,7 @@ 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. ## Shifting cultivation resolves to no destination diff --git a/docs/methodology.md b/docs/methodology.md index 1f73599..45deb0d 100644 --- a/docs/methodology.md +++ b/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. From dc2878eb796a6435473015eec04594955cc7d7c1 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Wed, 28 Oct 2026 16:41:17 -0700 Subject: [PATCH 06/21] fix quirk with four-band COG's using the fourth as an alpha channel --- jdluc/__tests__/geo_test.py | 29 +++++++++++++++++++++++++++++ jdluc/datasets/fdap.py | 1 + jdluc/geo.py | 1 + 3 files changed, 31 insertions(+) diff --git a/jdluc/__tests__/geo_test.py b/jdluc/__tests__/geo_test.py index d364377..fb5fad2 100644 --- a/jdluc/__tests__/geo_test.py +++ b/jdluc/__tests__/geo_test.py @@ -54,6 +54,35 @@ def test_validate_geotiff_holds_a_dataset_to_its_declared_dtype( validate_geotiff(dtype=declared, path_to_geotiff=path_to_geotiff) +@pytest.mark.parametrize( + ("photometric", "fails"), + ( + pytest.param(None, True, id="four-uint8-bands-default-to-rgba"), + pytest.param("MINISBLACK", False, id="minisblack-has-no-alpha"), + ), +) +def test_validate_geotiff_rejects_an_alpha_band( + tmp_path: pathlib.Path, photometric: str | None, fails: bool +) -> None: + path_to_geotiff = str(tmp_path / "geotiff.tif") + with rasterio.open( + path_to_geotiff, + "w", + count=4, + crs=4326, + driver="GTiff", + dtype="uint8", + height=4, + transform=rasterio.transform.from_origin(0, 10, 0.1, 0.1), + width=4, + **({} if photometric is None else {"photometric": photometric}), + ) as dataset: + dataset.write(numpy.zeros((4, 4, 4), dtype="uint8")) + + with pytest.raises(AssertionError) if fails else contextlib.nullcontext(): + validate_geotiff(dtype="uint8", path_to_geotiff=path_to_geotiff) + + @pytest.mark.parametrize( ("height", "width", "minimum_pixels", "expected"), ( diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index 83c9cbb..72a7e2e 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -218,6 +218,7 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: interleave="band", mode="w", nodata=NO_DATA, + photometric="MINISBLACK", tiled=True, transform=rasterio.transform.from_origin( columns.start * PIXEL_DEGREES, diff --git a/jdluc/geo.py b/jdluc/geo.py index 3badd3e..de6b36a 100644 --- a/jdluc/geo.py +++ b/jdluc/geo.py @@ -123,6 +123,7 @@ def validate_geotiff(dtype: str, path_to_geotiff: str) -> None: assert dataset.width > 0 assert dataset.transform.a > 0 assert dataset.transform.e < 0 + assert rasterio.enums.ColorInterp.alpha not in dataset.colorinterp def get_chunk_size( From 9471e1d39d786bf65bdabcd54fcad2b460d38c8f Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Wed, 28 Oct 2026 18:32:13 -0700 Subject: [PATCH 07/21] Log for each source one-degree tile to better show progress --- jdluc/datasets/fdap.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index 72a7e2e..f903811 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -229,7 +229,11 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: width=len(columns), ) as dataset, ): - for latitude, longitude in lat_lons: + for cell_idx, (latitude, longitude) in enumerate(lat_lons, start=1): + logger.info( + f"Thresholding source cell {cell_idx:d} of {len(lat_lons):d} " + f"({latitude=:d}, {longitude=:d}) onto {tile_id=:s}" + ) source_columns = get_pixel_range( start_degrees=longitude, stop_degrees=longitude + 1 ) From aecade2b66af6b05703bc58b8a250d833e6c0f8b Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Thu, 29 Oct 2026 10:00:55 -0700 Subject: [PATCH 08/21] update to final thresholds for coffee --- jdluc/datasets/fdap.py | 21 +++------------------ 1 file changed, 3 insertions(+), 18 deletions(-) diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index f903811..78126f9 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -162,10 +162,7 @@ def get_raster_dataset( assert set(iso_3166_year_to_threshold) == { (iso_3166, year) for iso_3166 in iso_3166s for year in YEARS } - assert all( - LOWEST_THRESHOLD <= threshold <= HIGHEST_THRESHOLD - for threshold in iso_3166_year_to_threshold.values() - ) + assert all(0 < threshold < 1 for threshold in iso_3166_year_to_threshold.values()) assert iso_3166s <= producer_iso_3166s assert set(excluded_admin_ids) == set(worldbank_jurisdictions.AdminLevel) assert ( @@ -469,7 +466,7 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("BRA", 2020): 0.658, ("BRA", 2022): 0.763, ("BRA", 2024): 0.593, - ("CIV", 2018): LOWEST_THRESHOLD, + ("CIV", 2018): 0.002, ("CIV", 2020): 0.005, ("CIV", 2022): 0.007, ("CIV", 2024): 0.008, @@ -497,14 +494,6 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("IND", 2020): 0.191, ("IND", 2022): 0.316, ("IND", 2024): 0.131, - ("KEN", 2018): 0.472, - ("KEN", 2020): 0.669, - ("KEN", 2022): HIGHEST_THRESHOLD, - ("KEN", 2024): 0.989, - ("MEX", 2018): 0.322, - ("MEX", 2020): 0.196, - ("MEX", 2022): 0.341, - ("MEX", 2024): 0.136, ("NIC", 2018): 0.317, ("NIC", 2020): 0.172, ("NIC", 2022): 0.223, @@ -513,12 +502,8 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("PER", 2020): 0.081, ("PER", 2022): 0.337, ("PER", 2024): 0.063, - ("PNG", 2018): 0.03, - ("PNG", 2020): 0.092, - ("PNG", 2022): 0.783, - ("PNG", 2024): 0.309, ("SLV", 2018): 0.005, - ("SLV", 2020): LOWEST_THRESHOLD, + ("SLV", 2020): 0.002, ("SLV", 2022): 0.024, ("SLV", 2024): 0.006, ("UGA", 2018): 0.751, From d1c4be2031b698f388771bb582ff24084ce74a6b Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Thu, 29 Oct 2026 11:06:56 -0700 Subject: [PATCH 09/21] update to final thresholds for cocoa --- jdluc/datasets/fdap.py | 54 ++++++------------------------------------ 1 file changed, 7 insertions(+), 47 deletions(-) diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index 78126f9..ca3d2bc 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -323,71 +323,31 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("BRA", 2018): 0.923, ("BRA", 2020): 0.946, ("BRA", 2022): 0.976, - ("BRA", 2024): 0.958, + ("BRA", 2024): 0.976, ("CIV", 2018): 0.972, ("CIV", 2020): 0.982, ("CIV", 2022): 0.969, ("CIV", 2024): 0.961, - ("CMR", 2018): 0.372, - ("CMR", 2020): 0.434, - ("CMR", 2022): 0.573, - ("CMR", 2024): 0.616, - ("COD", 2018): 0.628, - ("COD", 2020): 0.81, - ("COD", 2022): 0.91, - ("COD", 2024): 0.919, ("COL", 2018): 0.919, ("COL", 2020): 0.851, ("COL", 2022): 0.986, ("COL", 2024): 0.971, - ("DOM", 2018): 0.016, - ("DOM", 2020): LOWEST_THRESHOLD, - ("DOM", 2022): 0.004, - ("DOM", 2024): 0.017, - ("ECU", 2018): 0.16, - ("ECU", 2020): 0.292, + ("ECU", 2018): 0.862, + ("ECU", 2020): 0.862, ("ECU", 2022): 0.862, ("ECU", 2024): 0.812, ("GHA", 2018): 0.972, ("GHA", 2020): HIGHEST_THRESHOLD, ("GHA", 2022): HIGHEST_THRESHOLD, ("GHA", 2024): HIGHEST_THRESHOLD, - ("IDN", 2018): 0.257, - ("IDN", 2020): 0.357, - ("IDN", 2022): 0.355, - ("IDN", 2024): 0.418, - ("IND", 2018): 0.066, - ("IND", 2020): 0.043, - ("IND", 2022): 0.552, - ("IND", 2024): 0.078, - ("LBR", 2018): HIGHEST_THRESHOLD, - ("LBR", 2020): 0.988, - ("LBR", 2022): 0.986, - ("LBR", 2024): 0.984, - ("NGA", 2018): 0.415, - ("NGA", 2020): 0.319, - ("NGA", 2022): 0.544, - ("NGA", 2024): 0.533, + ("IDN", 2018): 0.540, + ("IDN", 2020): 0.540, + ("IDN", 2022): 0.540, + ("IDN", 2024): 0.540, ("PER", 2018): 0.558, ("PER", 2020): 0.866, ("PER", 2022): 0.91, ("PER", 2024): 0.812, - ("PNG", 2018): 0.579, - ("PNG", 2020): 0.589, - ("PNG", 2022): 0.779, - ("PNG", 2024): 0.568, - ("TZA", 2018): 0.249, - ("TZA", 2020): 0.329, - ("TZA", 2022): 0.58, - ("TZA", 2024): 0.431, - ("UGA", 2018): LOWEST_THRESHOLD, - ("UGA", 2020): LOWEST_THRESHOLD, - ("UGA", 2022): 0.014, - ("UGA", 2024): 0.028, - ("VEN", 2018): 0.849, - ("VEN", 2020): 0.466, - ("VEN", 2022): 0.958, - ("VEN", 2024): 0.819, }, # Countries with FAOSTAT area harvested of cocoa (item 661) in any year from 2018 to 2024 producer_iso_3166s={ From cc86a3fcc378c8726ec33d4b64f6e4430c4e2521 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Thu, 29 Oct 2026 11:51:20 -0700 Subject: [PATCH 10/21] Restrict producers to countries which devote >1% of their cropland to the crop --- jdluc/datasets/fdap.py | 82 ++++++------------------------------------ 1 file changed, 10 insertions(+), 72 deletions(-) diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index ca3d2bc..d826bcf 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -349,64 +349,38 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("PER", 2022): 0.91, ("PER", 2024): 0.812, }, - # Countries with FAOSTAT area harvested of cocoa (item 661) in any year from 2018 to 2024 + # FAOSTAT's producers: area harvested of cocoa (item 661) in some year from 2018 to + # 2024 and, for a country without thresholds, at least 1% of its cropland (item 6620) in one + # of those years. From the crops and land use bulk files: + # https://bulks-faostat.fao.org/production/Production_Crops_Livestock_E_All_Data_(Normalized).zip + # https://bulks-faostat.fao.org/production/Inputs_LandUse_E_All_Data_(Normalized).zip producer_iso_3166s={ - "AGO", - "BLZ", - "BOL", "BRA", - "CAF", "CIV", "CMR", - "COD", "COG", "COL", - "COM", "CRI", - "CUB", "DMA", "DOM", "ECU", - "FJI", "FSM", - "GAB", "GHA", - "GIN", "GNQ", "GRD", - "GTM", - "GUY", - "HND", - "HTI", "IDN", - "IND", - "JAM", "LBR", "LCA", - "LKA", - "MDG", - "MEX", - "MYS", "NGA", - "NIC", "PAN", "PER", - "PHL", "PNG", "SLB", "SLE", - "SLV", "STP", - "SUR", - "TGO", - "THA", - "TLS", "TTO", - "TZA", - "UGA", "VCT", "VEN", - "VNM", "VUT", "WSM", }, @@ -475,87 +449,51 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("VNM", 2022): 0.52, ("VNM", 2024): 0.296, }, - # Countries with FAOSTAT area harvested of coffee (item 656) in any year from 2018 to 2024 + # FAOSTAT's producers: area harvested of coffee (item 656) in some year from 2018 to + # 2024 and, for a country without thresholds, at least 1% of its cropland (item 6620) in one + # of those years. From the crops and land use bulk files: + # https://bulks-faostat.fao.org/production/Production_Crops_Livestock_E_All_Data_(Normalized).zip + # https://bulks-faostat.fao.org/production/Inputs_LandUse_E_All_Data_(Normalized).zip producer_iso_3166s={ - "AGO", - "BDI", - "BEN", - "BLZ", - "BOL", "BRA", "CAF", - "CHN", "CIV", "CMR", - "COD", "COG", - "COK", "COL", - "COM", - "CPV", "CRI", - "CUB", "DMA", "DOM", "ECU", "ETH", - "FJI", - "GAB", - "GHA", "GIN", "GNQ", "GTM", - "GUY", "HND", - "HTI", "IDN", "IND", "JAM", "KEN", - "KHM", "LAO", - "LBR", - "LKA", "MDG", "MEX", - "MMR", - "MOZ", - "MWI", - "MYS", - "NCL", - "NGA", "NIC", - "NPL", "PAN", "PER", "PHL", "PNG", "PRI", - "PRY", - "PYF", "RWA", - "SAU", - "SLE", "SLV", - "STP", - "SUR", "TGO", - "THA", "TLS", - "TON", "TTO", - "TWN", "TZA", "UGA", - "USA", "VCT", "VEN", "VNM", - "VUT", - "WSM", "YEM", - "ZMB", - "ZWE", }, version="v0", ) From e5d9e53608c09542731cb2fc3590b002effc4603 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Thu, 29 Oct 2026 12:58:19 -0700 Subject: [PATCH 11/21] remove flooring of tabulated thresholds --- jdluc/__tests__/datasets/fdap_test.py | 4 ++-- jdluc/datasets/fdap.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/jdluc/__tests__/datasets/fdap_test.py b/jdluc/__tests__/datasets/fdap_test.py index c1b4c34..18c39ac 100644 --- a/jdluc/__tests__/datasets/fdap_test.py +++ b/jdluc/__tests__/datasets/fdap_test.py @@ -95,13 +95,13 @@ def test_get_mask_leaves_out_exclusions_and_non_producers_even_at_1( @pytest.mark.parametrize( ("iso_3166", "year", "threshold"), [ - pytest.param("AAA", 2018, FDAP_FLOOR, id="fitted-below-the-floor"), + pytest.param("AAA", 2018, 0.2, id="fitted-below-the-floor"), pytest.param("AAA", 2020, 0.8, id="fitted-above-the-floor"), pytest.param( "AAA", 2022, fdap.HIGHEST_THRESHOLD, id="the-upper-end-of-the-sweep" ), pytest.param( - "AAA", 2024, FDAP_FLOOR, id="the-lower-end-of-the-sweep-below-the-floor" + "AAA", 2024, fdap.LOWEST_THRESHOLD, id="the-lower-end-of-the-sweep" ), pytest.param( "BBB", 2018, FDAP_FLOOR, id="a-producer-without-fitted-thresholds" diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index d826bcf..47aa07e 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -23,8 +23,8 @@ resolution, without resampling. Neighbouring files overlap by one pixel, which is taken from the file it falls in 2. Threshold: a pixel is in the mask when its probability is above its country's threshold for - that year, or above the commodity's floor if that is higher. A country without thresholds uses - the floor if FAOSTAT says it grows the commodity, and is left out otherwise + that year. A country without thresholds uses the commodity's floor if FAOSTAT says it grows the + commodity, and is left out otherwise 3. Exclude: pixels in the excluded areas are left out 4. Smooth: a pixel missing from the 2020 or 2022 mask is added back when it is in the masks of both the year before and the year after. 2018 and 2024 are not smoothed @@ -95,7 +95,7 @@ def get_mask( if national[iso_3166].intersects(box) ] fitted = [ - (national[iso_3166], max(threshold, floor)) + (national[iso_3166], threshold) for (iso_3166, the_year), threshold in sorted( iso_3166_year_to_threshold.items() ) From d0ca852c887e6ead4a6c39a55b8e2764004cb184 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Thu, 29 Oct 2026 14:30:13 -0700 Subject: [PATCH 12/21] Document the FAOSTAT producer population Remove unecessary assertion --- jdluc/datasets/fdap.py | 15 ++++----------- tools/fdap-producers.sql | 39 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 11 deletions(-) create mode 100644 tools/fdap-producers.sql diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index 47aa07e..9450f03 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -162,7 +162,6 @@ def get_raster_dataset( assert set(iso_3166_year_to_threshold) == { (iso_3166, year) for iso_3166 in iso_3166s for year in YEARS } - assert all(0 < threshold < 1 for threshold in iso_3166_year_to_threshold.values()) assert iso_3166s <= producer_iso_3166s assert set(excluded_admin_ids) == set(worldbank_jurisdictions.AdminLevel) assert ( @@ -349,11 +348,8 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("PER", 2022): 0.91, ("PER", 2024): 0.812, }, - # FAOSTAT's producers: area harvested of cocoa (item 661) in some year from 2018 to - # 2024 and, for a country without thresholds, at least 1% of its cropland (item 6620) in one - # of those years. From the crops and land use bulk files: - # https://bulks-faostat.fao.org/production/Production_Crops_Livestock_E_All_Data_(Normalized).zip - # https://bulks-faostat.fao.org/production/Inputs_LandUse_E_All_Data_(Normalized).zip + # FAOSTAT's producers that devote at least 1% of their cropland to cocoa, plus the countries + # with thresholds: `tools/fdap-producers.sql` with item = '661' producer_iso_3166s={ "BRA", "CIV", @@ -449,11 +445,8 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("VNM", 2022): 0.52, ("VNM", 2024): 0.296, }, - # FAOSTAT's producers: area harvested of coffee (item 656) in some year from 2018 to - # 2024 and, for a country without thresholds, at least 1% of its cropland (item 6620) in one - # of those years. From the crops and land use bulk files: - # https://bulks-faostat.fao.org/production/Production_Crops_Livestock_E_All_Data_(Normalized).zip - # https://bulks-faostat.fao.org/production/Inputs_LandUse_E_All_Data_(Normalized).zip + # FAOSTAT's producers that devote at least 1% of their cropland to coffee, plus the countries + # with thresholds: `tools/fdap-producers.sql` with item = '656' producer_iso_3166s={ "BRA", "CAF", diff --git a/tools/fdap-producers.sql b/tools/fdap-producers.sql new file mode 100644 index 0000000..1a51c4d --- /dev/null +++ b/tools/fdap-producers.sql @@ -0,0 +1,39 @@ +-- FAOSTAT producers of a crop that devote at least 1% of their cropland to it in some year from 2018 +-- to 2024. The producer sets in `jdluc/datasets/fdap.py` are these countries plus the countries with +-- thresholds. Area harvested (element 5312) and cropland (Land Use item 6620, element 5110, in +-- thousands of hectares) come from FAOSTAT's bulk files, and the ISO 3166 alpha-3 code for each M49 +-- code from a public list, since the bulk files carry only M49 codes. +-- +-- Run from the repo root, with the FAOSTAT item code: 661 for cocoa, 656 for coffee: +-- duckdb -cmd "SET VARIABLE item = '661';" -f tools/fdap-producers.sql +INSTALL httpfs; +LOAD httpfs; +INSTALL zipfs FROM community; +LOAD zipfs; + +WITH crop AS ( + SELECT "Area Code (M49)" AS m49, CAST("Year" AS INTEGER) AS year, CAST("Value" AS DOUBLE) AS area_hectares + FROM read_csv( + 'zip://https://bulks-faostat.fao.org/production/Production_Crops_Livestock_E_All_Data_(Normalized).zip/Production_Crops_Livestock_E_All_Data_(Normalized).csv', + all_varchar = true + ) + WHERE "Item Code" = getvariable('item') AND "Element Code" = '5312' AND "Flag" != 'M' AND "Value" != '' +), cropland AS ( + SELECT "Area Code (M49)" AS m49, CAST("Year" AS INTEGER) AS year, CAST("Value" AS DOUBLE) * 1000 AS cropland_hectares + FROM read_csv( + 'zip://https://bulks-faostat.fao.org/production/Inputs_LandUse_E_All_Data_(Normalized).zip/Inputs_LandUse_E_All_Data_(Normalized).csv', + all_varchar = true + ) + WHERE "Item Code" = '6620' AND "Element Code" = '5110' AND "Flag" != 'M' AND "Value" != '' +), iso AS ( + SELECT lpad("country-code", 3, '0') AS numeric, "alpha-3" AS iso_3166 + FROM read_csv('https://raw.githubusercontent.com/lukes/ISO-3166-Countries-with-Regional-Codes/master/all/all.csv', all_varchar = true) +) +SELECT iso_3166 +FROM crop +JOIN cropland USING (m49, year) +JOIN iso ON iso.numeric = ltrim(crop.m49, '''') +WHERE year BETWEEN 2018 AND 2024 AND area_hectares > 0 +GROUP BY iso_3166 +HAVING MAX(area_hectares / cropland_hectares) >= 0.01 +ORDER BY iso_3166; From 0f44a1dcb27b12e58e67c1525c2b3bf923cdc42b Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Thu, 29 Oct 2026 14:50:56 -0700 Subject: [PATCH 13/21] provisional exclusion lists --- jdluc/datasets/fdap.py | 34 ++++++++++++++++++++++++++++++++-- 1 file changed, 32 insertions(+), 2 deletions(-) diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index 9450f03..46306d8 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -314,7 +314,25 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: commodity=Commodity.COCOA, excluded_admin_ids={ worldbank_jurisdictions.AdminLevel.NATIONAL: set(), - worldbank_jurisdictions.AdminLevel.PROVINCIAL: set(), + worldbank_jurisdictions.AdminLevel.PROVINCIAL: { + # Semi-arid, and agronomically unsuitable for cocoa + "BRA002", # Alagoas + "BRA006", # Ceará + "BRA019", # Paraíba + "BRA021", # Pernambuco + "BRA022", # Piauí + "BRA024", # Rio Grande do Norte + # The detected signal is consistent with babaçu palm + "BRA010", # Maranhão + # Young oil palm and rubber + "IDN004", # Riau + "IDN005", # Jambi + "IDN009", # Kepulauan Bangka Belitung + "IDN022", # Kalimantan Selatan + "IDN024", # Kalimantan Utara + # No cocoa production + "ECU020", # Galápagos + }, worldbank_jurisdictions.AdminLevel.DISTRICT: set(), }, floor=0.54, @@ -387,7 +405,19 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: commodity=Commodity.COFFEE, excluded_admin_ids={ worldbank_jurisdictions.AdminLevel.NATIONAL: set(), - worldbank_jurisdictions.AdminLevel.PROVINCIAL: set(), + worldbank_jurisdictions.AdminLevel.PROVINCIAL: { + # Likely confused with other cropland or natural grassland + "BRA010", # Maranhão + "BRA011", # Mato Grosso + "BRA012", # Mato Grosso do Sul + "BRA018", # Pará + "BRA022", # Piauí + "BRA031", # Tocantins + # Confused with tea + "IDN012", # Jawa Barat + # Confused with natural ecosystems + "PER014", # Lambayeque + }, worldbank_jurisdictions.AdminLevel.DISTRICT: set(), }, floor=0.33, From d12f7a378bf2833ee182ecd8d573ccfaf50c1ba9 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Thu, 29 Oct 2026 15:15:19 -0700 Subject: [PATCH 14/21] GDAL setting to parallelize reads --- jdluc/datasets/fdap.py | 1 + 1 file changed, 1 insertion(+) diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index 46306d8..776f021 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -198,6 +198,7 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: rasterio.Env( CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif", GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", + GDAL_NUM_THREADS="ALL_CPUS", GS_USER_PROJECT=get_billing_project(), ), # Blocks no source writes are zeros when the file closes From 31ff96320031b12ea75b13f27b126274fa94f5ee Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Fri, 30 Oct 2026 10:27:22 -0700 Subject: [PATCH 15/21] Re-apply flooring of thresholds Pull threshold floors into constants and document them Various comment and doc polishes --- jdluc/__tests__/datasets/fdap_test.py | 6 +++--- jdluc/datasets/fdap.py | 27 ++++++++++++++++----------- 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/jdluc/__tests__/datasets/fdap_test.py b/jdluc/__tests__/datasets/fdap_test.py index 18c39ac..74873be 100644 --- a/jdluc/__tests__/datasets/fdap_test.py +++ b/jdluc/__tests__/datasets/fdap_test.py @@ -72,7 +72,7 @@ def get_fdap_mask(fill_probability: float, year: int) -> numpy.ndarray: ("AAA", 2018): 0.2, ("AAA", 2020): 0.8, ("AAA", 2022): fdap.HIGHEST_THRESHOLD, - ("AAA", 2024): fdap.LOWEST_THRESHOLD, + ("AAA", 2024): 0.1, }, probability=numpy.full( dtype="float32", fill_value=fill_probability, shape=FDAP_SHAPE @@ -95,13 +95,13 @@ def test_get_mask_leaves_out_exclusions_and_non_producers_even_at_1( @pytest.mark.parametrize( ("iso_3166", "year", "threshold"), [ - pytest.param("AAA", 2018, 0.2, id="fitted-below-the-floor"), + pytest.param("AAA", 2018, FDAP_FLOOR, id="fitted-below-the-floor"), pytest.param("AAA", 2020, 0.8, id="fitted-above-the-floor"), pytest.param( "AAA", 2022, fdap.HIGHEST_THRESHOLD, id="the-upper-end-of-the-sweep" ), pytest.param( - "AAA", 2024, fdap.LOWEST_THRESHOLD, id="the-lower-end-of-the-sweep" + "AAA", 2024, FDAP_FLOOR, id="the-lower-end-of-the-sweep-below-the-floor" ), pytest.param( "BBB", 2018, FDAP_FLOOR, id="a-producer-without-fitted-thresholds" diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index 776f021..cbee5ea 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -5,7 +5,7 @@ year: 2018, 2020, 2022, 2024 Google, for the Forest Data Partnership (2026). Commodity probabilities, release 2026a [Data set]. gs://earth-engine-public-requester-pays/forestdatapartnership/2026a -- the record ingested here, whose licence asks for the attribution "Produced by Google for the Forest Data Partnership" -TODO: add reference for the thresholds +TODO: add reference for the thresholds and exclusions https://github.com/google/forest-data-partnership @@ -23,8 +23,8 @@ resolution, without resampling. Neighbouring files overlap by one pixel, which is taken from the file it falls in 2. Threshold: a pixel is in the mask when its probability is above its country's threshold for - that year. A country without thresholds uses the commodity's floor if FAOSTAT says it grows the - commodity, and is left out otherwise + that year, or above the commodity's floor if that is higher. A country without thresholds uses + the floor if FAOSTAT says it grows the commodity, and is left out otherwise 3. Exclude: pixels in the excluded areas are left out 4. Smooth: a pixel missing from the 2020 or 2022 mask is added back when it is in the masks of both the year before and the year after. 2018 and 2024 are not smoothed @@ -54,7 +54,8 @@ YEARS: typing.Final[tuple[int, ...]] = (2018, 2020, 2022, 2024) # 10 m on the equator of the WGS 84 ellipsoid, whose semi-major axis is 6,378,137 m PIXEL_DEGREES = 10 / (math.pi * 6_378_137 / 180) -LOWEST_THRESHOLD = 0.0025 +# NB: The sweep for thresholds stopped here; a more precise value could be pursued but that would +# likely overfit within a region of the CDF that was quite steep HIGHEST_THRESHOLD = 0.99 NO_DATA = None @@ -95,7 +96,7 @@ def get_mask( if national[iso_3166].intersects(box) ] fitted = [ - (national[iso_3166], threshold) + (national[iso_3166], max(threshold, floor)) for (iso_3166, the_year), threshold in sorted( iso_3166_year_to_threshold.items() ) @@ -311,6 +312,8 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ) +# https://github.com/google/forest-data-partnership/tree/4d0fb3b8ec24f4788b9d46289cb19ad07f466ae0/models/model_2026a#model-selection-set +COCOA_PRECISION_RECALL_CROSSOVER = 0.54 COCOA_DATASET = get_raster_dataset( commodity=Commodity.COCOA, excluded_admin_ids={ @@ -336,7 +339,7 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: }, worldbank_jurisdictions.AdminLevel.DISTRICT: set(), }, - floor=0.54, + floor=COCOA_PRECISION_RECALL_CROSSOVER, iso_3166_year_to_threshold={ ("BRA", 2018): 0.923, ("BRA", 2020): 0.946, @@ -358,10 +361,10 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("GHA", 2020): HIGHEST_THRESHOLD, ("GHA", 2022): HIGHEST_THRESHOLD, ("GHA", 2024): HIGHEST_THRESHOLD, - ("IDN", 2018): 0.540, - ("IDN", 2020): 0.540, - ("IDN", 2022): 0.540, - ("IDN", 2024): 0.540, + ("IDN", 2018): COCOA_PRECISION_RECALL_CROSSOVER, + ("IDN", 2020): COCOA_PRECISION_RECALL_CROSSOVER, + ("IDN", 2022): COCOA_PRECISION_RECALL_CROSSOVER, + ("IDN", 2024): COCOA_PRECISION_RECALL_CROSSOVER, ("PER", 2018): 0.558, ("PER", 2020): 0.866, ("PER", 2022): 0.91, @@ -402,6 +405,8 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: version="v0", ) +# https://github.com/google/forest-data-partnership/tree/4d0fb3b8ec24f4788b9d46289cb19ad07f466ae0/models/model_2026a#model-selection-set +COFFEE_PRECISION_RECALL_CROSSOVER = 0.33 COFFEE_DATASET = get_raster_dataset( commodity=Commodity.COFFEE, excluded_admin_ids={ @@ -421,7 +426,7 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: }, worldbank_jurisdictions.AdminLevel.DISTRICT: set(), }, - floor=0.33, + floor=COFFEE_PRECISION_RECALL_CROSSOVER, iso_3166_year_to_threshold={ ("BRA", 2018): 0.554, ("BRA", 2020): 0.658, From f6a39753145ea29d9cd590422fa5ee6fd87be4f0 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Fri, 30 Oct 2026 10:44:22 -0700 Subject: [PATCH 16/21] Remove the 1% requirement Switch from FAOSTAT bulk to our ingested dataset --- jdluc/datasets/fdap.py | 82 ++++++++++++++++++++++++++++++++++++++-- tools/fdap-producers.sql | 39 ------------------- 2 files changed, 78 insertions(+), 43 deletions(-) delete mode 100644 tools/fdap-producers.sql diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index cbee5ea..a0d0f08 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -370,35 +370,67 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("PER", 2022): 0.91, ("PER", 2024): 0.812, }, - # FAOSTAT's producers that devote at least 1% of their cropland to cocoa, plus the countries - # with thresholds: `tools/fdap-producers.sql` with item = '661' + # FAOSTAT's producers of cocoa, from a local copy of the ingested crop table + # ($INGEST_ROOT/tabular/faostat/production-crops/v2/whole-world/world.parquet): + # duckdb -c "SELECT DISTINCT admin_id FROM 'world.parquet' WHERE commodity_name = 'COCO' + # AND year BETWEEN 2018 AND 2024 AND area_hectares > 0 ORDER BY admin_id" producer_iso_3166s={ + "AGO", + "BLZ", + "BOL", "BRA", + "CAF", "CIV", "CMR", + "COD", "COG", "COL", + "COM", "CRI", + "CUB", "DMA", "DOM", "ECU", + "FJI", "FSM", + "GAB", "GHA", + "GIN", "GNQ", "GRD", + "GTM", + "GUY", + "HND", + "HTI", "IDN", + "IND", + "JAM", "LBR", "LCA", + "LKA", + "MDG", + "MEX", + "MYS", "NGA", + "NIC", "PAN", "PER", + "PHL", "PNG", "SLB", "SLE", + "SLV", "STP", + "SUR", + "TGO", + "THA", + "TLS", "TTO", + "TZA", + "UGA", "VCT", "VEN", + "VNM", "VUT", "WSM", }, @@ -481,48 +513,90 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("VNM", 2022): 0.52, ("VNM", 2024): 0.296, }, - # FAOSTAT's producers that devote at least 1% of their cropland to coffee, plus the countries - # with thresholds: `tools/fdap-producers.sql` with item = '656' + # FAOSTAT's producers of coffee, from a local copy of the ingested crop table + # ($INGEST_ROOT/tabular/faostat/production-crops/v2/whole-world/world.parquet): + # duckdb -c "SELECT DISTINCT admin_id FROM 'world.parquet' WHERE commodity_name = 'COFF' + # AND year BETWEEN 2018 AND 2024 AND area_hectares > 0 ORDER BY admin_id" producer_iso_3166s={ + "AGO", + "BDI", + "BEN", + "BLZ", + "BOL", "BRA", "CAF", + "CHN", "CIV", "CMR", + "COD", "COG", + "COK", "COL", + "COM", + "CPV", "CRI", + "CUB", "DMA", "DOM", "ECU", "ETH", + "FJI", + "GAB", + "GHA", "GIN", "GNQ", "GTM", + "GUY", "HND", + "HTI", "IDN", "IND", "JAM", "KEN", + "KHM", "LAO", + "LBR", + "LKA", "MDG", "MEX", + "MMR", + "MOZ", + "MWI", + "MYS", + "NCL", + "NGA", "NIC", + "NPL", "PAN", "PER", "PHL", "PNG", "PRI", + "PRY", + "PYF", "RWA", + "SAU", + "SLE", "SLV", + "STP", + "SUR", "TGO", + "THA", "TLS", + "TON", "TTO", + "TWN", "TZA", "UGA", + "USA", "VCT", "VEN", "VNM", + "VUT", + "WSM", "YEM", + "ZMB", + "ZWE", }, version="v0", ) diff --git a/tools/fdap-producers.sql b/tools/fdap-producers.sql deleted file mode 100644 index 1a51c4d..0000000 --- a/tools/fdap-producers.sql +++ /dev/null @@ -1,39 +0,0 @@ --- FAOSTAT producers of a crop that devote at least 1% of their cropland to it in some year from 2018 --- to 2024. The producer sets in `jdluc/datasets/fdap.py` are these countries plus the countries with --- thresholds. Area harvested (element 5312) and cropland (Land Use item 6620, element 5110, in --- thousands of hectares) come from FAOSTAT's bulk files, and the ISO 3166 alpha-3 code for each M49 --- code from a public list, since the bulk files carry only M49 codes. --- --- Run from the repo root, with the FAOSTAT item code: 661 for cocoa, 656 for coffee: --- duckdb -cmd "SET VARIABLE item = '661';" -f tools/fdap-producers.sql -INSTALL httpfs; -LOAD httpfs; -INSTALL zipfs FROM community; -LOAD zipfs; - -WITH crop AS ( - SELECT "Area Code (M49)" AS m49, CAST("Year" AS INTEGER) AS year, CAST("Value" AS DOUBLE) AS area_hectares - FROM read_csv( - 'zip://https://bulks-faostat.fao.org/production/Production_Crops_Livestock_E_All_Data_(Normalized).zip/Production_Crops_Livestock_E_All_Data_(Normalized).csv', - all_varchar = true - ) - WHERE "Item Code" = getvariable('item') AND "Element Code" = '5312' AND "Flag" != 'M' AND "Value" != '' -), cropland AS ( - SELECT "Area Code (M49)" AS m49, CAST("Year" AS INTEGER) AS year, CAST("Value" AS DOUBLE) * 1000 AS cropland_hectares - FROM read_csv( - 'zip://https://bulks-faostat.fao.org/production/Inputs_LandUse_E_All_Data_(Normalized).zip/Inputs_LandUse_E_All_Data_(Normalized).csv', - all_varchar = true - ) - WHERE "Item Code" = '6620' AND "Element Code" = '5110' AND "Flag" != 'M' AND "Value" != '' -), iso AS ( - SELECT lpad("country-code", 3, '0') AS numeric, "alpha-3" AS iso_3166 - FROM read_csv('https://raw.githubusercontent.com/lukes/ISO-3166-Countries-with-Regional-Codes/master/all/all.csv', all_varchar = true) -) -SELECT iso_3166 -FROM crop -JOIN cropland USING (m49, year) -JOIN iso ON iso.numeric = ltrim(crop.m49, '''') -WHERE year BETWEEN 2018 AND 2024 AND area_hectares > 0 -GROUP BY iso_3166 -HAVING MAX(area_hectares / cropland_hectares) >= 0.01 -ORDER BY iso_3166; From 760aceab10377ab6ec7a614309efaaeb52a0bb68 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Fri, 30 Oct 2026 11:02:04 -0700 Subject: [PATCH 17/21] Skip a one-degree tile when it doesn't contain any producer countries --- jdluc/__tests__/datasets/fdap_test.py | 20 ++++++++++++++++ jdluc/datasets/fdap.py | 34 +++++++++++++++++++++++++-- 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/jdluc/__tests__/datasets/fdap_test.py b/jdluc/__tests__/datasets/fdap_test.py index 74873be..22bfbe7 100644 --- a/jdluc/__tests__/datasets/fdap_test.py +++ b/jdluc/__tests__/datasets/fdap_test.py @@ -138,3 +138,23 @@ def test_get_smoothed_masks_fills_one_year_gaps_in_the_inner_years( } ) assert tuple(int(smoothed[year][0, 0]) for year in fdap.YEARS) == expected + + +@pytest.mark.parametrize( + ("latitude", "longitude", "read"), + [ + pytest.param(0, 0, True, id="a-cell-inside-a-producer"), + pytest.param(0, 2, True, id="a-cell-touching-a-producer-at-its-edge"), + pytest.param(0, 3, False, id="a-cell-touching-only-a-non-producer"), + pytest.param(5, 5, False, id="a-cell-touching-no-country"), + ], +) +def test_touches_a_producer(latitude: int, longitude: int, read: bool) -> None: + assert ( + fdap.touches_a_producer( + lat_lon=(latitude, longitude), + national=FDAP_NATIONAL, + producer_iso_3166s={"AAA", "BBB"}, + ) + is read + ) diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index a0d0f08..633776f 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -31,6 +31,7 @@ """ import enum +import functools import logging import math import typing @@ -136,6 +137,19 @@ def get_smoothed_masks( } +def touches_a_producer( + lat_lon: tuple[int, int], + national: geopandas.GeoSeries, + producer_iso_3166s: set[str], +) -> bool: + """Whether the one-degree cell touches a producer; every other cell's masks are 0.""" + latitude, longitude = lat_lon + box = shapely.box(longitude, latitude, longitude + 1, latitude + 1) + return any( + national[iso_3166].intersects(box) for iso_3166 in sorted(producer_iso_3166s) + ) + + @utils.threadsafe_cache def get_billing_project() -> str: _, project = google.auth.default() @@ -191,10 +205,26 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: & (index["lat"] >= max_latitude - 10) & (index["lat"] < max_latitude) ] + all_lat_lons = set( + zip(sources["lat"].tolist(), sources["lng"].tolist(), strict=True) + ) + # The rest read 0 in every year, which the zero-filled file already holds lat_lons = sorted( - set(zip(sources["lat"].tolist(), sources["lng"].tolist(), strict=True)) + filter( + functools.partial( + touches_a_producer, + national=admin_level_to_geometries[ + worldbank_jurisdictions.AdminLevel.NATIONAL + ], + producer_iso_3166s=producer_iso_3166s, + ), + all_lat_lons, + ) + ) + logger.info( + f"Thresholding {len(lat_lons):d} of {len(all_lat_lons):d} source cells, the ones " + f"touching a producer, onto {tile_id=:s}" ) - logger.info(f"Thresholding {len(lat_lons):d} source cells onto {tile_id=:s}") with ( rasterio.Env( CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif", From 6a4042619aa1d2e16c635cfdaa5faa71c2802d63 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Fri, 30 Oct 2026 11:05:34 -0700 Subject: [PATCH 18/21] unify terminology --- jdluc/__tests__/datasets/fdap_test.py | 8 ++++---- jdluc/datasets/fdap.py | 13 +++++++------ 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/jdluc/__tests__/datasets/fdap_test.py b/jdluc/__tests__/datasets/fdap_test.py index 22bfbe7..b731ded 100644 --- a/jdluc/__tests__/datasets/fdap_test.py +++ b/jdluc/__tests__/datasets/fdap_test.py @@ -38,7 +38,7 @@ def test_get_pixel_range_reads_each_source_inside_its_file(start_degrees: int) - assert pixels.stop <= first_pixel + 11_133 -# Three countries side by side, one degree each: AAA has fitted thresholds, BBB is a producer +# Three countries side by side, one degree each: AAA has country thresholds, BBB is a producer # without them, and CCC is not a producer. AAA's western half is an excluded province. FDAP_NATIONAL = geopandas.GeoSeries( crs=4326, @@ -95,8 +95,8 @@ def test_get_mask_leaves_out_exclusions_and_non_producers_even_at_1( @pytest.mark.parametrize( ("iso_3166", "year", "threshold"), [ - pytest.param("AAA", 2018, FDAP_FLOOR, id="fitted-below-the-floor"), - pytest.param("AAA", 2020, 0.8, id="fitted-above-the-floor"), + pytest.param("AAA", 2018, FDAP_FLOOR, id="a-country-threshold-below-the-floor"), + pytest.param("AAA", 2020, 0.8, id="a-country-threshold-above-the-floor"), pytest.param( "AAA", 2022, fdap.HIGHEST_THRESHOLD, id="the-upper-end-of-the-sweep" ), @@ -104,7 +104,7 @@ def test_get_mask_leaves_out_exclusions_and_non_producers_even_at_1( "AAA", 2024, FDAP_FLOOR, id="the-lower-end-of-the-sweep-below-the-floor" ), pytest.param( - "BBB", 2018, FDAP_FLOOR, id="a-producer-without-fitted-thresholds" + "BBB", 2018, FDAP_FLOOR, id="a-producer-without-country-thresholds" ), ], ) diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index 633776f..deb24d6 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -23,9 +23,10 @@ resolution, without resampling. Neighbouring files overlap by one pixel, which is taken from the file it falls in 2. Threshold: a pixel is in the mask when its probability is above its country's threshold for - that year, or above the commodity's floor if that is higher. A country without thresholds uses - the floor if FAOSTAT says it grows the commodity, and is left out otherwise -3. Exclude: pixels in the excluded areas are left out + that year, or above the commodity's floor, the model's precision-recall crossover, if that is + higher. A country without thresholds uses the floor if FAOSTAT says it grows the commodity, and + is left out otherwise +3. Exclude: pixels in the excluded jurisdictions are left out 4. Smooth: a pixel missing from the 2020 or 2022 mask is added back when it is in the masks of both the year before and the year after. 2018 and 2024 are not smoothed """ @@ -91,12 +92,12 @@ def get_mask( *rasterio.transform.array_bounds(*probability.shape, transform=transform) ) national = admin_level_to_geometries[worldbank_jurisdictions.AdminLevel.NATIONAL] - defaults = [ + floors = [ (national[iso_3166], floor) for iso_3166 in sorted(producer_iso_3166s) if national[iso_3166].intersects(box) ] - fitted = [ + country_thresholds = [ (national[iso_3166], max(threshold, floor)) for (iso_3166, the_year), threshold in sorted( iso_3166_year_to_threshold.items() @@ -114,7 +115,7 @@ def get_mask( dtype="float32", fill=1, out_shape=probability.shape, - shapes=defaults + fitted + exclusions, + shapes=floors + country_thresholds + exclusions, transform=transform, ) return (probability > threshold).astype("uint8") From 4d2d2134f2f079e4abafa9e0b848f5479f1c1ef0 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Fri, 30 Oct 2026 18:54:40 -0700 Subject: [PATCH 19/21] fit FDaP thresholds for the producing countries where the floor would over-estimate their total area harvested --- jdluc/datasets/fdap.py | 121 ++++++++++++++- tools/fit-fdap-thresholds.py | 277 +++++++++++++++++++++++++++++++++++ 2 files changed, 397 insertions(+), 1 deletion(-) create mode 100644 tools/fit-fdap-thresholds.py diff --git a/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py index deb24d6..71154a0 100644 --- a/jdluc/datasets/fdap.py +++ b/jdluc/datasets/fdap.py @@ -372,6 +372,7 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: }, floor=COCOA_PRECISION_RECALL_CROSSOVER, iso_3166_year_to_threshold={ + # Fitted by Ad Astra ("BRA", 2018): 0.923, ("BRA", 2020): 0.946, ("BRA", 2022): 0.976, @@ -400,6 +401,64 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("PER", 2020): 0.866, ("PER", 2022): 0.91, ("PER", 2024): 0.812, + } + | { + # Fitted in this repo by `tools/fit-fdap-thresholds.py`, to FAOSTAT's area harvested as + # published now: producers where the floor maps several times FAOSTAT's area, and Côte + # d'Ivoire 2024, which Ad Astra fitted to an earlier FAOSTAT figure + ("CIV", 2024): 0.981, + ("COD", 2018): 0.63, + ("COD", 2020): 0.809, + ("COD", 2022): 0.912, + ("COD", 2024): 0.911, + ("FJI", 2018): 0.983, + ("FJI", 2020): HIGHEST_THRESHOLD, + ("FJI", 2022): HIGHEST_THRESHOLD, + ("FJI", 2024): 0.982, + ("GAB", 2018): 0.936, + ("GAB", 2020): 0.908, + ("GAB", 2022): 0.923, + ("GAB", 2024): 0.878, + ("GIN", 2018): 0.93, + ("GIN", 2020): 0.936, + ("GIN", 2022): 0.911, + ("GIN", 2024): 0.866, + ("GUY", 2018): 0.98, + ("GUY", 2020): 0.98, + ("GUY", 2022): HIGHEST_THRESHOLD, + ("GUY", 2024): HIGHEST_THRESHOLD, + ("LBR", 2018): HIGHEST_THRESHOLD, + ("LBR", 2020): 0.988, + ("LBR", 2022): 0.986, + ("LBR", 2024): 0.973, + ("PAN", 2018): 0.944, + ("PAN", 2020): 0.888, + ("PAN", 2022): 0.968, + ("PAN", 2024): 0.911, + ("SLE", 2018): 0.786, + ("SLE", 2020): 0.827, + ("SLE", 2022): 0.901, + ("SLE", 2024): 0.913, + ("SUR", 2018): HIGHEST_THRESHOLD, + ("SUR", 2020): HIGHEST_THRESHOLD, + ("SUR", 2022): HIGHEST_THRESHOLD, + ("SUR", 2024): HIGHEST_THRESHOLD, + ("TGO", 2018): 0.954, + ("TGO", 2020): 0.972, + ("TGO", 2022): 0.937, + ("TGO", 2024): 0.969, + ("VEN", 2018): 0.852, + ("VEN", 2020): 0.466, + ("VEN", 2022): 0.966, + ("VEN", 2024): 0.835, + ("VUT", 2018): 0.967, + ("VUT", 2020): 0.974, + ("VUT", 2022): 0.968, + ("VUT", 2024): 0.863, + ("WSM", 2018): 0.793, + ("WSM", 2020): 0.857, + ("WSM", 2022): 0.861, + ("WSM", 2024): 0.866, }, # FAOSTAT's producers of cocoa, from a local copy of the ingested crop table # ($INGEST_ROOT/tabular/faostat/production-crops/v2/whole-world/world.parquet): @@ -473,7 +532,10 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: COFFEE_DATASET = get_raster_dataset( commodity=Commodity.COFFEE, excluded_admin_ids={ - worldbank_jurisdictions.AdminLevel.NATIONAL: set(), + worldbank_jurisdictions.AdminLevel.NATIONAL: { + # Over 100x FAOSTAT's area even at `HIGHEST_THRESHOLD` + "MOZ", # Mozambique + }, worldbank_jurisdictions.AdminLevel.PROVINCIAL: { # Likely confused with other cropland or natural grassland "BRA010", # Maranhão @@ -491,6 +553,7 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: }, floor=COFFEE_PRECISION_RECALL_CROSSOVER, iso_3166_year_to_threshold={ + # Fitted by Ad Astra ("BRA", 2018): 0.554, ("BRA", 2020): 0.658, ("BRA", 2022): 0.763, @@ -543,6 +606,62 @@ def save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: ("VNM", 2020): 0.049, ("VNM", 2022): 0.52, ("VNM", 2024): 0.296, + } + | { + # Fitted in this repo by `tools/fit-fdap-thresholds.py`, to FAOSTAT's area harvested as + # published now, for producers where the floor maps several times FAOSTAT's area + ("AGO", 2018): 0.549, + ("AGO", 2020): 0.535, + ("AGO", 2022): 0.979, + ("AGO", 2024): 0.923, + ("BDI", 2018): HIGHEST_THRESHOLD, + ("BDI", 2020): 0.945, + ("BDI", 2022): HIGHEST_THRESHOLD, + ("BDI", 2024): HIGHEST_THRESHOLD, + ("BOL", 2018): 0.931, + ("BOL", 2020): 0.977, + ("BOL", 2022): 0.887, + ("BOL", 2024): 0.851, + ("COD", 2018): 0.921, + ("COD", 2020): 0.769, + ("COD", 2022): 0.974, + ("COD", 2024): 0.96, + ("ECU", 2018): 0.383, + ("ECU", 2020): 0.513, + ("ECU", 2022): 0.92, + ("ECU", 2024): 0.816, + ("HTI", 2018): 0.863, + ("HTI", 2020): 0.926, + ("HTI", 2022): 0.939, + ("HTI", 2024): 0.911, + ("KEN", 2018): 0.472, + ("KEN", 2020): 0.668, + ("KEN", 2022): HIGHEST_THRESHOLD, + ("KEN", 2024): 0.989, + ("MDG", 2018): 0.926, + ("MDG", 2020): 0.918, + ("MDG", 2022): 0.912, + ("MDG", 2024): 0.918, + ("MWI", 2018): HIGHEST_THRESHOLD, + ("MWI", 2020): HIGHEST_THRESHOLD, + ("MWI", 2022): HIGHEST_THRESHOLD, + ("MWI", 2024): HIGHEST_THRESHOLD, + ("PRY", 2018): HIGHEST_THRESHOLD, + ("PRY", 2020): HIGHEST_THRESHOLD, + ("PRY", 2022): HIGHEST_THRESHOLD, + ("PRY", 2024): HIGHEST_THRESHOLD, + ("RWA", 2018): 0.974, + ("RWA", 2020): 0.737, + ("RWA", 2022): 0.983, + ("RWA", 2024): 0.964, + ("TZA", 2018): 0.93, + ("TZA", 2020): 0.832, + ("TZA", 2022): 0.92, + ("TZA", 2024): 0.868, + ("ZWE", 2018): HIGHEST_THRESHOLD, + ("ZWE", 2020): HIGHEST_THRESHOLD, + ("ZWE", 2022): HIGHEST_THRESHOLD, + ("ZWE", 2024): HIGHEST_THRESHOLD, }, # FAOSTAT's producers of coffee, from a local copy of the ingested crop table # ($INGEST_ROOT/tabular/faostat/production-crops/v2/whole-world/world.parquet): diff --git a/tools/fit-fdap-thresholds.py b/tools/fit-fdap-thresholds.py new file mode 100644 index 0000000..a483459 --- /dev/null +++ b/tools/fit-fdap-thresholds.py @@ -0,0 +1,277 @@ +"""Fit FDaP thresholds for chosen countries to FAOSTAT area harvested. + +For every source file touching the countries, read the pixels whose centres fall in the file's own +degree, as ingest does, and add their hectares to one of 1,000 probability bins for the country and +year they fall in. Summed from the top, the bins give the area above every threshold. Each year's +threshold is where that area meets FAOSTAT's area harvested, interpolating log area against log +threshold between `LOWEST_THRESHOLD` and `fdap.HIGHEST_THRESHOLD`. Where the area can't reach +FAOSTAT's within that range, the threshold is the end of the range, marked OVER (the model maps +more than FAOSTAT even at the top) or UNDER. + +Prints each country and year's FAOSTAT area, its area at the commodity's floor and the fitted +threshold, and writes every curve to `--curves.csv` in the working directory. + +Run from the repo root. The source bucket is requester-pays and in us-central1, so run it on a VM +there; a country's files are a few hundred megabytes each: + + uv run python tools/fit-fdap-thresholds.py cocoa GIN SLE TGO + uv run python tools/fit-fdap-thresholds.py coffee BDI RWA --threads 16 +""" + +import argparse +import collections +import concurrent.futures +import logging +import math + +import geopandas +import numpy +import pandas +import rasterio +import rasterio.features +import rasterio.windows +import shapely +import xarray + +from jdluc import emit, storage +from jdluc.datasets import faostat, fdap, ifpri_mapspam, worldbank_jurisdictions + +logger = logging.getLogger(__name__) + +BINS = 1000 +# The bottom of the threshold sweep; its top is `fdap.HIGHEST_THRESHOLD` +LOWEST_THRESHOLD = 0.0025 +COMMODITY_TO_FLOOR = { + fdap.Commodity.COCOA: fdap.COCOA_PRECISION_RECALL_CROSSOVER, + fdap.Commodity.COFFEE: fdap.COFFEE_PRECISION_RECALL_CROSSOVER, +} +COMMODITY_TO_FAOSTAT_NAME = { + fdap.Commodity.COCOA: faostat.Crop.COCO.name, + fdap.Commodity.COFFEE: ifpri_mapspam.Crop2000.COFF.name, +} +# Ad Astra fitted the 2024 band to FAOSTAT 2023, before FAOSTAT published 2024, and so does this +YEAR_TO_FAOSTAT_YEAR = {2018: 2018, 2020: 2020, 2022: 2022, 2024: 2023} + + +def get_iso_3166_to_hectares_per_bin( + iso_3166s: list[str], + latitude: int, + longitude: int, + national: geopandas.GeoSeries, + uri: str, +) -> dict[str, numpy.ndarray]: + """Each country's hectares in each probability bin, within one source file's own degree.""" + columns = fdap.get_pixel_range(start_degrees=longitude, stop_degrees=longitude + 1) + rows = fdap.get_pixel_range(start_degrees=latitude, stop_degrees=latitude + 1) + with ( + rasterio.Env( + CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif", + GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", + GDAL_NUM_THREADS="ALL_CPUS", + GS_USER_PROJECT=fdap.get_billing_project(), + ), + rasterio.open(fp=storage.to_gdal_path(uri=uri)) as source, + ): + left = round(source.transform.c / fdap.PIXEL_DEGREES) + top = round(source.transform.f / fdap.PIXEL_DEGREES) + window = rasterio.windows.Window( + col_off=columns.start - left, + height=len(rows), + row_off=top - rows.stop, + width=len(columns), + ) + # NaN pixels compare false at ingest, so they count as 0 here too + probability = numpy.nan_to_num(source.read(1, window=window), nan=0.0) + transform = rasterio.windows.transform( + transform=source.transform, window=window + ) + labels = rasterio.features.rasterize( + dtype="uint8", + fill=0, + out_shape=probability.shape, + shapes=[ + (national[iso_3166], idx) for idx, iso_3166 in enumerate(iso_3166s, start=1) + ], + transform=transform, + ) + # Two columns, since the pixel width comes from the spacing of the x coordinates + row_hectares = emit.get_hectares_per_pixel( + darray=xarray.DataArray( + coords={ + "x": transform.c + (numpy.arange(2) + 0.5) * transform.a, + "y": transform.f + + (numpy.arange(probability.shape[0]) + 0.5) * transform.e, + }, + data=numpy.zeros(dtype="float32", shape=(probability.shape[0], 2)), + dims=("y", "x"), + ) + ).values[:, 0] + hectares_per_bin = numpy.zeros((len(iso_3166s) + 1) * BINS) + # In strips of rows, so the 64-bit bin and weight arrays stay small + for start in range(0, probability.shape[0], 1024): + strip = slice(start, start + 1024) + bins = numpy.clip((probability[strip] * BINS).astype(numpy.int64), 0, BINS - 1) + hectares_per_bin += numpy.bincount( + (labels[strip].astype(numpy.int64) * BINS + bins).ravel(), + minlength=hectares_per_bin.size, + weights=numpy.broadcast_to(row_hectares[strip, None], bins.shape).ravel(), + ) + return { + iso_3166: hectares_per_bin[idx * BINS : (idx + 1) * BINS] + for idx, iso_3166 in enumerate(iso_3166s, start=1) + } + + +def main() -> int: + logging.basicConfig( + format="%(asctime)s - %(levelname)s - %(message)s", level=logging.INFO + ) + parser = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) + parser.add_argument( + "commodity", choices=sorted(commodity.value for commodity in fdap.Commodity) + ) + parser.add_argument( + "iso_3166s", nargs="+", type=worldbank_jurisdictions.iso_3166_str + ) + parser.add_argument("--threads", default=8, type=int) + args = parser.parse_args() + commodity = fdap.Commodity(args.commodity) + + national = worldbank_jurisdictions.get_jurisdiction_for_admin_level( + admin_level=worldbank_jurisdictions.AdminLevel.NATIONAL + ).geometry + index = fdap.get_index() + sources = index[index["commodity"] == commodity] + lat_lon_to_iso_3166s = { + (latitude, longitude): [ + iso_3166 + for iso_3166 in args.iso_3166s + if national[iso_3166].intersects( + shapely.box(longitude, latitude, longitude + 1, latitude + 1) + ) + ] + for latitude, longitude in set( + zip(sources["lat"].tolist(), sources["lng"].tolist(), strict=True) + ) + } + tasks = [ + (latitude, longitude, year, uri) + for latitude, longitude, year, uri in zip( + sources["lat"].tolist(), + sources["lng"].tolist(), + sources["year"].tolist(), + sources["path"].tolist(), + strict=True, + ) + if lat_lon_to_iso_3166s[(latitude, longitude)] + ] + logger.info(f"Reading {len(tasks):d} source files for {args.iso_3166s}") + + iso_3166_year_to_hectares_per_bin: collections.defaultdict[ + tuple[str, int], numpy.ndarray + ] = collections.defaultdict(lambda: numpy.zeros(BINS)) + with concurrent.futures.ThreadPoolExecutor(max_workers=args.threads) as executor: + future_to_year = { + executor.submit( + get_iso_3166_to_hectares_per_bin, + iso_3166s=lat_lon_to_iso_3166s[(latitude, longitude)], + latitude=latitude, + longitude=longitude, + national=national, + uri=uri, + ): year + for latitude, longitude, year, uri in tasks + } + for n, future in enumerate( + concurrent.futures.as_completed(future_to_year), start=1 + ): + for iso_3166, hectares_per_bin in future.result().items(): + iso_3166_year_to_hectares_per_bin[ + (iso_3166, future_to_year[future]) + ] += hectares_per_bin + if n % 20 == 0: + logger.info(f"Read {n:d} of {len(tasks):d} source files") + + crops = faostat.load(dataset=faostat.CROP_DATASET).reset_index() + area_harvested = crops[ + crops["commodity_name"] == COMMODITY_TO_FAOSTAT_NAME[commodity] + ].set_index(["admin_id", "year"])["area_hectares"] + thresholds = numpy.arange(BINS) / BINS + in_range = (thresholds >= LOWEST_THRESHOLD) & (thresholds <= fdap.HIGHEST_THRESHOLD) + fits: list[dict[str, float | int | str]] = [] + curves: list[pandas.DataFrame] = [] + for (iso_3166, year), hectares_per_bin in sorted( + iso_3166_year_to_hectares_per_bin.items() + ): + # Each bin's lower edge, as a threshold, keeps that bin and every one above it + hectares_above = hectares_per_bin[::-1].cumsum()[::-1] + faostat_year = YEAR_TO_FAOSTAT_YEAR[year] + faostat_hectares = ( + float(area_harvested[(iso_3166, faostat_year)]) + if (iso_3166, faostat_year) in area_harvested.index + else 0.0 + ) + candidates, areas = thresholds[in_range], hectares_above[in_range] + if faostat_hectares == 0: + threshold, status = math.nan, "no FAOSTAT area" + elif areas[-1] > faostat_hectares: + threshold, status = fdap.HIGHEST_THRESHOLD, "OVER" + elif areas[0] < faostat_hectares: + threshold, status = LOWEST_THRESHOLD, "UNDER" + else: + # Area falls as the threshold rises, so interpolate over the reversed, positive part + positive = areas > 0 + threshold = float( + numpy.exp( + numpy.interp( + numpy.log(faostat_hectares), + numpy.log(areas[positive][::-1]), + numpy.log(candidates[positive][::-1]), + ) + ) + ) + status = "fitted" + floor_hectares = float( + hectares_above[round(COMMODITY_TO_FLOOR[commodity] * BINS)] + ) + fits.append( + { + "iso_3166": iso_3166, + "year": year, + "faostat_year": faostat_year, + "faostat_kha": faostat_hectares / 1e3, + "floor_kha": floor_hectares / 1e3, + "floor_ratio": ( + floor_hectares / faostat_hectares + if faostat_hectares > 0 + else math.nan + ), + "threshold": threshold, + "status": status, + } + ) + curves.append( + pandas.DataFrame( + { + "iso_3166": iso_3166, + "year": year, + "threshold": thresholds, + "hectares_above": hectares_above, + } + ) + ) + pandas.concat(curves).to_csv( + f"{commodity}-{'-'.join(args.iso_3166s)}-curves.csv", index=False + ) + print( + pandas.DataFrame(fits).to_string( + float_format=lambda v: f"{v:,.3f}", index=False + ) + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) From 937880bc0044658dc34c3bdf1d9f6709340108f5 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Fri, 30 Oct 2026 19:02:12 -0700 Subject: [PATCH 20/21] update docs --- docs/further_research.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/docs/further_research.md b/docs/further_research.md index 38d3ecd..2997774 100644 --- a/docs/further_research.md +++ b/docs/further_research.md @@ -31,6 +31,18 @@ Sampled over 1°×1° windows at 30 m on GLAD's land-cover map, not on GACED30, **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 **Issue**: GACED30 follows the FAO *Cropland* definition, which excludes land under shifting cultivation once its fallow exceeds four years. Forest cleared for a subsistence garden is therefore real conversion that no layer in the stack calls cropland: tree-cover loss dates it, and GACED30, Global Pasture Watch and the per-crop masks all decline it. Its carbon lands in `dropped-emissions` and is charged to nobody. This is the second group the woody-perennial entry above separates out — cassava, other roots, yams, sweet potato and other vegetables, together 43.0 Mha unbacked — and it is a different problem with a different fix: the crops are mapped, the definition excludes them. From cce2a1f0c4cfcf16b849e5fb7071dab66c613c39 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Fri, 30 Oct 2026 19:03:45 -0700 Subject: [PATCH 21/21] update faostat.parquet lock --- validation/data/sources.lock.json | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/validation/data/sources.lock.json b/validation/data/sources.lock.json index 34ddb5f..88b8319 100644 --- a/validation/data/sources.lock.json +++ b/validation/data/sources.lock.json @@ -1,10 +1,10 @@ { "faostat": { - "faostat_production.parquet": { - "bytes": 1377894, - "origin": "gs://cornerstone-ingest-us-central1/tabular/faostat/production-crops/v1/whole-world/world.parquet", - "path": "faostat_production.parquet", - "sha256": "bb608eb3af8a9bda1da884f0ec59d53813fb17461a0d867509c5fdcbcf0562fd" + "faostat.parquet": { + "bytes": 1416755, + "origin": "tabular/faostat/production-crops/v2/whole-world/world.parquet", + "path": "faostat.parquet", + "sha256": "22d11e597a4b3d611cedc88f75834bb546d6d97dd7b0793c111cc3de01bdae67" } }, "orbae": {