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..2997774 100644 --- a/docs/further_research.md +++ b/docs/further_research.md @@ -29,7 +29,19 @@ Sampled over 1°×1° windows at 30 m on GLAD's land-cover map, not on GACED30, **Potential impact:** Potentially the largest single gap for the global extension. The tropical acacia and oil palm emission factors (78.6 and 45.0 t CO₂e ha⁻¹ yr⁻¹) motivate the front-loaded shape of the entire peatland model, yet the pipeline can attribute close to zero peat emissions to those crops. Note that the usual reassurance — that per-crop factors are protected because the numerator only counts pixels that are cropland in 2020 — does not hold for woody perennials, where that same restriction is the mechanism of the undercount. The question is not whether the forest layer *sees* the clearing, which tree-cover loss does; it is that the clearing’s destination is never treated as converted. -**Potential improvement path:** A dedicated perennial-plantation layer is the route in, and the repo ingests one: Descals et al. (2024) global oil palm, whose 30 m pixels carry the planting year as well as the extent, so the peat pulse can be dated directly rather than inferred from a five-year GLAD span. Two properties bound what it can settle. Its planting year dates the plantation standing in 2021, so a replant reads as the replant year rather than as the year the land was first cleared — it says when the current crop was established, not when conversion happened, and pairing it with TCL or a GLAD span is what separates the two. And its accuracy is asymmetric: producer's and user's accuracy are 91.0 ± 2.5% and 91.8 ± 1.2% on industrial plantations against 71.4 ± 0.7% and 72.4 ± 1.8% on smallholders, so a smallholder frontier is the weak end of the layer as well as of GLAD. That has now been done for oil palm: the layer is admitted to the emissions numerator alongside the general cropland mask, and mapped plantation on peat is drained for the occupation term. It is read as extent at the assessment year and dates nothing, precisely because of the replant property above. One consequence is accepted rather than avoided: where the planting year precedes the tree-cover loss, the pixel still resolves to oil palm, so plantation-cycle canopy loss inside an existing estate is charged as conversion — 5.3-8.6% of the conversions the layer admits on Riau. Vetoing those would be dating with the layer, which it is not trusted to do. It also covers only oil palm, which is 352 Mt of the 1,071 Mt in the gap above; the remaining 719 Mt sits on crops with no dedicated extent layer here, cocoa, coconut and rubber among them, and coffee and tea with no open global map at all. That list is woody crops only, which is narrower than "perennial": sugarcane, banana and plantain are perennial but herbaceous, so GLAD's definition admits them. The classifier does not honor that admission where the canopy is dense — the screen above puts sugarcane at 64.2% coverage, banana at 33.9% and plantain at 28.6% — which brings those 71 (country, crop) pairs into scope alongside the woody crops. A general woody-cropland upgrade is also available: GLC_FCS30D (Zhang et al., 2024) carries a separable tree/shrub orchard cropland class, annual for 2000–2022, which would admit woody cropland to the numerator under a different carbon treatment than annual cropland — something GACED30's binary cropland class cannot support. Establish its per-class accuracy for that class first, since the published global overall accuracy of 80.88% spans 35 classes and says nothing about it, along with whether the 25.7% palm recall holds outside Riau and whether the 24% of that window it calls shrubland is immature palm. Adopting it also forces a regrid: its 30 m lattice does not nest into GLAD's 0.00025°, on which TCL, SDPT v2 and Global Pasture Watch are all native. +**Potential improvement path:** A dedicated perennial-plantation layer is the route in, and the repo ingests one: Descals et al. (2024) global oil palm, whose 30 m pixels carry the planting year as well as the extent, so the peat pulse can be dated directly rather than inferred from a five-year GLAD span. Two properties bound what it can settle. Its planting year dates the plantation standing in 2021, so a replant reads as the replant year rather than as the year the land was first cleared — it says when the current crop was established, not when conversion happened, and pairing it with TCL or a GLAD span is what separates the two. And its accuracy is asymmetric: producer's and user's accuracy are 91.0 ± 2.5% and 91.8 ± 1.2% on industrial plantations against 71.4 ± 0.7% and 72.4 ± 1.8% on smallholders, so a smallholder frontier is the weak end of the layer as well as of GLAD. That has now been done for oil palm: the layer is admitted to the emissions numerator alongside the general cropland mask, and mapped plantation on peat is drained for the occupation term. It is read as extent at the assessment year and dates nothing, precisely because of the replant property above. One consequence is accepted rather than avoided: where the planting year precedes the tree-cover loss, the pixel still resolves to oil palm, so plantation-cycle canopy loss inside an existing estate is charged as conversion — 5.3-8.6% of the conversions the layer admits on Riau. Vetoing those would be dating with the layer, which it is not trusted to do. It also covers only oil palm, which is 352 Mt of the 1,071 Mt in the gap above; the remaining 719 Mt sits on crops with no dedicated extent layer the emissions core reads, cocoa, coffee, coconut and rubber among them. For cocoa and coffee the repo ingests one, the Forest Data Partnership's 10 m masks for 2018–2024, which no stage reads yet; tea has no open global map at all. That list is woody crops only, which is narrower than "perennial": sugarcane, banana and plantain are perennial but herbaceous, so GLAD's definition admits them. The classifier does not honor that admission where the canopy is dense — the screen above puts sugarcane at 64.2% coverage, banana at 33.9% and plantain at 28.6% — which brings those 71 (country, crop) pairs into scope alongside the woody crops. A general woody-cropland upgrade is also available: GLC_FCS30D (Zhang et al., 2024) carries a separable tree/shrub orchard cropland class, annual for 2000–2022, which would admit woody cropland to the numerator under a different carbon treatment than annual cropland — something GACED30's binary cropland class cannot support. Establish its per-class accuracy for that class first, since the published global overall accuracy of 80.88% spans 35 classes and says nothing about it, along with whether the 25.7% palm recall holds outside Riau and whether the 24% of that window it calls shrubland is immature palm. Adopting it also forces a regrid: its 30 m lattice does not nest into GLAD's 0.00025°, on which TCL, SDPT v2 and Global Pasture Watch are all native. + +## FDaP masks match FAOSTAT's national area, not where it lies + +**Status**: open. No stage reads the masks yet, so this bounds what they can carry once one does. + +The masks threshold FDaP's per-pixel probabilities per country and year (`jdluc/datasets/fdap.py`). Where a country has a fitted threshold, it is set so the masked area meets FAOSTAT's area harvested; that matches the national total and says nothing about which pixels carry it. Other producers use the model's precision/recall crossover, so their area is not matched to FAOSTAT's at all, and countries with no FAOSTAT production mask nothing. The 2024 masks are fitted to FAOSTAT's 2023 area. + +Where even the highest threshold maps more than FAOSTAT's area harvested, the mask stays larger than FAOSTAT's figure: cocoa in Ghana (2020–2024), Liberia (2018), Fiji (2020 and 2022), Guyana (2022 and 2024) and Suriname, and coffee in Burundi (2018, 2022 and 2024), Kenya (2022), Malawi, Paraguay and Zimbabwe. Mozambique's coffee is excluded nationally instead, since its mask stays over 100 times FAOSTAT's area even at that threshold. + +**Potential impact:** Once the masks are a destination, an over-mapped country attributes more of its clearance to cocoa or coffee than FAOSTAT's area supports, and a producer whose mask maps less than FAOSTAT's area leaves the rest in `dropped-emissions`. + +**Potential improvement path:** `tools/fit-fdap-thresholds.py` refits any country against the current FAOSTAT release: the 2024 masks against FAOSTAT's 2024 area, and the thresholds fitted to earlier FAOSTAT figures. Where the pixels lie needs an independent check against field plots, which a fit to national totals cannot provide. ## Shifting cultivation resolves to no destination 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. 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/fdap_test.py b/jdluc/__tests__/datasets/fdap_test.py new file mode 100644 index 0000000..b731ded --- /dev/null +++ b/jdluc/__tests__/datasets/fdap_test.py @@ -0,0 +1,160 @@ +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 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, + 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): 0.1, + }, + 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="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" + ), + 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-country-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 + + +@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/__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/__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/__init__.py b/jdluc/datasets/__init__.py index cae9ac4..cb0e48c 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,8 @@ def _generate_next_value_( FAOSTAT_LIVESTOCK_PATTERNS = enum.auto() 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() @@ -65,6 +68,8 @@ 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.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/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/jdluc/datasets/fdap.py b/jdluc/datasets/fdap.py new file mode 100644 index 0000000..71154a0 --- /dev/null +++ b/jdluc/datasets/fdap.py @@ -0,0 +1,752 @@ +"""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 and exclusions + +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, 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 +""" + +import enum +import functools +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) +# 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 + + +class Commodity(enum.StrEnum): + COCOA = enum.auto() + COFFEE = 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] + floors = [ + (national[iso_3166], floor) + for iso_3166 in sorted(producer_iso_3166s) + if national[iso_3166].intersects(box) + ] + country_thresholds = [ + (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=floors + country_thresholds + 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) + } + + +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() + 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 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) + ] + 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( + 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}" + ) + with ( + 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 + 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, + photometric="MINISBLACK", + 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 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 + ) + 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, + ) + + +# 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={ + worldbank_jurisdictions.AdminLevel.NATIONAL: 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=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, + ("BRA", 2024): 0.976, + ("CIV", 2018): 0.972, + ("CIV", 2020): 0.982, + ("CIV", 2022): 0.969, + ("CIV", 2024): 0.961, + ("COL", 2018): 0.919, + ("COL", 2020): 0.851, + ("COL", 2022): 0.986, + ("COL", 2024): 0.971, + ("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): 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, + ("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): + # 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", + }, + 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={ + 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 + "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=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, + ("BRA", 2024): 0.593, + ("CIV", 2018): 0.002, + ("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, + ("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, + ("SLV", 2018): 0.005, + ("SLV", 2020): 0.002, + ("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, + } + | { + # 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): + # 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/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] 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( 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 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()) 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": { 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"})