From 1f43d087234f834b08a5a38414061698d822e979 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 13:24:06 -0800 Subject: [PATCH 01/16] recover dataset version comments --- jdluc/datasets/gnw_tcl.py | 1 + jdluc/datasets/usda_nass_cdl.py | 1 + 2 files changed, 2 insertions(+) diff --git a/jdluc/datasets/gnw_tcl.py b/jdluc/datasets/gnw_tcl.py index 2a6b2df..84dc66a 100644 --- a/jdluc/datasets/gnw_tcl.py +++ b/jdluc/datasets/gnw_tcl.py @@ -37,5 +37,6 @@ def _save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: product_name="tree-cover-loss", save_tile_id_to_local_path=_save_tile_id_to_local_path, source_name="gnw", + # v1: include 2025 version="v1", ) diff --git a/jdluc/datasets/usda_nass_cdl.py b/jdluc/datasets/usda_nass_cdl.py index 52442ee..2a47f6e 100644 --- a/jdluc/datasets/usda_nass_cdl.py +++ b/jdluc/datasets/usda_nass_cdl.py @@ -98,6 +98,7 @@ def _save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: product_name="cdl", save_tile_id_to_local_path=_save_tile_id_to_local_path, source_name="usda-nass", + # v1: reproject data so that when it doesn't cover a full tile it doesn't get stretched version="v1", ) From 7e7ccaec607d2390220829bf0c6a41e33a773bc0 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 13:27:31 -0800 Subject: [PATCH 02/16] Ingest 2020-2025 for USDA/NASS/CDL so it can cover all assessment years of interest --- jdluc/datasets/usda_nass_cdl.py | 134 +++++++++++++++++++------------- jdluc/jurisdictional_direct.py | 7 +- 2 files changed, 87 insertions(+), 54 deletions(-) diff --git a/jdluc/datasets/usda_nass_cdl.py b/jdluc/datasets/usda_nass_cdl.py index 2a47f6e..5d7b001 100644 --- a/jdluc/datasets/usda_nass_cdl.py +++ b/jdluc/datasets/usda_nass_cdl.py @@ -21,76 +21,105 @@ import tempfile import zipfile -import numpy +import rasterio import rasterio.enums import rasterio.transform -import xarray +import rasterio.vrt +import rasterio.windows -from jdluc import geo, tiling, utils +from jdluc import tiling, utils from jdluc.datasets import base logger = logging.getLogger(__name__) -YEAR = 2020 +YEARS = tuple(range(2020, 2026)) +BAND_NAMES = [f"{year=:d}" for year in YEARS] -# All tiles share one (expensive) download+reproject, so cache it +# A shorter warp reads the source more times; a taller one holds more of the tile in memory +DESTINATION_BLOCK_ROWS = 2048 + + +# All tiles share each year's (expensive) download, so cache it @utils.threadsafe_cache -def _get_dataarray() -> xarray.DataArray: - import rioxarray - - with tempfile.TemporaryDirectory() as local_dir: - path_to_zip = os.path.join(local_dir, "data.zip") - utils.save_remote_url_to_local_path( - local_path=path_to_zip, - params={}, - remote_url=f"https://www.nass.usda.gov/Research_and_Science/Cropland/Release/datasets/{YEAR:d}_30m_cdls.zip", - ) - logger.info(f"Extracting .tif from {path_to_zip=:s}") - path_to_geotiff: str | None = None - with zipfile.ZipFile(file=path_to_zip) as zf: - for name in zf.namelist(): - if name.endswith(f"{YEAR:d}_30m_cdls.tif"): - path_to_geotiff = zf.extract(name, local_dir) - break - assert path_to_geotiff is not None - logger.info(f"Loading {path_to_geotiff=:s} into xarray") - darray = rioxarray.open_rasterio( - path_to_geotiff, - chunks=geo.get_chunk_size(dtypes=[numpy.dtype("uint8")]), - ) - assert isinstance(darray, xarray.DataArray) - logger.info(f"Reprojecting from {darray.rio.crs.to_epsg():d} to 4326") - return ( - darray.isel(band=0) - .drop_vars("band") - .rio.reproject( - dst_crs=4326, - nodata=DATASET.no_data, - resampling=rasterio.enums.Resampling.nearest, - ) - ) +def get_path_to_geotiff_for_year(year: int) -> str: + """The path to `year`'s source GeoTIFF, extracted from its zip to a local directory. + + The zip's external overviews stay behind, so a warp can only read full resolution. + """ + # NB: deliberately not a TemporaryDirectory, so that the GeoTIFF outlives this call + local_dir = tempfile.mkdtemp() + path_to_zip = os.path.join(local_dir, "data.zip") + utils.save_remote_url_to_local_path( + local_path=path_to_zip, + params={}, + remote_url=f"https://www.nass.usda.gov/sites/default/files/Research_and_Science/Cropland/Release/datasets/{year:d}_30m_cdls.zip", + ) + logger.info(f"Extracting {year=:d}'s .tif from {path_to_zip=:s}") + basename = f"{year:d}_30m_cdls.tif" + with zipfile.ZipFile(file=path_to_zip) as zf: + zf.extract(basename, local_dir) + os.remove(path_to_zip) + return os.path.join(local_dir, basename) def _save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: - darray = _get_dataarray() max_lat, min_lon = tiling.get_lat_lon_for_tile_id(tile_id=tile_id) - tiled = darray.rio.reproject( - dst_crs=4326, - nodata=DATASET.no_data, - resampling=rasterio.enums.Resampling.nearest, - shape=(tiling.TileResolution.GLAD.x, tiling.TileResolution.GLAD.y), - transform=rasterio.transform.from_origin( - min_lon, max_lat, *tiling.TileResolution.GLAD.degrees_per_pixel - ), + transform = rasterio.transform.from_origin( + min_lon, max_lat, *tiling.TileResolution.GLAD.degrees_per_pixel ) - logger.info(f"Saving xarray to {local_path=:s} for {tile_id=:s}") - tiled.rio.to_raster(local_path, driver="GTiff", recalc_transform=False) + with rasterio.open( + BIGTIFF="YES", + NUM_THREADS="ALL_CPUS", + compress="deflate", + count=len(YEARS), + crs=4326, + driver="GTiff", + dtype=DATASET.dtype, + fp=local_path, + height=tiling.TileResolution.GLAD.y, + interleave="band", + mode="w", + nodata=DATASET.no_data, + tiled=True, + transform=transform, + width=tiling.TileResolution.GLAD.x, + ) as dataset: + for band_idx, year in enumerate(YEARS, start=1): + path_to_geotiff = get_path_to_geotiff_for_year(year=year) + logger.info(f"Warping {path_to_geotiff=:s} onto {tile_id=:s}") + with ( + rasterio.open(fp=path_to_geotiff) as source, + rasterio.vrt.WarpedVRT( + source, + crs=4326, + height=tiling.TileResolution.GLAD.y, + nodata=DATASET.no_data, + resampling=rasterio.enums.Resampling.nearest, + src_nodata=DATASET.no_data, + transform=transform, + width=tiling.TileResolution.GLAD.x, + ) as warped, + ): + for top in range( + 0, tiling.TileResolution.GLAD.y, DESTINATION_BLOCK_ROWS + ): + window = rasterio.windows.Window( + col_off=0, + height=min( + DESTINATION_BLOCK_ROWS, tiling.TileResolution.GLAD.y - top + ), + row_off=top, + width=tiling.TileResolution.GLAD.x, + ) + dataset.write( + warped.read(1, window=window), band_idx, window=window + ) DATASET = base.RasterDataset( - band_names=["crop-class"], + band_names=BAND_NAMES, band_type=base.BandType.CATEGORICAL, dtype="uint8", no_data=0, @@ -99,7 +128,8 @@ def _save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: save_tile_id_to_local_path=_save_tile_id_to_local_path, source_name="usda-nass", # v1: reproject data so that when it doesn't cover a full tile it doesn't get stretched - version="v1", + # v2: a band per year, 2020-2025, each warped straight from the source projection onto the tile + version="v2", ) diff --git a/jdluc/jurisdictional_direct.py b/jdluc/jurisdictional_direct.py index f7542b7..fbe3d23 100644 --- a/jdluc/jurisdictional_direct.py +++ b/jdluc/jurisdictional_direct.py @@ -38,6 +38,9 @@ class Crop(enum.Enum): DATASET_NAMES = (DatasetName.USDA_NASS_CDL,) +CROP_CLASS_BAND_NAME = usda_nass_cdl.DATASET.fully_qualified_band_names[ + usda_nass_cdl.YEARS.index(emit.ASSESSMENT_YEAR) +] EMISSIONS_AND_CROP_VARIABLE_NAMES = [ @@ -50,7 +53,7 @@ class Crop(enum.Enum): f"{emit.PASTURELAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha", *(f"{component!s}-emissions:tco2e-per-ha" for component in emit.EmissionComponent), # Crop - usda_nass_cdl.DATASET.fully_qualified_band_name, + CROP_CLASS_BAND_NAME, ] @@ -207,7 +210,7 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]: component: clipped[f"{component!s}-emissions:tco2e-per-ha"] for component in emit.EmissionComponent }, - crop_class=clipped[usda_nass_cdl.DATASET.fully_qualified_band_name], + crop_class=clipped[CROP_CLASS_BAND_NAME], cropland_occupation_per_hectare=clipped[ f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha" ], From 0c06fd990e8ec30bbc22430fd17f2e3a47cf96d7 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 14:59:34 -0800 Subject: [PATCH 03/16] Add --assessment-year to jdluc.emit --- jdluc/__tests__/emit_test.py | 127 +++++++++++++++++++----- jdluc/datasets/gnw_tcl.py | 1 + jdluc/emit.py | 138 ++++++++++++++++++--------- jdluc/jurisdictional_direct.py | 9 +- jdluc/jurisdictional_direct_proxy.py | 7 +- jdluc/statistical.py | 18 +++- tools/build-soy-brazil-figure.py | 4 +- 7 files changed, 223 insertions(+), 81 deletions(-) diff --git a/jdluc/__tests__/emit_test.py b/jdluc/__tests__/emit_test.py index b22548b..5ede0ef 100644 --- a/jdluc/__tests__/emit_test.py +++ b/jdluc/__tests__/emit_test.py @@ -20,12 +20,11 @@ from jdluc.datasets.liao_gaced30 import YEARS as CROPLAND_YEARS from jdluc.datasets.liao_gaced30 import Cropland from jdluc.emit import ( - ASSESSMENT_YEAR, CARBON_PER_BIOMASS_LIVE_WOOD, CLIMATE_ZONE_TO_SOC_RETENTION_FRACTION, CO2E_PER_CARBON, + DEFAULT_ASSESSMENT_YEAR, FROM_FOREST, - LOOKBACK_YEARS_RANGE, PEATLAND_EMISSIONS_ANNUAL_TCO2E_PER_HA, PEATLAND_EMISSIONS_PULSE_TCO2E_PER_HA, Conversion, @@ -45,10 +44,13 @@ get_last_departure_year, get_linear_discounted, get_linear_discounted_values, + get_lookback_years_range, get_mineral_soil_emissions, get_peatland_occupation_emissions, ) +LOOKBACK_YEARS_RANGE = get_lookback_years_range(assessment_year=DEFAULT_ASSESSMENT_YEAR) + def get_darray_for_data( data: collections.abc.Sequence[collections.abc.Sequence[int | float]], @@ -187,7 +189,7 @@ def mask(members: tuple[Conversion, ...]) -> xarray.DataArray: ), to_cropland=mask(to_cropland), to_pasture=mask(to_pasture), - year=get_darray_for_data(data=[[ASSESSMENT_YEAR] * width]), + year=get_darray_for_data(data=[[DEFAULT_ASSESSMENT_YEAR] * width]), ) @@ -226,7 +228,7 @@ def mask(value: bool) -> xarray.DataArray: from_rangeland=mask(False), to_cropland=mask(to_cropland), to_pasture=mask(False), - year=get_darray_for_data(data=[[ASSESSMENT_YEAR]]), + year=get_darray_for_data(data=[[DEFAULT_ASSESSMENT_YEAR]]), ) for from_forest, to_cropland in ((True, False), (False, True)): @@ -445,14 +447,32 @@ def test_get_dset_for_output_appends_each_band_its_units() -> None: } -def test_get_linear_discounted() -> None: - # Each year takes its own factor, from 0.25% in 2001 to 9.75% in 2020, and a pixel with no year - # takes none even where its value is missing - result = get_linear_discounted( - conversion_year=get_darray_for_data( - data=[[2001, 2005, 2006, 2011, 2020, 0, numpy.nan]] +@pytest.mark.parametrize( + ("assessment_year", "conversion_years", "values"), + ( + pytest.param( + 2020, + [2001, 2005, 2006, 2011, 2020, 0, numpy.nan], + [8, 8, 8, 8, 8, numpy.nan, numpy.nan], + id="a pixel with no year takes no factor even where its value is missing", + ), + pytest.param( + 2024, + [2005, 2009, 2010, 2015, 2024, 2004, numpy.nan], + [8, 8, 8, 8, 8, numpy.nan, numpy.nan], + id="the factors move with the assessment year, leaving its opening year out", ), - darray=get_darray_for_data(data=[[8, 8, 8, 8, 8, numpy.nan, numpy.nan]]), + ), +) +def test_get_linear_discounted( + assessment_year: int, conversion_years: list[float], values: list[float] +) -> None: + # Each year takes its own factor, from 0.25% nineteen years before the assessment year to 9.75% + # in it + result = get_linear_discounted( + assessment_year=assessment_year, + conversion_year=get_darray_for_data(data=[conversion_years]), + darray=get_darray_for_data(data=[values]), ) numpy.testing.assert_allclose(result.data, [[0.02, 0.18, 0.22, 0.42, 0.78, 0, 0]]) @@ -472,6 +492,7 @@ def test_get_linear_discounted_values_rejects_emissions_it_cannot_date( ) -> None: with pytest.raises(AssertionError): get_linear_discounted_values( + assessment_year=2020, conversion_years=numpy.array([2012, conversion_year]), values=numpy.array([1.0, 5.0]), ) @@ -481,6 +502,7 @@ def test_get_linear_discounted_rejects_inputs_on_different_grids() -> None: conversion_year = get_darray_for_data(data=[[2012, 2013]]) with pytest.raises(ValueError): get_linear_discounted( + assessment_year=2020, conversion_year=conversion_year, darray=get_darray_for_data(data=[[1, 1]]).assign_coords( x=conversion_year.x + 1 @@ -613,6 +635,7 @@ def test_get_last_departure_year(is_rangeland_by_year: list[bool], year: int) -> def get_dset_for_one_pixel( + assessment_year: int, grassland_by_year: list[float], is_cocoa: bool, is_coffee: bool, @@ -625,19 +648,21 @@ def get_dset_for_one_pixel( GRASSLAND.fully_qualified_band_names[ GRASSLAND_YEARS.index(grassland_year) ]: get_darray_for_data(data=[[grassland_by_year[offset]]]) - for offset, grassland_year in enumerate(LOOKBACK_YEARS_RANGE) + for offset, grassland_year in enumerate( + get_lookback_years_range(assessment_year=assessment_year) + ) } | { CROPLAND.fully_qualified_band_names[ - CROPLAND_YEARS.index(ASSESSMENT_YEAR) + CROPLAND_YEARS.index(assessment_year) ]: get_darray_for_data( data=[[Cropland.CROPLAND if is_cropland else Cropland.NOT_CROPLAND]] ), COCOA.fully_qualified_band_names[ - FDAP_YEARS.index(ASSESSMENT_YEAR) + FDAP_YEARS.index(assessment_year) ]: get_darray_for_data(data=[[int(is_cocoa)]]), COFFEE.fully_qualified_band_names[ - FDAP_YEARS.index(ASSESSMENT_YEAR) + FDAP_YEARS.index(assessment_year) ]: get_darray_for_data(data=[[int(is_coffee)]]), OIL_PALM.fully_qualified_band_name: get_darray_for_data( data=[[planting_year]] @@ -819,9 +844,9 @@ def get_dset_for_one_pixel( pytest.param( NO_GRASSLAND, True, - ASSESSMENT_YEAR, + DEFAULT_ASSESSMENT_YEAR, Conversion.FOREST_TO_CROPLAND, - ASSESSMENT_YEAR, + DEFAULT_ASSESSMENT_YEAR, id="loss in the assessment year is inside the lookback", ), pytest.param( @@ -874,14 +899,64 @@ def test_get_conversion_record( year: int, ) -> None: result = get_conversion_record( + assessment_year=DEFAULT_ASSESSMENT_YEAR, dset=get_dset_for_one_pixel( + assessment_year=DEFAULT_ASSESSMENT_YEAR, grassland_by_year=grassland_by_year, is_cocoa=False, is_coffee=False, is_cropland=is_cropland, loss_year=loss_year, planting_year=0, - ) + ), + ) + assert result.conversion.data == [[conversion]] + assert result.year.data == [[year]] + + +@pytest.mark.parametrize( + ("loss_year", "conversion", "year"), + ( + pytest.param( + 2004, + Conversion.NONE, + 0, + id="loss in the opening year, after the layer's base year, is outside the lookback", + ), + pytest.param( + 2005, + Conversion.FOREST_TO_CROPLAND, + 2005, + id="the year after the opening year is inside it", + ), + pytest.param( + 2024, + Conversion.FOREST_TO_CROPLAND, + 2024, + id="and so is the assessment year", + ), + pytest.param( + 2025, + Conversion.NONE, + 0, + id="loss after the assessment year is outside it", + ), + ), +) +def test_get_conversion_record_in_a_later_assessment_year( + loss_year: int, conversion: Conversion, year: int +) -> None: + result = get_conversion_record( + assessment_year=2024, + dset=get_dset_for_one_pixel( + assessment_year=2024, + grassland_by_year=NO_GRASSLAND, + is_cocoa=False, + is_coffee=False, + is_cropland=True, + loss_year=loss_year, + planting_year=0, + ), ) assert result.conversion.data == [[conversion]] assert result.year.data == [[year]] @@ -962,7 +1037,7 @@ def test_get_conversion_record( False, False, False, - ASSESSMENT_YEAR, + DEFAULT_ASSESSMENT_YEAR, DestinationDataset.DESCALS_OIL_PALM, True, False, @@ -973,7 +1048,7 @@ def test_get_conversion_record( False, False, False, - ASSESSMENT_YEAR + 1, + DEFAULT_ASSESSMENT_YEAR + 1, 0, False, False, @@ -1028,7 +1103,9 @@ def test_get_conversion_record_destination( ) -> None: """A destination keeps every bit that fired, but resolves to the earliest one.""" result = get_conversion_record( + assessment_year=DEFAULT_ASSESSMENT_YEAR, dset=get_dset_for_one_pixel( + assessment_year=DEFAULT_ASSESSMENT_YEAR, grassland_by_year=[ *NO_GRASSLAND[:-1], float(Grassland.CULTIVATED if is_cultivated else Grassland.OTHER), @@ -1038,7 +1115,7 @@ def test_get_conversion_record_destination( is_cropland=is_cropland, loss_year=2012, planting_year=planting_year, - ) + ), ) assert result.destination_dataset.data == [[destination_dataset]] assert result.to_cropland.data == [[to_cropland]] @@ -1148,14 +1225,16 @@ def test_the_halves_a_conversion_record_resolves( # The halves outlive the conversion they resolve into: one without the other names no # conversion, but each still drives something -- the dropped total and peat occupation. result = get_conversion_record( + assessment_year=DEFAULT_ASSESSMENT_YEAR, dset=get_dset_for_one_pixel( + assessment_year=DEFAULT_ASSESSMENT_YEAR, grassland_by_year=grassland_by_year, is_cocoa=False, is_coffee=False, is_cropland=is_cropland, loss_year=loss_year, planting_year=0, - ) + ), ) assert bool(result.from_forest) == from_forest assert bool(result.from_rangeland | result.from_pasture) == from_grassland @@ -1194,14 +1273,16 @@ def test_a_conversion_record_tells_pasture_from_rangeland_sources( ) -> None: # Pasture and rangeland sources carry different vegetation carbon, so the record keeps them apart result = get_conversion_record( + assessment_year=DEFAULT_ASSESSMENT_YEAR, dset=get_dset_for_one_pixel( + assessment_year=DEFAULT_ASSESSMENT_YEAR, grassland_by_year=grassland_by_year, is_cocoa=False, is_coffee=False, is_cropland=True, loss_year=loss_year, planting_year=0, - ) + ), ) assert bool(result.from_pasture) == from_pasture assert bool(result.from_rangeland) == from_rangeland diff --git a/jdluc/datasets/gnw_tcl.py b/jdluc/datasets/gnw_tcl.py index 84dc66a..3911643 100644 --- a/jdluc/datasets/gnw_tcl.py +++ b/jdluc/datasets/gnw_tcl.py @@ -15,6 +15,7 @@ from jdluc.datasets import base LOSS_YEAR_OFFSET = 2000 +LOSS_YEARS = tuple(range(LOSS_YEAR_OFFSET + 1, 2026)) def _save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: diff --git a/jdluc/emit.py b/jdluc/emit.py index ac75789..09bb8fc 100644 --- a/jdluc/emit.py +++ b/jdluc/emit.py @@ -17,11 +17,10 @@ Example invocations: uv run python jdluc/emit.py USA - uv run python jdluc/emit.py --backfill + uv run python jdluc/emit.py --assessment-year 2024 --backfill """ import argparse -import collections.abc import dataclasses import enum import logging @@ -228,19 +227,40 @@ def get_mineral_soil_emissions( LOOKBACK_YEARS = 20 -ASSESSMENT_YEAR = 2020 -LOOKBACK_YEARS_RANGE = tuple( - range(ASSESSMENT_YEAR - LOOKBACK_YEARS, ASSESSMENT_YEAR + 1) -) +# Bounded by the years the input datasets cover; see "Source datasets" in docs/methodology.md +ASSESSMENT_YEARS = (2020, 2022, 2024) +DEFAULT_ASSESSMENT_YEAR = 2020 +assert DEFAULT_ASSESSMENT_YEAR in ASSESSMENT_YEARS -# The LSR Guidance's linear amortization factor for a conversion in each year (Table 7.11): 9.75% in -# the assessment year, falling 0.5 points a year to 0.25% nineteen years before it -YEAR_TO_LINEAR_DISCOUNT_WEIGHT: dict[int, float] = { - year: (39 - 2 * (ASSESSMENT_YEAR - year)) / 400 - # NB: nothing can depart in the lookback's first year, which has no year before it to compare with - for year in LOOKBACK_YEARS_RANGE[1:] -} -assert math.isclose(sum(YEAR_TO_LINEAR_DISCOUNT_WEIGHT.values()), 1) + +def get_lookback_years_range(assessment_year: int) -> tuple[int, ...]: + assert assessment_year in ASSESSMENT_YEARS, f"{assessment_year=:d} is unsupported" + return tuple(range(assessment_year - LOOKBACK_YEARS, assessment_year + 1)) + + +def get_year_to_linear_discount_weight(assessment_year: int) -> dict[int, float]: + """The LSR Guidance's linear amortization factor for a conversion in each year (Table 7.11): + 9.75% in the assessment year, falling 0.5 points a year to 0.25% nineteen years before it. + """ + year_to_weight = { + year: (39 - 2 * (assessment_year - year)) / 400 + # NB: nothing can depart in the lookback's first year, which has no year before it to + # compare with + for year in get_lookback_years_range(assessment_year=assessment_year)[1:] + } + assert math.isclose(sum(year_to_weight.values()), 1) + return year_to_weight + + +assert all( + set(get_lookback_years_range(assessment_year=assessment_year)[1:]) + <= set(gnw_tcl.LOSS_YEARS) + and set(get_lookback_years_range(assessment_year=assessment_year)) + <= set(gpw_grassland.YEARS) + and assessment_year in liao_gaced30.YEARS + and assessment_year in fdap.YEARS + for assessment_year in ASSESSMENT_YEARS +) class Conversion(enum.IntEnum): @@ -289,45 +309,49 @@ class DestinationDataset(enum.IntFlag): GPW_GRASSLAND = enum.auto() -DESTINATION_DATASET_TO_PREDICATE: dict[ - DestinationDataset, collections.abc.Callable[[xarray.Dataset], xarray.DataArray] -] = { - DestinationDataset.DESCALS_OIL_PALM: lambda dset: ( +class DestinationPredicate(typing.Protocol): + def __call__( + self, assessment_year: int, dset: xarray.Dataset + ) -> xarray.DataArray: ... + + +DESTINATION_DATASET_TO_PREDICATE: dict[DestinationDataset, DestinationPredicate] = { + DestinationDataset.DESCALS_OIL_PALM: lambda assessment_year, dset: ( # NB: zero is "no oil palm" ( (planting_year := dset[descals_oil_palm.DATASET.fully_qualified_band_name]) > 0 ) - & (planting_year <= ASSESSMENT_YEAR) + & (planting_year <= assessment_year) ), - DestinationDataset.FDAP_COCOA: lambda dset: ( + DestinationDataset.FDAP_COCOA: lambda assessment_year, dset: ( dset[ fdap.COCOA_DATASET.fully_qualified_band_names[ - fdap.YEARS.index(ASSESSMENT_YEAR) + fdap.YEARS.index(assessment_year) ] ] == 1 ), - DestinationDataset.FDAP_COFFEE: lambda dset: ( + DestinationDataset.FDAP_COFFEE: lambda assessment_year, dset: ( dset[ fdap.COFFEE_DATASET.fully_qualified_band_names[ - fdap.YEARS.index(ASSESSMENT_YEAR) + fdap.YEARS.index(assessment_year) ] ] == 1 ), - DestinationDataset.LIAO_GACED30: lambda dset: ( + DestinationDataset.LIAO_GACED30: lambda assessment_year, dset: ( dset[ liao_gaced30.DATASET.fully_qualified_band_names[ - liao_gaced30.YEARS.index(ASSESSMENT_YEAR) + liao_gaced30.YEARS.index(assessment_year) ] ] == liao_gaced30.Cropland.CROPLAND ), - DestinationDataset.GPW_GRASSLAND: lambda dset: ( + DestinationDataset.GPW_GRASSLAND: lambda assessment_year, dset: ( dset[ gpw_grassland.DATASET.fully_qualified_band_names[ - gpw_grassland.YEARS.index(ASSESSMENT_YEAR) + gpw_grassland.YEARS.index(assessment_year) ] ] == gpw_grassland.Grassland.CULTIVATED @@ -355,7 +379,8 @@ class DestinationDataset(enum.IntFlag): def get_last_departure_year(is_source: xarray.DataArray) -> xarray.DataArray: """The last year a pixel left the class it held the year before, or 0 where it never did.""" assert "year" in is_source.coords - assert tuple(is_source.year.values) == LOOKBACK_YEARS_RANGE + # NB: the shift compares each year with the one before, so no year may be missing + assert (numpy.diff(is_source.year.values) == 1).all() departed = is_source.shift(fill_value=False, year=1) & ~is_source return departed.year.where(departed).max(dim="year").fillna(0) @@ -379,8 +404,12 @@ class ConversionRecord: year: xarray.DataArray -def get_conversion_record(dset: xarray.Dataset) -> ConversionRecord: +def get_conversion_record( + assessment_year: int, dset: xarray.Dataset +) -> ConversionRecord: """Which conversion fired on each pixel, and the year its source class ended.""" + lookback_years_range = get_lookback_years_range(assessment_year=assessment_year) + opening_year, *_ = lookback_years_range grassland = xarray.concat( [ dset[band_name] @@ -389,9 +418,9 @@ def get_conversion_record(dset: xarray.Dataset) -> ConversionRecord: gpw_grassland.DATASET.fully_qualified_band_names, strict=True, ) - if year in LOOKBACK_YEARS_RANGE + if year in lookback_years_range ], - dim=xarray.DataArray(list(LOOKBACK_YEARS_RANGE), dims="year", name="year"), + dim=xarray.DataArray(list(lookback_years_range), dims="year", name="year"), ) rangeland_departure_year = get_last_departure_year( @@ -411,17 +440,18 @@ def get_conversion_record(dset: xarray.Dataset) -> ConversionRecord: from_forest = ( # NB: Zero means no loss rather than a loss in the base year (year_of_loss > gnw_tcl.LOSS_YEAR_OFFSET) - # NB: the same bound as the line above while the lookback opens at the offset, so no test - # can tell them apart; they stop coinciding the moment the assessment year moves - & (year_of_loss > LOOKBACK_YEARS_RANGE[0]) - & (year_of_loss <= ASSESSMENT_YEAR) + # NB: the same bound as the line above for 2020, whose lookback opens at the offset, and + # the tighter one for every later year + & (year_of_loss > opening_year) + & (year_of_loss <= assessment_year) ) from_rangeland = ~from_forest & (rangeland_departure_year > pasture_departure_year) from_pasture = ~from_forest & (pasture_departure_year > rangeland_departure_year) # NB: destinations are flags, so we can superpose them to track intersections destination_dataset = sum( - predicate(dset) * numpy.uint8(destination_dataset) + predicate(assessment_year=assessment_year, dset=dset) + * numpy.uint8(destination_dataset) for destination_dataset, predicate in DESTINATION_DATASET_TO_PREDICATE.items() ) assert isinstance(destination_dataset, xarray.DataArray) @@ -555,16 +585,18 @@ def get_peatland_occupation_emissions( def get_linear_discounted_values( - conversion_years: numpy.ndarray, values: numpy.ndarray + conversion_years: numpy.ndarray, values: numpy.ndarray, *, assessment_year: int ) -> numpy.ndarray: """`values` times the discount weight of each one's conversion year, or zero where the year is zero (no conversion).""" indices = conversion_years.astype(numpy.int64) assert ( - (indices == conversion_years) & (indices >= 0) & (indices <= ASSESSMENT_YEAR) - ).all(), f"Conversion years must be whole years from 0 to {ASSESSMENT_YEAR:d}" - lookup = numpy.zeros(ASSESSMENT_YEAR + 1, dtype=numpy.float32) - for year, weight in YEAR_TO_LINEAR_DISCOUNT_WEIGHT.items(): + (indices == conversion_years) & (indices >= 0) & (indices <= assessment_year) + ).all(), f"Conversion years must be whole years from 0 to {assessment_year:d}" + lookup = numpy.zeros(assessment_year + 1, dtype=numpy.float32) + for year, weight in get_year_to_linear_discount_weight( + assessment_year=assessment_year + ).items(): lookup[year] = weight weights = lookup[indices] # NB: emissions without a weight -- no year, or one before the lookback -- would be zeroed @@ -576,7 +608,7 @@ def get_linear_discounted_values( def get_linear_discounted( - conversion_year: xarray.DataArray, darray: xarray.DataArray + assessment_year: int, conversion_year: xarray.DataArray, darray: xarray.DataArray ) -> xarray.DataArray: """`darray` times the discount weight of each pixel's conversion year, or zero without one.""" return xarray.apply_ufunc( @@ -586,6 +618,7 @@ def get_linear_discounted( dask="parallelized", # NB: arithmetic would align on an inner join, silently dropping pixels the grids don't share join="exact", + kwargs={"assessment_year": assessment_year}, output_dtypes=[numpy.result_type(darray.dtype, numpy.float32)], ).rename(darray.name) @@ -697,9 +730,12 @@ def merge_name_units(name: str, units: typing.Hashable | None) -> str: ) -@storage.cache_to_zarr(version=2) -def workflow(tile_id: str) -> xarray.Dataset: - logger.info(f"Running the land conversion and emissions worflow for {tile_id=:s}") +@storage.cache_to_zarr(version=0) +def workflow(assessment_year: int, tile_id: str) -> xarray.Dataset: + logger.info( + f"Running the land conversion and emissions worflow for {assessment_year=:d} and " + f"{tile_id=:s}" + ) dset = geo.exact_merge( harmonize.workflow( dataset_names=harmonize.Stack.LUC_AND_EMISSIONS.value, @@ -718,7 +754,9 @@ def workflow(tile_id: str) -> xarray.Dataset: ) logger.info("Resolving which conversion each pixel underwent, and when") - conversion_record: ConversionRecord = get_conversion_record(dset=dset) + conversion_record: ConversionRecord = get_conversion_record( + assessment_year=assessment_year, dset=dset + ) logger.info("Quantifying the carbon each source class holds") aboveground_biomass = dset[gnw_harris_agb.DATASET.fully_qualified_band_name] @@ -802,6 +840,12 @@ def main() -> int: nargs=argparse.ZERO_OR_MORE, type=worldbank_jurisdictions.iso_3166_str, ) + parser.add_argument( + "--assessment-year", + choices=ASSESSMENT_YEARS, + default=DEFAULT_ASSESSMENT_YEAR, + type=int, + ) parser.add_argument("--backfill", action="store_true", help="cover all GNW tiles") parser.add_argument( "--modulus", @@ -831,7 +875,7 @@ def main() -> int: ) ), ): - workflow(tile_id=tile_id) + workflow(assessment_year=int(args.assessment_year), tile_id=tile_id) return 0 diff --git a/jdluc/jurisdictional_direct.py b/jdluc/jurisdictional_direct.py index fbe3d23..d43554e 100644 --- a/jdluc/jurisdictional_direct.py +++ b/jdluc/jurisdictional_direct.py @@ -39,7 +39,7 @@ class Crop(enum.Enum): DATASET_NAMES = (DatasetName.USDA_NASS_CDL,) CROP_CLASS_BAND_NAME = usda_nass_cdl.DATASET.fully_qualified_band_names[ - usda_nass_cdl.YEARS.index(emit.ASSESSMENT_YEAR) + usda_nass_cdl.YEARS.index(emit.DEFAULT_ASSESSMENT_YEAR) ] @@ -158,7 +158,7 @@ def workflow( ), # NB: this should call the harmonize workflow with identical args and hit # the cache from the preceding call - emit.workflow(tile_id=tile_id), + emit.workflow(assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, tile_id=tile_id), # CDL harmonize.workflow( dataset_names=DATASET_NAMES, @@ -171,7 +171,9 @@ def workflow( merged = merged.assign( { f"{component!s}-emissions:tco2e-per-ha": emit.get_linear_discounted( - conversion_year=merged["conversion-year"], darray=emissions + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + conversion_year=merged["conversion-year"], + darray=emissions, ) for component, emissions in emit.get_component_to_emissions( dset=merged @@ -179,6 +181,7 @@ def workflow( } | { "dropped-emissions:tco2e-per-ha": emit.get_linear_discounted( + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, conversion_year=merged["conversion-year"], darray=merged["dropped-emissions:tco2e-per-ha"], ) diff --git a/jdluc/jurisdictional_direct_proxy.py b/jdluc/jurisdictional_direct_proxy.py index 930d71a..3537f6a 100644 --- a/jdluc/jurisdictional_direct_proxy.py +++ b/jdluc/jurisdictional_direct_proxy.py @@ -42,7 +42,9 @@ class the pixel became: the shares that divide it come from one year, so no span """ component_to_emissions = { component: emit.get_linear_discounted( - conversion_year=dset["conversion-year"], darray=emissions + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + conversion_year=dset["conversion-year"], + darray=emissions, ) for component, emissions in emit.get_component_to_emissions(dset=dset).items() } @@ -66,6 +68,7 @@ class the pixel became: the shares that divide it come from one year, so no span } | { "dropped-emissions:tco2e-per-ha": emit.get_linear_discounted( + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, conversion_year=dset["conversion-year"], darray=dset["dropped-emissions:tco2e-per-ha"], ) @@ -97,7 +100,7 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: ), # NB: this should call the harmonize workflow with identical args and hit # the cache from the preceding call - emit.workflow(tile_id=tile_id), + emit.workflow(assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, tile_id=tile_id), ) logger.info("Downsampling emissions from the GLAD grid to the CROPGRIDS grid") diff --git a/jdluc/statistical.py b/jdluc/statistical.py index b0d9651..057bab8 100644 --- a/jdluc/statistical.py +++ b/jdluc/statistical.py @@ -128,14 +128,16 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: ), # NB: this should call the harmonize workflow with identical args and hit # the cache from the preceding call - emit.workflow(tile_id=tile_id), + emit.workflow(assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, tile_id=tile_id), ) logger.info("Discounting emissions by component and splitting them by MapSPAM span") conversion_year = luc_and_emissions["conversion-year"] component_to_emissions = { component: emit.get_linear_discounted( - conversion_year=conversion_year, darray=emissions + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + conversion_year=conversion_year, + darray=emissions, ) for component, emissions in emit.get_component_to_emissions( dset=luc_and_emissions @@ -167,7 +169,9 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: logger.info("Discounting the source carbon no destination claimed") name = "dropped-emissions:tco2e-per-ha" luc_and_emissions[name] = emit.get_linear_discounted( - conversion_year=conversion_year, darray=luc_and_emissions[name] + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + conversion_year=conversion_year, + darray=luc_and_emissions[name], ) derived_variable_names.append(name) @@ -382,7 +386,9 @@ def get_commodity_to_peatland_conversion_share( MAPSPAM_SPANS = ((2000, 2005), (2005, 2010), (2010, 2020)) assert all( sum(before < year <= after for before, after in MAPSPAM_SPANS) == 1 - for year in emit.YEAR_TO_LINEAR_DISCOUNT_WEIGHT + for year in emit.get_year_to_linear_discount_weight( + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR + ) ) MAPSPAM_SNAPSHOT_YEARS = tuple( sorted({year for span in MAPSPAM_SPANS for year in span}) @@ -393,7 +399,9 @@ def get_commodity_to_peatland_conversion_share( SNAPSHOT_YEAR_TO_WEIGHT = { snapshot_year: sum( weight / 2 - for year, weight in emit.YEAR_TO_LINEAR_DISCOUNT_WEIGHT.items() + for year, weight in emit.get_year_to_linear_discount_weight( + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR + ).items() for before, after in MAPSPAM_SPANS if before < year <= after and snapshot_year in (before, after) ) diff --git a/tools/build-soy-brazil-figure.py b/tools/build-soy-brazil-figure.py index 9b60d32..8294101 100644 --- a/tools/build-soy-brazil-figure.py +++ b/tools/build-soy-brazil-figure.py @@ -82,7 +82,9 @@ Conversion = emit.Conversion # The lookback opens at the assessment year less twenty, and nothing can have departed in that # first year -- a departure needs the year before it to compare against -LOOKBACK_YEARS = emit.LOOKBACK_YEARS_RANGE[1:] +LOOKBACK_YEARS = emit.get_lookback_years_range( + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR +)[1:] # Hotter reads as more carbon released: forest sources red and orange, grassland tan and olive From f44e2dcefcf11c5a8b32e75cf900fe995d58052d Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 15:03:10 -0800 Subject: [PATCH 04/16] Add assessment_year to jurisdictional_direct.workflow --- jdluc/attribute.py | 8 +++++++- jdluc/jurisdictional_direct.py | 28 +++++++++++++++------------- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/jdluc/attribute.py b/jdluc/attribute.py index ac36074..3f4d298 100644 --- a/jdluc/attribute.py +++ b/jdluc/attribute.py @@ -25,11 +25,13 @@ import collections.abc import concurrent.futures import enum +import functools import logging import pandas from jdluc import ( + emit, jurisdictional_direct, jurisdictional_direct_proxy, statistical, @@ -64,9 +66,13 @@ def workflow( iso_3166s: collections.abc.Iterable[str], methodology: Methodology, ) -> pandas.DataFrame: + workflow_for_tile: collections.abc.Callable[..., pandas.DataFrame] match methodology: case Methodology.JURISDICTIONAL_DIRECT: - workflow_for_tile = jurisdictional_direct.workflow + workflow_for_tile = functools.partial( + jurisdictional_direct.workflow, + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + ) case Methodology.JURISDICTIONAL_DIRECT_PROXY: workflow_for_tile = jurisdictional_direct_proxy.workflow case Methodology.STATISTICAL: diff --git a/jdluc/jurisdictional_direct.py b/jdluc/jurisdictional_direct.py index d43554e..08b9a9e 100644 --- a/jdluc/jurisdictional_direct.py +++ b/jdluc/jurisdictional_direct.py @@ -38,12 +38,10 @@ class Crop(enum.Enum): DATASET_NAMES = (DatasetName.USDA_NASS_CDL,) -CROP_CLASS_BAND_NAME = usda_nass_cdl.DATASET.fully_qualified_band_names[ - usda_nass_cdl.YEARS.index(emit.DEFAULT_ASSESSMENT_YEAR) -] +assert set(emit.ASSESSMENT_YEARS) <= set(usda_nass_cdl.YEARS) -EMISSIONS_AND_CROP_VARIABLE_NAMES = [ +EMISSIONS_VARIABLE_NAMES = [ # Harmonized gnw_global_peatlands.DATASET.fully_qualified_band_name, # Emissions @@ -52,8 +50,6 @@ class Crop(enum.Enum): "hectares-per-pixel:ha", f"{emit.PASTURELAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha", *(f"{component!s}-emissions:tco2e-per-ha" for component in emit.EmissionComponent), - # Crop - CROP_CLASS_BAND_NAME, ] @@ -138,8 +134,9 @@ def get_commodity_name_to_totals( } -@storage.cache_to_parquet(version=2) +@storage.cache_to_parquet(version=0) def workflow( + assessment_year: int, commodity_names: tuple[str, ...], iso_3166: str, tile_id: str, @@ -147,7 +144,12 @@ def workflow( crops = tuple(Crop[commodity_name] for commodity_name in commodity_names) assert iso_3166 == "USA", "JD only supports USA today" - logger.info(f"Computing emissions for {crops=:} and {tile_id=:s}") + logger.info( + f"Computing emissions for {crops=:}, {assessment_year=:d} and {tile_id=:s}" + ) + crop_class_band_name = usda_nass_cdl.DATASET.fully_qualified_band_names[ + usda_nass_cdl.YEARS.index(assessment_year) + ] merged = geo.exact_merge( harmonize.workflow( dataset_names=harmonize.Stack.LUC_AND_EMISSIONS.value, @@ -158,7 +160,7 @@ def workflow( ), # NB: this should call the harmonize workflow with identical args and hit # the cache from the preceding call - emit.workflow(assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, tile_id=tile_id), + emit.workflow(assessment_year=assessment_year, tile_id=tile_id), # CDL harmonize.workflow( dataset_names=DATASET_NAMES, @@ -171,7 +173,7 @@ def workflow( merged = merged.assign( { f"{component!s}-emissions:tco2e-per-ha": emit.get_linear_discounted( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + assessment_year=assessment_year, conversion_year=merged["conversion-year"], darray=emissions, ) @@ -181,7 +183,7 @@ def workflow( } | { "dropped-emissions:tco2e-per-ha": emit.get_linear_discounted( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + assessment_year=assessment_year, conversion_year=merged["conversion-year"], darray=merged["dropped-emissions:tco2e-per-ha"], ) @@ -199,7 +201,7 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]: logger.info(f"Clipping to provincial geometry for {jurisdiction.id=:s}") try: clipped = geo.clip_dset( - dset=merged[EMISSIONS_AND_CROP_VARIABLE_NAMES], + dset=merged[[*EMISSIONS_VARIABLE_NAMES, crop_class_band_name]], geometry=jurisdiction.geometry, ) except NoDataInBounds as exc: @@ -213,7 +215,7 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]: component: clipped[f"{component!s}-emissions:tco2e-per-ha"] for component in emit.EmissionComponent }, - crop_class=clipped[CROP_CLASS_BAND_NAME], + crop_class=clipped[crop_class_band_name], cropland_occupation_per_hectare=clipped[ f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha" ], From 537d336717b916ac82f5b91d8addc63cd0f31613 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 15:05:46 -0800 Subject: [PATCH 05/16] Add assessment_year to jurisdictional_direct_proxy.workflow --- .../jurisdictional_direct_proxy_test.py | 10 +++++--- jdluc/attribute.py | 5 +++- jdluc/jurisdictional_direct_proxy.py | 23 ++++++++++++------- 3 files changed, 26 insertions(+), 12 deletions(-) diff --git a/jdluc/__tests__/jurisdictional_direct_proxy_test.py b/jdluc/__tests__/jurisdictional_direct_proxy_test.py index 45b0524..7312312 100644 --- a/jdluc/__tests__/jurisdictional_direct_proxy_test.py +++ b/jdluc/__tests__/jurisdictional_direct_proxy_test.py @@ -12,6 +12,7 @@ from jdluc.datasets.tang_cropgrids import Crop from jdluc.emit import ( CROPLAND, + DEFAULT_ASSESSMENT_YEAR, PASTURELAND, PEATLAND_EMISSIONS_ANNUAL_TCO2E_PER_HA, Conversion, @@ -122,9 +123,10 @@ def test_get_name_to_band_splits_each_component_by_land_class( conversion: Conversion, is_peat: bool, name_to_expected: dict[str, float] ) -> None: name_to_band = get_name_to_band( + assessment_year=DEFAULT_ASSESSMENT_YEAR, dset=get_dset_for_pixel( conversion=conversion, grassland=float(Grassland.OTHER), is_peat=is_peat - ) + ), ) # Every other land class and component is charged nothing assert { @@ -140,9 +142,10 @@ def test_get_name_to_band_charges_every_conversion_exactly_once( conversion: Conversion, is_peat: bool ) -> None: name_to_band = get_name_to_band( + assessment_year=DEFAULT_ASSESSMENT_YEAR, dset=get_dset_for_pixel( conversion=conversion, grassland=float(Grassland.OTHER), is_peat=is_peat - ) + ), ) charged = sum(name_to_band[name].item() for name in LAND_CLASS_COMPONENT_NAMES) assert charged == pytest.approx(0.0 if conversion == Conversion.NONE else 110.0) @@ -160,9 +163,10 @@ def test_get_name_to_band_charges_every_conversion_exactly_once( ) def test_get_name_to_band_pasture_fraction(grassland: float, expected: float) -> None: name_to_band = get_name_to_band( + assessment_year=DEFAULT_ASSESSMENT_YEAR, dset=get_dset_for_pixel( conversion=Conversion.NONE, grassland=grassland, is_peat=False - ) + ), ) assert name_to_band[f"pasture:fraction:{CROPGRIDS_YEAR:d}"].item() == expected diff --git a/jdluc/attribute.py b/jdluc/attribute.py index 3f4d298..9deffa1 100644 --- a/jdluc/attribute.py +++ b/jdluc/attribute.py @@ -74,7 +74,10 @@ def workflow( assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, ) case Methodology.JURISDICTIONAL_DIRECT_PROXY: - workflow_for_tile = jurisdictional_direct_proxy.workflow + workflow_for_tile = functools.partial( + jurisdictional_direct_proxy.workflow, + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + ) case Methodology.STATISTICAL: workflow_for_tile = statistical.workflow case _: diff --git a/jdluc/jurisdictional_direct_proxy.py b/jdluc/jurisdictional_direct_proxy.py index 3537f6a..c16c7b2 100644 --- a/jdluc/jurisdictional_direct_proxy.py +++ b/jdluc/jurisdictional_direct_proxy.py @@ -33,7 +33,9 @@ HERD_HEADS_COLUMN = statistical.get_herd_heads_column(year=tang_cropgrids.YEAR) -def get_name_to_band(dset: xarray.Dataset) -> dict[str, xarray.DataArray]: +def get_name_to_band( + assessment_year: int, dset: xarray.Dataset +) -> dict[str, xarray.DataArray]: """The 30 m bands this leg divides, each per hectare or a fraction of the pixel. Each emissions component is discounted by its pixel's conversion year and split by the land @@ -42,7 +44,7 @@ class the pixel became: the shares that divide it come from one year, so no span """ component_to_emissions = { component: emit.get_linear_discounted( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + assessment_year=assessment_year, conversion_year=dset["conversion-year"], darray=emissions, ) @@ -68,7 +70,7 @@ class the pixel became: the shares that divide it come from one year, so no span } | { "dropped-emissions:tco2e-per-ha": emit.get_linear_discounted( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + assessment_year=assessment_year, conversion_year=dset["conversion-year"], darray=dset["dropped-emissions:tco2e-per-ha"], ) @@ -88,7 +90,7 @@ class the pixel became: the shares that divide it come from one year, so no span @storage.cache_to_zarr(version=0) -def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: +def get_downscaled_luc_emissions(assessment_year: int, tile_id: str) -> xarray.Dataset: logger.info("Computing emissions on the GLAD grid") luc_and_emissions = geo.exact_merge( harmonize.workflow( @@ -100,7 +102,7 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: ), # NB: this should call the harmonize workflow with identical args and hit # the cache from the preceding call - emit.workflow(assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, tile_id=tile_id), + emit.workflow(assessment_year=assessment_year, tile_id=tile_id), ) logger.info("Downsampling emissions from the GLAD grid to the CROPGRIDS grid") @@ -118,7 +120,9 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: transform=grid.transform, width=grid.resolution.x, ) - for name, band in get_name_to_band(dset=luc_and_emissions).items() + for name, band in get_name_to_band( + assessment_year=assessment_year, dset=luc_and_emissions + ).items() } ) @@ -282,6 +286,7 @@ def get_commodity_name_to_totals( @storage.cache_to_parquet(version=0) def workflow( + assessment_year: int, commodity_names: tuple[str, ...], iso_3166: str, tile_id: str, @@ -297,9 +302,11 @@ def workflow( for commodity in commodities if isinstance(commodity, tang_cropgrids.Crop) ) - logger.info(f"Computing emissions for {commodities=:} and {tile_id=:s}") + logger.info( + f"Computing emissions for {commodities=:}, {assessment_year=:d} and {tile_id=:s}" + ) merged = geo.exact_merge( - get_downscaled_luc_emissions(tile_id=tile_id), + get_downscaled_luc_emissions(assessment_year=assessment_year, tile_id=tile_id), harmonize.workflow( dataset_names=DATASET_NAMES, ignore_missing_tiles=True, From 8524ecd844d4acd85d735c94b193f4356fa50df8 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 15:11:42 -0800 Subject: [PATCH 06/16] Add assessment_year to statistical.workflow --- jdluc/__tests__/statistical_test.py | 44 ++++++++++-- jdluc/__tests__/trace_test.py | 18 ++++- jdluc/attribute.py | 13 +--- jdluc/statistical.py | 103 +++++++++++++++++++--------- jdluc/trace.py | 9 ++- tools/measure-drift.py | 4 +- validation/capture.py | 4 +- 7 files changed, 140 insertions(+), 55 deletions(-) diff --git a/jdluc/__tests__/statistical_test.py b/jdluc/__tests__/statistical_test.py index bfed01f..90e2640 100644 --- a/jdluc/__tests__/statistical_test.py +++ b/jdluc/__tests__/statistical_test.py @@ -12,7 +12,6 @@ from jdluc.statistical import ( HERD_HEADS_COLUMNS, MAPSPAM_SPANS, - SNAPSHOT_YEAR_TO_WEIGHT, Commodity, Crop, Livestock, @@ -22,6 +21,8 @@ get_commodity_to_share, get_crop_to_area_share, get_livestock_to_grazing_share, + get_mapspam_span, + get_snapshot_year_to_weight, ) AREA = ifpri_mapspam.Quantity.PHYSICAL_AREA @@ -534,13 +535,43 @@ 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. -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 +@pytest.mark.parametrize( + ("assessment_year", "expected"), + ( + pytest.param( + 2020, + {2000: 0.03125, 2005: 0.125, 2010: 0.46875, 2020: 0.375}, + id="each year split between its span's snapshots", + ), + # 2005 alone splits between 2000 and 2005, and 2021-2024 give their 36% to 2020 whole + pytest.param( + 2024, + {2000: 0.00125, 2005: 0.045, 2010: 0.31875, 2020: 0.635}, + id="years after the last snapshot weigh on it alone", + ), + ), +) +def test_get_snapshot_year_to_weight( + assessment_year: int, expected: dict[int, float] +) -> None: + assert pytest.approx(expected) == get_snapshot_year_to_weight( + assessment_year=assessment_year ) +@pytest.mark.parametrize( + ("year", "span"), + ( + pytest.param(2001, (2000, 2005), id="spans are open below"), + pytest.param(2005, (2000, 2005), id="and closed above"), + pytest.param(2015, (2010, 2020), id="the decade without a 2015 snapshot"), + pytest.param(2024, (2010, 2020), id="years after the last snapshot"), + ), +) +def test_get_mapspam_span(year: int, span: tuple[int, int]) -> None: + assert get_mapspam_span(year=year) == span + + # 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({}) @@ -599,6 +630,9 @@ def get_shares( } ), occupation_shares=get_shares(shares=occupation_shares), + snapshot_year_to_weight=get_snapshot_year_to_weight( + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR + ), span_to_commodity_to_peatland_conversion_share=dict.fromkeys( MAPSPAM_SPANS, get_shares(shares=peatland_conversion_shares) ), diff --git a/jdluc/__tests__/trace_test.py b/jdluc/__tests__/trace_test.py index 8176125..a7ec0e2 100644 --- a/jdluc/__tests__/trace_test.py +++ b/jdluc/__tests__/trace_test.py @@ -359,6 +359,12 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: "herd_heads_2020", ], ).set_index(["admin_level", "admin_id", "commodity_name", "methodology"]) + + +# The 2020 assessment's snapshot weights +SNAPSHOT_YEAR_TO_WEIGHT = {2000: 0.03125, 2005: 0.125, 2010: 0.46875, 2020: 0.375} + + # 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( @@ -397,7 +403,9 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: def test_derive_statistical_production_kg_stays_indexed() -> None: result = derive_statistical_production_kg( - emissions=STATISTICAL_EMISSIONS, stocks=STOCKS + emissions=STATISTICAL_EMISSIONS, + snapshot_year_to_weight=SNAPSHOT_YEAR_TO_WEIGHT, + stocks=STOCKS, ) # regression guard: must return an indexed frame like its jurisdictional twin assert result.index.names == [ @@ -420,7 +428,9 @@ def test_derive_statistical_production_kg_prices_a_herd_at_each_snapshots_rate() None ): result = derive_statistical_production_kg( - emissions=STATISTICAL_EMISSIONS, stocks=STOCKS + emissions=STATISTICAL_EMISSIONS, + snapshot_year_to_weight=SNAPSHOT_YEAR_TO_WEIGHT, + stocks=STOCKS, ) # 100 heads at 10, 20, 40 and 80 kg, weighted 3.125%, 12.5%, 46.875% and 37.5% assert result.loc[ @@ -446,7 +456,9 @@ def test_methodologies_share_a_schema() -> None: ) statistical = attach_ratios( df=derive_statistical_production_kg( - emissions=STATISTICAL_EMISSIONS, stocks=STOCKS + emissions=STATISTICAL_EMISSIONS, + snapshot_year_to_weight=SNAPSHOT_YEAR_TO_WEIGHT, + stocks=STOCKS, ) ) assert set(jurisdictional.reset_index().columns) == set( diff --git a/jdluc/attribute.py b/jdluc/attribute.py index 9deffa1..aac4fd8 100644 --- a/jdluc/attribute.py +++ b/jdluc/attribute.py @@ -25,7 +25,6 @@ import collections.abc import concurrent.futures import enum -import functools import logging import pandas @@ -66,18 +65,11 @@ def workflow( iso_3166s: collections.abc.Iterable[str], methodology: Methodology, ) -> pandas.DataFrame: - workflow_for_tile: collections.abc.Callable[..., pandas.DataFrame] match methodology: case Methodology.JURISDICTIONAL_DIRECT: - workflow_for_tile = functools.partial( - jurisdictional_direct.workflow, - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, - ) + workflow_for_tile = jurisdictional_direct.workflow case Methodology.JURISDICTIONAL_DIRECT_PROXY: - workflow_for_tile = functools.partial( - jurisdictional_direct_proxy.workflow, - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, - ) + workflow_for_tile = jurisdictional_direct_proxy.workflow case Methodology.STATISTICAL: workflow_for_tile = statistical.workflow case _: @@ -97,6 +89,7 @@ def workflow( def workflow_for_tile_id(tile_id: str) -> list[pandas.DataFrame]: return [ workflow_for_tile( + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, commodity_names=commodity_names, iso_3166=iso_3166, tile_id=tile_id, diff --git a/jdluc/statistical.py b/jdluc/statistical.py index 057bab8..4ba43d9 100644 --- a/jdluc/statistical.py +++ b/jdluc/statistical.py @@ -115,8 +115,8 @@ class Livestock(enum.StrEnum): } -@storage.cache_to_zarr(version=2) -def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: +@storage.cache_to_zarr(version=0) +def get_downscaled_luc_emissions(assessment_year: int, tile_id: str) -> xarray.Dataset: logger.info("Computing emissions on the GLAD grid") luc_and_emissions = geo.exact_merge( harmonize.workflow( @@ -128,14 +128,14 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: ), # NB: this should call the harmonize workflow with identical args and hit # the cache from the preceding call - emit.workflow(assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, tile_id=tile_id), + emit.workflow(assessment_year=assessment_year, tile_id=tile_id), ) logger.info("Discounting emissions by component and splitting them by MapSPAM span") conversion_year = luc_and_emissions["conversion-year"] component_to_emissions = { component: emit.get_linear_discounted( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + assessment_year=assessment_year, conversion_year=conversion_year, darray=emissions, ) @@ -144,8 +144,18 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: ).items() } derived_variable_names: list[str] = [] + # NB: nothing departs in the lookback's opening year, so every other year can date a conversion + _, *conversion_years = emit.get_lookback_years_range( + assessment_year=assessment_year + ) for before, after in MAPSPAM_SPANS: - in_span = (conversion_year > before) & (conversion_year <= after) + in_span = conversion_year.isin( + [ + year + for year in conversion_years + if get_mapspam_span(year=year) == (before, after) + ] + ) for component in POOLED_COMPONENTS: name = f"{component!s}-emissions:tco2e-per-ha:{before:d}-{after:d}" luc_and_emissions[name] = component_to_emissions[component].where( @@ -169,7 +179,7 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: logger.info("Discounting the source carbon no destination claimed") name = "dropped-emissions:tco2e-per-ha" luc_and_emissions[name] = emit.get_linear_discounted( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + assessment_year=assessment_year, conversion_year=conversion_year, darray=luc_and_emissions[name], ) @@ -380,34 +390,54 @@ def get_commodity_to_peatland_conversion_share( } | get_livestock_to_grazing_share(dset=dset, year=after) -# A conversion takes the shares of the MapSPAM span holding its year, open below and closed above. -# MapSPAM has no 2015 snapshot, so every conversion from 2011 to 2020 takes the one expansion -# between 2010 and 2020. MAPSPAM_SPANS = ((2000, 2005), (2005, 2010), (2010, 2020)) -assert all( - sum(before < year <= after for before, after in MAPSPAM_SPANS) == 1 - for year in emit.get_year_to_linear_discount_weight( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR - ) -) MAPSPAM_SNAPSHOT_YEARS = tuple( sorted({year for span in MAPSPAM_SPANS for year in span}) ) assert set(MAPSPAM_SNAPSHOT_YEARS).issubset(gpw_livestock.YEARS) -# Each snapshot's weight in `get_discounted_snapshot_mean`: every conversion year's discount weight, -# split evenly between the two snapshots of the MapSPAM span holding that year -SNAPSHOT_YEAR_TO_WEIGHT = { - snapshot_year: sum( - weight / 2 - for year, weight in emit.get_year_to_linear_discount_weight( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR - ).items() + + +def get_mapspam_span(year: int) -> tuple[int, int]: + """The MapSPAM span whose shares a conversion in `year` takes: the one holding it, open below + and closed above. + + MapSPAM has no 2015 snapshot, so every conversion from 2011 to 2020 takes the one expansion + between 2010 and 2020, and it has none after 2020, so every later conversion takes it too. + """ + *_, last_snapshot_year = MAPSPAM_SNAPSHOT_YEARS + (span,) = ( + (before, after) for before, after in MAPSPAM_SPANS - if before < year <= after and snapshot_year in (before, after) + if before < min(year, last_snapshot_year) <= after ) - for snapshot_year in MAPSPAM_SNAPSHOT_YEARS -} -assert math.isclose(sum(SNAPSHOT_YEAR_TO_WEIGHT.values()), 1) + return span + + +def get_snapshot_year_to_weight(assessment_year: int) -> dict[int, float]: + """Each snapshot's weight in `get_discounted_snapshot_mean`: every conversion year's discount + weight, split evenly between the two snapshots of its MapSPAM span. A year after the last + snapshot gives its whole weight to that snapshot instead. + """ + *_, last_snapshot_year = MAPSPAM_SNAPSHOT_YEARS + year_to_weight = emit.get_year_to_linear_discount_weight( + assessment_year=assessment_year + ) + year_to_snapshot_years = { + year: (last_snapshot_year,) + if year > last_snapshot_year + else get_mapspam_span(year=year) + for year in year_to_weight + } + snapshot_year_to_weight = { + snapshot_year: sum( + weight / len(year_to_snapshot_years[year]) + for year, weight in year_to_weight.items() + if snapshot_year in year_to_snapshot_years[year] + ) + for snapshot_year in MAPSPAM_SNAPSHOT_YEARS + } + assert math.isclose(sum(snapshot_year_to_weight.values()), 1) + return snapshot_year_to_weight def get_herd_heads_column(year: int) -> str: @@ -420,12 +450,13 @@ def get_herd_heads_column(year: int) -> str: def get_discounted_snapshot_mean( + snapshot_year_to_weight: dict[int, float], year_to_value: dict[int, xarray.DataArray], ) -> xarray.DataArray: """A level -- area, production -- as a mean over the windows emissions are summed on.""" return sum( ( - SNAPSHOT_YEAR_TO_WEIGHT[year] * year_to_value[year] + snapshot_year_to_weight[year] * year_to_value[year] for year in MAPSPAM_SNAPSHOT_YEARS ), start=xarray.DataArray(numpy.float32(0)), @@ -436,6 +467,7 @@ def get_commodity_name_to_totals( commodities: tuple[Commodity, ...], dset: xarray.Dataset, occupation_shares: dict[Commodity, xarray.DataArray], + snapshot_year_to_weight: dict[int, float], span_to_commodity_to_peatland_conversion_share: dict[ tuple[int, int], dict[Commodity, xarray.DataArray] ], @@ -527,6 +559,7 @@ def get_commodity_name_to_totals( else: occupation = cropland_occupation * occupation_shares[commodity] totals["production_mt"] = get_discounted_snapshot_mean( + snapshot_year_to_weight=snapshot_year_to_weight, year_to_value={ year: ifpri_mapspam.get_canonical_quantity( canonical_crop_name=commodity.value, @@ -539,6 +572,7 @@ def get_commodity_name_to_totals( ) totals["commodity_hectares"] = get_discounted_snapshot_mean( + snapshot_year_to_weight=snapshot_year_to_weight, year_to_value={ year: get_commodity_hectares(commodity=commodity, dset=dset, year=year) for year in MAPSPAM_SNAPSHOT_YEARS @@ -662,8 +696,9 @@ def get_crop_to_area_share( } | dict.fromkeys(HERD_HEADS_COLUMNS, float) -@storage.cache_to_parquet(version=2) +@storage.cache_to_parquet(version=0) def workflow( + assessment_year: int, commodity_names: tuple[str, ...], iso_3166: str, tile_id: str, @@ -675,10 +710,15 @@ def workflow( for name in commodity_names ) crops = tuple(commodity for commodity in commodities if isinstance(commodity, Crop)) - logger.info(f"Computing emissions for {commodities=:} and {tile_id=:s}") + logger.info( + f"Computing emissions for {commodities=:}, {assessment_year=:d} and {tile_id=:s}" + ) + snapshot_year_to_weight = get_snapshot_year_to_weight( + assessment_year=assessment_year + ) merged = geo.exact_merge( # NB: this is deferred because it is expensive and would like to cache it - get_downscaled_luc_emissions(tile_id=tile_id), + get_downscaled_luc_emissions(assessment_year=assessment_year, tile_id=tile_id), # MAPSPAM harmonize.workflow( dataset_names=DATASET_NAMES, @@ -762,6 +802,7 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]: | get_livestock_to_grazing_share( dset=clipped, year=max(ifpri_mapspam.YEARS) ), + snapshot_year_to_weight=snapshot_year_to_weight, span_to_commodity_to_peatland_conversion_share=span_to_commodity_to_peatland_conversion_share, span_to_commodity_to_share=span_to_commodity_to_share, ) diff --git a/jdluc/trace.py b/jdluc/trace.py index fad0cb1..3debaac 100644 --- a/jdluc/trace.py +++ b/jdluc/trace.py @@ -117,7 +117,9 @@ def get_herd_production_kg( def derive_statistical_production_kg( - emissions: pandas.DataFrame, stocks: pandas.DataFrame + emissions: pandas.DataFrame, + snapshot_year_to_weight: dict[int, float], + 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 @@ -135,7 +137,7 @@ def derive_statistical_production_kg( stocks=stocks, years=(year,), ) - for year, weight in statistical.SNAPSHOT_YEAR_TO_WEIGHT.items() + for year, weight in snapshot_year_to_weight.items() ) return ret.drop(columns=["production_mt", *statistical.HERD_HEADS_COLUMNS]) @@ -340,6 +342,9 @@ def workflow( case attribute.Methodology.STATISTICAL: production = derive_statistical_production_kg( emissions=emissions, + snapshot_year_to_weight=statistical.get_snapshot_year_to_weight( + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR + ), stocks=faostat.load(dataset=faostat.LIVESTOCK_DATASET), ) case _: diff --git a/tools/measure-drift.py b/tools/measure-drift.py index 76ecf11..9800243 100644 --- a/tools/measure-drift.py +++ b/tools/measure-drift.py @@ -71,8 +71,8 @@ # a tropical forest frontier, and tropical peat. Every added country multiplies capture cost, so # widen this deliberately -- and use --iso to work on a subset meanwhile. HND and NIC are the # cheap exception: both fall entirely inside the one ten-degree tile 20N_090W, and the expensive -# layer -- statistical.get_downscaled_luc_emissions -- is keyed on tile_id alone, so whichever -# of the two runs second pays for little beyond its province clips. +# layer -- statistical.get_downscaled_luc_emissions -- is keyed on the tile and the assessment year +# alone, so whichever of the two runs second pays for little beyond its province clips. # # Central America earns five slots because it carries every 2000 crop group ifpri_mapspam still # decomposes, and the weight sits in different groups either side of a border: 44% of Honduran diff --git a/validation/capture.py b/validation/capture.py index 1365a18..a27b267 100644 --- a/validation/capture.py +++ b/validation/capture.py @@ -12,8 +12,8 @@ **Every commodity, not just the targets.** `attribute.get_commodity_names` fixes the commodity set per leg, so a capture covers every commodity the leg models rather than only the pairs `targets.json` names. It is cheaper than it looks: -`statistical.get_downscaled_luc_emissions` is keyed on `tile_id` alone, so the expensive -per-tile layer is shared and a crop adds only its share arithmetic. +`statistical.get_downscaled_luc_emissions` is keyed on the tile and the assessment year alone, so +the expensive per-tile layer is shared and a crop adds only its share arithmetic. **Merging is per (ISO, methodology), and rows keep their own `code_version`.** A capture scoped with `--isos` leaves other countries in place, so the artifact can span versions and a reader has to be From 7d31cadd231bcea036e862696299299fbc8bdbd1 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 15:16:43 -0800 Subject: [PATCH 07/16] Add --assessment-year to jdluc.attribute --- jdluc/attribute.py | 16 ++++++++++++---- jdluc/trace.py | 3 ++- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/jdluc/attribute.py b/jdluc/attribute.py index aac4fd8..e157c87 100644 --- a/jdluc/attribute.py +++ b/jdluc/attribute.py @@ -18,7 +18,7 @@ Example invocations: uv run python jdluc/attribute.py --methodology-name STATISTICAL USA - uv run python jdluc/attribute.py --methodology-name STATISTICAL --backfill + uv run python jdluc/attribute.py --assessment-year 2024 --methodology-name STATISTICAL --backfill """ import argparse @@ -60,8 +60,9 @@ def merge_dfs(*dfs: pandas.DataFrame) -> pandas.DataFrame: def workflow( - concurrency: int, + assessment_year: int, commodity_names: tuple[str, ...], + concurrency: int, iso_3166s: collections.abc.Iterable[str], methodology: Methodology, ) -> pandas.DataFrame: @@ -89,7 +90,7 @@ def workflow( def workflow_for_tile_id(tile_id: str) -> list[pandas.DataFrame]: return [ workflow_for_tile( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + assessment_year=assessment_year, commodity_names=commodity_names, iso_3166=iso_3166, tile_id=tile_id, @@ -167,6 +168,12 @@ def main() -> int: nargs=argparse.ZERO_OR_MORE, type=worldbank_jurisdictions.iso_3166_str, ) + parser.add_argument( + "--assessment-year", + choices=emit.ASSESSMENT_YEARS, + default=emit.DEFAULT_ASSESSMENT_YEAR, + type=int, + ) parser.add_argument( "--backfill", action="store_true", @@ -186,8 +193,9 @@ def main() -> int: methodology = Methodology[str(args.methodology_name)] df = workflow( - concurrency=int(args.concurrency), + assessment_year=int(args.assessment_year), commodity_names=get_commodity_names(methodology=methodology), + concurrency=int(args.concurrency), iso_3166s=( worldbank_jurisdictions.get_all_iso_3166s() if args.backfill diff --git a/jdluc/trace.py b/jdluc/trace.py index 3debaac..34c46a8 100644 --- a/jdluc/trace.py +++ b/jdluc/trace.py @@ -317,8 +317,9 @@ def workflow( methodology: attribute.Methodology, ) -> pandas.DataFrame: emissions = attribute.workflow( - concurrency=concurrency, + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, commodity_names=commodity_names, + concurrency=concurrency, iso_3166s=iso_3166s, methodology=methodology, ) From b0d63c590ae8a354e412a790991761bde4c51377 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 15:23:25 -0800 Subject: [PATCH 08/16] Add --assessment-year to jdluc.trace --- jdluc/__tests__/trace_test.py | 26 ++++++++++++--- jdluc/trace.py | 62 ++++++++++++++++++++++++----------- tools/export-graphs.py | 3 +- tools/measure-drift.py | 3 +- validation/capture.py | 5 +-- 5 files changed, 71 insertions(+), 28 deletions(-) diff --git a/jdluc/__tests__/trace_test.py b/jdluc/__tests__/trace_test.py index a7ec0e2..8042b0a 100644 --- a/jdluc/__tests__/trace_test.py +++ b/jdluc/__tests__/trace_test.py @@ -111,10 +111,22 @@ ).set_index(["admin_level", "admin_id", "jurisdiction_name", "commodity_name", "year"]) +def test_derive_jurisdictional_production_kg_needs_every_year_of_its_window() -> None: + # The yields stop at 2020, so a window reaching 2021 has a year to average that NASS never gave + with pytest.raises(AssertionError, match="NASS yields are missing a year"): + derive_jurisdictional_production_kg( + emissions=JURISDICTIONAL_EMISSIONS, + raw_yields=RAW_YIELDS, + years=(2018, 2019, 2020, 2021), + ) + + def test_derive_jurisdictional_production_kg_then_attach_ratios() -> None: result = attach_ratios( df=derive_jurisdictional_production_kg( - emissions=JURISDICTIONAL_EMISSIONS, raw_yields=RAW_YIELDS + emissions=JURISDICTIONAL_EMISSIONS, + raw_yields=RAW_YIELDS, + years=(2017, 2018, 2019, 2020), ) ) assert isinstance(result, pandas.DataFrame) @@ -451,7 +463,9 @@ def test_derive_statistical_production_kg_prices_a_herd_at_each_snapshots_rate() def test_methodologies_share_a_schema() -> None: jurisdictional = attach_ratios( df=derive_jurisdictional_production_kg( - emissions=JURISDICTIONAL_EMISSIONS, raw_yields=RAW_YIELDS + emissions=JURISDICTIONAL_EMISSIONS, + raw_yields=RAW_YIELDS, + years=(2017, 2018, 2019, 2020), ) ) statistical = attach_ratios( @@ -517,7 +531,9 @@ def test_with_commodity_ids_adds_together_the_names_sharing_an_item() -> None: 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 + emissions=JURISDICTIONAL_EMISSIONS, + raw_yields=RAW_YIELDS, + years=(2017, 2018, 2019, 2020), ), methodology=Methodology.JURISDICTIONAL_DIRECT, ) @@ -609,7 +625,9 @@ def test_derive_proxy_production_kg() -> None: emissions = emissions.set_index( ["admin_level", "admin_id", "commodity_name", "methodology"] ) - result = derive_proxy_production_kg(crops=crops, emissions=emissions, stocks=stocks) + result = derive_proxy_production_kg( + crops=crops, emissions=emissions, stocks=stocks, years=(2019, 2020, 2021) + ) production_kg = result["production_kg"].droplevel(["admin_level", "methodology"]) assert production_kg[("BRA001", "MAIZE")] == pytest.approx(120.0) assert production_kg[("BRA002", "MAIZE")] == pytest.approx(40.0) diff --git a/jdluc/trace.py b/jdluc/trace.py index 34c46a8..f2c087d 100644 --- a/jdluc/trace.py +++ b/jdluc/trace.py @@ -2,12 +2,14 @@ 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). + - JURISDICTIONAL_DIRECT: production = crop area x NASS QuickStats yield, averaged over the 4 + years ending at the assessment year. - 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. - - JURISDICTIONAL_DIRECT_PROXY: production = FAOSTAT national production split among provinces - by CROPGRIDS crop area for crops, and the attribute's 2020 GPW heads x FAOSTAT carcass weight - per head for livestock. + - JURISDICTIONAL_DIRECT_PROXY: production = FAOSTAT national production over the 3 years + centered on the assessment year, split among provinces by CROPGRIDS crop area for crops, and + the attribute's 2020 GPW heads x FAOSTAT carcass weight per head over the same years for + livestock. It then derives, identically for all three, each commodity's yield (kg per hectare) and emissions factor (kgCO2e per kg), and rolls the provincials up to national totals. Returns a cached @@ -20,7 +22,7 @@ Example invocations: uv run python jdluc/trace.py --methodology-name STATISTICAL USA - uv run python jdluc/trace.py --methodology-name STATISTICAL --backfill + uv run python jdluc/trace.py --assessment-year 2024 --methodology-name STATISTICAL --backfill """ import argparse @@ -50,9 +52,6 @@ logger = logging.getLogger(__name__) -NASS_YIELD_YEARS = (2017, 2018, 2019, 2020) -# The proxy leg's FAOSTAT years, centered on CROPGRIDS' -FAOSTAT_YEARS = tuple(range(tang_cropgrids.YEAR - 1, tang_cropgrids.YEAR + 2)) KG_PER_TONNE = 1000 CANONICAL_KEY = ("admin_level", "admin_id", "commodity_id", "methodology") ADDITIVE_COLUMNS = ( @@ -66,10 +65,13 @@ def derive_jurisdictional_production_kg( - emissions: pandas.DataFrame, raw_yields: pandas.DataFrame + emissions: pandas.DataFrame, raw_yields: pandas.DataFrame, years: tuple[int, ...] ) -> pandas.DataFrame: + """A crop's production is its area times NASS's state yield, averaged over `years`.""" + yield_years = raw_yields.index.get_level_values("year") + assert set(years) <= set(yield_years), f"NASS yields are missing a year of {years}" reduced_yields = ( - raw_yields[raw_yields.index.get_level_values("year").isin(NASS_YIELD_YEARS)] + raw_yields[yield_years.isin(years)] .groupby(level=["admin_id", "commodity_name"])["yield_kg_per_ha"] .mean() .reset_index() @@ -143,9 +145,12 @@ def derive_statistical_production_kg( def derive_proxy_production_kg( - crops: pandas.DataFrame, emissions: pandas.DataFrame, stocks: pandas.DataFrame + crops: pandas.DataFrame, + emissions: pandas.DataFrame, + stocks: pandas.DataFrame, + years: tuple[int, ...], ) -> pandas.DataFrame: - """A crop's production is FAOSTAT's national production of its item over `FAOSTAT_YEARS`, split + """A crop's production is FAOSTAT's national production of its item over `years`, split among provinces by CROPGRIDS crop area. A herd's is its heads times FAOSTAT's carcass weight per head. Rows with no FAOSTAT item have none. """ @@ -167,8 +172,14 @@ def derive_proxy_production_kg( | set(taxonomy.CO_PRODUCT_ITEM_CODE_TO_ITEM_CODE) | taxonomy.UNMAPPED_CROP_ITEM_CODES ) + assert set(years) <= set(table["year"]), ( + f"FAOSTAT production is missing a year of {years}" + ) + assert set(years) <= set(stocks.index.get_level_values("year")), ( + f"FAOSTAT stocks are missing a year of {years}" + ) national_production = ( - table[table["year"].isin(FAOSTAT_YEARS)] + table[table["year"].isin(years)] .groupby(["admin_id", "commodity_id"])["production_kg"] .mean() ) @@ -222,7 +233,7 @@ def derive_proxy_production_kg( heads=ret.loc[is_herd, jurisdictional_direct_proxy.HERD_HEADS_COLUMN], species=species, stocks=stocks, - years=FAOSTAT_YEARS, + years=years, ) return ret.drop(columns=jurisdictional_direct_proxy.HERD_HEADS_COLUMN) @@ -309,15 +320,16 @@ def safe_div(numer: float, denom: float) -> float: } # type: ignore -@storage.cache_to_parquet(ignored_args=["concurrency"], version=2) +@storage.cache_to_parquet(ignored_args=["concurrency"], version=0) def workflow( - concurrency: int, + assessment_year: int, commodity_names: tuple[str, ...], + concurrency: int, iso_3166s: tuple[str, ...], methodology: attribute.Methodology, ) -> pandas.DataFrame: emissions = attribute.workflow( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + assessment_year=assessment_year, commodity_names=commodity_names, concurrency=concurrency, iso_3166s=iso_3166s, @@ -332,19 +344,22 @@ def workflow( match methodology: case attribute.Methodology.JURISDICTIONAL_DIRECT: production = derive_jurisdictional_production_kg( - emissions=emissions, raw_yields=usda_nass_quickstats.load() + emissions=emissions, + raw_yields=usda_nass_quickstats.load(), + years=tuple(range(assessment_year - 3, assessment_year + 1)), ) case attribute.Methodology.JURISDICTIONAL_DIRECT_PROXY: production = derive_proxy_production_kg( crops=faostat.load(dataset=faostat.CROP_DATASET), emissions=emissions, stocks=faostat.load(dataset=faostat.LIVESTOCK_DATASET), + years=tuple(range(assessment_year - 1, assessment_year + 2)), ) case attribute.Methodology.STATISTICAL: production = derive_statistical_production_kg( emissions=emissions, snapshot_year_to_weight=statistical.get_snapshot_year_to_weight( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR + assessment_year=assessment_year ), stocks=faostat.load(dataset=faostat.LIVESTOCK_DATASET), ) @@ -407,6 +422,12 @@ def main() -> int: nargs=argparse.ZERO_OR_MORE, type=worldbank_jurisdictions.iso_3166_str, ) + parser.add_argument( + "--assessment-year", + choices=emit.ASSESSMENT_YEARS, + default=emit.DEFAULT_ASSESSMENT_YEAR, + type=int, + ) parser.add_argument( "--backfill", action="store_true", @@ -426,8 +447,9 @@ def main() -> int: methodology = attribute.Methodology[str(args.methodology_name)] df = workflow( - concurrency=int(args.concurrency), + assessment_year=int(args.assessment_year), commodity_names=attribute.get_commodity_names(methodology=methodology), + concurrency=int(args.concurrency), iso_3166s=tuple( sorted( worldbank_jurisdictions.get_all_iso_3166s() diff --git a/tools/export-graphs.py b/tools/export-graphs.py index 614d6e9..fd55f0d 100644 --- a/tools/export-graphs.py +++ b/tools/export-graphs.py @@ -11,7 +11,7 @@ import numpy import pandas -from jdluc import attribute, trace +from jdluc import attribute, emit, trace from jdluc.datasets import worldbank_jurisdictions logger = logging.getLogger(__name__) @@ -149,6 +149,7 @@ def iter_dfs( yield subset_to_commodity_ids( commodity_ids=commodity_ids, df=trace.workflow( + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, commodity_names=workflow_commodity_names, iso_3166s=(iso_3166,), methodology=methodology, diff --git a/tools/measure-drift.py b/tools/measure-drift.py index 9800243..c90cd08 100644 --- a/tools/measure-drift.py +++ b/tools/measure-drift.py @@ -61,7 +61,7 @@ import pandas -from jdluc import attribute, config, storage, taxonomy, trace +from jdluc import attribute, config, emit, storage, taxonomy, trace logger = logging.getLogger(__name__) @@ -928,6 +928,7 @@ def capture(directory: pathlib.Path, label: Side, targets: tuple[Target, ...]) - retired = {"skip_glad_crop_filter": False} parameters = inspect.signature(trace.workflow).parameters df = trace.workflow( + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, commodity_names=commodity_names, concurrency=attribute.DEFAULT_CONCURRENCY, iso_3166s=(target.iso_3166,), diff --git a/validation/capture.py b/validation/capture.py index a27b267..5b5c08f 100644 --- a/validation/capture.py +++ b/validation/capture.py @@ -39,7 +39,7 @@ import pandas -from jdluc import attribute, trace +from jdluc import attribute, emit, trace from jdluc.datasets import worldbank_jurisdictions from validation import pull, schema, targets @@ -166,8 +166,9 @@ def workflow(iso_3166s: tuple[str, ...], repo_root: pathlib.Path) -> pandas.Data ) frames.append( trace.workflow( - concurrency=attribute.DEFAULT_CONCURRENCY, + assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, commodity_names=commodity_names, + concurrency=attribute.DEFAULT_CONCURRENCY, iso_3166s=wanted, methodology=methodology, ) From 94a7c5dca1400ef460d11f18f3248d1c191f169e Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 15:25:00 -0800 Subject: [PATCH 09/16] change JD production window to 3-year window centered on assessment_year --- jdluc/trace.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/jdluc/trace.py b/jdluc/trace.py index f2c087d..e3ea295 100644 --- a/jdluc/trace.py +++ b/jdluc/trace.py @@ -2,8 +2,8 @@ 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, averaged over the 4 - years ending at the assessment year. + - JURISDICTIONAL_DIRECT: production = crop area x NASS QuickStats yield, averaged over the 3 + years centered on the assessment year. - 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. - JURISDICTIONAL_DIRECT_PROXY: production = FAOSTAT national production over the 3 years @@ -346,7 +346,7 @@ def workflow( production = derive_jurisdictional_production_kg( emissions=emissions, raw_yields=usda_nass_quickstats.load(), - years=tuple(range(assessment_year - 3, assessment_year + 1)), + years=tuple(range(assessment_year - 1, assessment_year + 2)), ) case attribute.Methodology.JURISDICTIONAL_DIRECT_PROXY: production = derive_proxy_production_kg( From 0d345c0ed0b11ce4dcd293ce8fef8330299eecf0 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 15:36:17 -0800 Subject: [PATCH 10/16] Add assessment_year to tabular outputs as a new index column --- jdluc/__tests__/trace_test.py | 187 ++++++++++++++++------------ jdluc/attribute.py | 21 +++- jdluc/trace.py | 24 ++-- tools/measure-drift.py | 18 ++- validation/__tests__/schema_test.py | 7 +- validation/capture.py | 33 +++-- validation/prepare.py | 9 +- validation/schema.py | 13 +- 8 files changed, 191 insertions(+), 121 deletions(-) diff --git a/jdluc/__tests__/trace_test.py b/jdluc/__tests__/trace_test.py index 8042b0a..208ef0a 100644 --- a/jdluc/__tests__/trace_test.py +++ b/jdluc/__tests__/trace_test.py @@ -25,61 +25,67 @@ [ # normal row ( + "USA008", "PROVINCIAL", + 2020, + 100.0, "MAIZE", + 200.0, + 120.0, "Delaware", "JURISDICTIONAL_DIRECT", - "USA008", - 100.0, 10.0, - 50.0, - 200.0, - 120.0, 30.0, + 50.0, ), # zero production (commodity_hectares == 0) ( + "USA008", "PROVINCIAL", + 2020, + 0.0, "SOYBEAN", + 80.0, + 50.0, "Delaware", "JURISDICTIONAL_DIRECT", - "USA008", - 0.0, 0.0, - 0.0, - 80.0, - 50.0, 30.0, + 0.0, ), # unmatched yield + zero emissions ( - "PROVINCIAL", - "WHEAT", - "Iowa", - "JURISDICTIONAL_DIRECT", "USA016", + "PROVINCIAL", + 2020, 50.0, + "WHEAT", 0.0, 0.0, + "Iowa", + "JURISDICTIONAL_DIRECT", 0.0, 0.0, 0.0, ), ], columns=[ + "admin_id", "admin_level", + "assessment_year", + "commodity_hectares", "commodity_name", + "emissions_mt", + "forest_emissions_mt", "jurisdiction_name", "methodology", - "admin_id", - "commodity_hectares", "peatland_commodity_hectares", - "peatland_occupation_emissions_mt", - "emissions_mt", - "forest_emissions_mt", "peatland_conversion_emissions_mt", + "peatland_occupation_emissions_mt", ], -).set_index(["admin_level", "admin_id", "commodity_name", "methodology"]) +).set_index( + ["admin_id", "admin_level", "assessment_year", "commodity_name", "methodology"] +) RAW_YIELDS = pandas.DataFrame.from_records( [ # MAIZE: 2016 must be excluded; mean(8, 10, 12, 10) over 2017-2020 == 10 @@ -133,20 +139,20 @@ def test_derive_jurisdictional_production_kg_then_attach_ratios() -> None: iter_result = result.iterrows() key, corn = next(iter_result) - assert key == ("PROVINCIAL", "USA008", "MAIZE", "JURISDICTIONAL_DIRECT") + assert key == ("USA008", "PROVINCIAL", 2020, "MAIZE", "JURISDICTIONAL_DIRECT") assert corn["yield_kg_per_ha"] == 10.0 # 4-year mean, 2016 excluded assert corn["production_kg"] == 1000.0 # 100 ha x 10 assert corn["emissions_factor_kgco2e_per_kg"] == 200.0 # 200 t x 1000 / 1000 kg # zero production -> EF guarded to NaN rather than inf key, soy = next(iter_result) - assert key == ("PROVINCIAL", "USA008", "SOYBEAN", "JURISDICTIONAL_DIRECT") + assert key == ("USA008", "PROVINCIAL", 2020, "SOYBEAN", "JURISDICTIONAL_DIRECT") assert soy["production_kg"] == 0.0 assert numpy.isnan(soy["emissions_factor_kgco2e_per_kg"]) # unmatched yield -> NaN yield, production and EF key, wheat = next(iter_result) - assert key == ("PROVINCIAL", "USA016", "WHEAT", "JURISDICTIONAL_DIRECT") + assert key == ("USA016", "PROVINCIAL", 2020, "WHEAT", "JURISDICTIONAL_DIRECT") assert numpy.isnan(wheat["yield_kg_per_ha"]) assert numpy.isnan(wheat["production_kg"]) assert numpy.isnan(wheat["emissions_factor_kgco2e_per_kg"]) @@ -169,6 +175,7 @@ def test_iter_national_from_provincials_single_province() -> None: ] ) (result,) = iter_national_from_provincials( + assessment_year=2020, iso_3166="ISO_A3", methodology=Methodology.STATISTICAL, national_name="NAME", @@ -177,6 +184,7 @@ def test_iter_national_from_provincials_single_province() -> None: assert result == { "admin_id": "ISO_A3", "admin_level": AdminLevel.NATIONAL.name, + "assessment_year": 2020, "commodity_hectares": 1, "commodity_id": 236, "emissions_factor_kgco2e_per_kg": 1000, @@ -235,6 +243,7 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: ] ) (maize, soya_beans) = iter_national_from_provincials( + assessment_year=2020, iso_3166="ISO_A3", methodology=Methodology.STATISTICAL, national_name="NAME", @@ -243,6 +252,7 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: assert maize == { "admin_id": "ISO_A3", "admin_level": AdminLevel.NATIONAL.name, + "assessment_year": 2020, "commodity_hectares": 2, "commodity_id": 56, "emissions_factor_kgco2e_per_kg": 1000, @@ -260,6 +270,7 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: assert soya_beans == { "admin_id": "ISO_A3", "admin_level": AdminLevel.NATIONAL.name, + "assessment_year": 2020, "commodity_hectares": 1, "commodity_id": 236, "emissions_factor_kgco2e_per_kg": 1000, @@ -279,98 +290,105 @@ def test_iter_national_from_provincials_multi_province_crop() -> None: STATISTICAL_EMISSIONS = pandas.DataFrame.from_records( [ ( - "PROVINCIAL", - "CORN", - "Delaware", - "STATISTICAL", "USA008", + "PROVINCIAL", + 2020, 100.0, - 10.0, - 50.0, + "CORN", 200.0, - 5.0, 120.0, - 30.0, 0.0, 0.0, 0.0, 0.0, + "Delaware", + "STATISTICAL", + 10.0, + 30.0, + 50.0, + 5.0, ), ( - "PROVINCIAL", - "BEEF_CATTLE", - "Acre", - "STATISTICAL", "BRA001", + "PROVINCIAL", + 2020, 100.0, - 0.0, - 0.0, + "BEEF_CATTLE", 200.0, - 0.0, 200.0, - 0.0, 100.0, 100.0, 100.0, 100.0, + "Acre", + "STATISTICAL", + 0.0, + 0.0, + 0.0, + 0.0, ), # A country FAOSTAT reports no herd for ( - "PROVINCIAL", - "BEEF_CATTLE", - "Niue", - "STATISTICAL", "NIU001", + "PROVINCIAL", + 2020, 100.0, - 0.0, - 0.0, + "BEEF_CATTLE", 200.0, - 0.0, 200.0, - 0.0, 100.0, 100.0, 100.0, 100.0, + "Niue", + "STATISTICAL", + 0.0, + 0.0, + 0.0, + 0.0, ), ( - "PROVINCIAL", - "BEEF_CATTLE", - "Concepción", - "STATISTICAL", "PRY001", + "PROVINCIAL", + 2020, 100.0, - 0.0, - 0.0, + "BEEF_CATTLE", 200.0, - 0.0, 200.0, - 0.0, 100.0, 100.0, 100.0, 100.0, + "Concepción", + "STATISTICAL", + 0.0, + 0.0, + 0.0, + 0.0, ), ], columns=[ - "admin_level", - "commodity_name", - "jurisdiction_name", - "methodology", "admin_id", + "admin_level", + "assessment_year", "commodity_hectares", - "peatland_commodity_hectares", - "peatland_occupation_emissions_mt", + "commodity_name", "emissions_mt", - "production_mt", "forest_emissions_mt", - "peatland_conversion_emissions_mt", "herd_heads_2000", "herd_heads_2005", "herd_heads_2010", "herd_heads_2020", + "jurisdiction_name", + "methodology", + "peatland_commodity_hectares", + "peatland_conversion_emissions_mt", + "peatland_occupation_emissions_mt", + "production_mt", ], -).set_index(["admin_level", "admin_id", "commodity_name", "methodology"]) +).set_index( + ["admin_id", "admin_level", "assessment_year", "commodity_name", "methodology"] +) # The 2020 assessment's snapshot weights @@ -421,15 +439,18 @@ def test_derive_statistical_production_kg_stays_indexed() -> None: ) # regression guard: must return an indexed frame like its jurisdictional twin assert result.index.names == [ - "admin_level", "admin_id", + "admin_level", + "assessment_year", "commodity_name", "methodology", ] assert "production_mt" not in result.columns assert not set(HERD_HEADS_COLUMNS) & set(result.columns) assert ( - result.loc[("PROVINCIAL", "USA008", "CORN", "STATISTICAL"), "production_kg"] + result.loc[ + ("USA008", "PROVINCIAL", 2020, "CORN", "STATISTICAL"), "production_kg" + ] == 5000.0 # 5 mt x 1000 ) # .copy() guard: attach_ratios mutates in place, so the source frame must be untouched @@ -446,17 +467,18 @@ def test_derive_statistical_production_kg_prices_a_herd_at_each_snapshots_rate() ) # 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" + ("BRA001", "PROVINCIAL", 2020, "BEEF_CATTLE", "STATISTICAL"), "production_kg" ] == pytest.approx(5156.25) assert ( result.loc[ - ("PROVINCIAL", "NIU001", "BEEF_CATTLE", "STATISTICAL"), "production_kg" + ("NIU001", "PROVINCIAL", 2020, "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" + ("PRY001", "PROVINCIAL", 2020, "BEEF_CATTLE", "STATISTICAL"), "production_kg" ] == pytest.approx(3000.0) @@ -488,6 +510,7 @@ def test_with_commodity_ids_adds_together_the_names_sharing_an_item() -> None: "BRA001", "ARABICA_COFFEE", "STATISTICAL", + 2020, "Acre", 1.0, 10.0, @@ -497,30 +520,34 @@ def test_with_commodity_ids_adds_together_the_names_sharing_an_item() -> None: "BRA001", "ROBUSTA_COFFEE", "STATISTICAL", + 2020, "Acre", 2.0, 20.0, ), - ("PROVINCIAL", "BRA001", "MAIZE", "STATISTICAL", "Acre", 4.0, 40.0), + ("PROVINCIAL", "BRA001", "MAIZE", "STATISTICAL", 2020, "Acre", 4.0, 40.0), ], columns=[ "admin_level", "admin_id", "commodity_name", "methodology", + "assessment_year", "jurisdiction_name", "commodity_hectares", "production_kg", ], - ).set_index(["admin_level", "admin_id", "commodity_name", "methodology"]) + ).set_index( + ["admin_id", "admin_level", "assessment_year", "commodity_name", "methodology"] + ) result = with_commodity_ids(df=emissions, methodology=Methodology.STATISTICAL) assert result.to_dict(orient="index") == { - ("PROVINCIAL", "BRA001", 56, "STATISTICAL"): { + ("BRA001", "PROVINCIAL", 2020, 56, "STATISTICAL"): { "commodity_hectares": 4.0, "jurisdiction_name": "Acre", "production_kg": 40.0, }, - ("PROVINCIAL", "BRA001", 656, "STATISTICAL"): { + ("BRA001", "PROVINCIAL", 2020, 656, "STATISTICAL"): { "commodity_hectares": 3.0, "jurisdiction_name": "Acre", "production_kg": 30.0, @@ -539,7 +566,7 @@ def test_with_commodity_ids_keeps_a_missing_yield_missing() -> None: ) assert pandas.isna( result.loc[ - ("PROVINCIAL", "USA016", 15, "JURISDICTIONAL_DIRECT"), + ("USA016", "PROVINCIAL", 2020, 15, "JURISDICTIONAL_DIRECT"), "production_kg", ] ) @@ -621,14 +648,20 @@ def test_derive_proxy_production_kg() -> None: ("BRA001", "DROPPED", 0.0, numpy.nan), ], columns=["admin_id", "commodity_name", "commodity_hectares", HERD_HEADS_COLUMN], - ).assign(admin_level="PROVINCIAL", methodology="JURISDICTIONAL_DIRECT_PROXY") + ).assign( + admin_level="PROVINCIAL", + assessment_year=2020, + methodology="JURISDICTIONAL_DIRECT_PROXY", + ) emissions = emissions.set_index( - ["admin_level", "admin_id", "commodity_name", "methodology"] + ["admin_id", "admin_level", "assessment_year", "commodity_name", "methodology"] ) result = derive_proxy_production_kg( crops=crops, emissions=emissions, stocks=stocks, years=(2019, 2020, 2021) ) - production_kg = result["production_kg"].droplevel(["admin_level", "methodology"]) + production_kg = result["production_kg"].droplevel( + ["admin_level", "methodology", "assessment_year"] + ) assert production_kg[("BRA001", "MAIZE")] == pytest.approx(120.0) assert production_kg[("BRA002", "MAIZE")] == pytest.approx(40.0) assert production_kg[("BRA001", "POPCORN")] == pytest.approx(40.0) diff --git a/jdluc/attribute.py b/jdluc/attribute.py index e157c87..7502176 100644 --- a/jdluc/attribute.py +++ b/jdluc/attribute.py @@ -9,9 +9,10 @@ by whether each pixel became cropland or pasture, and divides cropland's part among CROPGRIDS' crops by their 2020 crop-area shares and pasture's part among grazers. All three clip to provincial (World Bank admin-1) polygons and sum commodity area, peatland -commodity area, peatland-occupation emissions, and total emissions. Returns a pandas.DataFrame indexed by -(admin_level, admin_id, commodity_name, methodology); the per-(country, tile) sub-workflows are -cached and `--concurrency` bounds how many tiles are in flight at once. +commodity area, peatland-occupation emissions, and total emissions. Returns a pandas.DataFrame +indexed by `INDEX_NAMES`, (admin_id, admin_level, assessment_year, commodity_name, methodology); +the per-(country, tile) sub-workflows are cached and `--concurrency` bounds how many tiles are in +flight at once. The countries come from the positional ISO 3166 alpha-3 codes, or -- with `--backfill` -- from every country in the World Bank admin-0 layer. @@ -46,6 +47,15 @@ class Methodology(enum.IntEnum): STATISTICAL = enum.auto() +INDEX_NAMES = ( + "admin_id", + "admin_level", + "assessment_year", + "commodity_name", + "methodology", +) + + def merge_dfs(*dfs: pandas.DataFrame) -> pandas.DataFrame: combined = pandas.concat(list(dfs)) jurisdiction_name = ( @@ -127,8 +137,9 @@ def workflow_for_tile_id(tile_id: str) -> list[pandas.DataFrame]: ret = ( merge_dfs(all_jurisdiction_commodities, *dfs) - .assign(methodology=methodology.name) - .set_index("methodology", append=True) + .assign(assessment_year=assessment_year, methodology=methodology.name) + .set_index(["assessment_year", "methodology"], append=True) + .reorder_levels(list(INDEX_NAMES)) .sort_index() ) assert not ret.index.duplicated().any() diff --git a/jdluc/trace.py b/jdluc/trace.py index e3ea295..90fcf17 100644 --- a/jdluc/trace.py +++ b/jdluc/trace.py @@ -13,9 +13,9 @@ It then derives, identically for all three, 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_id, methodology), where -`commodity_id` is a FAOSTAT item code or a `taxonomy.InventedCommodityId`, and the `commodity_name` -column holds its name: FAOSTAT's, from its crop table, or `taxonomy`'s. +pandas.DataFrame indexed by `CANONICAL_KEY`, (admin_id, admin_level, assessment_year, commodity_id, +methodology), where `commodity_id` is a FAOSTAT item code or a `taxonomy.InventedCommodityId`, +and the `commodity_name` column holds its name: FAOSTAT's, from its crop table, or `taxonomy`'s. The countries come from the positional ISO 3166 alpha-3 codes, or -- with `--backfill` -- from every country in the World Bank admin-0 layer. @@ -53,7 +53,13 @@ KG_PER_TONNE = 1000 -CANONICAL_KEY = ("admin_level", "admin_id", "commodity_id", "methodology") +CANONICAL_KEY = ( + "admin_id", + "admin_level", + "assessment_year", + "commodity_id", + "methodology", +) ADDITIVE_COLUMNS = ( "commodity_hectares", "emissions_mt", @@ -291,6 +297,7 @@ def attach_ratios(df: pandas.DataFrame) -> pandas.DataFrame: def iter_national_from_provincials( + assessment_year: int, iso_3166: str, methodology: attribute.Methodology, national_name: str, @@ -317,6 +324,7 @@ def safe_div(numer: float, denom: float) -> float: numer=ret["emissions_mt"] * KG_PER_TONNE, denom=ret["production_kg"] ), "methodology": methodology.name, + "assessment_year": assessment_year, } # type: ignore @@ -335,12 +343,7 @@ def workflow( iso_3166s=iso_3166s, methodology=methodology, ) - assert emissions.index.names == [ - "admin_level", - "admin_id", - "commodity_name", - "methodology", - ] + assert emissions.index.names == list(attribute.INDEX_NAMES) match methodology: case attribute.Methodology.JURISDICTIONAL_DIRECT: production = derive_jurisdictional_production_kg( @@ -378,6 +381,7 @@ def workflow( national for iso_3166 in iso_3166s for national in iter_national_from_provincials( + assessment_year=assessment_year, iso_3166=iso_3166, methodology=methodology, national_name=str( diff --git a/tools/measure-drift.py b/tools/measure-drift.py index c90cd08..2cab418 100644 --- a/tools/measure-drift.py +++ b/tools/measure-drift.py @@ -100,9 +100,6 @@ PROVINCIAL = "PROVINCIAL" # The emissions columns are megatonnes while an emissions factor is kg CO2e per kg 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_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. @@ -320,7 +317,7 @@ def get_dirty_paths(cwd: pathlib.Path) -> list[str]: def load_capture( directory: pathlib.Path, target: Target ) -> tuple[pandas.DataFrame, tuple[str, ...]]: - """Read one side's parquet onto CANONICAL_KEY, whatever index it was written with. + """Read one side's parquet onto `trace.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. A side keyed by each leg's own commodity @@ -333,20 +330,21 @@ def load_capture( 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 "assessment_year" not in df.columns: + # NB: written before `trace` keyed on the assessment year, when every table was 2020's + df["assessment_year"] = 2020 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"] - ), + df=df.set_index(list(attribute.INDEX_NAMES)), methodology=target.methodology, ) ).reset_index() - missing = [name for name in CANONICAL_KEY if name not in df.columns] + missing = [name for name in trace.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?" ) - return df.set_index(list(CANONICAL_KEY)).sort_index(), written + return df.set_index(list(trace.CANONICAL_KEY)).sort_index(), written def check_alignment( @@ -372,7 +370,7 @@ def check_alignment( message=( f"{slug:s} has no shared rows, so NOTHING was compared -- the keys do not " f"line up. Written as {' + '.join(written[0])} against " - f"{' + '.join(written[1])}, re-keyed onto {' + '.join(CANONICAL_KEY)}" + f"{' + '.join(written[1])}, re-keyed onto {' + '.join(trace.CANONICAL_KEY)}" ), severity=Severity.BLOCKING, ) diff --git a/validation/__tests__/schema_test.py b/validation/__tests__/schema_test.py index 1c9ca8a..6d576ff 100644 --- a/validation/__tests__/schema_test.py +++ b/validation/__tests__/schema_test.py @@ -18,7 +18,7 @@ tools copy from the same source, so pinning each to jdluc keeps them consistent with each other. """ -from jdluc import attribute, trace +from jdluc import attribute, emit, trace from jdluc.datasets import worldbank_jurisdictions from validation import schema @@ -55,3 +55,8 @@ def test_kg_per_tonne_matches_trace() -> None: def test_canonical_key_matches_the_index_trace_writes() -> None: """Re-keying onto the wrong names would drop every row of a comparison.""" assert schema.CANONICAL_KEY == trace.CANONICAL_KEY + + +def test_assessment_year_is_one_the_pipeline_supports() -> None: + """A year the pipeline cannot produce would leave every comparison empty.""" + assert schema.ASSESSMENT_YEAR in emit.ASSESSMENT_YEARS diff --git a/validation/capture.py b/validation/capture.py index 5b5c08f..9e8352e 100644 --- a/validation/capture.py +++ b/validation/capture.py @@ -4,10 +4,10 @@ `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_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. +`trace.CANONICAL_KEY` -- `(admin_id, admin_level, assessment_year, 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. **Every commodity, not just the targets.** `attribute.get_commodity_names` fixes the commodity set per leg, so a capture covers every commodity the leg models rather than only the pairs @@ -15,8 +15,9 @@ `statistical.get_downscaled_luc_emissions` is keyed on the tile and the assessment year alone, so the expensive per-tile layer is shared and a crop adds only its share arithmetic. -**Merging is per (ISO, methodology), and rows keep their own `code_version`.** A capture scoped with -`--isos` leaves other countries in place, so the artifact can span versions and a reader has to be +**Merging is per (ISO, methodology, assessment year), and rows keep their own `code_version`.** A +capture scoped with `--isos` leaves other countries and years in place, so the artifact can span +versions and a reader has to be able to see that -- `report.get_provenance_findings` raises a DEFECT when it does. That is why provenance travels in the row rather than in a sidecar. @@ -126,25 +127,31 @@ def merge_efs( ) -> pandas.DataFrame: """Replace what was recomputed, keep everything else exactly as it was. - Keyed on (ISO, methodology) rather than on the whole index: a country recomputed for one leg - must lose all of that leg's rows, including any crop that no longer produces one, or a stale row - would survive beside its replacements and be summed with them. + Keyed on (ISO, methodology, assessment year) rather than on the whole index: a country + recomputed for one leg and year must lose all of those rows, including any crop that no longer + produces one, or a stale row would survive beside its replacements and be summed with them. """ if existing is None: return captured + assert existing.index.names == list(schema.CANONICAL_KEY), ( + f"{pull.EFS} is indexed by {existing.index.names}, not {schema.CANONICAL_KEY}; " + "delete it and capture again" + ) recomputed = { - (str(admin_id)[:3], str(methodology)) - for admin_id, methodology in zip( + (str(admin_id)[:3], str(methodology), int(assessment_year)) + for admin_id, methodology, assessment_year in zip( captured.index.get_level_values("admin_id"), captured.index.get_level_values("methodology"), + captured.index.get_level_values("assessment_year"), strict=True, ) } keep = [ - (str(admin_id)[:3], str(methodology)) not in recomputed - for admin_id, methodology in zip( + (str(admin_id)[:3], str(methodology), int(assessment_year)) not in recomputed + for admin_id, methodology, assessment_year in zip( existing.index.get_level_values("admin_id"), existing.index.get_level_values("methodology"), + existing.index.get_level_values("assessment_year"), strict=True, ) ] diff --git a/validation/prepare.py b/validation/prepare.py index 4d5561a..7f11711 100644 --- a/validation/prepare.py +++ b/validation/prepare.py @@ -746,14 +746,17 @@ def get_anchor_versions(measure: targets.Measure) -> dict[str, str]: def read_efs() -> pandas.DataFrame | None: - """The captured emissions factors, or None before a capture has run. + """The captured emissions factors for `schema.ASSESSMENT_YEAR`, or None before a capture has + run. Read here rather than through `validation.capture` so the reporting path never imports the - module that runs the pipeline. + module that runs the pipeline. Only the one year is kept, since the comparisons below pivot and + sum over the other levels and would add years together. """ if not pull.EFS.exists(): return None - return pandas.read_parquet(pull.EFS) + efs = pandas.read_parquet(pull.EFS) + return efs[efs.index.get_level_values("assessment_year") == schema.ASSESSMENT_YEAR] # The MapSPAM code each captured crop row is compared under. The capture keys on FAOSTAT items, and diff --git a/validation/schema.py b/validation/schema.py index 5c64a90..f3e37f0 100644 --- a/validation/schema.py +++ b/validation/schema.py @@ -33,8 +33,17 @@ RESULT_BEARING_PATHS = ("jdluc", "pyproject.toml", "uv.lock") # 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_id", "methodology") +# is a key level because both legs share one artifact, so it separates sLUC from jdLUC, and so is +# `assessment_year`, because captures of several years can share it too. +CANONICAL_KEY = ( + "admin_id", + "admin_level", + "assessment_year", + "commodity_id", + "methodology", +) +# The assessment year the comparisons read from efs.parquet +ASSESSMENT_YEAR = 2020 NATIONAL = "NATIONAL" PROVINCIAL = "PROVINCIAL" # The two legs, as `methodology` records them in efs.parquet. From 5ff13ca4ccb6e6dde99215b7241d771e88d6bdf5 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 15:39:11 -0800 Subject: [PATCH 11/16] Ensure tools and validation call with the correct arguments --- tools/export-graphs.py | 1 + tools/measure-drift.py | 18 ++++++++++++++---- validation/capture.py | 23 +++++++++++++++++++---- 3 files changed, 34 insertions(+), 8 deletions(-) diff --git a/tools/export-graphs.py b/tools/export-graphs.py index fd55f0d..dfef43f 100644 --- a/tools/export-graphs.py +++ b/tools/export-graphs.py @@ -151,6 +151,7 @@ def iter_dfs( df=trace.workflow( assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, commodity_names=workflow_commodity_names, + concurrency=attribute.DEFAULT_CONCURRENCY, iso_3166s=(iso_3166,), methodology=methodology, ), diff --git a/tools/measure-drift.py b/tools/measure-drift.py index 2cab418..bb66d6a 100644 --- a/tools/measure-drift.py +++ b/tools/measure-drift.py @@ -61,7 +61,7 @@ import pandas -from jdluc import attribute, config, emit, storage, taxonomy, trace +from jdluc import attribute, config, storage, taxonomy, trace logger = logging.getLogger(__name__) @@ -923,15 +923,25 @@ def capture(directory: pathlib.Path, label: Side, targets: tuple[Target, ...]) - f"({index:d}/{len(targets):d})" ) - retired = {"skip_glad_crop_filter": False} + # NB: a baseline runs this copy of the tool against its own code, so each of these goes + # only to a `trace.workflow` that takes it: one argument retired since, one introduced since + name_to_value = { + # NB: written out, since a baseline's `emit` may predate its default; drift is measured + # in 2020, which every table before the assessment year became a parameter holds + "assessment_year": 2020, + "skip_glad_crop_filter": False, + } parameters = inspect.signature(trace.workflow).parameters df = trace.workflow( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, commodity_names=commodity_names, concurrency=attribute.DEFAULT_CONCURRENCY, iso_3166s=(target.iso_3166,), methodology=target.methodology, - **{name: value for name, value in retired.items() if name in parameters}, + **{ + name: value + for name, value in name_to_value.items() + if name in parameters + }, ) path = directory / f"{target.slug:s}.parquet" df.to_parquet(path) diff --git a/validation/capture.py b/validation/capture.py index 9e8352e..e0a06a7 100644 --- a/validation/capture.py +++ b/validation/capture.py @@ -162,8 +162,12 @@ def merge_efs( return pandas.concat([existing[keep], captured]).sort_index() -def workflow(iso_3166s: tuple[str, ...], repo_root: pathlib.Path) -> pandas.DataFrame: - """Both legs over the given countries, stamped with the code that produced them.""" +def workflow( + assessment_year: int, iso_3166s: tuple[str, ...], repo_root: pathlib.Path +) -> pandas.DataFrame: + """Both legs over the given countries in `assessment_year`, stamped with the code that + produced them. + """ frames = [] for methodology, wanted in iter_methodologies(iso_3166s=iso_3166s): commodity_names = attribute.get_commodity_names(methodology=methodology) @@ -173,7 +177,7 @@ def workflow(iso_3166s: tuple[str, ...], repo_root: pathlib.Path) -> pandas.Data ) frames.append( trace.workflow( - assessment_year=emit.DEFAULT_ASSESSMENT_YEAR, + assessment_year=assessment_year, commodity_names=commodity_names, concurrency=attribute.DEFAULT_CONCURRENCY, iso_3166s=wanted, @@ -223,6 +227,13 @@ def main() -> int: level=logging.INFO, ) parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--assessment-year", + choices=emit.ASSESSMENT_YEARS, + default=schema.ASSESSMENT_YEAR, + help="the assessment year to capture; the comparisons read only schema.ASSESSMENT_YEAR", + type=int, + ) parser.add_argument( "--isos", default=(), @@ -249,7 +260,11 @@ def main() -> int: if args.dry_run: return 0 - captured = workflow(iso_3166s=iso_3166s, repo_root=args.repo_root) + captured = workflow( + assessment_year=int(args.assessment_year), + iso_3166s=iso_3166s, + repo_root=args.repo_root, + ) merged = merge_efs(captured=captured, existing=read_efs()) pull.CAPTURE.mkdir(exist_ok=True, parents=True) merged.to_parquet(pull.EFS) From 4e48b2d4af46c7df7932c4f9d62a49c9da0d77a6 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 15:44:05 -0800 Subject: [PATCH 12/16] update docs --- README.md | 2 +- docs/data.md | 6 +++--- docs/further_research.md | 6 +++--- docs/methodology.md | 29 ++++++++++++++++++++++------- docs/validation.md | 4 ++-- jdluc/emit.py | 2 +- 6 files changed, 32 insertions(+), 17 deletions(-) diff --git a/README.md b/README.md index e5e0def..b6e36c1 100644 --- a/README.md +++ b/README.md @@ -81,7 +81,7 @@ export GOOGLE_APPLICATION_CREDENTIALS="${HOME}/.config/gcloud/application_defaul ### Running the pipeline -Every stage is scoped by one or more **ISO 3166 alpha-3 country codes**, or by `--backfill` for a global run. Stage 4 additionally takes a `--methodology-name` (`STATISTICAL`, `JURISDICTIONAL_DIRECT` or `JURISDICTIONAL_DIRECT_PROXY`), and stage 5 takes the first two. The example below reproduces Honduras via the statistical leg. Honduras is the cheapest end-to-end run in the corpus: the unit of work is a 10° tile, and it falls inside a single one (`20N_090W`). Swap in `USA` for the jurisdictional-direct leg, which needs the US-only CDL and NASS inputs: +Every stage is scoped by one or more **ISO 3166 alpha-3 country codes**, or by `--backfill` for a global run. Stage 4 additionally takes a `--methodology-name` (`STATISTICAL`, `JURISDICTIONAL_DIRECT` or `JURISDICTIONAL_DIRECT_PROXY`), and stage 5 takes the first two. Stages 3–5 also take `--assessment-year` (2020, 2022 or 2024; 2020 by default), the year the emissions factors describe and the end of their 20-year lookback. Ingest and harmonize keep every year of each dataset, so they take no year and serve every assessment year alike. The example below reproduces Honduras via the statistical leg. Honduras is the cheapest end-to-end run in the corpus: the unit of work is a 10° tile, and it falls inside a single one (`20N_090W`). Swap in `USA` for the jurisdictional-direct leg, which needs the US-only CDL and NASS inputs: ```bash # 1. Ingest each source dataset for the countries (positional: dataset, then ISO codes). diff --git a/docs/data.md b/docs/data.md index 537226a..e97b123 100644 --- a/docs/data.md +++ b/docs/data.md @@ -57,7 +57,7 @@ harmonized = xarray.open_zarr("/.zarr", consolidated=Fa Written by the `emit` stage; a pipeline output, not a published artifact (see [Access](#access)). -The crop-agnostic, per-pixel LUC emissions produced by the `emit` stage, undiscounted. A conversion is dated to the year its source class ended, `conversion-year`, and each allocation leg discounts the emissions by that year (see [`methodology.md`](methodology.md)). +The crop-agnostic, per-pixel LUC emissions produced by the `emit` stage for one assessment year, undiscounted. A conversion is dated to the year its source class ended, `conversion-year`, within the 20 years ending at the assessment year, and each allocation leg discounts the emissions by that year (see [`methodology.md`](methodology.md)). | Variable | Type | Units | Description | | ---------------------------------------------- | ------- | --------------- | --------------------------------------------------------------------------------------------------------------------------- | @@ -82,7 +82,7 @@ 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_id, methodology)`, where `admin_level` is `PROVINCIAL` (World Bank admin-1) or `NATIONAL` (admin-0, summed from provincial rows), `methodology` is `JURISDICTIONAL_DIRECT`, `STATISTICAL` or `JURISDICTIONAL_DIRECT_PROXY`, 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. All three 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, as do the proxy leg's 16 fodder crops (-101 to -116), which FAOSTAT has no item for. [`coverage.md`](coverage.md#crops) lists each leg's crops with their codes. +The final per-(jurisdiction, commodity) table produced by the `trace` stage. Indexed by `(admin_id, admin_level, assessment_year, commodity_id, methodology)`, where `assessment_year` is the year the factors describe (2020, 2022 or 2024), `admin_level` is `PROVINCIAL` (World Bank admin-1) or `NATIONAL` (admin-0, summed from provincial rows), `methodology` is `JURISDICTIONAL_DIRECT`, `STATISTICAL` or `JURISDICTIONAL_DIRECT_PROXY`, 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. All three 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, as do the proxy leg's 16 fodder crops (-101 to -116), which FAOSTAT has no item for. [`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 every leg 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`, `UNATTRIBUTED` and the fodder crops have no production figure, so none of them has an emissions factor. See "Rows that name no commodity" in [`methodology.md`](methodology.md). @@ -109,4 +109,4 @@ import pandas emission_factors = pandas.read_parquet("emissions-factors.parquet") ``` -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, a CROPGRIDS crop for the proxy 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 the herd's heads: for the statistical leg `production_mt` and the heads at each MapSPAM snapshot, `herd_heads_2000` to `herd_heads_2020`, and for the proxy leg the heads in 2020, `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`). +The `attribute` stage writes an intermediate rollup parquet indexed by `(admin_id, admin_level, assessment_year, 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, a CROPGRIDS crop for the proxy 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 the herd's heads: for the statistical leg `production_mt` and the heads at each MapSPAM snapshot, `herd_heads_2000` to `herd_heads_2020`, and for the proxy leg the heads in 2020, `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/further_research.md b/docs/further_research.md index eed858a..d21490b 100644 --- a/docs/further_research.md +++ b/docs/further_research.md @@ -78,7 +78,7 @@ It is not a lag between clearing and grass: in Boquerón, 38–56% of forest cle ## sLUC grassland attribution (expansion-share over-attribution) -**Issue**: The statistical (sLUC) leg attributes each coarse cell's conversion emissions to crops by their share of crop *expansion* within the cell (`statistical.py`; see §3, Statistical). A post-refactor validation found this hands corn/soy/wheat substantially more grassland conversion than the per-pixel jurisdictional-direct (jdLUC) leg does on the *same* emissions layer (grassland EF: corn 0.032 vs 0.012, soy 0.126 vs 0.054 kg CO₂e/kg), driving the US corn+soy+wheat total to ~76 MtCO₂e — the highest of the satellite EF models and above jdLUC (61). Grassland (30.9 Mt) is the largest single US source and the biggest change from the refactor, yet it has **no reliable external anchor** since WRI is forest-only. The one weak reference — EPA's whole-cropland Grassland→Cropland, **10.6 Mt for 2020** — sits *below* the sLUC corn+soy+wheat *subset* (30.9): a three-crop subset nearly **tripling** the whole-cropland total. Both sides are annual, so this is not muddied by a temporal basis. The GHGP per-year weights integrate to exactly 1.0 across the 20-year window, so `YEAR_TO_LINEAR_DISCOUNT_WEIGHT` allocates each conversion's emissions to a *single sourcing year* rather than accumulating twenty of them, and the production denominator is a weighted mean over the same window — an annual numerator over an annual denominator. The one real mismatch is scope, and it runs against us: widening sLUC from three crops to all cropland can only raise 30.9. Whether expansion-share allocation over-attributes grassland — or the per-pixel CDL leg under-attributes it — is unresolved. +**Issue**: The statistical (sLUC) leg attributes each coarse cell's conversion emissions to crops by their share of crop *expansion* within the cell (`statistical.py`; see §3, Statistical). A post-refactor validation found this hands corn/soy/wheat substantially more grassland conversion than the per-pixel jurisdictional-direct (jdLUC) leg does on the *same* emissions layer (grassland EF: corn 0.032 vs 0.012, soy 0.126 vs 0.054 kg CO₂e/kg), driving the US corn+soy+wheat total to ~76 MtCO₂e — the highest of the satellite EF models and above jdLUC (61). Grassland (30.9 Mt) is the largest single US source and the biggest change from the refactor, yet it has **no reliable external anchor** since WRI is forest-only. The one weak reference — EPA's whole-cropland Grassland→Cropland, **10.6 Mt for 2020** — sits *below* the sLUC corn+soy+wheat *subset* (30.9): a three-crop subset nearly **tripling** the whole-cropland total. Both sides are annual, so this is not muddied by a temporal basis. The GHGP per-year weights integrate to exactly 1.0 across the 20-year window, so `emit.get_year_to_linear_discount_weight` allocates each conversion's emissions to a *single sourcing year* rather than accumulating twenty of them, and the production denominator is a weighted mean over the same window — an annual numerator over an annual denominator. The one real mismatch is scope, and it runs against us: widening sLUC from three crops to all cropland can only raise 30.9. Whether expansion-share allocation over-attributes grassland — or the per-pixel CDL leg under-attributes it — is unresolved. **Not the denominator of the expansion share.** That denominator is the obvious suspect and is not the cause: removing it leaves the US ratio of the statistical leg's total factor to the jurisdictional-direct leg's at 1.15–1.75×, and the same ratio for grassland at 0.97–2.34×. The figures quoted above reproduce from the current parquets (corn 0.0317 vs 0.0131, soy 0.1255 vs 0.0576), so the gap is stable and unexplained. All of these figures were measured with destination layers other than GACED30 and Global Pasture Watch, and need measuring again on the current pipeline. @@ -231,12 +231,12 @@ The 2010 snapshot ends the 2005–2010 span and starts the expansion both later **Issue**: The jurisdictional-direct leg divides emissions by `commodity_hectares x yield_kg_per_ha`, where the yield is a USDA NASS survey statistic and the hectares come from the CDL. Every one of the eleven crops `jurisdictional_direct.Crop` models carries a yield, so every one produces an emissions factor. Three properties of how that denominator is assembled are worth recording, because each is an assumption rather than a measurement. 1. **Unmatched (admin, crop) pairs are still only a log line.** NASS suppresses small and identifying cells under Title 7 U.S.C. and CIPSEA, and does not survey every crop in every state. Where the CDL sees a crop in a state that NASS does not publish a yield for, that state contributes hectares and emissions to the national rollup but no production, so the national factor is a ratio over the reporting states alone. The count of such rows is logged and then discarded. -2. **Cotton is converted between product forms by a single constant.** NASS reports cotton yield as ginned lint, while MapSPAM, FAOSTAT (item 328, not lint's 767) and WRI all carry seed cotton, and either mistake scales a comparison by the milling yield. `usda_nass_quickstats.Series.reported_fraction` divides the lint yield back up to the unginned crop at 0.36. That figure is what FAOSTAT's US seed-cotton yield implies against the NASS series over 2011–2020, and it agrees with the conventional ~35% gin turnout, but the same comparison implies ~0.403 over 2000–2010. The constant is therefore tied to the window `trace.NASS_YIELD_YEARS` currently averages, and a national constant hides whatever provincial variation in gin turnout exists. +2. **Cotton is converted between product forms by a single constant.** NASS reports cotton yield as ginned lint, while MapSPAM, FAOSTAT (item 328, not lint's 767) and WRI all carry seed cotton, and either mistake scales a comparison by the milling yield. `usda_nass_quickstats.Series.reported_fraction` divides the lint yield back up to the unginned crop at 0.36. That figure is what FAOSTAT's US seed-cotton yield implies against the NASS series over 2011–2020, and it agrees with the conventional ~35% gin turnout, but the same comparison implies ~0.403 over 2000–2010. The constant is therefore tied to the years the direct leg averages yields over, the 3 centered on the assessment year, and a national constant hides whatever provincial variation in gin turnout exists. 3. **Two series are not quite the crop the CDL scores.** NASS renamed its rolled-up dry-bean class in 2019 and publishes no `ALL CLASSES` to fall back on, so spanning a multi-year window takes both names; the pre-2019 aggregate includes chickpeas, which the CDL scores separately as `CHICK_PEAS`. Sugarcane publishes no `ALL UTILIZATION PRACTICES` row either, so its total is `SUGAR & SEED` and includes seed cane. **Potential impact:** Modest overall, and concentrated in one crop. Cotton is the pair whose assumption moves a published number materially: carried as lint rather than seed cotton its factor would read roughly 2.8x higher. The eight crops outside maize, soybean and wheat carry 9% of US jurisdictional-direct emissions between them, so an error in their denominators is bounded by that share. The residual coverage gap is concentrated in minor producing states for the smaller crops, so it biases a national factor by whatever those states' yields differ from the surveyed ones — likely small for maize and soybean, larger for crops NASS surveys in only a handful of states. The chickpea contamination affects half the current averaging window for one crop carrying 0.627 Mt. -**Potential improvement path:** Carry the unmatched count into the reported output rather than the log, so a factor computed over partial state coverage is visible to whoever reads it — the same treatment `validation.md` gives its other coverage measures. For cotton, derive the lint fraction per state and year from NASS's own cottonseed production rather than a single national constant, which would also make the pre-2011 drift moot. Re-measure the constant if `NASS_YIELD_YEARS` ever moves earlier. +**Potential improvement path:** Carry the unmatched count into the reported output rather than the log, so a factor computed over partial state coverage is visible to whoever reads it — the same treatment `validation.md` gives its other coverage measures. For cotton, derive the lint fraction per state and year from NASS's own cottonseed production rather than a single national constant, which would also make the pre-2011 drift moot. Re-measure the constant for an assessment year whose yield window falls before 2011. ## Higher resolution yield data diff --git a/docs/methodology.md b/docs/methodology.md index a4a9b29..1b1ef9d 100644 --- a/docs/methodology.md +++ b/docs/methodology.md @@ -91,7 +91,7 @@ Every input is ingested from its upstream publisher into cloud storage as tiled | **SoilGrids OCS (0–30 cm)** | Soil organic carbon stock | ISRIC SoilGrids (WCS) | raster | | **GNW Global Peatlands** | Binary peatland mask | GNW data-api (`gfw_peatlands` v20230315) | raster | | **IPCC climate zones** | Climate domain per pixel (drives IPCC factors) | Zenodo (doi:10.5281/zenodo.7303808) | raster | -| **USDA NASS CDL** | US per-pixel crop identity (jurisdictional-direct leg) | USDA Cropland Data Layer | raster | +| **USDA NASS CDL** | US per-pixel crop identity, 2020–2025 (jurisdictional-direct leg) | USDA Cropland Data Layer | raster | | **IFPRI MapSPAM** | Global per-crop physical area + production, 2000/2005/2010/2020 (statistical leg) | Harvard Dataverse | raster | | **USDA NASS QuickStats** | State-level crop yields (jurisdictional-direct production) | NASS QuickStats API | tabular | | **World Bank Official Boundaries** | Admin-0/1/2 jurisdiction polygons | World Bank | vector | @@ -110,10 +110,25 @@ GPW grassland is native to the GLAD tile grid: its published mosaics are 0.00025 Descals oil palm is one of three per-crop extent layers, beside FDaP's cocoa and coffee, and the only one that dates its own pixels: each 30 m pixel holds the year the plantation standing in 2021 was established, 0 where there is no oil palm. **It is read as extent, not as a date.** A pixel is an oil palm destination where its planting year is non-zero and no later than the assessment year; dating stays with tree-cover loss and the last grassland departure, as for every other conversion. The planting year does not date one, because a replant reads as the year it was replanted rather than as the year the land was first cleared — it bounds when the current crop was established, not when conversion happened. The upper bound is what tests the destination at the assessment year: a plantation established after 2020 is not yet evidence of what an earlier clearance was for, which excludes 2.2% of mapped palm globally. Extent accuracy is high on industrial plantations (producer's 91.0 ± 2.5%, user's 91.8 ± 1.2%) and much weaker on smallholders (71.4 ± 0.7% and 72.4 ± 1.8%); the planting year itself carries an RMSE of 2.65 years against field data. The publisher screens globally and then publishes only the 609 hundred-kilometre cells where oil palm was found, so 232 of the 280 ten-degree tiles touch no cell and are not ingested at all; a reader that opts into missing tiles sees them as zero, which is the right answer here, but it is the screen's recall that answers rather than a per-pixel classification. The companion 10 m industrial-versus-smallholder extent layer from the same record is not ingested: its palm footprint agrees with the planting-year layer to within 0.4%, so the only thing it adds is the two-way split. -The FDaP masks are the cocoa and coffee destinations between 24°S and 24°N, read at 2020; `jdluc/datasets/fdap.py` describes how ingest thresholds them. Harmonize subsamples the 10 m masks onto the GLAD tile grid by nearest, which keeps their area. +The FDaP masks are the cocoa and coffee destinations between 24°S and 24°N, read at the assessment year; `jdluc/datasets/fdap.py` describes how ingest thresholds them. Harmonize subsamples the 10 m masks onto the GLAD tile grid by nearest, which keeps their area. MapSPAM is the newest addition and the one that makes the statistical leg possible: it downscales sub-national crop statistics to a ~10 km (5 arc-minute) grid via a cross-entropy allocation, publishing physical area and production per crop for 2000, 2005, 2010, and 2020 (note: **no 2015 snapshot**). Its crop taxonomy is coarser in earlier years (the 2000 snapshot reports only 21, partly grouped, crops), which the statistical leg reconciles to a common per-crop taxonomy before use (see §3). +### Assessment years + +The pipeline produces emissions factors for the assessment years 2020, 2022 and 2024 (`emit.ASSESSMENT_YEARS`, chosen with `--assessment-year`, 2020 by default). The dated layers bound them: the 20-year lookback needs tree-cover loss and GPW grassland in every one of its years, the destinations need GACED30 and FDaP in the assessment year itself, and FDaP is published for even years only. The jurisdictional-direct leg also reads the CDL in the assessment year. + +Ingest and harmonize keep every year of each layer, so every assessment year reads the same harmonized tiles; only `emit` and the stages after it depend on the year. Some inputs have no snapshot after 2020, so a later assessment year holds them at their latest: + +| Input | Held at | Where it is read | +| ---------------- | -------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| Descals oil palm | The plantations standing in 2021 | The oil palm destination, which takes every plantation established by the assessment year | +| MapSPAM | 2020 | Statistical leg: conversions after 2020 take the 2010→2020 expansion shares, their production and area weigh on the 2020 snapshot, and occupation shares are 2020's | +| GPW livestock | 2020 | Statistical and proxy legs: grazing shares and herds | +| CROPGRIDS | 2020 | Proxy leg: crop-area shares, with the pasture fraction read in the same year | + +The production windows follow the year instead: the direct leg's NASS yields and the proxy leg's FAOSTAT production are each averaged over the 3 years centered on the assessment year, so those two legs' factors for 2024 need 2025's published statistics. + ### Three grids Every source is ingested in `EPSG:4326` and resampled onto one of three common grids, depending on which stage consumes them: @@ -174,7 +189,7 @@ For each emissive transition, emissions are the sum of carbon lost from vegetati ### Allocation to crop years (GHGP linear discounting) -All fluxes are first computed as if the conversion were instantaneous, then allocated over time using the 20-year linearly-discounted lookback prescribed by the GHGP Land Sector and Removals Standard (§7.2.1): the more recently a conversion occurred, the more heavily its emissions are weighted. A conversion is dated to the year its source class ended and takes that year's linear amortization factor from the LSR Guidance (Table 7.11; `YEAR_TO_LINEAR_DISCOUNT_WEIGHT` in `emit.py`): 9.75% for the assessment year, falling 0.5 points a year to 0.25% for 2001, summing to 100% over the window. `emit` writes the emissions undiscounted beside their year, and each leg discounts every pixel by its own year before clipping or downscaling it. +All fluxes are first computed as if the conversion were instantaneous, then allocated over time using the 20-year linearly-discounted lookback prescribed by the GHGP Land Sector and Removals Standard (§7.2.1): the more recently a conversion occurred, the more heavily its emissions are weighted. A conversion is dated to the year its source class ended and takes that year's linear amortization factor from the LSR Guidance (Table 7.11; `get_year_to_linear_discount_weight` in `emit.py`): 9.75% for the assessment year, falling 0.5 points a year to 0.25% nineteen years before it (2001 for the 2020 assessment), summing to 100% over the window. `emit` writes the emissions undiscounted beside their year, and each leg discounts every pixel by its own year before clipping or downscaling it. The per-pixel result is the discounted conversion emissions plus the current year's peatland occupation emissions, an annual flux that is not discounted, scaled by pixel area. This same temporal ramp is shared by both attribution methodologies. @@ -226,7 +241,7 @@ Concretely, for each (jurisdiction, crop): A crop's row covers only the pixels the crop map assigns to that crop, and it never includes emissions from peat drained under pasture, even where the crop map puts the crop on the same pixel. Everything else goes to `UNATTRIBUTED`: all emissions on pixels the crop map does not assign to a crop we model, and all emissions from peat drained under pasture. -Because direct attribution relies only on **current** production, production is later computed as crop area × NASS QuickStats yield (a multi-year state-level mean), converting NASS's reported yield per acre to kg/ha. NASS publishes that yield in four units — bushels, pounds, hundredweight and short tons per acre — and which one applies is a property of the commodity, so `usda_nass_quickstats.CropSeries` carries the unit and its weight per crop. The bushel weights are marketing bushels, from table 6 of USDA Agricultural Handbook 697 (ERS, June 1992), rather than the grading test weights of 7 CFR 810, which are a different quantity. Cotton is the one crop whose published yield is not the harvested crop: NASS reports ginned lint where MapSPAM and FAOSTAT carry seed cotton, so the lint yield is divided by a lint fraction of 0.36. This leg currently asserts `iso_3166 == "USA"`, since CDL is US-only. +Because direct attribution relies only on **current** production, production is later computed as crop area × NASS QuickStats yield (the state-level mean over the 3 years centered on the assessment year), converting NASS's reported yield per acre to kg/ha. NASS publishes that yield in four units — bushels, pounds, hundredweight and short tons per acre — and which one applies is a property of the commodity, so `usda_nass_quickstats.CropSeries` carries the unit and its weight per crop. The bushel weights are marketing bushels, from table 6 of USDA Agricultural Handbook 697 (ERS, June 1992), rather than the grading test weights of 7 CFR 810, which are a different quantity. Cotton is the one crop whose published yield is not the harvested crop: NASS reports ginned lint where MapSPAM and FAOSTAT carry seed cotton, so the lint yield is divided by a lint fraction of 0.36. This leg currently asserts `iso_3166 == "USA"`, since CDL is US-only. ### Statistical (`statistical.py`) @@ -235,7 +250,7 @@ When crop production is only available at low spatial resolution, spatial inters Concretely: 1. Discount each pixel's emissions by its conversion year and downsample them to the MapSPAM ~10 km grid — forest and grassland conversion emissions and the peat drainage pulse split by destination (cropland or pastureland), each in one band per MapSPAM span holding the conversion year, the cropland and pastureland peatland-occupation bands, and pasture extent at each MapSPAM snapshot. -2. For each MapSPAM span, compute each commodity's **expansion** and its **share** of total expansion in that cell. A crop's expansion is the positive change in its MapSPAM physical area; pastureland's is the positive change in the hectares Global Pasture Watch calls cultivated grassland, measured the same way. MapSPAM has no 2015 snapshot, so conversions from 2011 to 2020 all take the 2010→2020 expansion. The total is built by clipping each commodity's own movement, exactly as its numerator is, so the shares sum to **at most** one. +2. For each MapSPAM span, compute each commodity's **expansion** and its **share** of total expansion in that cell. A crop's expansion is the positive change in its MapSPAM physical area; pastureland's is the positive change in the hectares Global Pasture Watch calls cultivated grassland, measured the same way. MapSPAM has no 2015 snapshot and none after 2020, so conversions from 2011 onward all take the 2010→2020 expansion. The total is built by clipping each commodity's own movement, exactly as its numerator is, so the shares sum to **at most** one. 3. Attribute each span's discounted forest and grassland conversion emissions to commodities by that span's expansion share, and its peat drainage pulse by destination (below). 4. Compute each commodity's production denominator as a discount-weighted average over the lookback window (described below), rather than a single current-year snapshot. @@ -261,13 +276,13 @@ Where no reference year places the group in a pixel there is no basis for a spli **Tropical woody-perennial factors other than oil palm are lower bounds.** GACED30 admits woody crops in principle, but measured against independent references it recovers only a small fraction of them — 1.1% and 0.3% of mapped oil palm on two Riau windows — so a woody perennial with no layer of its own is one of the destinations bucket 4 above is blind to, and the conversion emissions reaching coconut and rubber are far below what external references report. Those crops' factors should be read as lower bounds rather than estimates until a per-crop destination layer supplies what the general cropland mask cannot. Oil palm has one, and coffee has FDaP's between 24°S and 24°N; `further_research.md` covers the mechanism for the rest, and separately a smaller allocation gap specific to oil palm. The asymmetry worth naming is that MapSPAM's expansion denominator *does* see these crops, so they dilute every other crop's share in the cells they grow in while collecting almost none of the emissions themselves. -**Windowed production denominator (a departure from WRI).** The emissions numerator is time-resolved: it links conversion emissions to each crop's *historical* expansion, weighted by recency. Dividing that by a single *current-year* production snapshot — as in WRI's methodology — would place the numerator and denominator on inconsistent time bases. A crop that expanded early in the window and then contracted would carry real conversion emissions against a shrunken (or zero) present-day production, yielding an inflated or undefined emissions factor. Instead, production is reduced over the same 20-year lookback window using the same linear temporal ramp applied to emissions: each conversion year's discount weight is split evenly between the two MapSPAM snapshots of the span holding it, and production is the snapshots' mean under those weights (`SNAPSHOT_YEAR_TO_WEIGHT` in `statistical.py`). This ties emission allocation and production to the same years with the same recency weighting. When production is flat across the window the result reduces exactly to the current-production value; it diverges only when production actually changed. This consistent linkage is a material improvement over attributing historical-expansion emissions to present-day production alone. +**Windowed production denominator (a departure from WRI).** The emissions numerator is time-resolved: it links conversion emissions to each crop's *historical* expansion, weighted by recency. Dividing that by a single *current-year* production snapshot — as in WRI's methodology — would place the numerator and denominator on inconsistent time bases. A crop that expanded early in the window and then contracted would carry real conversion emissions against a shrunken (or zero) present-day production, yielding an inflated or undefined emissions factor. Instead, production is reduced over the same 20-year lookback window using the same linear temporal ramp applied to emissions: each conversion year's discount weight is split evenly between the two MapSPAM snapshots of the span holding it, a year after 2020 gives its whole weight to the 2020 snapshot, and production is the snapshots' mean under those weights (`get_snapshot_year_to_weight` in `statistical.py`). This ties emission allocation and production to the same years with the same recency weighting. When production is flat across the window the result reduces exactly to the current-production value; it diverges only when production actually changed. This consistent linkage is a material improvement over attributing historical-expansion emissions to present-day production alone. ### Jurisdictional direct proxy (`jurisdictional_direct_proxy.py`) Each pixel's emissions are split by the land class it became, cropland or pasture, and averaged onto the CROPGRIDS grid. In each cell, each crop takes its share of the 2020 crop area of all 173 CROPGRIDS crops of what cropland is charged, and each livestock commodity takes its 2020 grazing share, as in the statistical leg, of what pasture is charged. CROPGRIDS has one year, so every span is divided by the same shares. A crop's hectares are its CROPGRIDS crop area, and `UNATTRIBUTED` is cropland's charge in cells with no CROPGRIDS crop. -Production is joined in `trace.py`, from FAOSTAT's 2019–2021 means, centered on CROPGRIDS' year. A crop's is its FAOSTAT item's national production, split among the country's provinces by their share of the country's CROPGRIDS crop area, over every crop that shares the item. `BEEF_CATTLE`'s is the cattle GPW maps in 2020 in cells holding any cultivated grassland, times FAOSTAT's national carcass weight per standing head, as in the statistical leg. FAOSTAT production where CROPGRIDS maps none of the item in that country has no province to go to; `trace.py` logs it and leaves it out. A crop's hectares are physical area while FAOSTAT counts every harvest, so `yield_kg_per_ha` reads above FAOSTAT's yield per harvested hectare for crops grown more than once a year, rice above all. +Production is joined in `trace.py`, from FAOSTAT's means over the 3 years centered on the assessment year (2019–2021 for 2020). The 2024 assessment needs FAOSTAT's 2025 production, and `trace.py` stops rather than average fewer years when a year of the window is missing. A crop's is its FAOSTAT item's national production, split among the country's provinces by their share of the country's CROPGRIDS crop area, over every crop that shares the item. `BEEF_CATTLE`'s is the cattle GPW maps in 2020 in cells holding any cultivated grassland, times FAOSTAT's national carcass weight per standing head, as in the statistical leg. FAOSTAT production where CROPGRIDS maps none of the item in that country has no province to go to; `trace.py` logs it and leaves it out. A crop's hectares are physical area while FAOSTAT counts every harvest, so `yield_kg_per_ha` reads above FAOSTAT's yield per harvested hectare for crops grown more than once a year, rice above all. ### From rollups to emissions factors (`trace.py`) diff --git a/docs/validation.md b/docs/validation.md index afa7632..d836446 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -41,7 +41,7 @@ The split is **by what the code is, not by what it computes**. `tools/` needs geopandas and reads a 93 MiB GeoPackage. Nothing imports it — the hyphenated filenames make that literal — so its cost stays off the reporting path. It runs rarely, and what it produces is a committed artifact rather than console output, because a generator whose result evaporates cannot be read by anything downstream. -`capture.py` is the only module that runs the pipeline itself, which is why it is run by hand rather than reached from `python -m validation`. Keeping it separate is what lets the report finish in seconds. +`capture.py` is the only module that runs the pipeline itself, which is why it is run by hand rather than reached from `python -m validation`. Keeping it separate is what lets the report finish in seconds. It captures one assessment year per run, `schema.ASSESSMENT_YEAR` unless `--assessment-year` names another, and `efs.parquet` keys every row on its year, so captures of several years sit side by side while the comparisons read `schema.ASSESSMENT_YEAR`'s rows alone. `report.py` reads only the frames it is handed, and computes nothing. That constraint is load-bearing: it forces every filter, rescale, factorization and cross-walk into `prepare.py`, where each is named and testable, and leaves a report section as a pivot plus a caption. @@ -141,7 +141,7 @@ Note this narrows what E3 is buying. It still gates the provincial factors `ORBA ### The reporting year is the window, not a vintage -WRI publishes `{LD,production,EF}_{2020..2024}`, and the year is not a data vintage: it names where the LSRS 20-year assessment window sits. Reporting year 2020 covers loss years 2001–2020, 2024 covers 2005–2024, and each loss year carries the linear discounting factor for its age — 9.75% for the reporting year itself, falling 0.5 pt/yr to 0.25% nineteen years back, summing to 100% (Fitts et al., 2025a, Table 4, which takes the schedule from the draft LSRG). Those are the factors in `emit.YEAR_TO_LINEAR_DISCOUNT_WEIGHT`, year for year. +WRI publishes `{LD,production,EF}_{2020..2024}`, and the year is not a data vintage: it names where the LSRS 20-year assessment window sits. Reporting year 2020 covers loss years 2001–2020, 2024 covers 2005–2024, and each loss year carries the linear discounting factor for its age — 9.75% for the reporting year itself, falling 0.5 pt/yr to 0.25% nineteen years back, summing to 100% (Fitts et al., 2025a, Table 4, which takes the schedule from the draft LSRG). Those are the factors `emit.get_year_to_linear_discount_weight` gives, year for year. **The window slides over real new loss data, so the five columns are a series.** Re-discounting a fixed 2001–2020 series would bound `LD_2024 / LD_2020` at 31/39 = 0.795 for any non-negative series; 62% of the 4,423 populated pairs exceed that bound, 864 rise monotonically and hold 49% of global `LD_2024`, and the global total is flat at 2,664 → 2,669 Mt where re-discounting would have to fall. The 19 pairs that decay to exactly zero are the kernel's fingerprint rather than a counterexample: `GEO LENT` and `Z07 OCER` both run 7 : 5 : 3 : 1 : 0, one loss year in 2004 and nothing since. The denominator slides too — `production_Y` is SPAM 2020 production rescaled by the FAOSTAT national ratio for year *Y* (Fitts et al., 2025a, step 5), which reproduces to a median 0.02% (2021), 0.05% (2022) and 0.24% (2023), then 6.2% for 2024, because WRI predicted 2024 with a per-(country, crop) random forest before FAOSTAT published it. It is not national production: over 1,133 pairs above 100 kt it is a median 12% of FAOSTAT's national figure — 52% for `BRA SOYB`, 5% for `USA MAIZ`, 93% for `CIV COCO` — which is the deforestation-linked denominator this section relies on, measured. diff --git a/jdluc/emit.py b/jdluc/emit.py index 09bb8fc..8a1ebc4 100644 --- a/jdluc/emit.py +++ b/jdluc/emit.py @@ -227,7 +227,7 @@ def get_mineral_soil_emissions( LOOKBACK_YEARS = 20 -# Bounded by the years the input datasets cover; see "Source datasets" in docs/methodology.md +# Bounded by the years the input datasets cover; see "Assessment years" in docs/methodology.md ASSESSMENT_YEARS = (2020, 2022, 2024) DEFAULT_ASSESSMENT_YEAR = 2020 assert DEFAULT_ASSESSMENT_YEAR in ASSESSMENT_YEARS From ee7ad8f5c91c5dfdaa3b4b9bd96a064fdf8bf86c Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 15:56:12 -0800 Subject: [PATCH 13/16] document datasets which are stale in 2024 --- docs/further_research.md | 8 ++++++++ docs/methodology.md | 4 ++-- jdluc/emit.py | 4 +++- jdluc/trace.py | 14 ++++++++++---- 4 files changed, 23 insertions(+), 7 deletions(-) diff --git a/docs/further_research.md b/docs/further_research.md index d21490b..9268a21 100644 --- a/docs/further_research.md +++ b/docs/further_research.md @@ -244,6 +244,14 @@ The 2010 snapshot ends the 2005–2010 span and starts the expansion both later **Potential impact:** This issue affects only the spatial distribution of emissions within each state, but could be significant at county-level. +## Oil palm planted after 2021 + +**Issue**: The Descals oil palm layer maps the plantations standing in 2021, so the 2022 and 2024 assessments cannot see palm planted after 2021. Forest cleared for those plantations takes cropland as its destination only where GACED30 calls the pixel cropland, and GACED30 recovers only a small fraction of tropical woody perennials (see `methodology.md`). Elsewhere the clearance resolves to no destination and is reported in `dropped-emissions`, so a later assessment year undercounts oil palm's conversion emissions. + +**Potential impact:** Not measured. It grows with the years between 2021 and the assessment year, and is largest where palm expanded fastest after 2021. + +**Potential improvement path:** Measure it first, by comparing the oil palm destination's converted area for 2022 with 2020's on a palm frontier tile. Then take a later release of the layer if one is published, or admit a second palm layer that maps after 2021. + ## Forest cleared for built-up land is charged to nobody **Superseded in part.** There is no built-up destination in the conversion matrix, so forest cleared for settlement now resolves to no destination and is reported in `dropped-emissions` rather than charged. That keeps it out of every published total. What remains open is that the carbon is real, is emitted, and is attributable to no commodity — bucket 5 in `methodology.md`. diff --git a/docs/methodology.md b/docs/methodology.md index 1b1ef9d..11caa3e 100644 --- a/docs/methodology.md +++ b/docs/methodology.md @@ -127,7 +127,7 @@ Ingest and harmonize keep every year of each layer, so every assessment year rea | GPW livestock | 2020 | Statistical and proxy legs: grazing shares and herds | | CROPGRIDS | 2020 | Proxy leg: crop-area shares, with the pasture fraction read in the same year | -The production windows follow the year instead: the direct leg's NASS yields and the proxy leg's FAOSTAT production are each averaged over the 3 years centered on the assessment year, so those two legs' factors for 2024 need 2025's published statistics. +The production windows follow the year instead: the direct leg's NASS yields and the proxy leg's FAOSTAT production are each averaged over the 3 years centered on the assessment year. FAOSTAT has not published 2025's production, so the proxy leg's 2024 factors take 2024's alone until it does. ### Three grids @@ -282,7 +282,7 @@ Where no reference year places the group in a pixel there is no basis for a spli Each pixel's emissions are split by the land class it became, cropland or pasture, and averaged onto the CROPGRIDS grid. In each cell, each crop takes its share of the 2020 crop area of all 173 CROPGRIDS crops of what cropland is charged, and each livestock commodity takes its 2020 grazing share, as in the statistical leg, of what pasture is charged. CROPGRIDS has one year, so every span is divided by the same shares. A crop's hectares are its CROPGRIDS crop area, and `UNATTRIBUTED` is cropland's charge in cells with no CROPGRIDS crop. -Production is joined in `trace.py`, from FAOSTAT's means over the 3 years centered on the assessment year (2019–2021 for 2020). The 2024 assessment needs FAOSTAT's 2025 production, and `trace.py` stops rather than average fewer years when a year of the window is missing. A crop's is its FAOSTAT item's national production, split among the country's provinces by their share of the country's CROPGRIDS crop area, over every crop that shares the item. `BEEF_CATTLE`'s is the cattle GPW maps in 2020 in cells holding any cultivated grassland, times FAOSTAT's national carcass weight per standing head, as in the statistical leg. FAOSTAT production where CROPGRIDS maps none of the item in that country has no province to go to; `trace.py` logs it and leaves it out. A crop's hectares are physical area while FAOSTAT counts every harvest, so `yield_kg_per_ha` reads above FAOSTAT's yield per harvested hectare for crops grown more than once a year, rice above all. +Production is joined in `trace.py`, from FAOSTAT's means over the 3 years centered on the assessment year (2019–2021 for 2020). Until FAOSTAT publishes 2025, the 2024 assessment takes 2024's production alone, and `trace.py` stops rather than average fewer years when a year of a window is missing. A crop's is its FAOSTAT item's national production, split among the country's provinces by their share of the country's CROPGRIDS crop area, over every crop that shares the item. `BEEF_CATTLE`'s is the cattle GPW maps in 2020 in cells holding any cultivated grassland, times FAOSTAT's national carcass weight per standing head, as in the statistical leg. FAOSTAT production where CROPGRIDS maps none of the item in that country has no province to go to; `trace.py` logs it and leaves it out. A crop's hectares are physical area while FAOSTAT counts every harvest, so `yield_kg_per_ha` reads above FAOSTAT's yield per harvested hectare for crops grown more than once a year, rice above all. ### From rollups to emissions factors (`trace.py`) diff --git a/jdluc/emit.py b/jdluc/emit.py index 8a1ebc4..ea35766 100644 --- a/jdluc/emit.py +++ b/jdluc/emit.py @@ -317,7 +317,9 @@ def __call__( DESTINATION_DATASET_TO_PREDICATE: dict[DestinationDataset, DestinationPredicate] = { DestinationDataset.DESCALS_OIL_PALM: lambda assessment_year, dset: ( - # NB: zero is "no oil palm" + # NB: zero is "no oil palm". The layer maps the plantations standing in 2021, so a later + # assessment year misses palm planted since; see "Oil palm planted after 2021" in + # docs/further_research.md ( (planting_year := dset[descals_oil_palm.DATASET.fully_qualified_band_name]) > 0 diff --git a/jdluc/trace.py b/jdluc/trace.py index 90fcf17..4d543bb 100644 --- a/jdluc/trace.py +++ b/jdluc/trace.py @@ -7,9 +7,9 @@ - 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. - JURISDICTIONAL_DIRECT_PROXY: production = FAOSTAT national production over the 3 years - centered on the assessment year, split among provinces by CROPGRIDS crop area for crops, and - the attribute's 2020 GPW heads x FAOSTAT carcass weight per head over the same years for - livestock. + centered on the assessment year (2024 alone until FAOSTAT publishes 2025), split among + provinces by CROPGRIDS crop area for crops, and the attribute's 2020 GPW heads x FAOSTAT + carcass weight per head over the same years for livestock. It then derives, identically for all three, each commodity's yield (kg per hectare) and emissions factor (kgCO2e per kg), and rolls the provincials up to national totals. Returns a cached @@ -356,7 +356,13 @@ def workflow( crops=faostat.load(dataset=faostat.CROP_DATASET), emissions=emissions, stocks=faostat.load(dataset=faostat.LIVESTOCK_DATASET), - years=tuple(range(assessment_year - 1, assessment_year + 2)), + # TODO: center 2024 like every other year once FAOSTAT publishes 2025's + # production; until then it takes its own year alone + years=( + (assessment_year,) + if assessment_year == 2024 + else tuple(range(assessment_year - 1, assessment_year + 2)) + ), ) case attribute.Methodology.STATISTICAL: production = derive_statistical_production_kg( From 52ea44051022eaeda5bf147b90b8255838af40fd Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 16:00:09 -0800 Subject: [PATCH 14/16] Inline some functions to shrink the statistical module --- jdluc/statistical.py | 69 +++++++++++++++++--------------------------- jdluc/trace.py | 2 +- 2 files changed, 28 insertions(+), 43 deletions(-) diff --git a/jdluc/statistical.py b/jdluc/statistical.py index 4ba43d9..0533d79 100644 --- a/jdluc/statistical.py +++ b/jdluc/statistical.py @@ -219,16 +219,6 @@ def get_downscaled_luc_emissions(assessment_year: int, tile_id: str) -> xarray.D ) -def get_pasture_hectares(dset: xarray.Dataset, year: int) -> xarray.DataArray: - darray = dset[f"pasture:fraction:{year:d}"] - return darray * emit.get_hectares_per_pixel(darray=darray) - - -def get_livestock_units_band_name(species: gpw_livestock.Species, year: int) -> str: - """The band carrying `species`' density in livestock units/ha in `year`.""" - return f"{species:s}:livestock-units-per-ha:{year:d}" - - def get_band_name_to_livestock_units( dset: xarray.Dataset, species_to_units_per_head: dict[gpw_livestock.Species, float] ) -> dict[str, xarray.DataArray]: @@ -238,7 +228,7 @@ def get_band_name_to_livestock_units( grazers; production counts heads, and reads the heads bands. """ return { - get_livestock_units_band_name(species=species, year=year): dset[ + f"{species:s}:livestock-units-per-ha:{year:d}": dset[ gpw_livestock.get_band_name(species=species, year=year) ] * species_to_units_per_head[species] @@ -255,7 +245,7 @@ def get_livestock_to_grazing_share( ratios between species' units per head matter, since a share divides each by the cell's total. """ species_to_units = { - species: dset[get_livestock_units_band_name(species=species, year=year)] + species: dset[f"{species:s}:livestock-units-per-ha:{year:d}"] for species in gpw_livestock.Species } zero = xarray.DataArray(numpy.float32(0)) @@ -275,8 +265,10 @@ def get_commodity_hectares( commodity: Commodity, dset: xarray.Dataset, year: int ) -> xarray.DataArray: if isinstance(commodity, Livestock): + pasture_fraction = dset[f"pasture:fraction:{year:d}"] return ( - get_pasture_hectares(dset=dset, year=year) + pasture_fraction + * emit.get_hectares_per_pixel(darray=pasture_fraction) * get_livestock_to_grazing_share(dset=dset, year=year)[commodity] ) else: @@ -341,9 +333,13 @@ def get_commodity_to_share( """Each commodity's share of everything a cell emitted in a span, whatever each pixel became.""" # Pastureland expands alongside the crops, measured net and clipped exactly as a crop's is: # GPW is annual, and its gross gain would charge pasture for churn MapSPAM cannot see + pasture_fraction_after = dset[f"pasture:fraction:{after:d}"] + pasture_fraction_before = dset[f"pasture:fraction:{before:d}"] pasture_expansion = ( - get_pasture_hectares(dset=dset, year=after) - - get_pasture_hectares(dset=dset, year=before) + pasture_fraction_after + * emit.get_hectares_per_pixel(darray=pasture_fraction_after) + - pasture_fraction_before + * emit.get_hectares_per_pixel(darray=pasture_fraction_before) ).clip(min=0) total_expansion = ( get_cropland_expansion(after=after, before=before, dset=dset) @@ -414,9 +410,9 @@ def get_mapspam_span(year: int) -> tuple[int, int]: def get_snapshot_year_to_weight(assessment_year: int) -> dict[int, float]: - """Each snapshot's weight in `get_discounted_snapshot_mean`: every conversion year's discount - weight, split evenly between the two snapshots of its MapSPAM span. A year after the last - snapshot gives its whole weight to that snapshot instead. + """Each snapshot's weight in the means of a crop's production and a commodity's area: every + conversion year's discount weight, split evenly between the two snapshots of its MapSPAM span. + A year after the last snapshot gives its whole weight to that snapshot instead. """ *_, last_snapshot_year = MAPSPAM_SNAPSHOT_YEARS year_to_weight = emit.get_year_to_linear_discount_weight( @@ -449,20 +445,6 @@ def get_herd_heads_column(year: int) -> str: ) -def get_discounted_snapshot_mean( - snapshot_year_to_weight: dict[int, float], - year_to_value: dict[int, xarray.DataArray], -) -> xarray.DataArray: - """A level -- area, production -- as a mean over the windows emissions are summed on.""" - return sum( - ( - snapshot_year_to_weight[year] * year_to_value[year] - for year in MAPSPAM_SNAPSHOT_YEARS - ), - start=xarray.DataArray(numpy.float32(0)), - ) - - def get_commodity_name_to_totals( commodities: tuple[Commodity, ...], dset: xarray.Dataset, @@ -558,25 +540,28 @@ def get_commodity_name_to_totals( } else: occupation = cropland_occupation * occupation_shares[commodity] - totals["production_mt"] = get_discounted_snapshot_mean( - snapshot_year_to_weight=snapshot_year_to_weight, - year_to_value={ - year: ifpri_mapspam.get_canonical_quantity( + # NB: a level, so a mean over the snapshots under the weights emissions are summed with + totals["production_mt"] = sum( + ( + snapshot_year_to_weight[year] + * ifpri_mapspam.get_canonical_quantity( canonical_crop_name=commodity.value, dset=dset, quantity=ifpri_mapspam.Quantity.PRODUCTION, year=year, ) for year in MAPSPAM_SNAPSHOT_YEARS - }, + ), + start=xarray.DataArray(numpy.float32(0)), ) - totals["commodity_hectares"] = get_discounted_snapshot_mean( - snapshot_year_to_weight=snapshot_year_to_weight, - year_to_value={ - year: get_commodity_hectares(commodity=commodity, dset=dset, year=year) + totals["commodity_hectares"] = sum( + ( + snapshot_year_to_weight[year] + * get_commodity_hectares(commodity=commodity, dset=dset, year=year) for year in MAPSPAM_SNAPSHOT_YEARS - }, + ), + start=xarray.DataArray(numpy.float32(0)), ) # NB: occupation is a flat annual rate over the area drained, so dividing it back out # recovers that area on `emit`'s 30 m intersection diff --git a/jdluc/trace.py b/jdluc/trace.py index 4d543bb..5378bfb 100644 --- a/jdluc/trace.py +++ b/jdluc/trace.py @@ -131,7 +131,7 @@ def derive_statistical_production_kg( ) -> 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. + over the window with the weights `statistical.get_snapshot_year_to_weight` gives a crop's. """ ret = emissions.copy() ret["production_kg"] = ret["production_mt"] * KG_PER_TONNE From 561fb109bda18c2cd6dc8394520f4aa2f744bd1c Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 16:07:23 -0800 Subject: [PATCH 15/16] Move livestock enums and functions from statistical to gpw_livestock --- .../jurisdictional_direct_proxy_test.py | 3 +- jdluc/__tests__/statistical_test.py | 15 ++- jdluc/attribute.py | 6 +- jdluc/datasets/gpw_livestock.py | 67 ++++++++++++ jdluc/jurisdictional_direct_proxy.py | 18 ++-- jdluc/statistical.py | 100 ++++-------------- jdluc/taxonomy.py | 14 +-- jdluc/trace.py | 6 +- 8 files changed, 119 insertions(+), 110 deletions(-) diff --git a/jdluc/__tests__/jurisdictional_direct_proxy_test.py b/jdluc/__tests__/jurisdictional_direct_proxy_test.py index 7312312..3f96bc3 100644 --- a/jdluc/__tests__/jurisdictional_direct_proxy_test.py +++ b/jdluc/__tests__/jurisdictional_direct_proxy_test.py @@ -6,7 +6,7 @@ from jdluc.datasets.gpw_grassland import DATASET as GRASSLAND from jdluc.datasets.gpw_grassland import YEARS as GRASSLAND_YEARS from jdluc.datasets.gpw_grassland import Grassland -from jdluc.datasets.gpw_livestock import Species, get_band_name +from jdluc.datasets.gpw_livestock import Livestock, Species, get_band_name from jdluc.datasets.tang_cropgrids import DATASET as CROPGRIDS from jdluc.datasets.tang_cropgrids import YEAR as CROPGRIDS_YEAR from jdluc.datasets.tang_cropgrids import Crop @@ -25,7 +25,6 @@ get_crop_to_area_share, get_name_to_band, ) -from jdluc.statistical import Livestock # A conversion in 2013 takes a factor of 0.0625, so 160 and 1600 discount to 10 and 100 CONVERSION_YEAR = 2013 diff --git a/jdluc/__tests__/statistical_test.py b/jdluc/__tests__/statistical_test.py index 90e2640..ffd8f8f 100644 --- a/jdluc/__tests__/statistical_test.py +++ b/jdluc/__tests__/statistical_test.py @@ -8,19 +8,21 @@ import xarray from jdluc import emit -from jdluc.datasets import gpw_livestock, ifpri_mapspam +from jdluc.datasets import faostat, gpw_livestock, ifpri_mapspam +from jdluc.datasets.gpw_livestock import ( + Livestock, + get_band_name_to_livestock_units, + get_livestock_to_grazing_share, +) from jdluc.statistical import ( HERD_HEADS_COLUMNS, MAPSPAM_SPANS, Commodity, Crop, - Livestock, - get_band_name_to_livestock_units, get_commodity_name_to_totals, get_commodity_to_peatland_conversion_share, get_commodity_to_share, get_crop_to_area_share, - get_livestock_to_grazing_share, get_mapspam_span, get_snapshot_year_to_weight, ) @@ -1061,3 +1063,8 @@ def test_the_beef_row_reports_its_herd( assert totals["commodity_hectares"] == pytest.approx( PIXELS * HECTARES_PER_CELL * pasture_fraction ) + + +def test_faostat_species_match_gpw_livestock() -> None: + """Every grazer GPW maps has FAOSTAT stocks and carcass weights to price it by.""" + assert {e.name for e in faostat.Species} == {e.name for e in gpw_livestock.Species} diff --git a/jdluc/attribute.py b/jdluc/attribute.py index 7502176..18cdc29 100644 --- a/jdluc/attribute.py +++ b/jdluc/attribute.py @@ -36,7 +36,7 @@ jurisdictional_direct_proxy, statistical, ) -from jdluc.datasets import tang_cropgrids, worldbank_jurisdictions +from jdluc.datasets import gpw_livestock, tang_cropgrids, worldbank_jurisdictions logger = logging.getLogger(__name__) @@ -152,11 +152,11 @@ def get_commodity_names(methodology: Methodology) -> tuple[str, ...]: return tuple(sorted(e.name for e in jurisdictional_direct.Crop)) case Methodology.JURISDICTIONAL_DIRECT_PROXY: return tuple( - sorted(e.name for e in (*tang_cropgrids.Crop, *statistical.Livestock)) + sorted(e.name for e in (*tang_cropgrids.Crop, *gpw_livestock.Livestock)) ) case Methodology.STATISTICAL: return tuple( - sorted(e.name for e in (*statistical.Crop, *statistical.Livestock)) + sorted(e.name for e in (*statistical.Crop, *gpw_livestock.Livestock)) ) case _: raise ValueError(methodology) diff --git a/jdluc/datasets/gpw_livestock.py b/jdluc/datasets/gpw_livestock.py index 7b3e8f5..b6212e2 100644 --- a/jdluc/datasets/gpw_livestock.py +++ b/jdluc/datasets/gpw_livestock.py @@ -34,6 +34,7 @@ import rasterio.enums import rasterio.transform import rasterio.warp +import xarray from jdluc import tiling, utils from jdluc.datasets import base @@ -149,3 +150,69 @@ def _save_tile_id_to_local_path(local_path: str, tile_id: str) -> None: def get_band_name(species: Species, year: int) -> str: """The fully-qualified band carrying `species`' density, in heads/ha, in `year`.""" return DATASET.fully_qualified_band_names[SPECIES_YEARS.index((species, year))] + + +@enum.unique +class Livestock(enum.StrEnum): + """The commodities dividing pasture's share of a cell; PASTURE keeps what the others don't.""" + + BEEF_CATTLE = enum.auto() + PASTURE = enum.auto() + + +LIVESTOCK_TO_SPECIES = { + Livestock.BEEF_CATTLE: (Species.CATTLE,), +} +assert set(LIVESTOCK_TO_SPECIES) == set(Livestock) - {Livestock.PASTURE} +# No grazer may be named by two commodities, or its livestock units would count toward both shares +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()) + + +def get_band_name_to_livestock_units( + dset: xarray.Dataset, species_to_units_per_head: dict[Species, float] +) -> dict[str, xarray.DataArray]: + """Each grazer's density in livestock units/ha, for assigning beside its density in heads. + + Livestock units weigh how much one head of a species grazes, so they divide pasture between + grazers; production counts heads, and reads the heads bands. + """ + return { + f"{species:s}:livestock-units-per-ha:{year:d}": dset[ + get_band_name(species=species, year=year) + ] + * species_to_units_per_head[species] + for species, year in SPECIES_YEARS + } + + +def get_livestock_to_grazing_share( + dset: xarray.Dataset, year: int +) -> dict[Livestock, xarray.DataArray]: + """Each commodity's share of a cell's pasture: its grazers' livestock units over all of them. + + Each livestock unit in a cell is taken to graze an equal slice of that cell's pasture. Only the + ratios between species' units per head matter, since a share divides each by the cell's total. + """ + species_to_units = { + species: dset[f"{species:s}:livestock-units-per-ha:{year:d}"] + for species in Species + } + zero = xarray.DataArray(numpy.float32(0)) + total = sum(species_to_units.values(), start=zero) + total = total.where(total > 0) + shares = { + # Share is zero where no grazer is mapped at all + livestock: ( + sum((species_to_units[one] for one in species), start=zero) / total + ).fillna(0) + for livestock, species in LIVESTOCK_TO_SPECIES.items() + } + return shares | {Livestock.PASTURE: 1 - sum(shares.values(), start=zero)} + + +def get_herd_heads_column(year: int) -> str: + return f"herd_heads_{year:d}" diff --git a/jdluc/jurisdictional_direct_proxy.py b/jdluc/jurisdictional_direct_proxy.py index c16c7b2..b4c84d3 100644 --- a/jdluc/jurisdictional_direct_proxy.py +++ b/jdluc/jurisdictional_direct_proxy.py @@ -8,7 +8,7 @@ import xarray from rioxarray.exceptions import NoDataInBounds -from jdluc import emit, geo, harmonize, statistical, storage, tiling, utils +from jdluc import emit, geo, harmonize, storage, tiling, utils from jdluc.datasets import ( DatasetName, faostat, @@ -24,13 +24,13 @@ assert tang_cropgrids.YEAR in gpw_grassland.YEARS assert tang_cropgrids.YEAR in gpw_livestock.YEARS -type Commodity = tang_cropgrids.Crop | statistical.Livestock +type Commodity = tang_cropgrids.Crop | gpw_livestock.Livestock assert not {e.name for e in tang_cropgrids.Crop} & { - e.name for e in statistical.Livestock + e.name for e in gpw_livestock.Livestock } DATASET_NAMES = (DatasetName.GPW_LIVESTOCK, DatasetName.TANG_CROPGRIDS_CROP_AREA) -HERD_HEADS_COLUMN = statistical.get_herd_heads_column(year=tang_cropgrids.YEAR) +HERD_HEADS_COLUMN = gpw_livestock.get_herd_heads_column(year=tang_cropgrids.YEAR) def get_name_to_band( @@ -189,7 +189,7 @@ def get_commodity_name_to_totals( for commodity in commodities: totals = commodity_to_totals[commodity] share = commodity_to_share[commodity] - if isinstance(commodity, statistical.Livestock): + if isinstance(commodity, gpw_livestock.Livestock): land_class = emit.PASTURELAND totals["commodity_hectares"] = ( dset[f"pasture:fraction:{tang_cropgrids.YEAR:d}"] * hectares * share @@ -222,7 +222,7 @@ def get_commodity_name_to_totals( gpw_livestock.get_band_name(species=species, year=tang_cropgrids.YEAR) ].where(dset[f"pasture:fraction:{tang_cropgrids.YEAR:d}"] > 0, other=0) * hectares - for livestock, (species,) in statistical.LIVESTOCK_TO_SPECIES.items() + for livestock, (species,) in gpw_livestock.LIVESTOCK_TO_SPECIES.items() if livestock in commodities } for livestock, heads in livestock_to_heads.items(): @@ -294,7 +294,7 @@ def workflow( commodities: tuple[Commodity, ...] = tuple( tang_cropgrids.Crop[name] if name in tang_cropgrids.Crop.__members__ - else statistical.Livestock[name] + else gpw_livestock.Livestock[name] for name in commodity_names ) crops = tuple( @@ -317,7 +317,7 @@ def workflow( ) # Each grazer's livestock units per head, as FAOSTAT counts them for this country merged = merged.assign( - statistical.get_band_name_to_livestock_units( + gpw_livestock.get_band_name_to_livestock_units( dset=merged, species_to_units_per_head={ gpw_livestock.Species[species.name]: units_per_head @@ -359,7 +359,7 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]: commodity_name_to_totals = get_commodity_name_to_totals( commodities=commodities, commodity_to_share=get_crop_to_area_share(crops=crops, dset=clipped) - | statistical.get_livestock_to_grazing_share( + | gpw_livestock.get_livestock_to_grazing_share( dset=clipped, year=tang_cropgrids.YEAR ), dset=clipped, diff --git a/jdluc/statistical.py b/jdluc/statistical.py index 0533d79..5d99e91 100644 --- a/jdluc/statistical.py +++ b/jdluc/statistical.py @@ -62,31 +62,9 @@ class Crop(enum.StrEnum): assert {e.value for e in Crop} == ifpri_mapspam.RECOVERABLE_CROP_NAMES -@enum.unique -class Livestock(enum.StrEnum): - """The commodities dividing pasture's share of a cell; PASTURE keeps what the others don't.""" - - BEEF_CATTLE = enum.auto() - PASTURE = enum.auto() - - -LIVESTOCK_TO_SPECIES = { - Livestock.BEEF_CATTLE: (gpw_livestock.Species.CATTLE,), -} -assert set(LIVESTOCK_TO_SPECIES) == set(Livestock) - {Livestock.PASTURE} -# No grazer may be named by two commodities, or its livestock units would count toward both shares -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} - - # `Crop` stays exactly MapSPAM's recoverable crops, because everything that divides by one looks # its name up there -type Commodity = Crop | Livestock +type Commodity = Crop | gpw_livestock.Livestock DATASET_NAMES = ( @@ -219,57 +197,17 @@ def get_downscaled_luc_emissions(assessment_year: int, tile_id: str) -> xarray.D ) -def get_band_name_to_livestock_units( - dset: xarray.Dataset, species_to_units_per_head: dict[gpw_livestock.Species, float] -) -> dict[str, xarray.DataArray]: - """Each grazer's density in livestock units/ha, for assigning beside its density in heads. - - Livestock units weigh how much one head of a species grazes, so they divide pasture between - grazers; production counts heads, and reads the heads bands. - """ - return { - f"{species:s}:livestock-units-per-ha:{year:d}": dset[ - gpw_livestock.get_band_name(species=species, year=year) - ] - * species_to_units_per_head[species] - for species, year in gpw_livestock.SPECIES_YEARS - } - - -def get_livestock_to_grazing_share( - dset: xarray.Dataset, year: int -) -> dict[Livestock, xarray.DataArray]: - """Each commodity's share of a cell's pasture: its grazers' livestock units over all of them. - - Each livestock unit in a cell is taken to graze an equal slice of that cell's pasture. Only the - ratios between species' units per head matter, since a share divides each by the cell's total. - """ - species_to_units = { - species: dset[f"{species:s}:livestock-units-per-ha:{year:d}"] - for species in gpw_livestock.Species - } - zero = xarray.DataArray(numpy.float32(0)) - total = sum(species_to_units.values(), start=zero) - total = total.where(total > 0) - shares = { - # Share is zero where no grazer is mapped at all - livestock: ( - sum((species_to_units[one] for one in species), start=zero) / total - ).fillna(0) - for livestock, species in LIVESTOCK_TO_SPECIES.items() - } - return shares | {Livestock.PASTURE: 1 - sum(shares.values(), start=zero)} - - def get_commodity_hectares( commodity: Commodity, dset: xarray.Dataset, year: int ) -> xarray.DataArray: - if isinstance(commodity, Livestock): + if isinstance(commodity, gpw_livestock.Livestock): pasture_fraction = dset[f"pasture:fraction:{year:d}"] return ( pasture_fraction * emit.get_hectares_per_pixel(darray=pasture_fraction) - * get_livestock_to_grazing_share(dset=dset, year=year)[commodity] + * gpw_livestock.get_livestock_to_grazing_share(dset=dset, year=year)[ + commodity + ] ) else: return ifpri_mapspam.get_canonical_quantity( @@ -359,7 +297,7 @@ def get_commodity_to_share( # livestock commodity's own expansion instead would count one grazer replacing another as # new pasture. livestock: pasture_share * grazing_share - for livestock, grazing_share in get_livestock_to_grazing_share( + for livestock, grazing_share in gpw_livestock.get_livestock_to_grazing_share( dset=dset, year=after ).items() } @@ -383,7 +321,7 @@ def get_commodity_to_peatland_conversion_share( / cropland_expansion ).fillna(0) for crop in crops - } | get_livestock_to_grazing_share(dset=dset, year=after) + } | gpw_livestock.get_livestock_to_grazing_share(dset=dset, year=after) MAPSPAM_SPANS = ((2000, 2005), (2005, 2010), (2010, 2020)) @@ -436,12 +374,8 @@ def get_snapshot_year_to_weight(assessment_year: int) -> dict[int, float]: return snapshot_year_to_weight -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 + gpw_livestock.get_herd_heads_column(year=year) for year in MAPSPAM_SNAPSHOT_YEARS ) @@ -484,7 +418,9 @@ def get_commodity_name_to_totals( for commodity in commodities: totals = commodity_to_totals[commodity] land_class = ( - emit.PASTURELAND if isinstance(commodity, Livestock) else emit.CROPLAND + emit.PASTURELAND + if isinstance(commodity, gpw_livestock.Livestock) + else emit.CROPLAND ) component_to_emissions = { component: sum( @@ -520,9 +456,9 @@ 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; a crop takes its area share of the cropland peat # band and reports MapSPAM's production - if isinstance(commodity, Livestock): + if isinstance(commodity, gpw_livestock.Livestock): occupation = pastureland_occupation * occupation_shares[commodity] - if commodity in LIVESTOCK_TO_SPECIES: + if commodity in gpw_livestock.LIVESTOCK_TO_SPECIES: # NB: GPW's density is heads per hectare of the whole cell, so times the cell's # hectares it is the cell's entire herd, all of it grazing the cell's pasture # however little of the cell that is. Scaling by the pasture fraction too would @@ -530,9 +466,9 @@ 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. - (species,) = LIVESTOCK_TO_SPECIES[commodity] + (species,) = gpw_livestock.LIVESTOCK_TO_SPECIES[commodity] totals |= { - get_herd_heads_column(year=year): dset[ + gpw_livestock.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 @@ -691,7 +627,7 @@ def workflow( from rioxarray.exceptions import NoDataInBounds commodities: tuple[Commodity, ...] = tuple( - Crop[name] if name in Crop.__members__ else Livestock[name] + Crop[name] if name in Crop.__members__ else gpw_livestock.Livestock[name] for name in commodity_names ) crops = tuple(commodity for commodity in commodities if isinstance(commodity, Crop)) @@ -716,7 +652,7 @@ def workflow( # Each grazer's livestock units per head, as FAOSTAT counts them for this country merged = merged.assign( - get_band_name_to_livestock_units( + gpw_livestock.get_band_name_to_livestock_units( dset=merged, species_to_units_per_head={ gpw_livestock.Species[species.name]: units_per_head @@ -784,7 +720,7 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]: dset=clipped, year=max(ifpri_mapspam.YEARS), ) - | get_livestock_to_grazing_share( + | gpw_livestock.get_livestock_to_grazing_share( dset=clipped, year=max(ifpri_mapspam.YEARS) ), snapshot_year_to_weight=snapshot_year_to_weight, diff --git a/jdluc/taxonomy.py b/jdluc/taxonomy.py index b79911f..ddef07b 100644 --- a/jdluc/taxonomy.py +++ b/jdluc/taxonomy.py @@ -18,7 +18,7 @@ import enum from jdluc import attribute, emit, jurisdictional_direct, statistical -from jdluc.datasets import faostat, tang_cropgrids +from jdluc.datasets import faostat, gpw_livestock, tang_cropgrids @enum.unique @@ -137,18 +137,18 @@ class InventedCommodityId(enum.IntEnum): statistical.Crop.SWEET_POTATO.name: 122, statistical.Crop.WHEAT.name: 15, statistical.Crop.YAM.name: 137, - statistical.Livestock.BEEF_CATTLE.name: faostat.SPECIES_TO_MEAT_ITEM_CODE[ + gpw_livestock.Livestock.BEEF_CATTLE.name: faostat.SPECIES_TO_MEAT_ITEM_CODE[ faostat.Species.CATTLE ], - statistical.Livestock.PASTURE.name: InventedCommodityId.PASTURE, + gpw_livestock.Livestock.PASTURE.name: InventedCommodityId.PASTURE, }, attribute.Methodology.JURISDICTIONAL_DIRECT_PROXY: { emit.NonCommodity.DROPPED.name: InventedCommodityId.DROPPED, emit.NonCommodity.UNATTRIBUTED.name: InventedCommodityId.UNATTRIBUTED, - statistical.Livestock.BEEF_CATTLE.name: faostat.SPECIES_TO_MEAT_ITEM_CODE[ + gpw_livestock.Livestock.BEEF_CATTLE.name: faostat.SPECIES_TO_MEAT_ITEM_CODE[ faostat.Species.CATTLE ], - statistical.Livestock.PASTURE.name: InventedCommodityId.PASTURE, + gpw_livestock.Livestock.PASTURE.name: InventedCommodityId.PASTURE, # The publisher's crop table (`CROPS_matching_v17.xlsx` in figshare's `CODES.zip`), through # the authors' 2020 FAOSTAT extract. Both kapok crops are 310, kapok fruit tang_cropgrids.Crop.ABACA.name: 809, @@ -352,7 +352,7 @@ class InventedCommodityId(enum.IntEnum): METHODOLOGY_TO_COMMODITY_NAME_TO_COMMODITY_ID[attribute.Methodology.STATISTICAL] ) == { *(crop.name for crop in statistical.Crop), - *(livestock.name for livestock in statistical.Livestock), + *(livestock.name for livestock in gpw_livestock.Livestock), *(non_commodity.name for non_commodity in emit.NonCommodity), } assert set( @@ -361,7 +361,7 @@ class InventedCommodityId(enum.IntEnum): ] ) == { *(crop.name for crop in tang_cropgrids.Crop), - *(livestock.name for livestock in statistical.Livestock), + *(livestock.name for livestock in gpw_livestock.Livestock), *(non_commodity.name for non_commodity in emit.NonCommodity), } assert len( diff --git a/jdluc/trace.py b/jdluc/trace.py index 5378bfb..c3c1bea 100644 --- a/jdluc/trace.py +++ b/jdluc/trace.py @@ -135,12 +135,12 @@ def derive_statistical_production_kg( """ ret = emissions.copy() ret["production_kg"] = ret["production_mt"] * KG_PER_TONNE - for livestock, (species,) in statistical.LIVESTOCK_TO_SPECIES.items(): + for livestock, (species,) in gpw_livestock.LIVESTOCK_TO_SPECIES.items(): is_herd = ret.index.get_level_values("commodity_name") == livestock.name ret.loc[is_herd, "production_kg"] = sum( weight * get_herd_production_kg( - heads=ret.loc[is_herd, statistical.get_herd_heads_column(year=year)], + heads=ret.loc[is_herd, gpw_livestock.get_herd_heads_column(year=year)], species=species, stocks=stocks, years=(year,), @@ -233,7 +233,7 @@ def derive_proxy_production_kg( ret = flat.drop(columns=["commodity_id", "iso_3166"]).set_index( emissions.index.names ) - for livestock, (species,) in statistical.LIVESTOCK_TO_SPECIES.items(): + for livestock, (species,) in gpw_livestock.LIVESTOCK_TO_SPECIES.items(): is_herd = ret.index.get_level_values("commodity_name") == livestock.name ret.loc[is_herd, "production_kg"] = get_herd_production_kg( heads=ret.loc[is_herd, jurisdictional_direct_proxy.HERD_HEADS_COLUMN], From d223e9c03c722fccc0ff1481dca17fbceab24447 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Tue, 10 Nov 2026 11:41:18 -0800 Subject: [PATCH 16/16] Revert JD production year-window to be 4 years ending at the assessment year Have JD proxy adopt this same window -- allowing us to remove the awkward coupling to FAOSTAT freshness another docs update --- docs/executive_summary.md | 2 +- docs/further_research.md | 2 +- docs/methodology.md | 6 +++--- jdluc/trace.py | 23 +++++++++-------------- 4 files changed, 14 insertions(+), 19 deletions(-) diff --git a/docs/executive_summary.md b/docs/executive_summary.md index 81f0aa2..d59f011 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_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`, `STATISTICAL` or `JURISDICTIONAL_DIRECT_PROXY`. Each row carries the headline emissions factor and the quantities behind it; all three legs 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_id, admin_level, assessment_year, commodity_id, methodology)`, where `admin_id` is the ISO 3166 alpha-3 code nationally and the World Bank admin-1 code provincially, `admin_level` is `NATIONAL` or `PROVINCIAL`, `assessment_year` is the year the factors describe (2020, 2022 or 2024), `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`, `STATISTICAL` or `JURISDICTIONAL_DIRECT_PROXY`. Each row carries the headline emissions factor and the quantities behind it; all three legs 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/further_research.md b/docs/further_research.md index 9268a21..9de5455 100644 --- a/docs/further_research.md +++ b/docs/further_research.md @@ -231,7 +231,7 @@ The 2010 snapshot ends the 2005–2010 span and starts the expansion both later **Issue**: The jurisdictional-direct leg divides emissions by `commodity_hectares x yield_kg_per_ha`, where the yield is a USDA NASS survey statistic and the hectares come from the CDL. Every one of the eleven crops `jurisdictional_direct.Crop` models carries a yield, so every one produces an emissions factor. Three properties of how that denominator is assembled are worth recording, because each is an assumption rather than a measurement. 1. **Unmatched (admin, crop) pairs are still only a log line.** NASS suppresses small and identifying cells under Title 7 U.S.C. and CIPSEA, and does not survey every crop in every state. Where the CDL sees a crop in a state that NASS does not publish a yield for, that state contributes hectares and emissions to the national rollup but no production, so the national factor is a ratio over the reporting states alone. The count of such rows is logged and then discarded. -2. **Cotton is converted between product forms by a single constant.** NASS reports cotton yield as ginned lint, while MapSPAM, FAOSTAT (item 328, not lint's 767) and WRI all carry seed cotton, and either mistake scales a comparison by the milling yield. `usda_nass_quickstats.Series.reported_fraction` divides the lint yield back up to the unginned crop at 0.36. That figure is what FAOSTAT's US seed-cotton yield implies against the NASS series over 2011–2020, and it agrees with the conventional ~35% gin turnout, but the same comparison implies ~0.403 over 2000–2010. The constant is therefore tied to the years the direct leg averages yields over, the 3 centered on the assessment year, and a national constant hides whatever provincial variation in gin turnout exists. +2. **Cotton is converted between product forms by a single constant.** NASS reports cotton yield as ginned lint, while MapSPAM, FAOSTAT (item 328, not lint's 767) and WRI all carry seed cotton, and either mistake scales a comparison by the milling yield. `usda_nass_quickstats.Series.reported_fraction` divides the lint yield back up to the unginned crop at 0.36. That figure is what FAOSTAT's US seed-cotton yield implies against the NASS series over 2011–2020, and it agrees with the conventional ~35% gin turnout, but the same comparison implies ~0.403 over 2000–2010. The constant is therefore tied to the years the direct leg averages yields over, the 4 ending at the assessment year, and a national constant hides whatever provincial variation in gin turnout exists. 3. **Two series are not quite the crop the CDL scores.** NASS renamed its rolled-up dry-bean class in 2019 and publishes no `ALL CLASSES` to fall back on, so spanning a multi-year window takes both names; the pre-2019 aggregate includes chickpeas, which the CDL scores separately as `CHICK_PEAS`. Sugarcane publishes no `ALL UTILIZATION PRACTICES` row either, so its total is `SUGAR & SEED` and includes seed cane. **Potential impact:** Modest overall, and concentrated in one crop. Cotton is the pair whose assumption moves a published number materially: carried as lint rather than seed cotton its factor would read roughly 2.8x higher. The eight crops outside maize, soybean and wheat carry 9% of US jurisdictional-direct emissions between them, so an error in their denominators is bounded by that share. The residual coverage gap is concentrated in minor producing states for the smaller crops, so it biases a national factor by whatever those states' yields differ from the surveyed ones — likely small for maize and soybean, larger for crops NASS surveys in only a handful of states. The chickpea contamination affects half the current averaging window for one crop carrying 0.627 Mt. diff --git a/docs/methodology.md b/docs/methodology.md index 11caa3e..2a776f3 100644 --- a/docs/methodology.md +++ b/docs/methodology.md @@ -127,7 +127,7 @@ Ingest and harmonize keep every year of each layer, so every assessment year rea | GPW livestock | 2020 | Statistical and proxy legs: grazing shares and herds | | CROPGRIDS | 2020 | Proxy leg: crop-area shares, with the pasture fraction read in the same year | -The production windows follow the year instead: the direct leg's NASS yields and the proxy leg's FAOSTAT production are each averaged over the 3 years centered on the assessment year. FAOSTAT has not published 2025's production, so the proxy leg's 2024 factors take 2024's alone until it does. +The production windows follow the year instead: the direct leg's NASS yields and the proxy leg's FAOSTAT production are both averaged over the 4 years ending at the assessment year (2017–2020 for 2020). ### Three grids @@ -241,7 +241,7 @@ Concretely, for each (jurisdiction, crop): A crop's row covers only the pixels the crop map assigns to that crop, and it never includes emissions from peat drained under pasture, even where the crop map puts the crop on the same pixel. Everything else goes to `UNATTRIBUTED`: all emissions on pixels the crop map does not assign to a crop we model, and all emissions from peat drained under pasture. -Because direct attribution relies only on **current** production, production is later computed as crop area × NASS QuickStats yield (the state-level mean over the 3 years centered on the assessment year), converting NASS's reported yield per acre to kg/ha. NASS publishes that yield in four units — bushels, pounds, hundredweight and short tons per acre — and which one applies is a property of the commodity, so `usda_nass_quickstats.CropSeries` carries the unit and its weight per crop. The bushel weights are marketing bushels, from table 6 of USDA Agricultural Handbook 697 (ERS, June 1992), rather than the grading test weights of 7 CFR 810, which are a different quantity. Cotton is the one crop whose published yield is not the harvested crop: NASS reports ginned lint where MapSPAM and FAOSTAT carry seed cotton, so the lint yield is divided by a lint fraction of 0.36. This leg currently asserts `iso_3166 == "USA"`, since CDL is US-only. +Because direct attribution relies only on **current** production, production is later computed as crop area × NASS QuickStats yield (the state-level mean over the 4 years ending at the assessment year), converting NASS's reported yield per acre to kg/ha. NASS publishes that yield in four units — bushels, pounds, hundredweight and short tons per acre — and which one applies is a property of the commodity, so `usda_nass_quickstats.CropSeries` carries the unit and its weight per crop. The bushel weights are marketing bushels, from table 6 of USDA Agricultural Handbook 697 (ERS, June 1992), rather than the grading test weights of 7 CFR 810, which are a different quantity. Cotton is the one crop whose published yield is not the harvested crop: NASS reports ginned lint where MapSPAM and FAOSTAT carry seed cotton, so the lint yield is divided by a lint fraction of 0.36. This leg currently asserts `iso_3166 == "USA"`, since CDL is US-only. ### Statistical (`statistical.py`) @@ -282,7 +282,7 @@ Where no reference year places the group in a pixel there is no basis for a spli Each pixel's emissions are split by the land class it became, cropland or pasture, and averaged onto the CROPGRIDS grid. In each cell, each crop takes its share of the 2020 crop area of all 173 CROPGRIDS crops of what cropland is charged, and each livestock commodity takes its 2020 grazing share, as in the statistical leg, of what pasture is charged. CROPGRIDS has one year, so every span is divided by the same shares. A crop's hectares are its CROPGRIDS crop area, and `UNATTRIBUTED` is cropland's charge in cells with no CROPGRIDS crop. -Production is joined in `trace.py`, from FAOSTAT's means over the 3 years centered on the assessment year (2019–2021 for 2020). Until FAOSTAT publishes 2025, the 2024 assessment takes 2024's production alone, and `trace.py` stops rather than average fewer years when a year of a window is missing. A crop's is its FAOSTAT item's national production, split among the country's provinces by their share of the country's CROPGRIDS crop area, over every crop that shares the item. `BEEF_CATTLE`'s is the cattle GPW maps in 2020 in cells holding any cultivated grassland, times FAOSTAT's national carcass weight per standing head, as in the statistical leg. FAOSTAT production where CROPGRIDS maps none of the item in that country has no province to go to; `trace.py` logs it and leaves it out. A crop's hectares are physical area while FAOSTAT counts every harvest, so `yield_kg_per_ha` reads above FAOSTAT's yield per harvested hectare for crops grown more than once a year, rice above all. +Production is joined in `trace.py`, from FAOSTAT's means over the 4 years ending at the assessment year (2017–2020 for 2020), the same years the direct leg averages NASS yields over. `trace.py` stops rather than average fewer years when a year of the window is missing. A crop's is its FAOSTAT item's national production, split among the country's provinces by their share of the country's CROPGRIDS crop area, over every crop that shares the item. `BEEF_CATTLE`'s is the cattle GPW maps in 2020 in cells holding any cultivated grassland, times FAOSTAT's national carcass weight per standing head, as in the statistical leg. FAOSTAT production where CROPGRIDS maps none of the item in that country has no province to go to; `trace.py` logs it and leaves it out. A crop's hectares are physical area while FAOSTAT counts every harvest, so `yield_kg_per_ha` reads above FAOSTAT's yield per harvested hectare for crops grown more than once a year, rice above all. ### From rollups to emissions factors (`trace.py`) diff --git a/jdluc/trace.py b/jdluc/trace.py index c3c1bea..89a349d 100644 --- a/jdluc/trace.py +++ b/jdluc/trace.py @@ -2,14 +2,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, averaged over the 3 - years centered on the assessment year. + - JURISDICTIONAL_DIRECT: production = crop area x NASS QuickStats yield, averaged over the 4 + years ending at the assessment year. - 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. - - JURISDICTIONAL_DIRECT_PROXY: production = FAOSTAT national production over the 3 years - centered on the assessment year (2024 alone until FAOSTAT publishes 2025), split among - provinces by CROPGRIDS crop area for crops, and the attribute's 2020 GPW heads x FAOSTAT - carcass weight per head over the same years for livestock. + - JURISDICTIONAL_DIRECT_PROXY: production = FAOSTAT national production over the same 4 years, + split among provinces by CROPGRIDS crop area for crops, and the attribute's 2020 GPW heads x + FAOSTAT carcass weight per head over the same years for livestock. It then derives, identically for all three, each commodity's yield (kg per hectare) and emissions factor (kgCO2e per kg), and rolls the provincials up to national totals. Returns a cached @@ -344,25 +343,21 @@ def workflow( methodology=methodology, ) assert emissions.index.names == list(attribute.INDEX_NAMES) + # The years both direct legs average their production statistics over + production_years = tuple(range(assessment_year - 3, assessment_year + 1)) match methodology: case attribute.Methodology.JURISDICTIONAL_DIRECT: production = derive_jurisdictional_production_kg( emissions=emissions, raw_yields=usda_nass_quickstats.load(), - years=tuple(range(assessment_year - 1, assessment_year + 2)), + years=production_years, ) case attribute.Methodology.JURISDICTIONAL_DIRECT_PROXY: production = derive_proxy_production_kg( crops=faostat.load(dataset=faostat.CROP_DATASET), emissions=emissions, stocks=faostat.load(dataset=faostat.LIVESTOCK_DATASET), - # TODO: center 2024 like every other year once FAOSTAT publishes 2025's - # production; until then it takes its own year alone - years=( - (assessment_year,) - if assessment_year == 2024 - else tuple(range(assessment_year - 1, assessment_year + 2)) - ), + years=production_years, ) case attribute.Methodology.STATISTICAL: production = derive_statistical_production_kg(