diff --git a/docs/coverage.md b/docs/coverage.md index 8e94802..b50a554 100644 --- a/docs/coverage.md +++ b/docs/coverage.md @@ -240,66 +240,66 @@ Taiwan is its own country, `TWN`, following CEDA's convention. The World Bank la ### sLUC (statistical) -The statistical pipeline covers the following 32 crops, identified by IFPRI MapSPAM crop code. +The statistical pipeline covers the following 32 crops, identified by IFPRI MapSPAM crop code. "Published as" is the `commodity_id` and `commodity_name` each crop's rows carry in the emissions-factors table (see [`data.md`](data.md#3-emissions-factors-parquet)): a FAOSTAT item, or a negative code of our own where FAOSTAT has none. The two coffees share 656 and the two millets share 79, so each pair is published as one row. -| Crop | MapSPAM code | -| -------------- | ------------ | -| ARABICA_COFFEE | ACOF | -| BANANA | BANA | -| BARLEY | BARL | -| BEAN | BEAN | -| CASSAVA | CASS | -| CHICKPEA | CHIC | -| COCONUT | CNUT | -| COTTON | COTT | -| COWPEA | COWP | -| GROUNDNUT | GROU | -| LENTIL | LENT | -| MAIZE | MAIZ | -| OILPALM | OILP | -| OTHER_OILCROPS | OOIL | -| OTHER_PULSES | OPUL | -| PEARL_MILLET | PMIL | -| PIGEONPEA | PIGE | -| PLANTAIN | PLNT | -| POTATO | POTA | -| RAPESEED | RAPE | -| RICE | RICE | -| ROBUSTA_COFFEE | RCOF | -| SESAME_SEED | SESA | -| SMALL_MILLET | SMIL | -| SORGHUM | SORG | -| SOYBEAN | SOYB | -| SUGARBEET | SUGB | -| SUGARCANE | SUGC | -| SUNFLOWER | SUNF | -| SWEET_POTATO | SWPO | -| WHEAT | WHEA | -| YAM | YAMS | +| Crop | MapSPAM code | Published as | +| -------------- | ------------ | ----------------------------------- | +| ARABICA_COFFEE | ACOF | 656 (Coffee, green) | +| BANANA | BANA | 486 (Bananas) | +| BARLEY | BARL | 44 (Barley) | +| BEAN | BEAN | 176 (Beans, dry) | +| CASSAVA | CASS | 125 (Cassava, fresh) | +| CHICKPEA | CHIC | 191 (Chick peas, dry) | +| COCONUT | CNUT | 249 (Coconuts, in shell) | +| COTTON | COTT | 328 (Seed cotton, unginned) | +| COWPEA | COWP | 195 (Cow peas, dry) | +| GROUNDNUT | GROU | 242 (Groundnuts, excluding shelled) | +| LENTIL | LENT | 201 (Lentils, dry) | +| MAIZE | MAIZ | 56 (Maize (corn)) | +| OILPALM | OILP | 254 (Oil palm fruit) | +| OTHER_OILCROPS | OOIL | -201 (OTHER_OILCROPS) | +| OTHER_PULSES | OPUL | -202 (OTHER_PULSES) | +| PEARL_MILLET | PMIL | 79 (Millet) | +| PIGEONPEA | PIGE | 197 (Pigeon peas, dry) | +| PLANTAIN | PLNT | 489 (Plantains and cooking bananas) | +| POTATO | POTA | 116 (Potatoes) | +| RAPESEED | RAPE | 270 (Rape or colza seed) | +| RICE | RICE | 27 (Rice) | +| ROBUSTA_COFFEE | RCOF | 656 (Coffee, green) | +| SESAME_SEED | SESA | 289 (Sesame seed) | +| SMALL_MILLET | SMIL | 79 (Millet) | +| SORGHUM | SORG | 83 (Sorghum) | +| SOYBEAN | SOYB | 236 (Soya beans) | +| SUGARBEET | SUGB | 157 (Sugar beet) | +| SUGARCANE | SUGC | 156 (Sugar cane) | +| SUNFLOWER | SUNF | 267 (Sunflower seed) | +| SWEET_POTATO | SWPO | 122 (Sweet potatoes) | +| WHEAT | WHEA | 15 (Wheat) | +| YAM | YAMS | 137 (Yams) | ### sLUC livestock (statistical) -The statistical pipeline also divides pasture's share of each cell among livestock commodities, by livestock units from Global Pasture Watch's headcounts. +The statistical pipeline also divides pasture's share of each cell among livestock commodities, by livestock units from Global Pasture Watch's headcounts. "Published as" is as for the crops. -| Commodity | Grazers | -| ----------- | ---------------------------------------------------------------------------------- | -| BEEF_CATTLE | Cattle | -| PASTURE | The residual: buffalo, goats, horses and sheep, and pasture no grazer is mapped on | +| Commodity | Grazers | Published as | +| ----------- | ---------------------------------------------------------------------------------- | ---------------------------------------------------- | +| BEEF_CATTLE | Cattle | 867 (Meat of cattle with the bone, fresh or chilled) | +| PASTURE | The residual: buffalo, goats, horses and sheep, and pasture no grazer is mapped on | -3 (PASTURE) | ### jdLUC (jurisdictional-direct) -The jurisdictional-direct pipeline (USA only) covers the following 11 crops, mapped to USDA NASS CDL crop classes. +The jurisdictional-direct pipeline (USA only) covers the following 11 crops, mapped to USDA NASS CDL crop classes. "Published as" is as for the statistical crops, and each code is the one the statistical crop of the same name publishes under. -| Crop | USDA NASS CDL class(es) | -| --------- | --------------------------------------- | -| BARLEY | BARLEY | -| BEAN | DRY_BEANS | -| COTTON | COTTON | -| MAIZE | CORN | -| POTATO | POTATOES | -| RICE | RICE | -| SORGHUM | SORGHUM | -| SOYBEAN | SOYBEANS | -| SUGARBEET | SUGARBEETS | -| SUGARCANE | SUGARCANE | -| WHEAT | DURUM_WHEAT, SPRING_WHEAT, WINTER_WHEAT | +| Crop | USDA NASS CDL class(es) | Published as | +| --------- | --------------------------------------- | --------------------------- | +| BARLEY | BARLEY | 44 (Barley) | +| BEAN | DRY_BEANS | 176 (Beans, dry) | +| COTTON | COTTON | 328 (Seed cotton, unginned) | +| MAIZE | CORN | 56 (Maize (corn)) | +| POTATO | POTATOES | 116 (Potatoes) | +| RICE | RICE | 27 (Rice) | +| SORGHUM | SORGHUM | 83 (Sorghum) | +| SOYBEAN | SOYBEANS | 236 (Soya beans) | +| SUGARBEET | SUGARBEETS | 157 (Sugar beet) | +| SUGARCANE | SUGARCANE | 156 (Sugar cane) | +| WHEAT | DURUM_WHEAT, SPRING_WHEAT, WINTER_WHEAT | 15 (Wheat) | diff --git a/docs/data.md b/docs/data.md index 39e6751..9e8ae16 100644 --- a/docs/data.md +++ b/docs/data.md @@ -83,11 +83,14 @@ emissions = xarray.open_zarr("/.zarr", consolidated=Fal Deposited per data version as `emissions-factors.parquet` (see [Access](#access)). -The final per-(jurisdiction, commodity) table produced by the `trace` stage. Indexed by `(admin_level, admin_id, commodity_name, methodology)` — all strings — where `admin_level` is `PROVINCIAL` (World Bank admin-1) or `NATIONAL` (admin-0, summed from provincial rows) and `methodology` is `JURISDICTIONAL_DIRECT` or `STATISTICAL`. Four `commodity_name` values are not crops. `BEEF_CATTLE` is a livestock commodity: it carries its share of pasture's emissions and area, and production like a crop's. `PASTURE` keeps the rest of pasture's share, with emissions and area but no production. `DROPPED` and `UNATTRIBUTED` appear in both legs and are not commodities at all. `UNATTRIBUTED` holds emissions that count toward the jurisdiction's total but that none of the commodities we model takes, so the commodity rows plus `UNATTRIBUTED` add up to that total. `DROPPED` holds carbon released where, as far as our maps show, the land became neither cropland nor pasture; it does not count toward the total and comes on top of it. `PASTURE`, `DROPPED` and `UNATTRIBUTED` have no production figure, so none of them has an emissions factor. See "Rows that name no commodity" in [`methodology.md`](methodology.md). +The final per-(jurisdiction, commodity) table produced by the `trace` stage. Indexed by `(admin_level, admin_id, commodity_id, methodology)`, where `admin_level` is `PROVINCIAL` (World Bank admin-1) or `NATIONAL` (admin-0, summed from provincial rows), `methodology` is `JURISDICTIONAL_DIRECT` or `STATISTICAL`, and `commodity_id` is an integer: the FAOSTAT item code where FAOSTAT has an item, and a negative code of our own where it does not. `commodity_name` holds FAOSTAT's item name (`Maize (corn)`, `Soya beans`), or our own name for a negative code. Both legs use the same codes, so the same crop has the same `commodity_id` in each. Where MapSPAM splits one FAOSTAT item in two, the statistical leg publishes one row: arabica and robusta coffee add together as 656 (`Coffee, green`), and pearl and small millet as 79 (`Millet`). The statistical leg's `OTHER_OILCROPS` (-201) and `OTHER_PULSES` (-202) are MapSPAM groups of several FAOSTAT items, so they take codes of our own. [`coverage.md`](coverage.md#crops) lists each leg's crops with their codes. + +Four rows are not crops. Beef, 867 (`Meat of cattle with the bone, fresh or chilled`), is a livestock commodity: it carries its share of pasture's emissions and area, and production like a crop's. `PASTURE` (-3) keeps the rest of pasture's share, with emissions and area but no production. `DROPPED` (-1) and `UNATTRIBUTED` (-2) appear in both legs and are not commodities at all. `UNATTRIBUTED` holds emissions that count toward the jurisdiction's total but that none of the commodities we model takes, so the commodity rows plus `UNATTRIBUTED` add up to that total. `DROPPED` holds carbon released where, as far as our maps show, the land became neither cropland nor pasture; it does not count toward the total and comes on top of it. `PASTURE`, `DROPPED` and `UNATTRIBUTED` have no production figure, so none of them has an emissions factor. See "Rows that name no commodity" in [`methodology.md`](methodology.md). | Column | Type | Units | Description | | ---------------------------------- | ------- | ------------ | -------------------------------------------------------------------------------------------------------------------------------- | | `jurisdiction_name` | string | — | Display name (e.g. `United States of America \| Iowa`) | +| `commodity_name` | string | — | FAOSTAT's item name for a positive `commodity_id`, our own name for a negative one | | `commodity_hectares` | float64 | ha | Area the commodity occupies (land occupation), on production's discounted 2000–2020 window | | `peatland_commodity_hectares` | float64 | ha | Commodity area on peatland | | `emissions_mt` | float64 | t CO₂e | Total allocated LUC emissions | @@ -107,4 +110,4 @@ import pandas emission_factors = pandas.read_parquet("emissions-factors.parquet") ``` -The `attribute` stage writes an intermediate rollup parquet with the same index and the non-derived subset of these columns (`commodity_hectares`, `peatland_commodity_hectares`, `emissions_mt`, `peatland_occupation_emissions_mt`, `forest_emissions_mt`, `grassland_emissions_mt`, `peatland_conversion_emissions_mt`, and — statistical only — `production_mt`); `trace` joins production and appends the ratio columns (`production_kg`, `yield_kg_per_ha`, `emissions_factor_kgco2e_per_kg`). +The `attribute` stage writes an intermediate rollup parquet indexed by `(admin_level, admin_id, commodity_name, methodology)`, where `commodity_name` is the leg's own name for the commodity (a CDL crop for the direct leg, a MapSPAM crop for the statistical leg, `BEEF_CATTLE` for beef), with the non-derived subset of these columns (`commodity_hectares`, `peatland_commodity_hectares`, `emissions_mt`, `peatland_occupation_emissions_mt`, `forest_emissions_mt`, `grassland_emissions_mt`, `peatland_conversion_emissions_mt`, and — statistical only — `production_mt` and the herd's heads at each MapSPAM snapshot, `herd_heads_2000` to `herd_heads_2020`); `trace` maps those names onto `commodity_id`, adding together the rows that share one, joins production and appends the ratio columns (`production_kg`, `yield_kg_per_ha`, `emissions_factor_kgco2e_per_kg`). diff --git a/docs/executive_summary.md b/docs/executive_summary.md index 4d3e2ea..3d77760 100644 --- a/docs/executive_summary.md +++ b/docs/executive_summary.md @@ -19,7 +19,7 @@ This document is aimed at readers assessing how the methodology maps to the LSRS ## Data output -The published deliverable is a single per-(jurisdiction, commodity) table — the **emissions-factors parquet** — indexed by `(admin_level, admin_id, commodity_name, methodology)`, where `admin_level` is `NATIONAL` or `PROVINCIAL`, `admin_id` is the ISO 3166 alpha-3 code nationally and the World Bank admin-1 code provincially, and `methodology` is `JURISDICTIONAL_DIRECT` or `STATISTICAL`. Each row carries the headline emissions factor and the quantities behind it; both attribution legs (jdLUC and sLUC) emit the same columns. The full schema is in [`data.md`](data.md#3-emissions-factors-parquet). +The published deliverable is a single per-(jurisdiction, commodity) table — the **emissions-factors parquet** — indexed by `(admin_level, admin_id, commodity_id, methodology)`, where `admin_level` is `NATIONAL` or `PROVINCIAL`, `admin_id` is the ISO 3166 alpha-3 code nationally and the World Bank admin-1 code provincially, `commodity_id` is the FAOSTAT item code (or a negative code of our own where FAOSTAT has none) with its name in `commodity_name`, and `methodology` is `JURISDICTIONAL_DIRECT` or `STATISTICAL`. Each row carries the headline emissions factor and the quantities behind it; both attribution legs (jdLUC and sLUC) emit the same columns. The full schema is in [`data.md`](data.md#3-emissions-factors-parquet). | Column | Units | What it is | | -------------------------------- | ------------ | ---------------------------------------------------------------------------------------------------------- | diff --git a/docs/methodology.md b/docs/methodology.md index 4901272..0826c07 100644 --- a/docs/methodology.md +++ b/docs/methodology.md @@ -269,13 +269,13 @@ Each pixel's emissions are split by the land class it became, cropland or pastur ### From rollups to emissions factors (`trace.py`) -`trace.py` takes the attribution rollup and derives the emissions factor identically for the jurisdictional-direct and statistical methodologies — the only methodology-specific step is where production comes from (NASS yield × area for direct; for statistical, MapSPAM production for crops and GPW heads × FAOSTAT carcass weight per head for livestock). The proxy leg has no production yet: +`trace.py` takes the attribution rollup and derives the emissions factor identically for the jurisdictional-direct and statistical methodologies — the only methodology-specific step is where production comes from (NASS yield × area for direct; for statistical, MapSPAM production for crops and GPW heads × FAOSTAT carcass weight per head for livestock). It then maps each leg's own commodity names onto FAOSTAT item codes (`jdluc/taxonomy.py`), so the legs share one vocabulary, and adds together the rows that share a code, such as MapSPAM's arabica and robusta coffee. The proxy leg has no production yet, and no codes: ``` emissions_factor_kgco2e_per_kg = Σ emissions / Σ production ``` -Provincial (admin-1) rows are then summed up to national (admin-0) totals. The final output is a `pandas.DataFrame` indexed by `(admin_level, admin_id, commodity_name, methodology)` with crop hectares, peatland crop hectares, emissions (with the forest, grassland, peatland-conversion and peatland-occupation components broken out), production, and the emissions factor; see `data.md` for the full column schema. The full stage-by-stage pipeline that produces this table, and how its outputs are cached, is described in `architecture.md`. +Provincial (admin-1) rows are then summed up to national (admin-0) totals. The final output is a `pandas.DataFrame` indexed by `(admin_level, admin_id, commodity_id, methodology)` with the commodity's name, crop hectares, peatland crop hectares, emissions (with the forest, grassland, peatland-conversion and peatland-occupation components broken out), production, and the emissions factor; see `data.md` for the full column schema. The full stage-by-stage pipeline that produces this table, and how its outputs are cached, is described in `architecture.md`. ______________________________________________________________________ diff --git a/docs/validation.md b/docs/validation.md index 1534563..a870e01 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -53,7 +53,7 @@ The split is **by what the code is, not by what it computes**. **One artifact does not reproduce from the lock: the capture.** `efs.parquet` is a re-keyed extract of the pipeline's own parquet outputs, so it is neither committed here nor re-fetchable; producing it means re-running `validation.capture`, which takes hours of pipeline. Everything upstream of a capture reproduces in seconds — eligibility, yields, the anchor tables, and `ORBAE_OVER_WRI`, whose two sides are both external. Only the FAOSTAT half of that needs more than the lock: `prepare.read_faostat_production` reads the ingested parquet rather than a second copy of the upstream archive, so a checkout that has never run the ingest fetches it once and caches it beside the pulled anchors, with `pull.record_digest` pinning what it got. Every measure with a term of ours does not. A report carrying those numbers is reproducible only alongside the run that produced them, and `code_version` and `source_version` on every row name which run that is. -FAOSTAT **is** an ingested dataset, because it is a plausible pipeline input: national tabular production statistics keyed on jdluc's own `admin_id` and `commodity_name`, the global analogue of `usda_nass_quickstats`. WRI, EPA and Orbae stay outside `jdluc/`. Every pipeline module imports `jdluc.datasets`, so putting the yardstick there would let `statistical.py` consume the number it is measured against, with no contract preventing it. For livestock it does: `statistical.py` reads FAOSTAT's livestock table to price `BEEF_CATTLE`'s production, which therefore reproduces FAOSTAT's by construction, less the herd GPW places off mapped pasture. Livestock production has no independent anchor here; the crop table stays outside the attribution path. +FAOSTAT **is** an ingested dataset, because it is a plausible pipeline input: national tabular production statistics keyed on jdluc's own `admin_id` and `commodity_name`, the global analogue of `usda_nass_quickstats`. WRI, EPA and Orbae stay outside `jdluc/`. Every pipeline module imports `jdluc.datasets`, so putting the yardstick there would let `statistical.py` consume the number it is measured against, with no contract preventing it. For livestock it does: `trace.py` reads FAOSTAT's livestock table to price `BEEF_CATTLE`'s production, which therefore reproduces FAOSTAT's by construction, less the herd GPW places off mapped pasture. Livestock production has no independent anchor here; the crop table stays outside the attribution path. ## 3. Choosing what to validate diff --git a/jdluc/__tests__/statistical_test.py b/jdluc/__tests__/statistical_test.py index f2be593..bbf1438 100644 --- a/jdluc/__tests__/statistical_test.py +++ b/jdluc/__tests__/statistical_test.py @@ -10,7 +10,8 @@ from jdluc import emit from jdluc.datasets import gpw_livestock, ifpri_mapspam from jdluc.statistical import ( - MAPSPAM_SNAPSHOT_YEARS, + HERD_HEADS_COLUMNS, + SNAPSHOT_YEAR_TO_WEIGHT, SPAN_TO_MAPSPAM_SPAN, Commodity, Crop, @@ -20,6 +21,7 @@ get_commodity_to_peatland_conversion_share, get_commodity_to_share, get_crop_to_area_share, + get_discounted_snapshot_mean, get_livestock_to_grazing_share, ) @@ -533,12 +535,31 @@ def test_get_canonical_quantity_conserves_the_group_total(group_name: str) -> No # those: a test that redoes the weighting can agree with a wrong implementation. -# Defaults for `get_name_to_totals`, read-only so every call can safely share them. A test passes -# only the shares and yields it cares about, and everything else is zero. +def test_snapshot_year_to_weight() -> None: + assert ( + pytest.approx({2000: 0.03125, 2005: 0.125, 2010: 0.46875, 2020: 0.375}) + == SNAPSHOT_YEAR_TO_WEIGHT + ) + + +def test_snapshot_year_to_weight_reduces_as_the_discounted_snapshot_mean_does() -> None: + year_to_value = {2000: 1.0, 2005: 10.0, 2010: 100.0, 2020: 1000.0} + assert get_value( + darray=get_discounted_snapshot_mean( + year_to_value={ + year: xarray.DataArray(value) for year, value in year_to_value.items() + } + ) + ) == pytest.approx( + sum( + weight * year_to_value[year] + for year, weight in SNAPSHOT_YEAR_TO_WEIGHT.items() + ) + ) + + +# The default for `get_name_to_totals`'s shares, read-only so every call can share it NO_SHARES: collections.abc.Mapping[Commodity, float] = types.MappingProxyType({}) -NO_YIELD: collections.abc.Mapping[int, float] = types.MappingProxyType( - dict.fromkeys(MAPSPAM_SNAPSHOT_YEARS, 0.0) -) def get_name_to_totals( @@ -552,12 +573,11 @@ def get_name_to_totals( pasture_occupation: float = 0.0, pasture_pulse: float = 0.0, peatland_conversion_shares: collections.abc.Mapping[Commodity, float] = NO_SHARES, - year_to_kg_per_head: collections.abc.Mapping[int, float] = NO_YIELD, ) -> dict[str, dict[str, float]]: """`get_commodity_name_to_totals` over WHEAT and the livestock rows on `get_dset`'s 2x2 grid. - Every band, share and rate a caller does not set is zero, so each test moves one thing. The - shares are the same in every span. + Every band and share a caller does not set is zero, so each test moves one thing. The shares + are the same in every span. """ commodities = (Crop.WHEAT, *Livestock) @@ -600,9 +620,6 @@ def get_shares( span_to_commodity_to_share=dict.fromkeys( SPAN_TO_MAPSPAM_SPAN.values(), get_shares(shares=commodity_to_share) ), - species_to_year_to_kg_per_head={ - species: dict(year_to_kg_per_head) for species in gpw_livestock.Species - }, ) @@ -987,37 +1004,26 @@ def test_get_livestock_to_grazing_share( ), ) @pytest.mark.parametrize( - ("pasture_fraction", "year_to_kg_per_head", "expected_kg_per_ha"), + ("pasture_fraction", "expected_heads_per_ha"), ( pytest.param( 0.5, - dict.fromkeys(MAPSPAM_SNAPSHOT_YEARS, 40.0), - 4.0, + 0.1, id="the-whole-herd-counts-where-the-cell-holds-any-pasture-not-a-fraction-weighted-part", ), pytest.param( 0.0, - dict.fromkeys(MAPSPAM_SNAPSHOT_YEARS, 40.0), 0.0, - id="a-herd-in-cells-without-pasture-earns-no-production", - ), - pytest.param( - # 0.1 heads/ha x 40 kg in 2020 alone, at 2020's 37.5% of the window - 0.5, - {2000: 0.0, 2005: 0.0, 2010: 0.0, 2020: 40.0}, - 1.5, - id="each-snapshot-pairs-its-own-heads-with-its-own-rate", + id="a-herd-in-cells-without-pasture-counts-no-heads", ), ), ) -def test_the_beef_row_reports_its_herds_carcass_weight( +def test_the_beef_row_reports_its_herd( pasture_fraction: float, - year_to_kg_per_head: dict[int, float], - expected_kg_per_ha: float, + expected_heads_per_ha: float, units_per_head: dict[gpw_livestock.Species, float], ) -> None: - # Production counts heads, so however a head of cattle is weighed in livestock units, the herd's - # carcass weight is the same + # Livestock units weigh grazing, not the herd totals = get_name_to_totals( dset=get_dset( density=get_density(heads_per_ha={CATTLE: 0.1}), @@ -1025,11 +1031,13 @@ def test_the_beef_row_reports_its_herds_carcass_weight( units_per_head=units_per_head, values={}, ), - year_to_kg_per_head=year_to_kg_per_head, )[Livestock.BEEF_CATTLE.name] - assert totals["production_mt"] == pytest.approx( - PIXELS * HECTARES_PER_CELL * expected_kg_per_ha / 1_000 + assert {column: totals[column] for column in HERD_HEADS_COLUMNS} == pytest.approx( + dict.fromkeys( + HERD_HEADS_COLUMNS, PIXELS * HECTARES_PER_CELL * expected_heads_per_ha + ) ) + assert math.isnan(totals["production_mt"]) # Cattle alone graze the cell, so the beef row holds all its pasture assert totals["commodity_hectares"] == pytest.approx( PIXELS * HECTARES_PER_CELL * pasture_fraction diff --git a/jdluc/__tests__/taxonomy_test.py b/jdluc/__tests__/taxonomy_test.py new file mode 100644 index 0000000..f2907ce --- /dev/null +++ b/jdluc/__tests__/taxonomy_test.py @@ -0,0 +1,49 @@ +import pytest + +from jdluc import jurisdictional_direct, statistical +from jdluc.attribute import Methodology +from jdluc.datasets import faostat, ifpri_mapspam +from jdluc.taxonomy import METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID + +JURISDICTIONAL_DIRECT_NAME_TO_COMMODITY_ID = ( + METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[Methodology.JURISDICTIONAL_DIRECT] +) +STATISTICAL_NAME_TO_COMMODITY_ID = METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[ + Methodology.STATISTICAL +] + + +# The statistical leg's crop of the same name already has its FAOSTAT item, chosen independently +@pytest.mark.parametrize( + "crop", list(jurisdictional_direct.Crop), ids=lambda crop: crop.name.lower() +) +def test_jurisdictional_direct_codes_match_statistical_items( + crop: jurisdictional_direct.Crop, +) -> None: + statistical_crop = statistical.Crop[crop.name] + item_code = faostat.CROP_TO_ITEM_CODE[faostat.Crop[statistical_crop.value]] + assert JURISDICTIONAL_DIRECT_NAME_TO_COMMODITY_ID[crop.name] == item_code + + +@pytest.mark.parametrize( + "crop", + [crop for crop in statistical.Crop if crop.value in faostat.Crop.__members__], + ids=lambda crop: crop.name.lower(), +) +def test_statistical_codes_match_faostat_items(crop: statistical.Crop) -> None: + item_code = faostat.CROP_TO_ITEM_CODE[faostat.Crop[crop.value]] + assert STATISTICAL_NAME_TO_COMMODITY_ID[crop.name] == item_code + + +def test_both_coffees_are_faostat_coffee() -> None: + item_code = faostat.COMBINED_CROP_NAME_TO_ITEM_CODE[ + ifpri_mapspam.Crop2000.COFF.name + ] + assert ( + STATISTICAL_NAME_TO_COMMODITY_ID[statistical.Crop.ARABICA_COFFEE.name] + == item_code + ) + assert ( + STATISTICAL_NAME_TO_COMMODITY_ID[statistical.Crop.ROBUSTA_COFFEE.name] + == item_code + ) diff --git a/jdluc/__tests__/trace_test.py b/jdluc/__tests__/trace_test.py index aa02b59..b30d47e 100644 --- a/jdluc/__tests__/trace_test.py +++ b/jdluc/__tests__/trace_test.py @@ -1,13 +1,17 @@ import numpy import pandas +import pytest from jdluc.attribute import Methodology from jdluc.datasets.worldbank_jurisdictions import AdminLevel +from jdluc.statistical import HERD_HEADS_COLUMNS +from jdluc.taxonomy import CommodityId from jdluc.trace import ( attach_ratios, derive_jurisdictional_production_kg, derive_statistical_production_kg, iter_national_from_provincials, + with_commodity_ids, ) JURISDICTIONAL_EMISSIONS = pandas.DataFrame.from_records( @@ -134,7 +138,7 @@ def test_iter_national_from_provincials_single_province() -> None: [ { "commodity_hectares": 1, - "commodity_name": "CROP_1", + "commodity_id": CommodityId.SOYA_BEANS.value, "emissions_mt": 1, "forest_emissions_mt": 1, "grassland_emissions_mt": 1, @@ -155,7 +159,8 @@ def test_iter_national_from_provincials_single_province() -> None: "admin_id": "ISO_A3", "admin_level": AdminLevel.NATIONAL.name, "commodity_hectares": 1, - "commodity_name": "CROP_1", + "commodity_id": CommodityId.SOYA_BEANS.value, + "commodity_name": "Soya beans", "emissions_factor_kgco2e_per_kg": 1000, "emissions_mt": 1, "forest_emissions_mt": 1, @@ -176,7 +181,7 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: { "admin_id": "ADM_0", "commodity_hectares": 1, - "commodity_name": "CROP_0", + "commodity_id": CommodityId.MAIZE.value, "emissions_mt": 1, "forest_emissions_mt": 1, "grassland_emissions_mt": 1, @@ -188,7 +193,7 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: { "admin_id": "ADM_0", "commodity_hectares": 1, - "commodity_name": "CROP_1", + "commodity_id": CommodityId.SOYA_BEANS.value, "emissions_mt": 1, "forest_emissions_mt": 1, "grassland_emissions_mt": 1, @@ -200,7 +205,7 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: { "admin_id": "ADM_1", "commodity_hectares": 1, - "commodity_name": "CROP_0", + "commodity_id": CommodityId.MAIZE.value, "emissions_mt": 1, "forest_emissions_mt": 1, "grassland_emissions_mt": 1, @@ -211,17 +216,18 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: }, ] ) - (crop_0, crop_1) = iter_national_from_provincials( + (maize, soya_beans) = iter_national_from_provincials( iso_3166="ISO_A3", methodology=Methodology.STATISTICAL, national_name="NAME", provincials=provincials, ) - assert crop_0 == { + assert maize == { "admin_id": "ISO_A3", "admin_level": AdminLevel.NATIONAL.name, "commodity_hectares": 2, - "commodity_name": "CROP_0", + "commodity_id": CommodityId.MAIZE.value, + "commodity_name": "Maize (corn)", "emissions_factor_kgco2e_per_kg": 1000, "emissions_mt": 2, "forest_emissions_mt": 2, @@ -234,11 +240,12 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: "production_kg": 2, "yield_kg_per_ha": 1, } - assert crop_1 == { + assert soya_beans == { "admin_id": "ISO_A3", "admin_level": AdminLevel.NATIONAL.name, "commodity_hectares": 1, - "commodity_name": "CROP_1", + "commodity_id": CommodityId.SOYA_BEANS.value, + "commodity_name": "Soya beans", "emissions_factor_kgco2e_per_kg": 1000, "emissions_mt": 1, "forest_emissions_mt": 1, @@ -268,6 +275,65 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: 5.0, 120.0, 30.0, + 0.0, + 0.0, + 0.0, + 0.0, + ), + ( + "PROVINCIAL", + "BEEF_CATTLE", + "Acre", + "STATISTICAL", + "BRA001", + 100.0, + 0.0, + 0.0, + 200.0, + 0.0, + 200.0, + 0.0, + 100.0, + 100.0, + 100.0, + 100.0, + ), + # A country FAOSTAT reports no herd for + ( + "PROVINCIAL", + "BEEF_CATTLE", + "Niue", + "STATISTICAL", + "NIU001", + 100.0, + 0.0, + 0.0, + 200.0, + 0.0, + 200.0, + 0.0, + 100.0, + 100.0, + 100.0, + 100.0, + ), + ( + "PROVINCIAL", + "BEEF_CATTLE", + "Concepción", + "STATISTICAL", + "PRY001", + 100.0, + 0.0, + 0.0, + 200.0, + 0.0, + 200.0, + 0.0, + 100.0, + 100.0, + 100.0, + 100.0, ), ], columns=[ @@ -283,12 +349,52 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: "production_mt", "forest_emissions_mt", "peatland_conversion_emissions_mt", + "herd_heads_2000", + "herd_heads_2005", + "herd_heads_2010", + "herd_heads_2020", ], ).set_index(["admin_level", "admin_id", "commodity_name", "methodology"]) +# Brazil's cattle rate doubles each snapshot, and its sheep rate catches a herd priced as another +# species. Paraguay reports meat for its herd only in 2020. +STOCKS = pandas.DataFrame.from_records( + [ + *( + ("NATIONAL", "BRA", "Brazil", "CATTLE", year, 1000.0, kg) + for year, kg in ( + (2000, 10_000.0), + (2005, 20_000.0), + (2010, 40_000.0), + (2020, 80_000.0), + ) + ), + ("NATIONAL", "BRA", "Brazil", "SHEEP", 2020, 1.0, 1_000_000.0), + *( + ("NATIONAL", "PRY", "Paraguay", "CATTLE", year, 1000.0, kg) + for year, kg in ( + (2000, numpy.nan), + (2005, numpy.nan), + (2010, numpy.nan), + (2020, 80_000.0), + ) + ), + ], + columns=[ + "admin_level", + "admin_id", + "jurisdiction_name", + "commodity_name", + "year", + "stocks_head", + "production_kg", + ], +).set_index(["admin_level", "admin_id", "jurisdiction_name", "commodity_name", "year"]) def test_derive_statistical_production_kg_stays_indexed() -> None: - result = derive_statistical_production_kg(emissions=STATISTICAL_EMISSIONS) + result = derive_statistical_production_kg( + emissions=STATISTICAL_EMISSIONS, stocks=STOCKS + ) # regression guard: must return an indexed frame like its jurisdictional twin assert result.index.names == [ "admin_level", @@ -297,11 +403,37 @@ def test_derive_statistical_production_kg_stays_indexed() -> None: "methodology", ] assert "production_mt" not in result.columns - assert result.iloc[0]["production_kg"] == 5000.0 # 5 mt x 1000 + assert not set(HERD_HEADS_COLUMNS) & set(result.columns) + assert ( + result.loc[("PROVINCIAL", "USA008", "CORN", "STATISTICAL"), "production_kg"] + == 5000.0 # 5 mt x 1000 + ) # .copy() guard: attach_ratios mutates in place, so the source frame must be untouched assert "production_kg" not in STATISTICAL_EMISSIONS.columns +def test_derive_statistical_production_kg_prices_a_herd_at_each_snapshots_rate() -> ( + None +): + result = derive_statistical_production_kg( + emissions=STATISTICAL_EMISSIONS, stocks=STOCKS + ) + # 100 heads at 10, 20, 40 and 80 kg, weighted 3.125%, 12.5%, 46.875% and 37.5% + assert result.loc[ + ("PROVINCIAL", "BRA001", "BEEF_CATTLE", "STATISTICAL"), "production_kg" + ] == pytest.approx(5156.25) + assert ( + result.loc[ + ("PROVINCIAL", "NIU001", "BEEF_CATTLE", "STATISTICAL"), "production_kg" + ] + == 0.0 + ) + # 100 heads at 80 kg, weighted 37.5%, in 2020 alone + assert result.loc[ + ("PROVINCIAL", "PRY001", "BEEF_CATTLE", "STATISTICAL"), "production_kg" + ] == pytest.approx(3000.0) + + def test_methodologies_share_a_schema() -> None: jurisdictional = attach_ratios( df=derive_jurisdictional_production_kg( @@ -309,8 +441,75 @@ def test_methodologies_share_a_schema() -> None: ) ) statistical = attach_ratios( - df=derive_statistical_production_kg(emissions=STATISTICAL_EMISSIONS) + df=derive_statistical_production_kg( + emissions=STATISTICAL_EMISSIONS, stocks=STOCKS + ) ) assert set(jurisdictional.reset_index().columns) == set( statistical.reset_index().columns ) + + +def test_with_commodity_ids_adds_together_the_names_sharing_an_item() -> None: + emissions = pandas.DataFrame.from_records( + [ + ( + "PROVINCIAL", + "BRA001", + "ARABICA_COFFEE", + "STATISTICAL", + "Acre", + 1.0, + 10.0, + ), + ( + "PROVINCIAL", + "BRA001", + "ROBUSTA_COFFEE", + "STATISTICAL", + "Acre", + 2.0, + 20.0, + ), + ("PROVINCIAL", "BRA001", "MAIZE", "STATISTICAL", "Acre", 4.0, 40.0), + ], + columns=[ + "admin_level", + "admin_id", + "commodity_name", + "methodology", + "jurisdiction_name", + "commodity_hectares", + "production_kg", + ], + ).set_index(["admin_level", "admin_id", "commodity_name", "methodology"]) + result = with_commodity_ids(df=emissions, methodology=Methodology.STATISTICAL) + assert result.to_dict(orient="index") == { + ("PROVINCIAL", "BRA001", CommodityId.MAIZE, "STATISTICAL"): { + "commodity_hectares": 4.0, + "commodity_name": "Maize (corn)", + "jurisdiction_name": "Acre", + "production_kg": 40.0, + }, + ("PROVINCIAL", "BRA001", CommodityId.COFFEE_GREEN, "STATISTICAL"): { + "commodity_hectares": 3.0, + "commodity_name": "Coffee, green", + "jurisdiction_name": "Acre", + "production_kg": 30.0, + }, + } + + +def test_with_commodity_ids_keeps_a_missing_yield_missing() -> None: + result = with_commodity_ids( + df=derive_jurisdictional_production_kg( + emissions=JURISDICTIONAL_EMISSIONS, raw_yields=RAW_YIELDS + ), + methodology=Methodology.JURISDICTIONAL_DIRECT, + ) + assert pandas.isna( + result.loc[ + ("PROVINCIAL", "USA016", CommodityId.WHEAT, "JURISDICTIONAL_DIRECT"), + "production_kg", + ] + ) diff --git a/jdluc/statistical.py b/jdluc/statistical.py index 514333e..f16e21e 100644 --- a/jdluc/statistical.py +++ b/jdluc/statistical.py @@ -2,6 +2,7 @@ import collections.abc import enum import logging +import math import numpy import pandas @@ -77,6 +78,8 @@ class Livestock(enum.StrEnum): assert sum(len(species) for species in LIVESTOCK_TO_SPECIES.values()) == len( set(LIVESTOCK_TO_SPECIES.values()) ) +# `trace` prices a commodity's herd at one species' carcass weight per head +assert all(len(species) == 1 for species in LIVESTOCK_TO_SPECIES.values()) assert {e.name for e in faostat.Species} == {e.name for e in gpw_livestock.Species} @@ -368,6 +371,26 @@ def get_commodity_to_peatland_conversion_share( sorted({year for span in SPAN_TO_MAPSPAM_SPAN.values() for year in span}) ) assert set(MAPSPAM_SNAPSHOT_YEARS).issubset(gpw_livestock.YEARS) +# Each snapshot's weight in `get_discounted_snapshot_mean` +SNAPSHOT_YEAR_TO_WEIGHT = { + year: sum( + weight / 2 + for span, weight in emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT.items() + if year in SPAN_TO_MAPSPAM_SPAN[span] + ) + / sum(emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT.values()) + for year in MAPSPAM_SNAPSHOT_YEARS +} +assert math.isclose(sum(SNAPSHOT_YEAR_TO_WEIGHT.values()), 1) + + +def get_herd_heads_column(year: int) -> str: + return f"herd_heads_{year:d}" + + +HERD_HEADS_COLUMNS = tuple( + get_herd_heads_column(year=year) for year in MAPSPAM_SNAPSHOT_YEARS +) def get_discounted_snapshot_mean( @@ -390,9 +413,9 @@ def get_commodity_name_to_totals( emit.SpanType, dict[Commodity, xarray.DataArray] ], span_to_commodity_to_share: dict[emit.SpanType, dict[Commodity, xarray.DataArray]], - species_to_year_to_kg_per_head: dict[gpw_livestock.Species, dict[int, float]], ) -> dict[str, dict[str, float]]: - """Each commodity's emissions, area and production, summed over a jurisdiction's cells. + """Each commodity's emissions and area, a crop's production and a herd's heads, summed over a + jurisdiction's cells. Both share mappings are keyed by MapSPAM span. Each emissions span uses the shares of the MapSPAM span that contains it (`SPAN_TO_MAPSPAM_SPAN`). @@ -454,8 +477,8 @@ def get_commodity_name_to_totals( } # Peat occupation and production: a livestock commodity takes its grazing share of the - # pasture peat band and reports its herd's carcass weight; a crop takes its area share of - # the cropland peat band and reports MapSPAM's production + # pasture peat band and reports its herd; a crop takes its area share of the cropland peat + # band and reports MapSPAM's production if isinstance(commodity, Livestock): occupation = pastureland_occupation * occupation_shares[commodity] if commodity in LIVESTOCK_TO_SPECIES: @@ -466,22 +489,14 @@ def get_commodity_name_to_totals( # rangeland -- produce nothing. Counting production where the herd stands, as a # crop's is counted where it grows, keeps it on the same animals the grazing split # weighs. - year_to_kg_per_ha = { - year: sum( - ( - dset[gpw_livestock.get_band_name(species=one, year=year)] - * species_to_year_to_kg_per_head[one][year] - for one in sorted(LIVESTOCK_TO_SPECIES[commodity]) - ), - start=xarray.DataArray(numpy.float32(0)), - ).where(dset[f"pasture:fraction:{year:d}"] > 0, other=0) + (species,) = LIVESTOCK_TO_SPECIES[commodity] + totals |= { + get_herd_heads_column(year=year): dset[ + gpw_livestock.get_band_name(species=species, year=year) + ].where(dset[f"pasture:fraction:{year:d}"] > 0, other=0) + * hectares for year in MAPSPAM_SNAPSHOT_YEARS } - totals["production_mt"] = ( - get_discounted_snapshot_mean(year_to_value=year_to_kg_per_ha) - * hectares - / faostat.KG_PER_TONNE - ) else: occupation = cropland_occupation * occupation_shares[commodity] totals["production_mt"] = get_discounted_snapshot_mean( @@ -565,7 +580,9 @@ def get_commodity_name_to_totals( # NB: the residual pasture row names no herd, so it carries no production and `trace` declines # to publish an emission factor for it rather than dividing by zero return { - name: {"production_mt": numpy.nan} | totals + name: {"production_mt": numpy.nan} + | dict.fromkeys(HERD_HEADS_COLUMNS, numpy.nan) + | totals for name, totals in utils.get_sum_totals( enum_to_name_to_darray=commodity_to_totals ).items() @@ -617,7 +634,7 @@ def get_crop_to_area_share( "peatland_conversion_emissions_mt": float, "peatland_occupation_emissions_mt": float, "production_mt": float, -} +} | dict.fromkeys(HERD_HEADS_COLUMNS, float) @storage.cache_to_parquet(version=2) @@ -647,7 +664,6 @@ def workflow( ), ) - stocks = faostat.load(dataset=faostat.LIVESTOCK_DATASET) # Each grazer's livestock units per head, as FAOSTAT counts them for this country merged = merged.assign( get_band_name_to_livestock_units( @@ -659,7 +675,7 @@ def workflow( livestock_units=faostat.load( dataset=faostat.LIVESTOCK_PATTERNS_DATASET ), - stocks=stocks, + stocks=faostat.load(dataset=faostat.LIVESTOCK_DATASET), ).items() }, ) @@ -668,26 +684,6 @@ def workflow( # OOM risk. Doing so speeds up the many downstream reads for each jurisdiction merged = merged.load() - # National carcass weight per standing head. A year FAOSTAT reports no meat for gets a rate of - # 0, which is what `attribute` would make of a missing production anyway. - livestock = stocks.reset_index() - livestock = livestock[ - (livestock["admin_id"] == iso_3166) - & livestock["year"].isin(MAPSPAM_SNAPSHOT_YEARS) - ] - kg_per_head = (livestock["production_kg"] / livestock["stocks_head"]).fillna(0) - species_to_year_to_kg_per_head = { - species: dict.fromkeys(MAPSPAM_SNAPSHOT_YEARS, 0.0) - for species in gpw_livestock.Species - } - for name, year, value in zip( - livestock["commodity_name"].astype(str), - livestock["year"].astype(int), - kg_per_head.astype(float), - strict=True, - ): - species_to_year_to_kg_per_head[gpw_livestock.Species[name]][year] = value - def it() -> collections.abc.Iterator[dict[str, float | str]]: for ( jurisdiction @@ -743,7 +739,6 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]: ), span_to_commodity_to_peatland_conversion_share=span_to_commodity_to_peatland_conversion_share, span_to_commodity_to_share=span_to_commodity_to_share, - species_to_year_to_kg_per_head=species_to_year_to_kg_per_head, ) for commodity_name, totals in commodity_name_to_totals.items(): yield totals | { diff --git a/jdluc/taxonomy.py b/jdluc/taxonomy.py new file mode 100644 index 0000000..767055f --- /dev/null +++ b/jdluc/taxonomy.py @@ -0,0 +1,191 @@ +"""The commodity codes that key `trace`'s output, and each leg's commodities mapped onto them. + +`attribute` names each leg's rows by its own map: CDL crops for the jurisdictional direct leg, +MapSPAM crops for the statistical leg and CROPGRIDS crops for the proxy leg. +`METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID` maps those names onto codes, one entry per leg, and +`trace` adds together the rows that share a code. A code is a FAOSTAT item code where FAOSTAT has an +item, and a negative code of our own where it does not: + +- -1 and -2: DROPPED and UNATTRIBUTED, which every leg has +- -3: PASTURE, the pasture no grazer takes +- -201 and -202: OTHER_OILCROPS and OTHER_PULSES, MapSPAM groups of several FAOSTAT items +""" + +import enum + +from jdluc import attribute, emit, jurisdictional_direct, statistical + + +@enum.unique +class CommodityId(enum.IntEnum): + WHEAT = 15 + RICE = 27 + BARLEY = 44 + MAIZE = 56 + MILLET = 79 + SORGHUM = 83 + POTATOES = 116 + SWEET_POTATOES = 122 + CASSAVA = 125 + YAMS = 137 + SUGAR_CANE = 156 + SUGAR_BEET = 157 + BEANS_DRY = 176 + CHICK_PEAS = 191 + COW_PEAS = 195 + PIGEON_PEAS = 197 + LENTILS = 201 + SOYA_BEANS = 236 + GROUNDNUTS = 242 + COCONUTS = 249 + OIL_PALM_FRUIT = 254 + SUNFLOWER_SEED = 267 + RAPE_OR_COLZA_SEED = 270 + SESAME_SEED = 289 + SEED_COTTON = 328 + BANANAS = 486 + PLANTAINS = 489 + COFFEE_GREEN = 656 + CATTLE_MEAT = 867 + DROPPED = -1 + UNATTRIBUTED = -2 + PASTURE = -3 + OTHER_OILCROPS = -201 + OTHER_PULSES = -202 + + +# FAOSTAT's item names, as its bulk downloads spell them +COMMODITY_ID_TO_NAME = { + CommodityId.WHEAT: "Wheat", + CommodityId.RICE: "Rice", + CommodityId.BARLEY: "Barley", + CommodityId.MAIZE: "Maize (corn)", + CommodityId.MILLET: "Millet", + CommodityId.SORGHUM: "Sorghum", + CommodityId.POTATOES: "Potatoes", + CommodityId.SWEET_POTATOES: "Sweet potatoes", + CommodityId.CASSAVA: "Cassava, fresh", + CommodityId.YAMS: "Yams", + CommodityId.SUGAR_CANE: "Sugar cane", + CommodityId.SUGAR_BEET: "Sugar beet", + CommodityId.BEANS_DRY: "Beans, dry", + CommodityId.CHICK_PEAS: "Chick peas, dry", + CommodityId.COW_PEAS: "Cow peas, dry", + CommodityId.PIGEON_PEAS: "Pigeon peas, dry", + CommodityId.LENTILS: "Lentils, dry", + CommodityId.SOYA_BEANS: "Soya beans", + CommodityId.GROUNDNUTS: "Groundnuts, excluding shelled", + CommodityId.COCONUTS: "Coconuts, in shell", + CommodityId.OIL_PALM_FRUIT: "Oil palm fruit", + CommodityId.SUNFLOWER_SEED: "Sunflower seed", + CommodityId.RAPE_OR_COLZA_SEED: "Rape or colza seed", + CommodityId.SESAME_SEED: "Sesame seed", + CommodityId.SEED_COTTON: "Seed cotton, unginned", + CommodityId.BANANAS: "Bananas", + CommodityId.PLANTAINS: "Plantains and cooking bananas", + CommodityId.COFFEE_GREEN: "Coffee, green", + CommodityId.CATTLE_MEAT: "Meat of cattle with the bone, fresh or chilled", + CommodityId.DROPPED: "DROPPED", + CommodityId.UNATTRIBUTED: "UNATTRIBUTED", + CommodityId.PASTURE: "PASTURE", + CommodityId.OTHER_OILCROPS: "OTHER_OILCROPS", + CommodityId.OTHER_PULSES: "OTHER_PULSES", +} +assert set(COMMODITY_ID_TO_NAME) == set(CommodityId) +assert len(set(COMMODITY_ID_TO_NAME.values())) == len(COMMODITY_ID_TO_NAME) + + +# COTTON is seed cotton: `usda_nass_quickstats` converts NASS's ginned lint to it, as MapSPAM +# carries it. MapSPAM splits coffee into arabica and robusta and millet into pearl and small, where +# FAOSTAT has one item for each +METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID = { + attribute.Methodology.JURISDICTIONAL_DIRECT: { + emit.NonCommodity.DROPPED.name: CommodityId.DROPPED, + emit.NonCommodity.UNATTRIBUTED.name: CommodityId.UNATTRIBUTED, + jurisdictional_direct.Crop.BARLEY.name: CommodityId.BARLEY, + jurisdictional_direct.Crop.BEAN.name: CommodityId.BEANS_DRY, + jurisdictional_direct.Crop.COTTON.name: CommodityId.SEED_COTTON, + jurisdictional_direct.Crop.MAIZE.name: CommodityId.MAIZE, + jurisdictional_direct.Crop.POTATO.name: CommodityId.POTATOES, + jurisdictional_direct.Crop.RICE.name: CommodityId.RICE, + jurisdictional_direct.Crop.SORGHUM.name: CommodityId.SORGHUM, + jurisdictional_direct.Crop.SOYBEAN.name: CommodityId.SOYA_BEANS, + jurisdictional_direct.Crop.SUGARBEET.name: CommodityId.SUGAR_BEET, + jurisdictional_direct.Crop.SUGARCANE.name: CommodityId.SUGAR_CANE, + jurisdictional_direct.Crop.WHEAT.name: CommodityId.WHEAT, + }, + attribute.Methodology.STATISTICAL: { + emit.NonCommodity.DROPPED.name: CommodityId.DROPPED, + emit.NonCommodity.UNATTRIBUTED.name: CommodityId.UNATTRIBUTED, + statistical.Crop.ARABICA_COFFEE.name: CommodityId.COFFEE_GREEN, + statistical.Crop.BANANA.name: CommodityId.BANANAS, + statistical.Crop.BARLEY.name: CommodityId.BARLEY, + statistical.Crop.BEAN.name: CommodityId.BEANS_DRY, + statistical.Crop.CASSAVA.name: CommodityId.CASSAVA, + statistical.Crop.CHICKPEA.name: CommodityId.CHICK_PEAS, + statistical.Crop.COCONUT.name: CommodityId.COCONUTS, + statistical.Crop.COTTON.name: CommodityId.SEED_COTTON, + statistical.Crop.COWPEA.name: CommodityId.COW_PEAS, + statistical.Crop.GROUNDNUT.name: CommodityId.GROUNDNUTS, + statistical.Crop.LENTIL.name: CommodityId.LENTILS, + statistical.Crop.MAIZE.name: CommodityId.MAIZE, + statistical.Crop.OILPALM.name: CommodityId.OIL_PALM_FRUIT, + statistical.Crop.OTHER_OILCROPS.name: CommodityId.OTHER_OILCROPS, + statistical.Crop.OTHER_PULSES.name: CommodityId.OTHER_PULSES, + statistical.Crop.PEARL_MILLET.name: CommodityId.MILLET, + statistical.Crop.PIGEONPEA.name: CommodityId.PIGEON_PEAS, + statistical.Crop.PLANTAIN.name: CommodityId.PLANTAINS, + statistical.Crop.POTATO.name: CommodityId.POTATOES, + statistical.Crop.RAPESEED.name: CommodityId.RAPE_OR_COLZA_SEED, + statistical.Crop.RICE.name: CommodityId.RICE, + statistical.Crop.ROBUSTA_COFFEE.name: CommodityId.COFFEE_GREEN, + statistical.Crop.SESAME_SEED.name: CommodityId.SESAME_SEED, + statistical.Crop.SMALL_MILLET.name: CommodityId.MILLET, + statistical.Crop.SORGHUM.name: CommodityId.SORGHUM, + statistical.Crop.SOYBEAN.name: CommodityId.SOYA_BEANS, + statistical.Crop.SUGARBEET.name: CommodityId.SUGAR_BEET, + statistical.Crop.SUGARCANE.name: CommodityId.SUGAR_CANE, + statistical.Crop.SUNFLOWER.name: CommodityId.SUNFLOWER_SEED, + statistical.Crop.SWEET_POTATO.name: CommodityId.SWEET_POTATOES, + statistical.Crop.WHEAT.name: CommodityId.WHEAT, + statistical.Crop.YAM.name: CommodityId.YAMS, + statistical.Livestock.BEEF_CATTLE.name: CommodityId.CATTLE_MEAT, + statistical.Livestock.PASTURE.name: CommodityId.PASTURE, + }, +} +assert all( + set(commodity_name_to_commodity_id) + == { + *attribute.get_commodity_names(methodology=methodology), + *(non_commodity.name for non_commodity in emit.NonCommodity), + } + for methodology, commodity_name_to_commodity_id in ( + METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID.items() + ) +) +assert set( + METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[ + attribute.Methodology.JURISDICTIONAL_DIRECT + ] +) == { + *(crop.name for crop in jurisdictional_direct.Crop), + *(non_commodity.name for non_commodity in emit.NonCommodity), +} +assert set( + METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[attribute.Methodology.STATISTICAL] +) == { + *(crop.name for crop in statistical.Crop), + *(livestock.name for livestock in statistical.Livestock), + *(non_commodity.name for non_commodity in emit.NonCommodity), +} +assert len( + set( + METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[ + attribute.Methodology.JURISDICTIONAL_DIRECT + ].values() + ) +) == len( + METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[ + attribute.Methodology.JURISDICTIONAL_DIRECT + ] +) diff --git a/jdluc/trace.py b/jdluc/trace.py index 56a8747..de1ad19 100644 --- a/jdluc/trace.py +++ b/jdluc/trace.py @@ -3,12 +3,13 @@ Starting from the per-(admin, commodity) attribution rollup (`attribute.workflow`), determines each commodity's total production — the only step that depends on methodology: - JURISDICTIONAL_DIRECT: production = crop area x NASS QuickStats yield (4-year mean). - - STATISTICAL: production = MAPSPAM production for crops, and GPW heads x FAOSTAT carcass - weight per head for livestock, carried straight through the attribute. + - STATISTICAL: production = MAPSPAM production for crops, carried straight through the + attribute, and the attribute's GPW heads x FAOSTAT carcass weight per head for livestock. It then derives, identically for both, each commodity's yield (kg per hectare) and emissions factor (kgCO2e per kg), and rolls the provincials up to national totals. Returns a cached -pandas.DataFrame indexed by (admin_level, admin_id, commodity_name, methodology). +pandas.DataFrame indexed by (admin_level, admin_id, commodity_id, methodology), where +`commodity_id` is a `taxonomy.CommodityId` and the `commodity_name` column holds its name. The countries come from the positional ISO 3166 alpha-3 codes, or -- with `--backfill` -- from every country in the World Bank admin-0 layer. @@ -21,18 +22,19 @@ import argparse import collections.abc import logging +import typing import pandas -from jdluc import attribute, emit, storage -from jdluc.datasets import usda_nass_quickstats, worldbank_jurisdictions +from jdluc import attribute, emit, statistical, storage, taxonomy +from jdluc.datasets import faostat, usda_nass_quickstats, worldbank_jurisdictions logger = logging.getLogger(__name__) NASS_YIELD_YEARS = (2017, 2018, 2019, 2020) KG_PER_TONNE = 1000 -CANONICAL_KEY = ("admin_level", "admin_id", "commodity_name", "methodology") +CANONICAL_KEY = ("admin_level", "admin_id", "commodity_id", "methodology") ADDITIVE_COLUMNS = ( "commodity_hectares", "emissions_mt", @@ -61,15 +63,78 @@ def derive_jurisdictional_production_kg( merged["production_kg"] = merged["commodity_hectares"] * merged["yield_kg_per_ha"] return ( merged.drop(columns="yield_kg_per_ha") - .set_index(list(CANONICAL_KEY)) + .set_index(emissions.index.names) .sort_index() ) -def derive_statistical_production_kg(emissions: pandas.DataFrame) -> pandas.DataFrame: +def derive_statistical_production_kg( + emissions: pandas.DataFrame, stocks: pandas.DataFrame +) -> pandas.DataFrame: + """A crop's production is MapSPAM's, carried through `attribute`. A herd's is its heads at + each snapshot times FAOSTAT's national carcass weight per standing head that year, reduced + over the window with the weights `statistical.get_discounted_snapshot_mean` gives a crop's. + """ + livestock = stocks.reset_index() + # No meat in a year, or no herd in FAOSTAT, prices at 0 + kg_per_head = ( + (livestock["production_kg"] / livestock["stocks_head"]) + .fillna(0) + .set_axis( + pandas.MultiIndex.from_frame( + livestock[["admin_id", "commodity_name", "year"]] + ) + ) + ) ret = emissions.copy() ret["production_kg"] = ret["production_mt"] * KG_PER_TONNE - return ret.drop(columns="production_mt") + for livestock_commodity, (species,) in statistical.LIVESTOCK_TO_SPECIES.items(): + is_herd = ( + ret.index.get_level_values("commodity_name") == livestock_commodity.name + ) + iso_3166s = ret.index.get_level_values("admin_id")[is_herd].str[:3] + ret.loc[is_herd, "production_kg"] = sum( + weight + * ret.loc[is_herd, statistical.get_herd_heads_column(year=year)] + * kg_per_head.reindex( + fill_value=0, + index=pandas.MultiIndex.from_arrays( + [ + iso_3166s, + [species.name] * len(iso_3166s), + [year] * len(iso_3166s), + ] + ), + ).to_numpy() + for year, weight in statistical.SNAPSHOT_YEAR_TO_WEIGHT.items() + ) + return ret.drop(columns=["production_mt", *statistical.HERD_HEADS_COLUMNS]) + + +def with_commodity_ids( + df: pandas.DataFrame, methodology: attribute.Methodology +) -> pandas.DataFrame: + """Re-key rows named by the leg's own map onto `taxonomy.CommodityId`, adding together the rows + that share a code, such as arabica and robusta coffee. + """ + commodity_name_to_commodity_id = ( + taxonomy.METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[methodology] + ) + flat = df.reset_index() + flat["commodity_id"] = [ + int(commodity_name_to_commodity_id[commodity_name]) + for commodity_name in flat["commodity_name"] + ] + grouped = flat.drop(columns="commodity_name").groupby(by=list(CANONICAL_KEY)) + # min_count keeps a row with no production at all, such as one NASS has no yield for, at NaN + ret = grouped.sum(min_count=1, numeric_only=True).assign( + jurisdiction_name=grouped["jurisdiction_name"].first() + ) + ret["commodity_name"] = [ + taxonomy.COMMODITY_ID_TO_NAME[taxonomy.CommodityId(commodity_id)] + for commodity_id in ret.index.get_level_values("commodity_id") + ] + return ret def attach_ratios(df: pandas.DataFrame) -> pandas.DataFrame: @@ -92,7 +157,8 @@ def iter_national_from_provincials( national_name: str, provincials: pandas.DataFrame, ) -> collections.abc.Iterator[dict[str, str | float]]: - for commodity_name, group in sorted(provincials.groupby(by="commodity_name")): + for key, group in sorted(provincials.groupby(by="commodity_id")): + commodity_id = taxonomy.CommodityId(typing.cast(int, key)) def safe_div(numer: float, denom: float) -> float: return float("nan") if denom == 0 else numer / denom @@ -103,7 +169,8 @@ def safe_div(numer: float, denom: float) -> float: } yield ret | { "admin_level": worldbank_jurisdictions.AdminLevel.NATIONAL.name, - "commodity_name": str(commodity_name), + "commodity_id": int(commodity_id), + "commodity_name": taxonomy.COMMODITY_ID_TO_NAME[commodity_id], "jurisdiction_name": national_name, "admin_id": iso_3166, "yield_kg_per_ha": safe_div( @@ -129,15 +196,25 @@ def workflow( iso_3166s=iso_3166s, methodology=methodology, ) - assert tuple(emissions.index.names) == CANONICAL_KEY - emissions_and_yields = ( + assert emissions.index.names == [ + "admin_level", + "admin_id", + "commodity_name", + "methodology", + ] + production = ( derive_jurisdictional_production_kg( emissions=emissions, raw_yields=usda_nass_quickstats.load() ) if methodology == attribute.Methodology.JURISDICTIONAL_DIRECT - else derive_statistical_production_kg(emissions=emissions) + else derive_statistical_production_kg( + emissions=emissions, + stocks=faostat.load(dataset=faostat.LIVESTOCK_DATASET), + ) + ) + provincials = attach_ratios( + df=with_commodity_ids(df=production, methodology=methodology) ) - provincials = attach_ratios(df=emissions_and_yields) assert all( admin_level_name == worldbank_jurisdictions.AdminLevel.PROVINCIAL.name for admin_level_name in map( diff --git a/pyproject.toml b/pyproject.toml index 6e37fee..6073cbb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -76,6 +76,7 @@ name = "ETL" type = "layers" layers = [ "jdluc.trace", + "jdluc.taxonomy", "jdluc.attribute", "jdluc.jurisdictional_direct | jdluc.statistical", "jdluc.emit", @@ -95,6 +96,7 @@ source_modules = [ ] forbidden_modules = [ "jdluc.attribute", + "jdluc.taxonomy", "jdluc.emit", "jdluc.harmonize", "jdluc.ingest", @@ -112,6 +114,7 @@ source_modules = [ ] forbidden_modules = [ "jdluc.attribute", + "jdluc.taxonomy", "jdluc.emit", "jdluc.harmonize", "jdluc.ingest", diff --git a/tools/export-graphs.py b/tools/export-graphs.py index 60a0e78..fcaf9cf 100644 --- a/tools/export-graphs.py +++ b/tools/export-graphs.py @@ -11,8 +11,8 @@ import numpy import pandas -from jdluc import attribute, statistical, trace -from jdluc.datasets import ifpri_mapspam, worldbank_jurisdictions +from jdluc import attribute, taxonomy, trace +from jdluc.datasets import worldbank_jurisdictions logger = logging.getLogger(__name__) @@ -38,11 +38,6 @@ def from_dataframe(cls, df: pandas.DataFrame) -> typing.Self: production_kg=float(srs.production_kg), ) - @property - def crop_description(self) -> str: - crop_name_2005 = statistical.Crop[self.commodity_name].value - return ifpri_mapspam.Crop2005[crop_name_2005].value - @property def payload(self) -> dict[str, float | str]: emissions_kgco2e = self.emissions_factor_kgco2e_per_kg * self.production_kg @@ -82,30 +77,56 @@ def csv(self) -> str: "Brazil | Name Unknown", } -METHODOLOGY_TO_ISO_3166_TO_CROP_NAMES = { - attribute.Methodology.JURISDICTIONAL_DIRECT: {"USA": ("MAIZE", "SOYBEAN", "WHEAT")}, +METHODOLOGY_TO_ISO_3166_TO_COMMODITY_IDS = { + attribute.Methodology.JURISDICTIONAL_DIRECT: { + "USA": ( + taxonomy.CommodityId.MAIZE, + taxonomy.CommodityId.SOYA_BEANS, + taxonomy.CommodityId.WHEAT, + ) + }, attribute.Methodology.STATISTICAL: { - "BRA": ("MAIZE", "SOYBEAN", "OILPALM"), - "BOL": ("SOYBEAN", "PLANTAIN", "BANANA"), - "ARG": ("SOYBEAN",), - "CAN": ("SOYBEAN",), - "IDN": ("OILPALM",), - "PNG": ("OILPALM",), - "GHA": ("RICE", "WHEAT"), - "VNM": ("RICE", "SORGHUM"), - "MYS": ("RICE", "SORGHUM", "OILPALM"), - "CUB": ("SUGARCANE", "RICE", "BANANA"), - "ZAF": ("RICE", "WHEAT", "SORGHUM"), + "BRA": ( + taxonomy.CommodityId.MAIZE, + taxonomy.CommodityId.SOYA_BEANS, + taxonomy.CommodityId.OIL_PALM_FRUIT, + ), + "BOL": ( + taxonomy.CommodityId.SOYA_BEANS, + taxonomy.CommodityId.PLANTAINS, + taxonomy.CommodityId.BANANAS, + ), + "ARG": (taxonomy.CommodityId.SOYA_BEANS,), + "CAN": (taxonomy.CommodityId.SOYA_BEANS,), + "IDN": (taxonomy.CommodityId.OIL_PALM_FRUIT,), + "PNG": (taxonomy.CommodityId.OIL_PALM_FRUIT,), + "GHA": (taxonomy.CommodityId.RICE, taxonomy.CommodityId.WHEAT), + "VNM": (taxonomy.CommodityId.RICE, taxonomy.CommodityId.SORGHUM), + "MYS": ( + taxonomy.CommodityId.RICE, + taxonomy.CommodityId.SORGHUM, + taxonomy.CommodityId.OIL_PALM_FRUIT, + ), + "CUB": ( + taxonomy.CommodityId.SUGAR_CANE, + taxonomy.CommodityId.RICE, + taxonomy.CommodityId.BANANAS, + ), + "ZAF": ( + taxonomy.CommodityId.RICE, + taxonomy.CommodityId.WHEAT, + taxonomy.CommodityId.SORGHUM, + ), }, } -def subset_to_crop_names( - df: pandas.DataFrame, crop_names: tuple[str, ...] +def subset_to_commodity_ids( + commodity_ids: tuple[taxonomy.CommodityId, ...], df: pandas.DataFrame ) -> pandas.DataFrame: - ret = df[df.index.get_level_values("commodity_name").isin(crop_names)] + ret = df[df.index.get_level_values("commodity_id").isin(commodity_ids)] missing = sorted( - set(crop_names) - set(map(str, ret.index.get_level_values("commodity_name"))) + set(commodity_ids) - set(map(int, ret.index.get_level_values("commodity_id"))) ) if missing: logger.warning(f"Workflow returned no rows for {missing=}") @@ -113,18 +134,18 @@ def subset_to_crop_names( def iter_dfs( - iso_3166_to_crop_names: dict[str, tuple[str, ...]], + iso_3166_to_commodity_ids: dict[str, tuple[taxonomy.CommodityId, ...]], methodology: attribute.Methodology, ) -> collections.abc.Iterator[pandas.DataFrame]: # The full commodity list, so this reuses the caches an ordinary trace.py run populates workflow_commodity_names = attribute.get_commodity_names(methodology=methodology) - for iso_3166, crop_names in iso_3166_to_crop_names.items(): + for iso_3166, commodity_ids in iso_3166_to_commodity_ids.items(): logger.info( f"Tracing {methodology.name:s}/{iso_3166:s} with " - f"{len(workflow_commodity_names):d} commodities, tabulating {len(crop_names):d}" + f"{len(workflow_commodity_names):d} commodities, tabulating {len(commodity_ids):d}" ) - yield subset_to_crop_names( - crop_names=crop_names, + yield subset_to_commodity_ids( + commodity_ids=commodity_ids, df=trace.workflow( commodity_names=workflow_commodity_names, iso_3166s=(iso_3166,), @@ -136,15 +157,13 @@ def iter_dfs( def iter_path_record( df: pandas.DataFrame, ) -> collections.abc.Generator[tuple[tuple[str, ...], Record]]: - for commodity_name, crop_df in df.reset_index(drop=False).groupby( - by="commodity_name" - ): + for commodity_id, crop_df in df.reset_index(drop=False).groupby(by="commodity_id"): nationals = crop_df[ crop_df.admin_level == worldbank_jurisdictions.AdminLevel.NATIONAL.name ] for national_name, national_df in nationals.groupby(by="jurisdiction_name"): record = Record.from_dataframe(df=national_df) - yield (record.crop_description, national_name), record + yield (record.commodity_name, national_name), record provincials = crop_df[ ( crop_df.admin_level @@ -157,13 +176,13 @@ def iter_path_record( ): if provincial_name in PROVINCE_NAME_BLOCKLIST: logger.warning( - f"Skipping {commodity_name=:s}/{provincial_name=:s} because it is blocklisted" + f"Skipping {commodity_id=:d}/{provincial_name=:s} because it is blocklisted" ) continue else: record = Record.from_dataframe(df=provincial_df) yield ( - (record.crop_description, *provincial_name.split(" | ")), + (record.commodity_name, *provincial_name.split(" | ")), record, ) @@ -176,7 +195,7 @@ def write_data_to_graph( with zipfile.ZipFile(buffer, compression=zipfile.ZIP_DEFLATED, mode="w") as zf: for df in dfs: for path_tokens, record in iter_path_record(df=df): - arcname = "/".join((*path_tokens, record.crop_description + ".csv")) + arcname = "/".join((*path_tokens, record.commodity_name + ".csv")) zf.writestr(arcname, record.csv) with tempfile.NamedTemporaryFile(delete=False, mode="wb", suffix=".zip") as fp: @@ -194,11 +213,11 @@ def main() -> int: dfs=[ pandas.concat( iter_dfs( - iso_3166_to_crop_names=iso_3166_to_crop_names, + iso_3166_to_commodity_ids=iso_3166_to_commodity_ids, methodology=methodology, ) ) - for methodology, iso_3166_to_crop_names in METHODOLOGY_TO_ISO_3166_TO_CROP_NAMES.items() + for methodology, iso_3166_to_commodity_ids in METHODOLOGY_TO_ISO_3166_TO_COMMODITY_IDS.items() ] ) diff --git a/tools/measure-drift.py b/tools/measure-drift.py index 4add01a..e43fc9f 100644 --- a/tools/measure-drift.py +++ b/tools/measure-drift.py @@ -13,8 +13,9 @@ cache keys are in the run log. `--isolated` sidesteps it with a private SCRATCH_ROOT per side, at the price of the full pipeline rather than a cache read. -Every commodity in the capture gets checked -- Target.commodity_names is the same -attribute.get_commodity_names call that drives the capture, so the two cannot diverge. The cache +Every commodity in the capture gets checked -- Target.commodity_ids is the same +attribute.get_commodity_names call that drives the capture, mapped onto `taxonomy` codes, so the +two cannot diverge. The cache key includes commodity_names, so asking for a subset would also cold-miss every layer an ordinary `uv run python jdluc/trace.py USA` populated; captures are one country at a time over the full per-methodology list. Scoring the same list means each sum drift is a true all-crop total, and @@ -60,7 +61,7 @@ import pandas -from jdluc import attribute, config, storage, trace +from jdluc import attribute, config, storage, taxonomy, trace logger = logging.getLogger(__name__) @@ -101,7 +102,7 @@ KG_PER_TONNE = 1000 # Two sides can be written with different index levels, so both are re-keyed onto this. admin_id # rather than jurisdiction_name: a machine identifier, where a display string could be renamed. -CANONICAL_KEY = ("admin_level", "admin_id", "commodity_name", "methodology") +CANONICAL_KEY = ("admin_level", "admin_id", "commodity_id", "methodology") # Columns this repo has renamed, old name to current. A baseline older than a rename writes the # old name, and this tool grades both sides with one vocabulary -- without which the two columns # would fail the "in both frames" guard and drop out of every check below unremarked. @@ -222,11 +223,23 @@ class Target: iso_3166: str @property - def commodity_names(self) -> tuple[str, ...]: + def commodity_ids(self) -> tuple[int, ...]: # The same list the capture hands the pipeline, so what gets checked can never be less # than what ran. A curated subset scores a partial sum, which moves whenever a crop # trades with a sibling outside it -- indistinguishable from mass being created. - return attribute.get_commodity_names(methodology=self.methodology) + commodity_name_to_commodity_id = ( + taxonomy.METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[self.methodology] + ) + return tuple( + sorted( + { + int(commodity_name_to_commodity_id[commodity_name]) + for commodity_name in attribute.get_commodity_names( + methodology=self.methodology + ) + } + ) + ) @property def slug(self) -> str: @@ -305,18 +318,30 @@ def get_dirty_paths(cwd: pathlib.Path) -> list[str]: def load_capture( - directory: pathlib.Path, slug: str + directory: pathlib.Path, target: Target ) -> tuple[pandas.DataFrame, tuple[str, ...]]: """Read one side's parquet onto CANONICAL_KEY, whatever index it was written with. Whichever of admin_id and jurisdiction_name is not in the index is a column, so resetting and - re-keying realigns two sides written differently. Also returns the index each was *written* - with, since a schema difference and a value disagreement look identical once rows fail to line - up. + re-keying realigns two sides written differently. A side keyed by each leg's own commodity + names, from before `trace` keyed on `commodity_id`, is re-keyed the way `trace` does it, with + its ratios derived again over the rows that were added together. Also returns the index each + was *written* with, since a schema difference and a value disagreement look identical once rows + fail to line up. """ + slug = target.slug raw = pandas.read_parquet(directory / f"{slug:s}.parquet") written = tuple(str(name) for name in raw.index.names if name is not None) df = raw.reset_index().rename(columns=RENAMED_COLUMNS) + if "commodity_id" not in df.columns: + df = trace.attach_ratios( + df=trace.with_commodity_ids( + df=df.set_index( + ["admin_level", "admin_id", "commodity_name", "methodology"] + ), + methodology=target.methodology, + ) + ).reset_index() missing = [name for name in CANONICAL_KEY if name not in df.columns] assert not missing, ( f"{slug:s} in {directory} lacks {missing}; is it the schema you think it is?" @@ -326,7 +351,7 @@ def load_capture( def check_alignment( baseline: pandas.DataFrame, - commodity_names: tuple[str, ...], + commodity_ids: tuple[int, ...], head: pandas.DataFrame, slug: str, written: tuple[tuple[str, ...], tuple[str, ...]], @@ -367,7 +392,7 @@ def check_alignment( ) # A target that silently stopped existing would otherwise be reported as nothing at all if absent := sorted( - set(commodity_names) - set(df.index.get_level_values("commodity_name")) + set(commodity_ids) - set(df.index.get_level_values("commodity_id")) ): findings.append( Finding( @@ -419,15 +444,15 @@ def print_worst_rows( for value, column, row, before_value, after_value in sorted(records, reverse=True)[ :limit ]: - admin_id, commodity_name = row[1], row[2] + admin_id, commodity_id = row[1], row[2] print( - f" {commodity_name!s:<16}{admin_id!s:<9}{column:<34}" + f" {commodity_id!s:<16}{admin_id!s:<9}{column:<34}" f"{before_value:>13.5g}{after_value:>13.5g}{value:>11.3e}" ) def check_excluded_crops( - baseline: pandas.DataFrame, crops: list[str], head: pandas.DataFrame, rtol: float + baseline: pandas.DataFrame, crops: list[int], head: pandas.DataFrame, rtol: float ) -> list[Finding]: """Say what the hardcoded crop filter left unchecked, and whether it moved. @@ -441,11 +466,11 @@ def get_outside(df: pandas.DataFrame) -> pandas.DataFrame: # Per frame, because the two sides need not hold the same rows: a change that adds a row # -- pastureland becoming a commodity of its own -- makes one side longer than the other, # and a mask built on the baseline cannot index the head - return df[~df.index.get_level_values("commodity_name").isin(crops)] + return df[~df.index.get_level_values("commodity_id").isin(crops)] left, right = get_outside(baseline), get_outside(head) shared = left.index.intersection(right.index) - excluded = sorted(set(left.index.get_level_values("commodity_name"))) + excluded = sorted(set(left.index.get_level_values("commodity_id"))) if not len(shared) or not excluded: return [] worst, where, flips, blind = 0.0, "", 0, 0 @@ -461,7 +486,7 @@ def get_outside(df: pandas.DataFrame) -> pandas.DataFrame: blind += 1 elif float(relative.max()) > worst: row = typing.cast(tuple[str, ...], relative.idxmax()) - worst, where = float(relative.max()), f"{column:s} {row[2]:s}/{row[1]:s}" + worst, where = float(relative.max()), f"{column:s} {row[2]!s}/{row[1]:s}" # Three distinct outcomes, and the first two used to print alike: something moved, nothing was # comparable, and everything was comparable and sat exactly still measured = ( @@ -482,7 +507,7 @@ def get_outside(df: pandas.DataFrame) -> pandas.DataFrame: Finding( message=( f"{len(excluded):d} crop(s) outside the checked set moved where nothing " - f"checked them: {measured:s}, plus {flips:d} NaN flip(s). Target.commodity_names " + f"checked them: {measured:s}, plus {flips:d} NaN flip(s). Target.commodity_ids " "should cover every crop the capture holds, so this means they diverged" ), severity=Severity.ADVISORY, @@ -517,9 +542,9 @@ def check_value_drift( Filtered to the hardcoded crops, which is what that subset is for; alignment ran over the whole table. """ - crops = list(target.commodity_names) - left_all = baseline[baseline.index.get_level_values("commodity_name").isin(crops)] - right_all = head[head.index.get_level_values("commodity_name").isin(crops)] + crops = list(target.commodity_ids) + left_all = baseline[baseline.index.get_level_values("commodity_id").isin(crops)] + right_all = head[head.index.get_level_values("commodity_id").isin(crops)] shared = left_all.index.intersection(right_all.index) columns = sorted( set(left_all.select_dtypes("number").columns) @@ -611,15 +636,13 @@ def check_sum_drift( def nationals(df: pandas.DataFrame) -> pandas.DataFrame: return df[ (df.index.get_level_values("admin_level") == NATIONAL) - & df.index.get_level_values("commodity_name").isin( - list(target.commodity_names) - ) + & df.index.get_level_values("commodity_id").isin(list(target.commodity_ids)) ] left, right = nationals(baseline), nationals(head) print( f"\n {target.slug:s}: {len(left):d} vs {len(right):d} national row(s) over " - f"{len(target.commodity_names):d} hardcoded crop(s)" + f"{len(target.commodity_ids):d} hardcoded crop(s)" ) print(f" {'column':<34}{'baseline':>15}{'head':>15}{'delta':>14}{'rel':>11}") totals = [ @@ -688,8 +711,8 @@ def worst_relative(left: pandas.Series, right: pandas.Series) -> float: broken: dict[Invariant, str] = {} levels = df.index.get_level_values("admin_level") - national = df[levels == NATIONAL].reset_index().set_index("commodity_name") - provincial = df[levels == PROVINCIAL].reset_index().groupby("commodity_name") + national = df[levels == NATIONAL].reset_index().set_index("commodity_id") + provincial = df[levels == PROVINCIAL].reset_index().groupby("commodity_id") columns = [column for column in trace.ADDITIVE_COLUMNS if column in df.columns] rollup = { @@ -774,7 +797,7 @@ def compare_invariants( class LegAgreement: """How far apart the two attribution legs price one crop, at one commit.""" - commodity_name: str + commodity_id: int statistical: float jurisdictional_direct: float @@ -798,17 +821,17 @@ def factors(df: pandas.DataFrame) -> pandas.Series: national = df[df.index.get_level_values("admin_level") == NATIONAL] priced = national[national["production_kg"] > 0] factor = priced["emissions_mt"] * KG_PER_TONNE / priced["production_kg"] - return factor.groupby(level="commodity_name").first() + return factor.groupby(level="commodity_id").first() left, right = factors(statistical), factors(jurisdictional_direct) return [ LegAgreement( - commodity_name=str(commodity_name), - jurisdictional_direct=float(right[commodity_name]), - statistical=float(left[commodity_name]), + commodity_id=int(commodity_id), + jurisdictional_direct=float(right[commodity_id]), + statistical=float(left[commodity_id]), ) - for commodity_name in sorted(set(left.index) & set(right.index)) - if right[commodity_name] + for commodity_id in sorted(set(left.index) & set(right.index)) + if right[commodity_id] ] @@ -826,15 +849,15 @@ def compare_leg_agreement( if not head: print(" no crop is priced by both legs, so there is nothing to compare") return [] - before = {agreement.commodity_name: agreement for agreement in baseline} + before = {agreement.commodity_id: agreement for agreement in baseline} print( f" {'crop':<16}{'jurisdictional':>16}{'statistical':>14}{'ratio':>8}" f"{'was':>14}{'ratio':>8}" ) for agreement in head: - was = before.get(agreement.commodity_name) + was = before.get(agreement.commodity_id) print( - f" {agreement.commodity_name:<16}{agreement.jurisdictional_direct:>16.4f}" + f" {agreement.commodity_id:<16}{agreement.jurisdictional_direct:>16.4f}" f"{agreement.statistical:>14.4f}{agreement.ratio:>8.2f}" + ( f"{was.statistical:>14.4f}{was.ratio:>8.2f}" @@ -1062,7 +1085,7 @@ def compare_target( Alignment first, because everything after it assumes the rows line up. """ pair = { - side: load_capture(directory=out / str(side), slug=target.slug) + side: load_capture(directory=out / str(side), target=target) for side in (Side.BASELINE, Side.HEAD) } left, right = pair[Side.BASELINE][0], pair[Side.HEAD][0] @@ -1070,7 +1093,7 @@ def compare_target( print(f"\n=== {target.slug:s}: alignment ===") findings = check_alignment( baseline=left, - commodity_names=target.commodity_names, + commodity_ids=target.commodity_ids, head=right, slug=target.slug, written=(pair[Side.BASELINE][1], pair[Side.HEAD][1]), @@ -1240,7 +1263,7 @@ def main() -> int: side: { methodology: load_capture( directory=out / str(side), - slug=Target(iso_3166=iso_3166, methodology=methodology).slug, + target=Target(iso_3166=iso_3166, methodology=methodology), )[0] for methodology in methodologies } diff --git a/validation/__tests__/crop_codes_test.py b/validation/__tests__/crop_codes_test.py new file mode 100644 index 0000000..909c6df --- /dev/null +++ b/validation/__tests__/crop_codes_test.py @@ -0,0 +1,82 @@ +"""Checks on how captured rows, keyed by FAOSTAT item, are put back into MapSPAM codes.""" + +import pandas +import pytest + +from jdluc import statistical +from jdluc.attribute import Methodology +from jdluc.datasets import ifpri_mapspam +from jdluc.taxonomy import METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID, CommodityId +from validation import prepare + +STATISTICAL_NAME_TO_COMMODITY_ID = METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[ + Methodology.STATISTICAL +] + + +@pytest.mark.parametrize( + "crop", + [ + crop + for crop in statistical.Crop + if crop.value not in prepare.SPLIT_CROP_CODE_TO_COMBINED_CROP_CODE + ], + ids=lambda crop: crop.name.lower(), +) +def test_a_crop_with_its_own_item_keeps_its_own_code(crop: statistical.Crop) -> None: + commodity_id = STATISTICAL_NAME_TO_COMMODITY_ID[crop.name] + assert prepare.COMMODITY_ID_TO_CROP_CODE[commodity_id] == crop.value + + +@pytest.mark.parametrize( + ("crop", "combined_crop_code"), + [ + pytest.param( + statistical.Crop.ARABICA_COFFEE, + ifpri_mapspam.Crop2000.COFF.name, + id="arabica-coffee", + ), + pytest.param( + statistical.Crop.ROBUSTA_COFFEE, + ifpri_mapspam.Crop2000.COFF.name, + id="robusta-coffee", + ), + pytest.param( + statistical.Crop.PEARL_MILLET, + ifpri_mapspam.Crop2000.MILL.name, + id="pearl-millet", + ), + pytest.param( + statistical.Crop.SMALL_MILLET, + ifpri_mapspam.Crop2000.MILL.name, + id="small-millet", + ), + ], +) +def test_a_split_crop_is_compared_under_its_combined_code( + crop: statistical.Crop, combined_crop_code: str +) -> None: + commodity_id = STATISTICAL_NAME_TO_COMMODITY_ID[crop.name] + assert prepare.COMMODITY_ID_TO_CROP_CODE[commodity_id] == combined_crop_code + assert ( + prepare.SPLIT_CROP_CODE_TO_COMBINED_CROP_CODE[crop.value] == combined_crop_code + ) + + +def test_with_crop_codes_keeps_only_crop_rows() -> None: + frame = pandas.DataFrame( + { + "commodity_id": [ + CommodityId.MAIZE, + CommodityId.COFFEE_GREEN, + CommodityId.CATTLE_MEAT, + CommodityId.PASTURE, + CommodityId.DROPPED, + CommodityId.UNATTRIBUTED, + ] + } + ).astype({"commodity_id": int}) + assert prepare.with_crop_codes(frame=frame)["crop_code"].tolist() == [ + statistical.Crop.MAIZE.value, + ifpri_mapspam.Crop2000.COFF.name, + ] diff --git a/validation/capture.py b/validation/capture.py index 663970b..1365a18 100644 --- a/validation/capture.py +++ b/validation/capture.py @@ -4,7 +4,7 @@ `python -m validation`, and its output is cached rather than recomputed per report. **One capture, both legs.** `efs.parquet` holds sLUC and jdLUC together, indexed on -`trace.CANONICAL_KEY` -- `(admin_level, admin_id, commodity_name, methodology)` -- so the US +`trace.CANONICAL_KEY` -- `(admin_level, admin_id, commodity_id, methodology)` -- so the US sLUC-versus-jdLUC head-to-head is a filter on `methodology` within one table rather than a join between two files. A join there would invite exactly the alignment failure the comparison exists to rule out. diff --git a/validation/prepare.py b/validation/prepare.py index 708b5db..6448de0 100644 --- a/validation/prepare.py +++ b/validation/prepare.py @@ -30,7 +30,7 @@ import iso3166 import pandas -from jdluc import statistical +from jdluc import attribute, statistical, taxonomy from jdluc.datasets import faostat, ifpri_mapspam from validation import pull, schema, targets @@ -756,8 +756,58 @@ def read_efs() -> pandas.DataFrame | None: return pandas.read_parquet(pull.EFS) -def keep_crop_rows(frame: pandas.DataFrame) -> pandas.DataFrame: - """Only the crop rows, kept before anything indexes `statistical.Crop`. +# The MapSPAM code each captured crop row is compared under. The capture keys on FAOSTAT items, and +# FAOSTAT has one item each for coffee and millet where MapSPAM 2005 onward has two, so those two +# take MapSPAM 2000's combined codes. +COMMODITY_ID_TO_CROP_CODE = { + taxonomy.CommodityId.BANANAS: statistical.Crop.BANANA.value, + taxonomy.CommodityId.BARLEY: statistical.Crop.BARLEY.value, + taxonomy.CommodityId.BEANS_DRY: statistical.Crop.BEAN.value, + taxonomy.CommodityId.CASSAVA: statistical.Crop.CASSAVA.value, + taxonomy.CommodityId.CHICK_PEAS: statistical.Crop.CHICKPEA.value, + taxonomy.CommodityId.COCONUTS: statistical.Crop.COCONUT.value, + taxonomy.CommodityId.COFFEE_GREEN: ifpri_mapspam.Crop2000.COFF.name, + taxonomy.CommodityId.COW_PEAS: statistical.Crop.COWPEA.value, + taxonomy.CommodityId.GROUNDNUTS: statistical.Crop.GROUNDNUT.value, + taxonomy.CommodityId.LENTILS: statistical.Crop.LENTIL.value, + taxonomy.CommodityId.MAIZE: statistical.Crop.MAIZE.value, + taxonomy.CommodityId.MILLET: ifpri_mapspam.Crop2000.MILL.name, + taxonomy.CommodityId.OIL_PALM_FRUIT: statistical.Crop.OILPALM.value, + taxonomy.CommodityId.OTHER_OILCROPS: statistical.Crop.OTHER_OILCROPS.value, + taxonomy.CommodityId.OTHER_PULSES: statistical.Crop.OTHER_PULSES.value, + taxonomy.CommodityId.PIGEON_PEAS: statistical.Crop.PIGEONPEA.value, + taxonomy.CommodityId.PLANTAINS: statistical.Crop.PLANTAIN.value, + taxonomy.CommodityId.POTATOES: statistical.Crop.POTATO.value, + taxonomy.CommodityId.RAPE_OR_COLZA_SEED: statistical.Crop.RAPESEED.value, + taxonomy.CommodityId.RICE: statistical.Crop.RICE.value, + taxonomy.CommodityId.SEED_COTTON: statistical.Crop.COTTON.value, + taxonomy.CommodityId.SESAME_SEED: statistical.Crop.SESAME_SEED.value, + taxonomy.CommodityId.SORGHUM: statistical.Crop.SORGHUM.value, + taxonomy.CommodityId.SOYA_BEANS: statistical.Crop.SOYBEAN.value, + taxonomy.CommodityId.SUGAR_BEET: statistical.Crop.SUGARBEET.value, + taxonomy.CommodityId.SUGAR_CANE: statistical.Crop.SUGARCANE.value, + taxonomy.CommodityId.SUNFLOWER_SEED: statistical.Crop.SUNFLOWER.value, + taxonomy.CommodityId.SWEET_POTATOES: statistical.Crop.SWEET_POTATO.value, + taxonomy.CommodityId.WHEAT: statistical.Crop.WHEAT.value, + taxonomy.CommodityId.YAMS: statistical.Crop.YAM.value, +} +assert set(COMMODITY_ID_TO_CROP_CODE) == { + taxonomy.METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[ + attribute.Methodology.STATISTICAL + ][crop.name] + for crop in statistical.Crop +} +# WRI publishes the two coffees and the two millets apart, as MapSPAM 2005 onward does +SPLIT_CROP_CODE_TO_COMBINED_CROP_CODE = { + statistical.Crop.ARABICA_COFFEE.value: ifpri_mapspam.Crop2000.COFF.name, + statistical.Crop.PEARL_MILLET.value: ifpri_mapspam.Crop2000.MILL.name, + statistical.Crop.ROBUSTA_COFFEE.value: ifpri_mapspam.Crop2000.COFF.name, + statistical.Crop.SMALL_MILLET.value: ifpri_mapspam.Crop2000.MILL.name, +} + + +def with_crop_codes(frame: pandas.DataFrame) -> pandas.DataFrame: + """Only the crop rows, each with the MapSPAM code it is compared under in `crop_code`. Some rows are not crops: the livestock rows, and the dropped and unattributed rows (see `emit.NonCommodity`). None has a MapSPAM code to key a join on, and neither WRI nor Orbae @@ -766,7 +816,12 @@ def keep_crop_rows(frame: pandas.DataFrame) -> pandas.DataFrame: stated as what is kept rather than what is removed, so a row added later is excluded by default rather than by being remembered here. """ - return frame[frame["commodity_name"].isin({crop.name for crop in statistical.Crop})] + crops = frame[frame["commodity_id"].isin(COMMODITY_ID_TO_CROP_CODE)].copy() + crops["crop_code"] = [ + COMMODITY_ID_TO_CROP_CODE[taxonomy.CommodityId(commodity_id)] + for commodity_id in crops["commodity_id"] + ] + return crops def iter_sluc_jdluc_comparisons() -> collections.abc.Iterator[dict[str, object]]: @@ -786,12 +841,14 @@ def iter_sluc_jdluc_comparisons() -> collections.abc.Iterator[dict[str, object]] efs = read_efs() if efs is None: return - national = efs[ - efs.index.get_level_values("admin_level") == schema.NATIONAL - ].reset_index() + national = with_crop_codes( + frame=efs[ + efs.index.get_level_values("admin_level") == schema.NATIONAL + ].reset_index() + ) by_methodology = national.pivot_table( columns="methodology", - index=["admin_id", "commodity_name"], + index=["admin_id", "crop_code"], values="emissions_factor_kgco2e_per_kg", ) # NB: the jurisdictional-direct leg is US-only, so a capture scoped to countries it does not @@ -816,7 +873,7 @@ def iter_sluc_jdluc_comparisons() -> collections.abc.Iterator[dict[str, object]] yield { "iso_3166": str(row["admin_id"]), # The MapSPAM code, because `get_comparisons` renames back to the canonical name once. - "crop_name": statistical.Crop[str(row["commodity_name"])].value, + "crop_name": str(row["crop_code"]), "emission_pool": schema.EmissionPool.TOTAL, "measure": targets.Measure.SLUC_OVER_JDLUC, "statistic": schema.Statistic.RATIO, @@ -865,21 +922,25 @@ def iter_sluc_wri_comparisons() -> collections.abc.Iterator[dict[str, object]]: efs = read_efs() if efs is None: return - deforestation = { - (row["iso_3166"], row["crop_name"]): float(row["deforestation_tonnes"]) - for row in read_wri_national_emissions(gas_scope="CO2e").to_dict("records") - if row["year"] == REFERENCE_YEAR and row["deforestation_tonnes"] > 0 - } + wri = read_wri_national_emissions(gas_scope="CO2e") + wri = wri[(wri["year"] == REFERENCE_YEAR) & (wri["deforestation_tonnes"] > 0)] + # Tonnes add, so WRI's two coffees and two millets sum to the one FAOSTAT item each that our + # rows carry + deforestation = ( + wri.assign( + crop_name=wri["crop_name"].replace(SPLIT_CROP_CODE_TO_COMBINED_CROP_CODE) + ) + .groupby(["iso_3166", "crop_name"])["deforestation_tonnes"] + .sum() + .to_dict() + ) levels = efs.index.get_level_values - national = keep_crop_rows( + national = with_crop_codes( frame=efs[ (levels("admin_level") == schema.NATIONAL) & (levels("methodology") == schema.STATISTICAL) ].reset_index() - ).copy() - national["crop_code"] = [ - statistical.Crop[name].value for name in national["commodity_name"] - ] + ) for row in sorted( national.to_dict("records"), key=lambda record: (record["admin_id"], record["crop_code"]), @@ -948,12 +1009,9 @@ def iter_orbae_capture_comparisons() -> collections.abc.Iterator[dict[str, objec orbae = orbae.dropna(subset=["admin_id"]) levels = efs.index.get_level_values - provincial = keep_crop_rows( + provincial = with_crop_codes( frame=efs[levels("admin_level") == schema.PROVINCIAL].reset_index() - ).copy() - provincial["crop_code"] = [ - statistical.Crop[name].value for name in provincial["commodity_name"] - ] + ) # `admin_id` leads with the ISO, which is the only place the country appears on our rows. provincial["iso_3166"] = provincial["admin_id"].str[:3] provincial["luc_only_kg_per_kg"] = ( diff --git a/validation/schema.py b/validation/schema.py index bf57abb..5c64a90 100644 --- a/validation/schema.py +++ b/validation/schema.py @@ -34,7 +34,7 @@ # efs.parquet's own index, and what any other capture is re-keyed onto. admin_id rather than # jurisdiction_name: a machine identifier, where a display string could be renamed. `methodology` # is a key level because both legs share one artifact, so it separates sLUC from jdLUC. -CANONICAL_KEY = ("admin_level", "admin_id", "commodity_name", "methodology") +CANONICAL_KEY = ("admin_level", "admin_id", "commodity_id", "methodology") NATIONAL = "NATIONAL" PROVINCIAL = "PROVINCIAL" # The two legs, as `methodology` records them in efs.parquet.