From 85da16071d7c7d845bcdc1fb94ef40dbf064b0cb Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 10:15:30 -0800 Subject: [PATCH 1/7] Have emit provide undiscounted emissions and their year Have downstreams perform the discounting --- jdluc/__tests__/emit_test.py | 66 ++++++---- .../jurisdictional_direct_proxy_test.py | 23 ++-- jdluc/__tests__/jurisdictional_direct_test.py | 21 ++-- jdluc/emit.py | 117 ++++++++---------- jdluc/jurisdictional_direct.py | 33 ++--- jdluc/jurisdictional_direct_proxy.py | 10 +- jdluc/statistical.py | 9 +- tools/build-soy-brazil-figure.py | 24 ++-- 8 files changed, 166 insertions(+), 137 deletions(-) diff --git a/jdluc/__tests__/emit_test.py b/jdluc/__tests__/emit_test.py index ac73a48..d177b82 100644 --- a/jdluc/__tests__/emit_test.py +++ b/jdluc/__tests__/emit_test.py @@ -44,6 +44,7 @@ get_grassland_carbon, get_hectares_per_pixel, get_last_departure_year, + get_linear_discounted, get_linear_discounted_total, get_mineral_soil_emissions, get_peatland_occupation_emissions, @@ -456,20 +457,18 @@ def test_get_span_to_charge(conversion_year: int, spans: set[SpanType]) -> None: def test_get_dset_for_output_appends_each_band_its_units() -> None: # Every emitted band is keyed on this name, so a change here strands the zarr caches keyed on - # the old one and KeyErrors the two legs that read it. The units come off the DataArray and - # slot in after the first colon, which is what puts the span last. + # the old one and KeyErrors the legs that read it. The units come off the DataArray and slot in + # after the first colon. result = get_dset_for_output( name_to_darray={ "conversion": get_darray_for_data(data=[[1]]).rename(None), - "soil-emissions:2015-2020": get_darray_for_data(data=[[1]]).rename( - "tco2e-per-ha" - ), + "soil-emissions": get_darray_for_data(data=[[1]]).rename("tco2e-per-ha"), "hectares-per-pixel": get_darray_for_data(data=[[1]]).rename("ha"), } ) assert set(result) == { "conversion", - "soil-emissions:tco2e-per-ha:2015-2020", + "soil-emissions:tco2e-per-ha", "hectares-per-pixel:ha", } @@ -486,6 +485,17 @@ def test_get_linear_discounted_total() -> None: assert numpy.array_equal(result.data, [[0.625]]) +def test_get_linear_discounted() -> None: + # A pixel with no conversion, or no year, takes no weight + result = get_linear_discounted( + conversion_year=get_darray_for_data( + data=[[2001, 2010, 2011, 2020, 0, numpy.nan]] + ), + darray=get_darray_for_data(data=[[8] * 6]), + ) + numpy.testing.assert_allclose(result.data, [[0.1, 0.3, 0.5, 0.7, 0, 0]]) + + def test_get_hectares_per_pixel() -> None: result = get_hectares_per_pixel(darray=get_darray_for_data(data=[[0] * 2] * 2)) assert result.name == "ha" @@ -524,19 +534,16 @@ def test_get_linear_discounted_total_over_the_weight_total_is_a_mean() -> None: def get_dset_for_pixel(conversion: Conversion, is_peat: float) -> xarray.Dataset: - """One pixel, charged the same emissions in every span.""" + """One pixel, converted in 2012.""" return xarray.Dataset( { "conversion": get_darray_for_data(data=[[conversion]]), + "conversion-year": get_darray_for_data(data=[[2012]]), PEATLANDS.fully_qualified_band_name: get_darray_for_data(data=[[is_peat]]), - } - | { - name: get_darray_for_data(data=[[value]]) - for before, after in SPAN_TO_LINEAR_DISCOUNT_WEIGHT - for name, value in ( - (f"vegetation-emissions:tco2e-per-ha:{before:d}-{after:d}", VEGETATION), - (f"soil-emissions:tco2e-per-ha:{before:d}-{after:d}", SOIL), - ) + "soil-emissions:tco2e-per-ha": get_darray_for_data(data=[[SOIL]]), + "vegetation-emissions:tco2e-per-ha": get_darray_for_data( + data=[[VEGETATION]] + ), } ) @@ -581,22 +588,39 @@ def get_dset_for_pixel(conversion: Conversion, is_peat: float) -> xarray.Dataset def test_get_span_to_component_to_emissions( conversion: Conversion, is_peat: float, source_to_expected: dict[str, float] ) -> None: - # `emissions_mt` is computed from the span totals rather than by summing the component columns, - # so the columns decompose it only if the sources claim every unit exactly once. Each case - # names the sources that fire and must account for the whole pixel, which makes a pool - # dropped, double-counted, or gated on the wrong conversion fail here. + # The legs sum the components into `emissions_mt`, so a pool dropped, double-counted or gated + # on the wrong conversion fails here. assert sum(source_to_expected.values()) == VEGETATION + SOIL span_to_component_to_emissions = get_span_to_component_to_emissions( dset=get_dset_for_pixel(conversion=conversion, is_peat=is_peat) ) assert set(span_to_component_to_emissions) == set(SPAN_TO_LINEAR_DISCOUNT_WEIGHT) - for component_to_emissions in span_to_component_to_emissions.values(): + for span, component_to_emissions in span_to_component_to_emissions.items(): assert set(component_to_emissions) == set(EmissionComponent) assert { component: float(darray.sum()) for component, darray in component_to_emissions.items() if float(darray.sum()) - } == source_to_expected + } == (source_to_expected if span == (2010, 2015) else {}) + + +def test_get_span_to_component_to_emissions_charges_each_pixel_to_its_year() -> None: + dset = xarray.Dataset( + { + "conversion": get_darray_for_data( + data=[[Conversion.FOREST_TO_CROPLAND] * 3] + ), + "conversion-year": get_darray_for_data(data=[[2003, 2012, 2020]]), + PEATLANDS.fully_qualified_band_name: get_darray_for_data(data=[[0] * 3]), + "soil-emissions:tco2e-per-ha": get_darray_for_data(data=[[10, 20, 30]]), + "vegetation-emissions:tco2e-per-ha": get_darray_for_data(data=[[1, 2, 3]]), + } + ) + result = get_span_to_component_to_emissions(dset=dset) + assert { + span: float(component_to_emissions[EmissionComponent.FOREST].sum()) + for span, component_to_emissions in result.items() + } == {(2000, 2005): 11, (2005, 2010): 0, (2010, 2015): 22, (2015, 2020): 33} @pytest.mark.parametrize( diff --git a/jdluc/__tests__/jurisdictional_direct_proxy_test.py b/jdluc/__tests__/jurisdictional_direct_proxy_test.py index f3bfe51..827af07 100644 --- a/jdluc/__tests__/jurisdictional_direct_proxy_test.py +++ b/jdluc/__tests__/jurisdictional_direct_proxy_test.py @@ -14,7 +14,6 @@ CROPLAND, PASTURELAND, PEATLAND_EMISSIONS_ANNUAL_TCO2E_PER_HA, - SPAN_TO_LINEAR_DISCOUNT_WEIGHT, Conversion, EmissionComponent, NonCommodity, @@ -27,9 +26,10 @@ ) from jdluc.statistical import Livestock -# The discount weights add up to 0.2, so 50 and 500 per span discount to 10 and 100 -VEGETATION_PER_SPAN = 50.0 -SOIL_PER_SPAN = 500.0 +# A conversion in 2012 takes the 2010-2015 weight of 0.0625, so 160 and 1600 discount to 10 and 100 +CONVERSION_YEAR = 2012 +VEGETATION = 160.0 +SOIL = 1600.0 LAND_CLASS_COMPONENT_NAMES = [ f"{land_class:s}-{component!s}-emissions:tco2e-per-ha" for land_class in (CROPLAND, PASTURELAND) @@ -48,9 +48,10 @@ def get_darray_for_cells(data: list[list[float]]) -> xarray.DataArray: def get_dset_for_pixel( conversion: Conversion, grassland: float, is_peat: bool ) -> xarray.Dataset: - """One pixel, charged the same emissions in every span.""" + """One pixel, converted in `CONVERSION_YEAR`.""" values = { "conversion": float(conversion), + "conversion-year": float(CONVERSION_YEAR), f"{CROPLAND:s}-peatland-occupation:tco2e-per-ha": 0.0, "dropped-emissions:tco2e-per-ha": 0.0, GRASSLAND.fully_qualified_band_names[ @@ -58,16 +59,8 @@ def get_dset_for_pixel( ]: grassland, f"{PASTURELAND:s}-peatland-occupation:tco2e-per-ha": 0.0, PEATLANDS.fully_qualified_band_name: float(is_peat), - } | { - name: value - for before, after in SPAN_TO_LINEAR_DISCOUNT_WEIGHT - for name, value in ( - (f"soil-emissions:tco2e-per-ha:{before:d}-{after:d}", SOIL_PER_SPAN), - ( - f"vegetation-emissions:tco2e-per-ha:{before:d}-{after:d}", - VEGETATION_PER_SPAN, - ), - ) + "soil-emissions:tco2e-per-ha": SOIL, + "vegetation-emissions:tco2e-per-ha": VEGETATION, } return xarray.Dataset( {name: get_darray_for_cells(data=[[value]]) for name, value in values.items()} diff --git a/jdluc/__tests__/jurisdictional_direct_test.py b/jdluc/__tests__/jurisdictional_direct_test.py index a7cc1f9..3d4212b 100644 --- a/jdluc/__tests__/jurisdictional_direct_test.py +++ b/jdluc/__tests__/jurisdictional_direct_test.py @@ -34,16 +34,16 @@ def test_get_commodity_name_to_totals() -> None: }, # crop pixels: (0,0), (0,1), (1,1); of those, peatland at (0,0) and (1,1). crop_class=get_darray_for_data([[1, 1, 0], [0, 1, 0]]), + cropland_occupation_per_hectare=get_darray_for_data( + data=[[5, 5, 5], [5, 5, 5]] + ), crops=(Crop.MAIZE,), dropped_per_hectare=get_darray_for_data([[1, 0, 1], [0, 1, 0]]), - # 99s sit on non-crop pixels — they must be excluded by the mask. - emissions_per_hectare=get_darray_for_data([[1, 2, 99], [99, 4, 99]]), hectares_per_pixel=get_darray_for_data([[10, 10, 10], [10, 10, 10]]), is_peatland=get_darray_for_data([[True, False, False], [False, True, False]]), pastureland_occupation_per_hectare=get_darray_for_data( data=[[0, 0, 0], [0, 0, 0]] ), - peatland_occupation_per_hectare=get_darray_for_data([[5, 5, 5], [5, 5, 5]]), ) assert set(result) == { Crop.MAIZE.name, @@ -56,7 +56,7 @@ def test_get_commodity_name_to_totals() -> None: } # (1+1+1) x 10 assert result[Crop.MAIZE.name] == { "commodity_hectares": 30, # 3 crop pixels x 10 ha - "emissions_mt": 70, # (1+2+4) x 10; non-crop 99s excluded + "emissions_mt": 340, # (9+11+14) x 10; components plus occupation "forest_emissions_mt": 100, # (2+3+5) x 10; non-crop 99s excluded "grassland_emissions_mt": 30, # (1+1+1) x 10; non-crop 99s excluded "peatland_commodity_hectares": 20, # only (0,0) and (1,1) are crop AND peatland @@ -71,13 +71,12 @@ def test_get_commodity_name_to_totals_with_no_crop_pixels() -> None: emit.EmissionComponent, get_darray_for_data([[1, 2], [3, 4]]) ), crop_class=get_darray_for_data([[0, 0], [0, 0]]), + cropland_occupation_per_hectare=get_darray_for_data(data=[[5, 5], [5, 5]]), crops=(Crop.SOYBEAN, Crop.WHEAT), dropped_per_hectare=get_darray_for_data([[0, 0, 0], [0, 0, 0]]), - emissions_per_hectare=get_darray_for_data([[1, 2], [3, 4]]), hectares_per_pixel=get_darray_for_data([[10, 10], [10, 10]]), is_peatland=get_darray_for_data([[True, True], [True, True]]), pastureland_occupation_per_hectare=get_darray_for_data(data=[[0, 0], [0, 0]]), - peatland_occupation_per_hectare=get_darray_for_data([[5, 5], [5, 5]]), ) assert result[emit.NonCommodity.DROPPED.name] == {"emissions_mt": 0} zeros = { @@ -98,8 +97,9 @@ def test_get_commodity_name_to_totals_with_no_crop_pixels() -> None: def test_the_crop_rows_and_the_unattributed_row_sum_to_the_charged_total() -> None: # Pixel by pixel: maize, a class no modeled crop names, and maize again on land that became # pasture over drained peat. The middle pixel goes to UNATTRIBUTED whole, and from the last - # only the pasture peat does, since no crop row takes peat drained under pasture. - emissions_per_hectare = get_darray_for_data(data=[[10, 20, 30]]) + # only the pasture peat does, since no crop row takes peat drained under pasture. Each pixel's + # charge is its components plus both occupation bands. + charged_per_hectare = get_darray_for_data(data=[[10, 20, 30]]) result = get_commodity_name_to_totals( component_to_per_hectare={ emit.EmissionComponent.FOREST: get_darray_for_data(data=[[5, 10, 12]]), @@ -109,13 +109,12 @@ def test_the_crop_rows_and_the_unattributed_row_sum_to_the_charged_total() -> No ), }, crop_class=get_darray_for_data(data=[[1, 0, 1]]), + cropland_occupation_per_hectare=get_darray_for_data(data=[[2, 0, 0]]), crops=(Crop.MAIZE,), dropped_per_hectare=get_darray_for_data(data=[[7, 7, 7]]), - emissions_per_hectare=emissions_per_hectare, hectares_per_pixel=get_darray_for_data(data=[[1, 1, 1]]), is_peatland=get_darray_for_data(data=[[True, True, True]]), pastureland_occupation_per_hectare=get_darray_for_data(data=[[0, 4, 6]]), - peatland_occupation_per_hectare=get_darray_for_data(data=[[2, 0, 0]]), ) assert result[Crop.MAIZE.name]["emissions_mt"] == 34 # 10 + (30 - 6) assert result[emit.NonCommodity.UNATTRIBUTED.name] == { @@ -127,7 +126,7 @@ def test_the_crop_rows_and_the_unattributed_row_sum_to_the_charged_total() -> No } assert result[Crop.MAIZE.name]["emissions_mt"] + result[ emit.NonCommodity.UNATTRIBUTED.name - ]["emissions_mt"] == float(emissions_per_hectare.sum()) + ]["emissions_mt"] == float(charged_per_hectare.sum()) # DROPPED was never charged, so it is no part of that total assert result[emit.NonCommodity.DROPPED.name] == {"emissions_mt": 21} diff --git a/jdluc/emit.py b/jdluc/emit.py index 6a0c112..f9f9b11 100644 --- a/jdluc/emit.py +++ b/jdluc/emit.py @@ -1,13 +1,13 @@ """Compute per-pixel land-conversion emissions from the harmonized dataset. Resolves which of five conversions each pixel underwent and when, releases the carbon pools that -conversion carries, charges them to the five-year span holding that year, adds ongoing -peatland-occupation emissions split by destination, and applies the GHGP 20-year linear discount. +conversion carries, and adds ongoing peatland-occupation emissions split by destination. The +emissions are undiscounted; each allocation leg discounts them by `conversion-year`. Returns a cached xarray.Dataset: the conversion, the year its source class ended and the datasets -that claimed its destination; vegetation and soil emissions per span; the two occupation bands; -the discounted per-hectare total; the source carbon no destination claimed; and a -hectares-per-pixel band for downstream area-scaling. +that claimed its destination; the vegetation and soil emissions of that conversion; the two +occupation bands; the source carbon no destination claimed; and a hectares-per-pixel band for +downstream area-scaling. Source carbon that no destination claimed is charged to nobody and reported on its own as `dropped-emissions`. @@ -221,7 +221,7 @@ def get_mineral_soil_emissions( output_dtypes=[numpy.float32], ) # NB: SoilGrids has genuine gaps (water, rock, ice) that harmonize turns into NaN. - # Treat a missing stock as zero rather than letting NaN reach emissions-per-hectare, + # Treat a missing stock as zero rather than letting NaN reach a leg's emissions total, # where it would silently discard the pixel's vegetation emissions too. * soil_organic_carbon.fillna(0) ).rename("tco2e-per-ha") @@ -588,6 +588,15 @@ def get_linear_discounted_total( return typing.cast(xarray.DataArray, ret) +def get_linear_discounted( + conversion_year: xarray.DataArray, darray: xarray.DataArray +) -> xarray.DataArray: + """`darray` times the discount weight of the span holding each pixel's conversion year.""" + return get_linear_discounted_total( + span_to_value=get_span_to_charge(conversion_year=conversion_year, darray=darray) + ) + + def get_hectares_per_pixel(darray: xarray.DataArray) -> xarray.DataArray: # https://en.wikipedia.org/wiki/Earth%27s_circumference equator_meters_per_longitude_degrees = 40_075_017 / 360 @@ -646,9 +655,7 @@ def get_span_to_component_to_emissions( A span's emissions are vegetation plus soil, and every unit of both is claimed exactly once: peat takes all the soil it sits under whatever the source class, and the source class takes - the soil that is left. So the three sum back to the span's own total: the statistical leg - builds `emissions_mt` by adding them up, and the jurisdictional-direct leg, which reads - `emissions_mt` from `emissions-per-hectare`, can check them against it. + the soil that is left. So the three sum back to the span's own total. There is no residual, because `get_conversion_emissions` charges only where a conversion fired and every conversion names a source class. Carbon that reaches no conversion is real, @@ -658,18 +665,32 @@ def get_span_to_component_to_emissions( is_peat = dset[gnw_global_peatlands.DATASET.fully_qualified_band_name] == 1 from_forest = conversion.isin(FROM_FOREST) from_grassland = conversion.isin(FROM_GRASSLAND) + span_to_vegetation = get_span_to_charge( + conversion_year=dset["conversion-year"], + darray=dset["vegetation-emissions:tco2e-per-ha"], + ) + span_to_soil = get_span_to_charge( + conversion_year=dset["conversion-year"], + darray=dset["soil-emissions:tco2e-per-ha"], + ) span_to_component_to_emissions: dict[ SpanType, dict[EmissionComponent, xarray.DataArray] ] = {} - for before, after in SPAN_TO_LINEAR_DISCOUNT_WEIGHT: - vegetation = dset[f"vegetation-emissions:tco2e-per-ha:{before:d}-{after:d}"] - soil = dset[f"soil-emissions:tco2e-per-ha:{before:d}-{after:d}"] - span_to_component_to_emissions[(before, after)] = { - EmissionComponent.FOREST: vegetation.where(from_forest, other=0) - + soil.where(from_forest & ~is_peat, other=0), - EmissionComponent.GRASSLAND: vegetation.where(from_grassland, other=0) - + soil.where(from_grassland & ~is_peat, other=0), - EmissionComponent.PEATLAND_CONVERSION: soil.where(is_peat, other=0), + for span in SPAN_TO_LINEAR_DISCOUNT_WEIGHT: + vegetation = span_to_vegetation[span] + soil = span_to_soil[span] + span_to_component_to_emissions[span] = { + EmissionComponent.FOREST: ( + vegetation.where(from_forest, other=0) + + soil.where(from_forest & ~is_peat, other=0) + ).rename("tco2e-per-ha"), + EmissionComponent.GRASSLAND: ( + vegetation.where(from_grassland, other=0) + + soil.where(from_grassland & ~is_peat, other=0) + ).rename("tco2e-per-ha"), + EmissionComponent.PEATLAND_CONVERSION: soil.where(is_peat, other=0).rename( + "tco2e-per-ha" + ), } assert all( set(component_to_emissions) == set(EmissionComponent) @@ -751,22 +772,8 @@ def workflow(tile_id: str) -> xarray.Dataset: is_peatland=is_peatland, soil_organic_carbon=dset[soilgrids_ocs.DATASET.fully_qualified_band_name], ) - span_to_vegetation_emissions = get_span_to_charge( - conversion_year=conversion_record.year, darray=conversion_emissions.vegetation - ) - span_to_soil_emissions = get_span_to_charge( - conversion_year=conversion_record.year, darray=conversion_emissions.soil - ) - - logger.info("Summing emissions from vegetation and soil") - span_to_emissions: dict[SpanType, xarray.DataArray] = { - span: ( - span_to_vegetation_emissions[span] + span_to_soil_emissions[span] - ).rename("tco2e-per-ha") - for span in span_to_vegetation_emissions - } - logger.info("Quantifying and adding peatland occupation emissions") + logger.info("Quantifying peatland occupation emissions") # NB: split by destination on the 30 m grid so each allocation leg can take the peat drained # under its own land class whole, rather than assuming peat is spread evenly across a cell cropland_occupation_emissions = get_peatland_occupation_emissions( @@ -775,49 +782,29 @@ def workflow(tile_id: str) -> xarray.Dataset: pastureland_occupation_emissions = get_peatland_occupation_emissions( destination=conversion_record.to_pasture, is_peatland=is_peatland ) - emissions_per_hectare: xarray.DataArray = ( - get_linear_discounted_total(span_to_value=span_to_emissions) - + cropland_occupation_emissions - + pastureland_occupation_emissions - ) logger.info("Quantifying the source carbon no destination claimed") - dropped_emissions = get_linear_discounted_total( - span_to_value=get_span_to_charge( - conversion_year=conversion_record.year, - darray=get_dropped_emissions( - forest_carbon=forest_carbon, - from_forest=conversion_record.from_forest, - from_pasture=conversion_record.from_pasture, - from_rangeland=conversion_record.from_rangeland, - grassland_carbon=grassland_carbon, - has_destination=conversion_record.to_cropland - | conversion_record.to_pasture, - ), - ) - ).rename("tco2e-per-ha") + dropped_emissions = get_dropped_emissions( + forest_carbon=forest_carbon, + from_forest=conversion_record.from_forest, + from_pasture=conversion_record.from_pasture, + from_rangeland=conversion_record.from_rangeland, + grassland_carbon=grassland_carbon, + has_destination=conversion_record.to_cropland | conversion_record.to_pasture, + ) return get_dset_for_output( name_to_darray={ "conversion": conversion_record.conversion, "conversion-year": conversion_record.year, - "destination-dataset": conversion_record.destination_dataset, - } - | { - f"soil-emissions:{before:d}-{after:d}": darray - for (before, after), darray in span_to_soil_emissions.items() - } - | { - f"vegetation-emissions:{before:d}-{after:d}": darray - for (before, after), darray in span_to_vegetation_emissions.items() - } - | { f"{CROPLAND:s}-{PEATLAND_OCCUPATION:s}": cropland_occupation_emissions, + "destination-dataset": conversion_record.destination_dataset, # NB: charged to nobody, so it is reported beside the total rather than inside it "dropped-emissions": dropped_emissions, - "emissions-per-hectare": emissions_per_hectare, - "hectares-per-pixel": get_hectares_per_pixel(darray=emissions_per_hectare), + "hectares-per-pixel": get_hectares_per_pixel(darray=dropped_emissions), f"{PASTURELAND:s}-{PEATLAND_OCCUPATION:s}": pastureland_occupation_emissions, + "soil-emissions": conversion_emissions.soil, + "vegetation-emissions": conversion_emissions.vegetation, } ) diff --git a/jdluc/jurisdictional_direct.py b/jdluc/jurisdictional_direct.py index f2d81a3..3a7edf2 100644 --- a/jdluc/jurisdictional_direct.py +++ b/jdluc/jurisdictional_direct.py @@ -46,7 +46,6 @@ class Crop(enum.Enum): # Emissions f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha", "dropped-emissions:tco2e-per-ha", - "emissions-per-hectare:tco2e-per-ha", "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), @@ -58,28 +57,27 @@ class Crop(enum.Enum): def get_commodity_name_to_totals( component_to_per_hectare: dict[emit.EmissionComponent, xarray.DataArray], crop_class: xarray.DataArray, + cropland_occupation_per_hectare: xarray.DataArray, crops: tuple[Crop, ...], dropped_per_hectare: xarray.DataArray, - emissions_per_hectare: xarray.DataArray, hectares_per_pixel: xarray.DataArray, is_peatland: xarray.DataArray, pastureland_occupation_per_hectare: xarray.DataArray, - peatland_occupation_per_hectare: xarray.DataArray, ) -> dict[str, dict[str, float]]: - peatland_occupation_emissions = peatland_occupation_per_hectare * hectares_per_pixel + cropland_occupation_emissions = cropland_occupation_per_hectare * hectares_per_pixel pastureland_occupation_emissions = ( pastureland_occupation_per_hectare * hectares_per_pixel ) - # NB: this leg divides emissions among crops alone, so no crop row takes the peat drained under - # pasture, even on a pixel the crop map calls that crop - emissions_mt = ( - emissions_per_hectare - pastureland_occupation_per_hectare - ) * hectares_per_pixel component_to_emissions = { component: per_hectare * hectares_per_pixel for component, per_hectare in component_to_per_hectare.items() } + # NB: this leg divides emissions among crops alone, so no crop row takes the peat drained under + # pasture, even on a pixel the crop map calls that crop + emissions_mt = sum( + component_to_emissions.values(), start=cropland_occupation_emissions + ) crop_to_totals: dict[Crop | emit.NonCommodity, dict[str, xarray.DataArray]] = ( collections.defaultdict(dict) @@ -94,7 +92,7 @@ def get_commodity_name_to_totals( is_peatland ) crop_to_totals[crop]["peatland_occupation_emissions_mt"] = ( - peatland_occupation_emissions.where(crop_mask) + cropland_occupation_emissions.where(crop_mask) ) crop_to_totals[crop]["emissions_mt"] = emissions_mt.where(crop_mask) for component, emissions in component_to_emissions.items(): @@ -110,7 +108,7 @@ def get_commodity_name_to_totals( # NB: zero rather than NaN on crop pixels, so the pasture peat drained under them is kept "emissions_mt": emissions_mt.where(is_unattributed, other=0) + pastureland_occupation_emissions, - "peatland_occupation_emissions_mt": peatland_occupation_emissions.where( + "peatland_occupation_emissions_mt": cropland_occupation_emissions.where( is_unattributed, other=0 ) + pastureland_occupation_emissions, @@ -179,6 +177,12 @@ def workflow( ) for component in emit.EmissionComponent } + | { + "dropped-emissions:tco2e-per-ha": emit.get_linear_discounted( + conversion_year=merged["conversion-year"], + darray=merged["dropped-emissions:tco2e-per-ha"], + ) + } ) def it() -> collections.abc.Iterator[dict[str, float | str]]: @@ -207,9 +211,11 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]: for component in emit.EmissionComponent }, crop_class=clipped[usda_nass_cdl.DATASET.fully_qualified_band_name], + cropland_occupation_per_hectare=clipped[ + f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha" + ], crops=crops, dropped_per_hectare=clipped["dropped-emissions:tco2e-per-ha"], - emissions_per_hectare=clipped["emissions-per-hectare:tco2e-per-ha"], hectares_per_pixel=clipped["hectares-per-pixel:ha"], # NB: no data has to read as not peat rather than truth-testing to peat is_peatland=clipped[ @@ -219,9 +225,6 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]: pastureland_occupation_per_hectare=clipped[ f"{emit.PASTURELAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha" ], - peatland_occupation_per_hectare=clipped[ - f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha" - ], ) for commodity_name, totals in commodity_name_to_totals.items(): yield totals | { diff --git a/jdluc/jurisdictional_direct_proxy.py b/jdluc/jurisdictional_direct_proxy.py index 56dbe3d..6807241 100644 --- a/jdluc/jurisdictional_direct_proxy.py +++ b/jdluc/jurisdictional_direct_proxy.py @@ -37,7 +37,8 @@ def get_name_to_band(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 split by the land class its pixel became, and discounted over - the spans: the shares that divide it come from one year, so no span needs its own band. + the spans: the shares that divide it come from one year, so no span needs its own band. The + source carbon no destination claimed is discounted the same way. """ span_to_component_to_emissions = emit.get_span_to_component_to_emissions(dset=dset) return ( @@ -60,10 +61,15 @@ def get_name_to_band(dset: xarray.Dataset) -> dict[str, xarray.DataArray]: name: dset[name] for name in ( f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha", - "dropped-emissions:tco2e-per-ha", f"{emit.PASTURELAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha", ) } + | { + "dropped-emissions:tco2e-per-ha": emit.get_linear_discounted( + conversion_year=dset["conversion-year"], + darray=dset["dropped-emissions:tco2e-per-ha"], + ) + } | { # NB: a mask here and an averaging downscale after, so this arrives as a cell fraction f"pasture:fraction:{tang_cropgrids.YEAR:d}": ( diff --git a/jdluc/statistical.py b/jdluc/statistical.py index f16e21e..b2d25ae 100644 --- a/jdluc/statistical.py +++ b/jdluc/statistical.py @@ -104,7 +104,6 @@ class Livestock(enum.StrEnum): EMIT_VARIABLE_NAMES = [ f"{emit.CROPLAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha", - "dropped-emissions:tco2e-per-ha", f"{emit.PASTURELAND:s}-{emit.PEATLAND_OCCUPATION:s}:tco2e-per-ha", ] @@ -153,6 +152,14 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: ].where(luc_and_emissions["conversion"].isin(conversions), other=0) derived_variable_names.append(name) + 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=luc_and_emissions["conversion-year"], + darray=luc_and_emissions[name], + ) + derived_variable_names.append(name) + logger.info("Measuring pasture extent at each MAPSPAM snapshot") for year in MAPSPAM_SNAPSHOT_YEARS: name = f"pasture:fraction:{year:d}" diff --git a/tools/build-soy-brazil-figure.py b/tools/build-soy-brazil-figure.py index faded56..b8946ef 100644 --- a/tools/build-soy-brazil-figure.py +++ b/tools/build-soy-brazil-figure.py @@ -23,8 +23,8 @@ ``--overview-resampling mode``: the default averages, the decimated read below pulls from those overviews, and one converted pixel among a hundred averages to a fraction, which survives the "is it NONE" test and lands in the first class -- painting the whole frame. - ``emissions-per-hectare:tco2e-per-ha`` and the two soybean bands are continuous, so they take - the default. Three runs per tile:: + ``vegetation-emissions:tco2e-per-ha``, ``soil-emissions:tco2e-per-ha`` and the two soybean + bands are continuous, so they take the default. Three runs per tile:: DUMP="uv run python tools/dump-zarr-to-cogs.py" $DUMP --overview-resampling mode cogs/00N_050W @@ -37,7 +37,8 @@ `.tif`, so nothing has to be renamed:: cd cogs - for band in conversion conversion-year emissions-per-hectare:tco2e-per-ha \ + for band in conversion conversion-year \ + vegetation-emissions:tco2e-per-ha soil-emissions:tco2e-per-ha \ ifpri:mapspam-physical-area-2000:soybean:ha \ ifpri:mapspam-physical-area-2020:soybean:ha; do gdalbuildvrt "$band.tif.vrt" 00N_050W/"$band".tif 10S_050W/"$band".tif @@ -52,7 +53,8 @@ `emit` publishes no land class per year, so the categorical panel is the `conversion` band -- which of the five conversions fired on each pixel -- and the panel beside it is `conversion-year`, the year that pixel's source class ended. Both the conversion codes and the -lookback are imported rather than mirrored, so a change to either cannot leave this behind. +lookback are imported rather than mirrored, so a change to either cannot leave this behind. The +emissions panel is vegetation plus soil, undiscounted. """ import argparse @@ -97,7 +99,9 @@ } # Fully transparent RGBA, for the pixels no conversion reached TRANSPARENT = (0.0, 0.0, 0.0, 0.0) -EMISSIONS_LABEL = "Land-use-change emissions (tonnes CO₂-equivalent per hectare)" +EMISSIONS_LABEL = ( + "Land-use-change emissions before discounting (tonnes CO₂-equivalent per hectare)" +) SOYBEAN_LABEL = "Soybean area (hectares per 10-kilometre cell)" CONVERSION_LABEL = { @@ -335,9 +339,15 @@ def out(name: str) -> str: # --- read layers (windowed, decimated) --- conversion = read_window(cog("conversion.tif"), bbox, args.width) conversion_year = read_window(cog("conversion-year.tif"), bbox, args.width) - emissions = read_window( - cog("emissions-per-hectare:tco2e-per-ha.tif"), bbox, args.width + vegetation = read_window( + bbox=bbox, + path=cog(name="vegetation-emissions:tco2e-per-ha.tif"), + width=args.width, ) + soil = read_window( + bbox=bbox, path=cog(name="soil-emissions:tco2e-per-ha.tif"), width=args.width + ) + emissions = Layer(data=vegetation.data + soil.data, extent=vegetation.extent) # MapSPAM is ~10 km, so give it a coarser read at the same extent soy_pair = { y: read_window( From 3e737455b31a1abd2daa861130ac6c628b88540e Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 10:26:49 -0800 Subject: [PATCH 2/7] apply discounting per year instead of per span --- jdluc/__tests__/emit_test.py | 112 +++---------------- jdluc/__tests__/statistical_test.py | 73 +++++-------- jdluc/emit.py | 128 ++++++++-------------- jdluc/jurisdictional_direct.py | 13 +-- jdluc/jurisdictional_direct_proxy.py | 24 ++-- jdluc/statistical.py | 158 ++++++++++++++------------- tools/build-soy-brazil-figure.py | 60 ++++------ 7 files changed, 209 insertions(+), 359 deletions(-) diff --git a/jdluc/__tests__/emit_test.py b/jdluc/__tests__/emit_test.py index d177b82..dd4e43c 100644 --- a/jdluc/__tests__/emit_test.py +++ b/jdluc/__tests__/emit_test.py @@ -28,14 +28,13 @@ LOOKBACK_YEARS_RANGE, PEATLAND_EMISSIONS_ANNUAL_TCO2E_PER_HA, PEATLAND_EMISSIONS_PULSE_TCO2E_PER_HA, - SPAN_TO_LINEAR_DISCOUNT_WEIGHT, Conversion, ConversionEmissions, ConversionRecord, DestinationDataset, EmissionComponent, - SpanType, get_belowground_carbon, + get_component_to_emissions, get_conversion_emissions, get_conversion_record, get_dead_organic_matter_carbon, @@ -45,11 +44,8 @@ get_hectares_per_pixel, get_last_departure_year, get_linear_discounted, - get_linear_discounted_total, get_mineral_soil_emissions, get_peatland_occupation_emissions, - get_span_to_charge, - get_span_to_component_to_emissions, ) @@ -430,31 +426,6 @@ def test_get_dropped_emissions() -> None: ) -@pytest.mark.parametrize( - ("conversion_year", "spans"), - ( - pytest.param(2001, {(2000, 2005)}, id="the first year a conversion can carry"), - pytest.param( - 2005, {(2000, 2005)}, id="a shared boundary belongs to the span it closes" - ), - pytest.param(2006, {(2005, 2010)}, id="and not to the one it opens"), - pytest.param( - ASSESSMENT_YEAR, {(2015, 2020)}, id="the last year a conversion can carry" - ), - pytest.param(0, set(), id="no conversion is charged to no span"), - ), -) -def test_get_span_to_charge(conversion_year: int, spans: set[SpanType]) -> None: - # Spans share their boundary years, so a charge landing in two of them would be counted twice - # and one landing in none would vanish without a trace. - result = get_span_to_charge( - conversion_year=get_darray_for_data(data=[[conversion_year]]), - darray=get_darray_for_data(data=[[1]]), - ) - assert set(result) == set(SPAN_TO_LINEAR_DISCOUNT_WEIGHT) - assert {span for span, darray in result.items() if float(darray.sum())} == spans - - def test_get_dset_for_output_appends_each_band_its_units() -> None: # Every emitted band is keyed on this name, so a change here strands the zarr caches keyed on # the old one and KeyErrors the legs that read it. The units come off the DataArray and slot in @@ -473,27 +444,16 @@ def test_get_dset_for_output_appends_each_band_its_units() -> None: } -def test_get_linear_discounted_total() -> None: - result = get_linear_discounted_total( - span_to_value={ - (2000, 2005): get_darray_for_data(data=[[1]]), - (2005, 2010): get_darray_for_data(data=[[2]]), - (2010, 2015): get_darray_for_data(data=[[3]]), - (2015, 2020): get_darray_for_data(data=[[4]]), - } - ) - assert numpy.array_equal(result.data, [[0.625]]) - - def test_get_linear_discounted() -> None: - # A pixel with no conversion, or no year, takes no weight + # A boundary year takes the weight of the span it closes, 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, 2010, 2011, 2020, 0, numpy.nan]] + data=[[2001, 2005, 2006, 2011, 2020, 0, numpy.nan]] ), - darray=get_darray_for_data(data=[[8] * 6]), + darray=get_darray_for_data(data=[[8, 8, 8, 8, 8, numpy.nan, numpy.nan]]), ) - numpy.testing.assert_allclose(result.data, [[0.1, 0.3, 0.5, 0.7, 0, 0]]) + numpy.testing.assert_allclose(result.data, [[0.1, 0.1, 0.3, 0.5, 0.7, 0, 0]]) def test_get_hectares_per_pixel() -> None: @@ -507,38 +467,15 @@ def test_get_hectares_per_pixel() -> None: ) -def test_get_linear_discounted_total_imposes_no_quantity() -> None: - # Area and production run through this too, so a name of its own would be a lie for two of - # the three callers. The caller names the result. - result = get_linear_discounted_total( - span_to_value=dict.fromkeys( - SPAN_TO_LINEAR_DISCOUNT_WEIGHT, get_darray_for_data(data=[[1]]) - ) - ) - assert result.name is None - - -def test_get_linear_discounted_total_over_the_weight_total_is_a_mean() -> None: - # How `statistical` puts area and production on the emissions' window: the same reduction, - # normalised. A quantity that never moves must survive it unchanged. - result = get_linear_discounted_total( - span_to_value=dict.fromkeys( - SPAN_TO_LINEAR_DISCOUNT_WEIGHT, get_darray_for_data(data=[[7]]) - ) - ) / sum(SPAN_TO_LINEAR_DISCOUNT_WEIGHT.values()) - numpy.testing.assert_allclose(result.data, [[7]]) - - VEGETATION = 10.0 SOIL = 100.0 def get_dset_for_pixel(conversion: Conversion, is_peat: float) -> xarray.Dataset: - """One pixel, converted in 2012.""" + """One pixel, with every emissions band set.""" return xarray.Dataset( { "conversion": get_darray_for_data(data=[[conversion]]), - "conversion-year": get_darray_for_data(data=[[2012]]), PEATLANDS.fully_qualified_band_name: get_darray_for_data(data=[[is_peat]]), "soil-emissions:tco2e-per-ha": get_darray_for_data(data=[[SOIL]]), "vegetation-emissions:tco2e-per-ha": get_darray_for_data( @@ -585,42 +522,21 @@ def get_dset_for_pixel(conversion: Conversion, is_peat: float) -> xarray.Dataset ), ), ) -def test_get_span_to_component_to_emissions( +def test_get_component_to_emissions( conversion: Conversion, is_peat: float, source_to_expected: dict[str, float] ) -> None: # The legs sum the components into `emissions_mt`, so a pool dropped, double-counted or gated # on the wrong conversion fails here. assert sum(source_to_expected.values()) == VEGETATION + SOIL - span_to_component_to_emissions = get_span_to_component_to_emissions( + component_to_emissions = get_component_to_emissions( dset=get_dset_for_pixel(conversion=conversion, is_peat=is_peat) ) - assert set(span_to_component_to_emissions) == set(SPAN_TO_LINEAR_DISCOUNT_WEIGHT) - for span, component_to_emissions in span_to_component_to_emissions.items(): - assert set(component_to_emissions) == set(EmissionComponent) - assert { - component: float(darray.sum()) - for component, darray in component_to_emissions.items() - if float(darray.sum()) - } == (source_to_expected if span == (2010, 2015) else {}) - - -def test_get_span_to_component_to_emissions_charges_each_pixel_to_its_year() -> None: - dset = xarray.Dataset( - { - "conversion": get_darray_for_data( - data=[[Conversion.FOREST_TO_CROPLAND] * 3] - ), - "conversion-year": get_darray_for_data(data=[[2003, 2012, 2020]]), - PEATLANDS.fully_qualified_band_name: get_darray_for_data(data=[[0] * 3]), - "soil-emissions:tco2e-per-ha": get_darray_for_data(data=[[10, 20, 30]]), - "vegetation-emissions:tco2e-per-ha": get_darray_for_data(data=[[1, 2, 3]]), - } - ) - result = get_span_to_component_to_emissions(dset=dset) + assert set(component_to_emissions) == set(EmissionComponent) assert { - span: float(component_to_emissions[EmissionComponent.FOREST].sum()) - for span, component_to_emissions in result.items() - } == {(2000, 2005): 11, (2005, 2010): 0, (2010, 2015): 22, (2015, 2020): 33} + component: float(darray.sum()) + for component, darray in component_to_emissions.items() + if float(darray.sum()) + } == source_to_expected @pytest.mark.parametrize( diff --git a/jdluc/__tests__/statistical_test.py b/jdluc/__tests__/statistical_test.py index bbf1438..bfed01f 100644 --- a/jdluc/__tests__/statistical_test.py +++ b/jdluc/__tests__/statistical_test.py @@ -11,8 +11,8 @@ from jdluc.datasets import gpw_livestock, ifpri_mapspam from jdluc.statistical import ( HERD_HEADS_COLUMNS, + MAPSPAM_SPANS, SNAPSHOT_YEAR_TO_WEIGHT, - SPAN_TO_MAPSPAM_SPAN, Commodity, Crop, Livestock, @@ -21,7 +21,6 @@ get_commodity_to_peatland_conversion_share, get_commodity_to_share, get_crop_to_area_share, - get_discounted_snapshot_mean, get_livestock_to_grazing_share, ) @@ -542,22 +541,6 @@ def test_snapshot_year_to_weight() -> None: ) -def test_snapshot_year_to_weight_reduces_as_the_discounted_snapshot_mean_does() -> None: - year_to_value = {2000: 1.0, 2005: 10.0, 2010: 100.0, 2020: 1000.0} - assert get_value( - darray=get_discounted_snapshot_mean( - year_to_value={ - year: xarray.DataArray(value) for year, value in year_to_value.items() - } - ) - ) == pytest.approx( - sum( - weight * year_to_value[year] - for year, weight in SNAPSHOT_YEAR_TO_WEIGHT.items() - ) - ) - - # The default for `get_name_to_totals`'s shares, read-only so every call can share it NO_SHARES: collections.abc.Mapping[Commodity, float] = types.MappingProxyType({}) @@ -576,8 +559,9 @@ def get_name_to_totals( ) -> dict[str, dict[str, float]]: """`get_commodity_name_to_totals` over WHEAT and the livestock rows on `get_dset`'s 2x2 grid. - Every band and share a caller does not set is zero, so each test moves one thing. The shares - are the same in every span. + Every band and share a caller does not set is zero, so each test moves one thing. Each + emissions band is discounted and set in the 2010-2020 span alone, and the shares are the same + in every span. """ commodities = (Crop.WHEAT, *Livestock) @@ -603,22 +587,23 @@ def get_shares( + pasture_occupation, } | { - f"{name:s}:tco2e-per-ha:{before:d}-{after:d}": grid + emissions + f"{name:s}-emissions:tco2e-per-ha:{before:d}-{after:d}": grid + + (emissions if (before, after) == (2010, 2020) else 0.0) for name, emissions in { "forest": forest_emissions, "grassland": grassland_emissions, f"{emit.CROPLAND:s}-peatland-conversion": cropland_pulse, f"{emit.PASTURELAND:s}-peatland-conversion": pasture_pulse, }.items() - for before, after in emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT + for before, after in MAPSPAM_SPANS } ), occupation_shares=get_shares(shares=occupation_shares), span_to_commodity_to_peatland_conversion_share=dict.fromkeys( - SPAN_TO_MAPSPAM_SPAN.values(), get_shares(shares=peatland_conversion_shares) + MAPSPAM_SPANS, get_shares(shares=peatland_conversion_shares) ), span_to_commodity_to_share=dict.fromkeys( - SPAN_TO_MAPSPAM_SPAN.values(), get_shares(shares=commodity_to_share) + MAPSPAM_SPANS, get_shares(shares=commodity_to_share) ), ) @@ -818,8 +803,8 @@ def test_get_commodity_to_peatland_conversion_share( UNATTRIBUTED = emit.NonCommodity.UNATTRIBUTED.name -# Expected values are per hectare of the cell. Every band and share is the same in each span, so a -# conversion component is band x share x 0.2, the four spans' discount weights summed. +# Expected values are per hectare of the cell. Each emissions band is already discounted and the +# shares are the same in every span, so a conversion component is band x share. @pytest.mark.parametrize( ("pasture_fraction", "kwargs", "expected"), ( @@ -858,9 +843,9 @@ def test_get_commodity_to_peatland_conversion_share( }, }, { - (WHEAT, "peatland_conversion_emissions_mt"): 0.1, - (BEEF, "peatland_conversion_emissions_mt"): 0.6, - (PASTURE, "peatland_conversion_emissions_mt"): 0.4, + (WHEAT, "peatland_conversion_emissions_mt"): 0.5, + (BEEF, "peatland_conversion_emissions_mt"): 3.0, + (PASTURE, "peatland_conversion_emissions_mt"): 2.0, }, id="crops-divide-the-cropland-pulse-and-livestock-the-pasture-pulse-by-their-own-shares", ), @@ -870,7 +855,7 @@ def test_get_commodity_to_peatland_conversion_share( "commodity_to_share": {Livestock.PASTURE: 0.25}, "forest_emissions": 10.0, }, - {(PASTURE, "emissions_mt"): 0.5}, + {(PASTURE, "emissions_mt"): 2.5}, id="the-pasture-row-takes-its-pooled-share-of-every-conversion-not-only-pasture-ones", ), pytest.param( @@ -885,11 +870,11 @@ def test_get_commodity_to_peatland_conversion_share( "peatland_conversion_shares": {Crop.WHEAT: 0.25}, }, { - (WHEAT, "forest_emissions_mt"): 0.1, - (WHEAT, "grassland_emissions_mt"): 1.0, - (WHEAT, "peatland_conversion_emissions_mt"): 5.0, + (WHEAT, "forest_emissions_mt"): 0.5, + (WHEAT, "grassland_emissions_mt"): 5.0, + (WHEAT, "peatland_conversion_emissions_mt"): 25.0, (WHEAT, "peatland_occupation_emissions_mt"): 100.0, - (WHEAT, "emissions_mt"): 106.1, + (WHEAT, "emissions_mt"): 130.5, }, id="emissions-mt-is-the-components-plus-occupation-each-by-its-own-share", ), @@ -904,19 +889,19 @@ def test_get_commodity_to_peatland_conversion_share( "occupation_shares": {Crop.WHEAT: 0.1}, "peatland_conversion_shares": {Crop.WHEAT: 0.25}, }, - # Charged: forest 0.2, grassland 2.0, peat conversion 20 and occupation 1000, so - # 1022.2 in all; wheat's row and the unattributed row split each between them + # Charged: forest 1, grassland 10, peat conversion 100 and occupation 1000, so 1111 in + # all; wheat's row and the unattributed row split each between them { - (WHEAT, "forest_emissions_mt"): 0.1, - (WHEAT, "grassland_emissions_mt"): 1.0, - (WHEAT, "peatland_conversion_emissions_mt"): 5.0, + (WHEAT, "forest_emissions_mt"): 0.5, + (WHEAT, "grassland_emissions_mt"): 5.0, + (WHEAT, "peatland_conversion_emissions_mt"): 25.0, (WHEAT, "peatland_occupation_emissions_mt"): 100.0, - (WHEAT, "emissions_mt"): 106.1, - (UNATTRIBUTED, "forest_emissions_mt"): 0.1, - (UNATTRIBUTED, "grassland_emissions_mt"): 1.0, - (UNATTRIBUTED, "peatland_conversion_emissions_mt"): 15.0, + (WHEAT, "emissions_mt"): 130.5, + (UNATTRIBUTED, "forest_emissions_mt"): 0.5, + (UNATTRIBUTED, "grassland_emissions_mt"): 5.0, + (UNATTRIBUTED, "peatland_conversion_emissions_mt"): 75.0, (UNATTRIBUTED, "peatland_occupation_emissions_mt"): 900.0, - (UNATTRIBUTED, "emissions_mt"): 916.1, + (UNATTRIBUTED, "emissions_mt"): 980.5, (UNATTRIBUTED, "production_mt"): math.nan, }, id="the-unattributed-row-is-what-the-shares-leave-of-the-charge-so-the-rows-sum-to-it", diff --git a/jdluc/emit.py b/jdluc/emit.py index f9f9b11..2c083d0 100644 --- a/jdluc/emit.py +++ b/jdluc/emit.py @@ -227,28 +227,24 @@ def get_mineral_soil_emissions( ).rename("tco2e-per-ha") -SpanType = tuple[int, int] -SPAN_TO_LINEAR_DISCOUNT_WEIGHT: dict[SpanType, float] = { - (2000, 2005): 0.0125, - (2005, 2010): 0.0375, - (2010, 2015): 0.0625, - (2015, 2020): 0.0875, -} -assert all((after - before) == 5 for (before, after) in SPAN_TO_LINEAR_DISCOUNT_WEIGHT) -assert math.isclose(sum(SPAN_TO_LINEAR_DISCOUNT_WEIGHT.values()), 0.2) LOOKBACK_YEARS = 20 -assert ( - max(after for _, after in SPAN_TO_LINEAR_DISCOUNT_WEIGHT) - - min(before for before, _ in SPAN_TO_LINEAR_DISCOUNT_WEIGHT) - == LOOKBACK_YEARS -) - - ASSESSMENT_YEAR = 2020 LOOKBACK_YEARS_RANGE = tuple( range(ASSESSMENT_YEAR - LOOKBACK_YEARS, ASSESSMENT_YEAR + 1) ) +# The GHGP linear discount weight by conversion year: the mean of the LSR Guidance's annual factors +# (Table 7.11) over the five-year span holding that year +YEAR_TO_LINEAR_DISCOUNT_WEIGHT: dict[int, float] = ( + dict.fromkeys(range(2001, 2006), 0.0125) + | dict.fromkeys(range(2006, 2011), 0.0375) + | dict.fromkeys(range(2011, 2016), 0.0625) + | dict.fromkeys(range(2016, 2021), 0.0875) +) +# NB: nothing can depart in the lookback's first year, which has no year before it to compare with +assert tuple(YEAR_TO_LINEAR_DISCOUNT_WEIGHT) == LOOKBACK_YEARS_RANGE[1:] +assert math.isclose(sum(YEAR_TO_LINEAR_DISCOUNT_WEIGHT.values()), 1) + class Conversion(enum.IntEnum): @staticmethod @@ -561,40 +557,21 @@ def get_peatland_occupation_emissions( ).rename("tco2e-per-ha") -def get_span_to_charge( - conversion_year: xarray.DataArray, darray: xarray.DataArray -) -> dict[SpanType, xarray.DataArray]: - return { - (before, after): darray.where( - (conversion_year > before) & (conversion_year <= after), other=0 - ) - for before, after in SPAN_TO_LINEAR_DISCOUNT_WEIGHT - } - - -def get_linear_discounted_total( - span_to_value: dict[SpanType, xarray.DataArray], -) -> xarray.DataArray: - """`span_to_value` reduced by the linear discount each span carries. - - Quantity-neutral: a flux -- emissions -- takes the weighted sum as it stands, while a level - -- area, production -- divides the result by the weights' total, making the same reduction - a weighted mean over the same windows. - """ - ret = sum( - value * SPAN_TO_LINEAR_DISCOUNT_WEIGHT[span] - for span, value in span_to_value.items() - ) - return typing.cast(xarray.DataArray, ret) - - def get_linear_discounted( conversion_year: xarray.DataArray, darray: xarray.DataArray ) -> xarray.DataArray: - """`darray` times the discount weight of the span holding each pixel's conversion year.""" - return get_linear_discounted_total( - span_to_value=get_span_to_charge(conversion_year=conversion_year, darray=darray) + """`darray` times the discount weight of each pixel's conversion year, or zero without one.""" + lookup = numpy.zeros(ASSESSMENT_YEAR + 1, dtype=numpy.float32) + for year, weight in YEAR_TO_LINEAR_DISCOUNT_WEIGHT.items(): + lookup[year] = weight + weights = xarray.apply_ufunc( + lookup.__getitem__, + conversion_year.fillna(0).astype(numpy.int64), + dask="parallelized", + output_dtypes=[numpy.float32], ) + # NB: NaN times a zero weight is still NaN, so a pixel with no year is zeroed, not weighed + return (darray * weights).where(weights > 0, other=0).rename(darray.name) def get_hectares_per_pixel(darray: xarray.DataArray) -> xarray.DataArray: @@ -648,14 +625,14 @@ def column(self) -> str: return f"{self.name.lower():s}_emissions_mt" -def get_span_to_component_to_emissions( +def get_component_to_emissions( dset: xarray.Dataset, -) -> dict[SpanType, dict[EmissionComponent, xarray.DataArray]]: - """Split each span's emissions into the three components that claim them. +) -> dict[EmissionComponent, xarray.DataArray]: + """Split each pixel's emissions into the three components that claim them, undiscounted. - A span's emissions are vegetation plus soil, and every unit of both is claimed exactly once: + A pixel's emissions are vegetation plus soil, and every unit of both is claimed exactly once: peat takes all the soil it sits under whatever the source class, and the source class takes - the soil that is left. So the three sum back to the span's own total. + the soil that is left. So the three sum back to the pixel's own total. There is no residual, because `get_conversion_emissions` charges only where a conversion fired and every conversion names a source class. Carbon that reaches no conversion is real, @@ -665,38 +642,23 @@ def get_span_to_component_to_emissions( is_peat = dset[gnw_global_peatlands.DATASET.fully_qualified_band_name] == 1 from_forest = conversion.isin(FROM_FOREST) from_grassland = conversion.isin(FROM_GRASSLAND) - span_to_vegetation = get_span_to_charge( - conversion_year=dset["conversion-year"], - darray=dset["vegetation-emissions:tco2e-per-ha"], - ) - span_to_soil = get_span_to_charge( - conversion_year=dset["conversion-year"], - darray=dset["soil-emissions:tco2e-per-ha"], - ) - span_to_component_to_emissions: dict[ - SpanType, dict[EmissionComponent, xarray.DataArray] - ] = {} - for span in SPAN_TO_LINEAR_DISCOUNT_WEIGHT: - vegetation = span_to_vegetation[span] - soil = span_to_soil[span] - span_to_component_to_emissions[span] = { - EmissionComponent.FOREST: ( - vegetation.where(from_forest, other=0) - + soil.where(from_forest & ~is_peat, other=0) - ).rename("tco2e-per-ha"), - EmissionComponent.GRASSLAND: ( - vegetation.where(from_grassland, other=0) - + soil.where(from_grassland & ~is_peat, other=0) - ).rename("tco2e-per-ha"), - EmissionComponent.PEATLAND_CONVERSION: soil.where(is_peat, other=0).rename( - "tco2e-per-ha" - ), - } - assert all( - set(component_to_emissions) == set(EmissionComponent) - for component_to_emissions in span_to_component_to_emissions.values() - ) - return span_to_component_to_emissions + vegetation = dset["vegetation-emissions:tco2e-per-ha"] + soil = dset["soil-emissions:tco2e-per-ha"] + component_to_emissions = { + EmissionComponent.FOREST: ( + vegetation.where(from_forest, other=0) + + soil.where(from_forest & ~is_peat, other=0) + ).rename("tco2e-per-ha"), + EmissionComponent.GRASSLAND: ( + vegetation.where(from_grassland, other=0) + + soil.where(from_grassland & ~is_peat, other=0) + ).rename("tco2e-per-ha"), + EmissionComponent.PEATLAND_CONVERSION: soil.where(is_peat, other=0).rename( + "tco2e-per-ha" + ), + } + assert set(component_to_emissions) == set(EmissionComponent) + return component_to_emissions def get_dset_for_output(name_to_darray: dict[str, xarray.DataArray]) -> xarray.Dataset: diff --git a/jdluc/jurisdictional_direct.py b/jdluc/jurisdictional_direct.py index 3a7edf2..f7542b7 100644 --- a/jdluc/jurisdictional_direct.py +++ b/jdluc/jurisdictional_direct.py @@ -167,15 +167,12 @@ def workflow( ) merged = merged.assign( { - f"{component!s}-emissions:tco2e-per-ha": emit.get_linear_discounted_total( - span_to_value={ - span: component_to_emissions[component] - for span, component_to_emissions in emit.get_span_to_component_to_emissions( - dset=merged - ).items() - } + f"{component!s}-emissions:tco2e-per-ha": emit.get_linear_discounted( + conversion_year=merged["conversion-year"], darray=emissions ) - for component in emit.EmissionComponent + for component, emissions in emit.get_component_to_emissions( + dset=merged + ).items() } | { "dropped-emissions:tco2e-per-ha": emit.get_linear_discounted( diff --git a/jdluc/jurisdictional_direct_proxy.py b/jdluc/jurisdictional_direct_proxy.py index 6807241..930d71a 100644 --- a/jdluc/jurisdictional_direct_proxy.py +++ b/jdluc/jurisdictional_direct_proxy.py @@ -36,21 +36,21 @@ def get_name_to_band(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 split by the land class its pixel became, and discounted over - the spans: the shares that divide it come from one year, so no span needs its own band. The - source carbon no destination claimed is discounted the same way. + Each emissions component is discounted by its pixel's conversion year and split by the land + class the pixel became: the shares that divide it come from one year, so no span needs its own + band. The source carbon no destination claimed is discounted the same way. """ - span_to_component_to_emissions = emit.get_span_to_component_to_emissions(dset=dset) + component_to_emissions = { + component: emit.get_linear_discounted( + conversion_year=dset["conversion-year"], darray=emissions + ) + for component, emissions in emit.get_component_to_emissions(dset=dset).items() + } return ( { - f"{land_class:s}-{component!s}-emissions:tco2e-per-ha": emit.get_linear_discounted_total( - span_to_value={ - span: component_to_emissions[component].where( - dset["conversion"].isin(conversions), other=0 - ) - for span, component_to_emissions in span_to_component_to_emissions.items() - } - ) + f"{land_class:s}-{component!s}-emissions:tco2e-per-ha": component_to_emissions[ + component + ].where(dset["conversion"].isin(conversions), other=0) for land_class, conversions in ( (emit.CROPLAND, emit.TO_CROPLAND), (emit.PASTURELAND, emit.TO_PASTURE), diff --git a/jdluc/statistical.py b/jdluc/statistical.py index b2d25ae..b0d9651 100644 --- a/jdluc/statistical.py +++ b/jdluc/statistical.py @@ -131,14 +131,24 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: emit.workflow(tile_id=tile_id), ) - logger.info("Splitting emissions by component per span") + 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 + ) + for component, emissions in emit.get_component_to_emissions( + dset=luc_and_emissions + ).items() + } derived_variable_names: list[str] = [] - for (before, after), component_to_darray in emit.get_span_to_component_to_emissions( - dset=luc_and_emissions - ).items(): + for before, after in MAPSPAM_SPANS: + in_span = (conversion_year > before) & (conversion_year <= after) for component in POOLED_COMPONENTS: - name = f"{component!s}:tco2e-per-ha:{before:d}-{after:d}" - luc_and_emissions[name] = component_to_darray[component] + name = f"{component!s}-emissions:tco2e-per-ha:{before:d}-{after:d}" + luc_and_emissions[name] = component_to_emissions[component].where( + in_span, other=0 + ) derived_variable_names.append(name) # NB: `emit` puts soil emissions only on converted pixels, and every conversion ends in # cropland or pasture, so the two bands add up to the whole pulse @@ -146,17 +156,18 @@ def get_downscaled_luc_emissions(tile_id: str) -> xarray.Dataset: (emit.CROPLAND, emit.TO_CROPLAND), (emit.PASTURELAND, emit.TO_PASTURE), ): - name = f"{land_class:s}-{emit.EmissionComponent.PEATLAND_CONVERSION!s}:tco2e-per-ha:{before:d}-{after:d}" - luc_and_emissions[name] = component_to_darray[ + name = f"{land_class:s}-{emit.EmissionComponent.PEATLAND_CONVERSION!s}-emissions:tco2e-per-ha:{before:d}-{after:d}" + luc_and_emissions[name] = component_to_emissions[ emit.EmissionComponent.PEATLAND_CONVERSION - ].where(luc_and_emissions["conversion"].isin(conversions), other=0) + ].where( + in_span & luc_and_emissions["conversion"].isin(conversions), other=0 + ) derived_variable_names.append(name) 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=luc_and_emissions["conversion-year"], - darray=luc_and_emissions[name], + conversion_year=conversion_year, darray=luc_and_emissions[name] ) derived_variable_names.append(name) @@ -365,28 +376,28 @@ def get_commodity_to_peatland_conversion_share( } | get_livestock_to_grazing_share(dset=dset, year=after) -SPAN_TO_MAPSPAM_SPAN = { - (2000, 2005): (2000, 2005), - (2005, 2010): (2005, 2010), - # NB: because MAPSPAM has no 2015 snapshot, we assume crop expansion is constant between - # 2010 and 2020 -- and then use that one crop expansion for both spans - (2010, 2015): (2010, 2020), - (2015, 2020): (2010, 2020), -} -assert set(SPAN_TO_MAPSPAM_SPAN) == set(emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT) +# 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.YEAR_TO_LINEAR_DISCOUNT_WEIGHT +) MAPSPAM_SNAPSHOT_YEARS = tuple( - sorted({year for span in SPAN_TO_MAPSPAM_SPAN.values() for year in span}) + 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` +# 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 = { - year: sum( + snapshot_year: sum( weight / 2 - for span, weight in emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT.items() - if year in SPAN_TO_MAPSPAM_SPAN[span] + for year, weight in emit.YEAR_TO_LINEAR_DISCOUNT_WEIGHT.items() + for before, after in MAPSPAM_SPANS + if before < year <= after and snapshot_year in (before, after) ) - / sum(emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT.values()) - for year in MAPSPAM_SNAPSHOT_YEARS + for snapshot_year in MAPSPAM_SNAPSHOT_YEARS } assert math.isclose(sum(SNAPSHOT_YEAR_TO_WEIGHT.values()), 1) @@ -404,12 +415,13 @@ def get_discounted_snapshot_mean( year_to_value: dict[int, xarray.DataArray], ) -> xarray.DataArray: """A level -- area, production -- as a mean over the windows emissions are summed on.""" - return emit.get_linear_discounted_total( - span_to_value={ - span: (year_to_value[before] + year_to_value[after]) / 2 - for span, (before, after) in SPAN_TO_MAPSPAM_SPAN.items() - } - ) / sum(emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT.values()) + 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( @@ -417,15 +429,17 @@ def get_commodity_name_to_totals( dset: xarray.Dataset, occupation_shares: dict[Commodity, xarray.DataArray], span_to_commodity_to_peatland_conversion_share: dict[ - emit.SpanType, dict[Commodity, xarray.DataArray] + tuple[int, int], dict[Commodity, xarray.DataArray] + ], + span_to_commodity_to_share: dict[ + tuple[int, int], dict[Commodity, xarray.DataArray] ], - span_to_commodity_to_share: dict[emit.SpanType, dict[Commodity, xarray.DataArray]], ) -> dict[str, dict[str, float]]: """Each commodity's emissions and area, a crop's production and a herd's heads, summed over a jurisdiction's cells. - Both share mappings are keyed by MapSPAM span. Each emissions span uses the shares of the - MapSPAM span that contains it (`SPAN_TO_MAPSPAM_SPAN`). + Both share mappings are keyed by MapSPAM span, and each divides that span's discounted + emissions bands. Beside the commodities sit `DROPPED` and `UNATTRIBUTED` (see `emit.NonCommodity`). Here `UNATTRIBUTED` is what the shares leave of the charge: conversion and occupation by the @@ -451,31 +465,29 @@ def get_commodity_name_to_totals( emit.PASTURELAND if isinstance(commodity, Livestock) else emit.CROPLAND ) component_to_emissions = { - component: emit.get_linear_discounted_total( - span_to_value={ - (before, after): dset[ - f"{component!s}:tco2e-per-ha:{before:d}-{after:d}" - ] + component: sum( + ( + dset[f"{component!s}-emissions:tco2e-per-ha:{before:d}-{after:d}"] * hectares - * span_to_commodity_to_share[SPAN_TO_MAPSPAM_SPAN[(before, after)]][ - commodity - ] - for (before, after) in emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT - } + * span_to_commodity_to_share[(before, after)][commodity] + for before, after in MAPSPAM_SPANS + ), + start=xarray.zeros_like(hectares), ) for component in POOLED_COMPONENTS } | { - emit.EmissionComponent.PEATLAND_CONVERSION: emit.get_linear_discounted_total( - span_to_value={ - (before, after): dset[ - f"{land_class:s}-{emit.EmissionComponent.PEATLAND_CONVERSION!s}:tco2e-per-ha:{before:d}-{after:d}" + emit.EmissionComponent.PEATLAND_CONVERSION: sum( + ( + dset[ + f"{land_class:s}-{emit.EmissionComponent.PEATLAND_CONVERSION!s}-emissions:tco2e-per-ha:{before:d}-{after:d}" ] * hectares - * span_to_commodity_to_peatland_conversion_share[ - SPAN_TO_MAPSPAM_SPAN[(before, after)] - ][commodity] - for (before, after) in emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT - } + * span_to_commodity_to_peatland_conversion_share[(before, after)][ + commodity + ] + for before, after in MAPSPAM_SPANS + ), + start=xarray.zeros_like(hectares), ) } totals |= { @@ -538,28 +550,26 @@ def get_commodity_name_to_totals( # The whole charge each component carries, as a commodity holding every share would take it component_to_charged = { - component: emit.get_linear_discounted_total( - span_to_value={ - (before, after): dset[ - f"{component!s}:tco2e-per-ha:{before:d}-{after:d}" - ] + component: sum( + ( + dset[f"{component!s}-emissions:tco2e-per-ha:{before:d}-{after:d}"] * hectares - for (before, after) in emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT - } + for before, after in MAPSPAM_SPANS + ), + start=xarray.zeros_like(hectares), ) for component in POOLED_COMPONENTS } | { - emit.EmissionComponent.PEATLAND_CONVERSION: emit.get_linear_discounted_total( - span_to_value={ - (before, after): sum( - dset[ - f"{land_class:s}-{emit.EmissionComponent.PEATLAND_CONVERSION!s}:tco2e-per-ha:{before:d}-{after:d}" - ] - for land_class in (emit.CROPLAND, emit.PASTURELAND) - ) + emit.EmissionComponent.PEATLAND_CONVERSION: sum( + ( + dset[ + f"{land_class:s}-{emit.EmissionComponent.PEATLAND_CONVERSION!s}-emissions:tco2e-per-ha:{before:d}-{after:d}" + ] * hectares - for (before, after) in emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT - } + for land_class in (emit.CROPLAND, emit.PASTURELAND) + for before, after in MAPSPAM_SPANS + ), + start=xarray.zeros_like(hectares), ) } unattributed = commodity_to_totals[emit.NonCommodity.UNATTRIBUTED] @@ -712,7 +722,7 @@ def it() -> collections.abc.Iterator[dict[str, float | str]]: crops=crops, dset=clipped, ) - for (before, after) in dict.fromkeys(SPAN_TO_MAPSPAM_SPAN.values()) + for (before, after) in MAPSPAM_SPANS } span_to_commodity_to_peatland_conversion_share = { (before, after): get_commodity_to_peatland_conversion_share( diff --git a/tools/build-soy-brazil-figure.py b/tools/build-soy-brazil-figure.py index b8946ef..9b60d32 100644 --- a/tools/build-soy-brazil-figure.py +++ b/tools/build-soy-brazil-figure.py @@ -83,10 +83,6 @@ # 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:] -# A conversion is charged to the span holding its year, and each span carries its own discount -# weight, so colouring by span says what a continuous ramp could not: how heavily that pixel -# counts. Sorted oldest first, and shaded light to dark with recency. -SPANS = sorted(emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT) # Hotter reads as more carbon released: forest sources red and orange, grassland tan and olive @@ -199,47 +195,27 @@ def draw_conversion(ax: plt.Axes, layer: Layer, title: str) -> None: _style_geo_axes(ax, title, layer.extent) -def span_artists() -> tuple[ - matplotlib.colors.ListedColormap, matplotlib.colors.BoundaryNorm, list[Line2D] -]: - """Discrete colormap keyed to the span index, plus legend handles reading oldest first.""" - colours = matplotlib.colormaps["magma"](numpy.linspace(0.75, 0.2, len(SPANS))) +def draw_conversion_year( + ax: plt.Axes, layer: Layer, title: str +) -> matplotlib.image.AxesImage: + """Each converted pixel's year, one colour per year, darker as the discount weighs it more.""" + colours = matplotlib.colormaps["magma"]( + numpy.linspace(0.75, 0.2, len(LOOKBACK_YEARS)) + ) cmap = matplotlib.colors.ListedColormap(colours).with_extremes(bad=TRANSPARENT) norm = matplotlib.colors.BoundaryNorm( - [index - 0.5 for index in range(len(SPANS) + 1)], cmap.N + [*(year - 0.5 for year in LOOKBACK_YEARS), LOOKBACK_YEARS[-1] + 0.5], cmap.N ) - handles = [ - Line2D( - [], - [], - label=f"{before + 1:d}-{after:d}", - linestyle="", - marker="s", - markeredgecolor="none", - markerfacecolor=colours[index], - markersize=10, - ) - for index, (before, after) in enumerate(SPANS) - ] - return cmap, norm, handles - - -def draw_conversion_span(ax: plt.Axes, layer: Layer, title: str) -> None: - """Which five-year span each converted pixel was charged to; zero means it never converted.""" - cmap, norm, _ = span_artists() - data = numpy.full(layer.data.shape, numpy.nan) - for index, (before, after) in enumerate(SPANS): - # Open below and closed above, exactly as `emit.get_span_to_charge` bins them - data = numpy.where((layer.data > before) & (layer.data <= after), index, data) - ax.imshow( - data, + im = ax.imshow( + numpy.where(numpy.isin(layer.data, LOOKBACK_YEARS), layer.data, numpy.nan), cmap=cmap, extent=layer.extent, interpolation="nearest", norm=norm, origin="upper", ) - _style_geo_axes(ax, title, layer.extent) + _style_geo_axes(ax=ax, extent=layer.extent, title=title) + return im def draw_emissions( @@ -375,11 +351,15 @@ def out(name: str) -> str: ax_sy0, ax_sy1 = fig.add_subplot(left[1, 0]), fig.add_subplot(left[1, 1]) ax_em = fig.add_subplot(outer[0, 1]) - draw_conversion_span(ax_lc0, conversion_year, "") - _, _, span_handles = span_artists() - ax_lc0.legend( - fontsize=7, framealpha=0.85, handles=span_handles, loc="lower left", ncol=2 + im_year = draw_conversion_year(ax=ax_lc0, layer=conversion_year, title="") + # NB: inset where a legend would sit, so this panel keeps the size of the one beside it + year_bar = fig.colorbar( + im_year, + cax=ax_lc0.inset_axes((0.05, 0.1, 0.5, 0.04)), + orientation="horizontal", + ticks=(LOOKBACK_YEARS[0], *LOOKBACK_YEARS[4::5]), ) + year_bar.ax.tick_params(labelsize=7) draw_conversion(ax_lc1, conversion, "") ax_lc1.legend(fontsize=7, framealpha=0.85, handles=legend_handles, loc="lower left") From e1f26e5b47be7edb05a0968845b54a12362b200f Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 10:32:43 -0800 Subject: [PATCH 3/7] update docs --- docs/data.md | 27 +++++++++++++-------------- docs/executive_summary.md | 2 +- docs/further_research.md | 2 +- docs/methodology.md | 18 +++++++++--------- docs/validation.md | 2 +- 5 files changed, 25 insertions(+), 26 deletions(-) diff --git a/docs/data.md b/docs/data.md index e860fed..537226a 100644 --- a/docs/data.md +++ b/docs/data.md @@ -57,20 +57,19 @@ 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, per-span LUC emissions produced by the `emit` stage. A conversion is dated to the year its source class ended, and charged to the span containing that year; the spans are `2000-2005`, `2005-2010`, `2010-2015` and `2015-2020`. - -| Variable | Type | Units | Description | -| ---------------------------------------------------- | ------- | --------------- | -------------------------------------------------------------------------------------------------------- | -| `conversion` | float32 | Conversion code | Which of the five conversions fired on the pixel, or none | -| `conversion-year` | float32 | year | The year the pixel's source class ended — its tree-cover loss year, or its last departure from grassland | -| `destination-dataset` | float32 | Dataset bitmask | Which layers claimed the pixel at the assessment year — one bit each, so a contested pixel keeps both | -| `vegetation-emissions:tco2e-per-ha:{before}-{after}` | float32 | t CO₂e/ha | Per-span vegetation-carbon loss (above-ground, below-ground, dead organic matter, grassland) | -| `soil-emissions:tco2e-per-ha:{before}-{after}` | float32 | t CO₂e/ha | Per-span soil-carbon loss (mineral stock change + peatland drainage pulse) | -| `cropland-peatland-occupation:tco2e-per-ha` | float32 | t CO₂e/ha | Current-year annual peatland-occupation emissions under cropland | -| `pastureland-peatland-occupation:tco2e-per-ha` | float32 | t CO₂e/ha | The same under pasture, split from cropland on the 30 m grid | -| `dropped-emissions:tco2e-per-ha` | float32 | t CO₂e/ha | Source carbon no destination claimed, reported beside the total rather than inside it | -| `emissions-per-hectare:tco2e-per-ha` | float32 | t CO₂e/ha | 20-year linearly-discounted LUC emissions + peatland occupation — the headline per-hectare layer | -| `hectares-per-pixel:ha` | float32 | ha | Pixel area (varies with latitude) | +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)). + +| Variable | Type | Units | Description | +| ---------------------------------------------- | ------- | --------------- | --------------------------------------------------------------------------------------------------------------------------- | +| `conversion` | float32 | Conversion code | Which of the five conversions fired on the pixel, or none | +| `conversion-year` | float32 | year | The year the pixel's source class ended — its tree-cover loss year, or its last departure from grassland | +| `destination-dataset` | float32 | Dataset bitmask | Which layers claimed the pixel at the assessment year — one bit each, so a contested pixel keeps both | +| `vegetation-emissions:tco2e-per-ha` | float32 | t CO₂e/ha | Vegetation-carbon loss of the pixel's conversion (above-ground, below-ground, dead organic matter, grassland), undiscounted | +| `soil-emissions:tco2e-per-ha` | float32 | t CO₂e/ha | Soil-carbon loss of the pixel's conversion (mineral stock change + peatland drainage pulse), undiscounted | +| `cropland-peatland-occupation:tco2e-per-ha` | float32 | t CO₂e/ha | Current-year annual peatland-occupation emissions under cropland | +| `pastureland-peatland-occupation:tco2e-per-ha` | float32 | t CO₂e/ha | The same under pasture, split from cropland on the 30 m grid | +| `dropped-emissions:tco2e-per-ha` | float32 | t CO₂e/ha | Source carbon no destination claimed, undiscounted; reported beside the charged emissions rather than inside them | +| `hectares-per-pixel:ha` | float32 | ha | Pixel area (varies with latitude) | ```python import xarray diff --git a/docs/executive_summary.md b/docs/executive_summary.md index b420603..32592ff 100644 --- a/docs/executive_summary.md +++ b/docs/executive_summary.md @@ -9,7 +9,7 @@ This document is aimed at readers assessing how the methodology maps to the LSRS | LSRS requirement | How this methodology addresses it | Status | | ----------------------------------------------------------- | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | -------------------------------------------------- | | **Land-use-change emissions in scope** | Detects conversions out of forest, rangeland and pasture per ~30 m pixel from annual layers spanning 2000–2020, and quantifies the carbon lost. | Met | -| **20-year assessment period, linearly discounted (§7.2.1)** | All fluxes are allocated over a 20-year linearly-discounted lookback; more recent conversions are weighted more heavily. A conversion is dated to the year its source class ended and charged to the five-year span holding that year, which carries one aggregated discount weight. | Met | +| **20-year assessment period, linearly discounted (§7.2.1)** | All fluxes are allocated over a 20-year linearly-discounted lookback; more recent conversions are weighted more heavily. A conversion is dated to the year its source class ended, and each year carries the mean discount weight of the five-year span holding it. | Met | | **Carbon pools (vegetation + soil)** | Includes key carbon pools described on the LSRS as follows: Above-ground biomass (Harris 2021), below-ground/root biomass (Huang 2021), dead organic matter (CDM AR-TOOL-12, forest only), mineral soil organic carbon (IPCC 2019 Tier 1 stock-change, 0–30 cm) and complements with peatland emissions (newly drained + long term drained). | Met | | **Jurisdictional direct vs. statistical land use change** | Two legs share one emissions core: jdLUC via spatial intersection with high-resolution crop maps (USDA CDL), and sLUC via each crop's share of local cropland expansion (IFPRI MapSPAM). | Met | | **Allocation to commodities** | Product expansion approach is used for sLUC calculation; pasture's share is divided among livestock commodities by grazing pressure. | Met | diff --git a/docs/further_research.md b/docs/further_research.md index ca6f5e1..eed858a 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 `SPAN_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 spans — 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 `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. diff --git a/docs/methodology.md b/docs/methodology.md index 8399524..81b9f96 100644 --- a/docs/methodology.md +++ b/docs/methodology.md @@ -128,7 +128,7 @@ ______________________________________________________________________ ## 2. Quantifying emissions -The emissions core (`emit.py`) computes **per-pixel, per-span LUC emissions** on the GLAD tile grid, independent of any crop. This single layer feeds both attribution legs. +The emissions core (`emit.py`) computes **per-pixel LUC emissions**, dated by conversion year, on the GLAD tile grid, independent of any crop. This single layer feeds both attribution legs. ### Cataloguing transitions @@ -174,9 +174,9 @@ 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 charged to the five-year span containing that year, and the GHGP's per-year weights are aggregated into one weight per span as the unbiased mean of that span's candidate conversion years (`SPAN_TO_LINEAR_DISCOUNT_WEIGHT` in `emit.py`). The span stays the accounting unit even though tree-cover loss dates annually, because it is what the MapSPAM snapshots behind the statistical leg can speak to. +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 each year carries the mean of the GHGP's per-year weights over the five-year span holding it (`YEAR_TO_LINEAR_DISCOUNT_WEIGHT` in `emit.py`): 1.25% for 2001–2005, 3.75% for 2006–2010, 6.25% for 2011–2015 and 8.75% for 2016–2020, 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 a discounted sum of span conversion emissions plus the current year's peatland occupation emissions, scaled by pixel area. This same temporal ramp is shared by both attribution methodologies. +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. ### What the emissions layer contains @@ -191,7 +191,7 @@ Every unit of gross land-carbon decrease inside the lookback falls into one of s | 5 | **Dropped: emissive, attributable to nobody** | Forest → built-up, and forest → natural rangeland | | 6 | **Out of frame** | Clearance before the lookback opens; degradation short of clearing the stand; regrowth, which the LSRS treats as a removal rather than a reduction | -Buckets 3 to 5 are published per pixel as `dropped-emissions`, deliberately outside `emissions-per-hectare` so that no total silently contains carbon charged to nobody. In the attribution tables, buckets 3 to 5 are the `DROPPED` row and bucket 2 is the `UNATTRIBUTED` row (see "Rows that name no commodity" in §3). +Buckets 3 to 5 are published per pixel as `dropped-emissions`, deliberately apart from the vegetation and soil emissions bands so that no total silently contains carbon charged to nobody. In the attribution tables, buckets 3 to 5 are the `DROPPED` row and bucket 2 is the `UNATTRIBUTED` row (see "Rows that name no commodity" in §3). The claim this supports, restricted to conversions that are both emissive and attributable to a commodity: **complete for annual-cropland, oil palm and pasture destinations, and blind to every other woody and aquatic commodity destination**. Oil palm is the one woody perennial with a dedicated destination layer; cocoa, coconut, rubber, coffee, tea and plantation timber have none here and remain in bucket 4. Buckets 3 and 4 are conformance gaps rather than scoping preferences — LSRS Requirement 11 names tree crops and plantation forest explicitly. Bucket 5 is a legitimate scope choice, since neither destination is attributable to a commodity, but both are still reported. @@ -207,7 +207,7 @@ All three legs quantify and allocate LUC emissions to a crop from a specific reg ### Rows that name no commodity -Each jurisdiction's table, in every leg, has two rows that are not commodities: `DROPPED` and `UNATTRIBUTED`. Neither has a production figure, so neither has an emissions factor, and neither changes any crop's factor. They are easy to mix up. The difference is whether their emissions count toward the jurisdiction's emissions total, which is the `emissions-per-hectare` layer added up over the jurisdiction: +Each jurisdiction's table, in every leg, has two rows that are not commodities: `DROPPED` and `UNATTRIBUTED`. Neither has a production figure, so neither has an emissions factor, and neither changes any crop's factor. They are easy to mix up. The difference is whether their emissions count toward the jurisdiction's emissions total, which is the jurisdiction's discounted conversion emissions plus its peatland occupation emissions: - **`DROPPED` does not count toward the total.** It is carbon released by clearing forest or grassland where, as far as our maps show, the land did not become cropland or pasture: for example, it was built on, or it became a tree crop we have no map for (§2, buckets 3 to 5). Its row has only a total, `emissions_mt`. - **`UNATTRIBUTED` counts toward the total.** It is emissions from land that did become cropland or pasture, but that none of the commodities we model takes, mostly land used for crops or grazing we don't model (§2, bucket 2). Its row breaks the emissions down the same way a commodity row does. @@ -234,9 +234,9 @@ When crop production is only available at low spatial resolution, spatial inters Concretely: -1. Downsample the per-pixel emissions to the MapSPAM ~10 km grid — per-span forest and grassland conversion emissions, the per-span peat drainage pulse split by destination (cropland or pastureland), 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. Because MapSPAM has no 2015 snapshot, the 2010→2015 and 2015→2020 spans both use 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 forest and grassland conversion emissions to commodities by that span's expansion share, and its peat drainage pulse by destination (below), weighting by the same GHGP span discount weights. +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. +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. **Pastureland is an expansion term, divided among its grazers.** Grazing land is inside the boundary the LSRS draws, and forest cleared for pasture emits whether or not a crop ever follows, so pasture takes its share of the cell alongside the crops. This matches the external references, which divide by cropland plus cultivated pasture rather than by cropland alone. The consequence for every crop is that its share is diluted by whatever pasture expanded beside it — in pasture-led frontiers that is most of the cell. @@ -261,7 +261,7 @@ 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 span's production is the mean of its two MapSPAM snapshot years, weighted by that span's discount weight and normalized by the sum of the weights. 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, 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. ### Jurisdictional direct proxy (`jurisdictional_direct_proxy.py`) diff --git a/docs/validation.md b/docs/validation.md index 47e7d98..c8bf266 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -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 weights in `emit.SPAN_TO_LINEAR_DISCOUNT_WEIGHT`: WRI's five-year block sums are exactly 5× ours, because ours multiplies a span *total* by a span-average factor where WRI multiplies annual emissions by annual factors. +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). `emit.YEAR_TO_LINEAR_DISCOUNT_WEIGHT` gives each loss year the mean of those factors over its five-year block, so its block sums equal WRI's exactly and only the slope within a block differs. **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. From 9f780444cd05228599704259e784378602fd191b Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 10:44:58 -0800 Subject: [PATCH 4/7] fix measure-drift API --- tools/measure-drift.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tools/measure-drift.py b/tools/measure-drift.py index e43fc9f..76ecf11 100644 --- a/tools/measure-drift.py +++ b/tools/measure-drift.py @@ -929,6 +929,7 @@ def capture(directory: pathlib.Path, label: Side, targets: tuple[Target, ...]) - parameters = inspect.signature(trace.workflow).parameters df = trace.workflow( 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}, From 01357d4390c646e017066da4b65df08729bce69d Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 11:45:47 -0800 Subject: [PATCH 5/7] Switch weights from plateaus within the 5-year spans, to annual variation --- jdluc/__tests__/emit_test.py | 6 +++--- .../jurisdictional_direct_proxy_test.py | 4 ++-- jdluc/emit.py | 17 +++++++---------- 3 files changed, 12 insertions(+), 15 deletions(-) diff --git a/jdluc/__tests__/emit_test.py b/jdluc/__tests__/emit_test.py index dd4e43c..bc0a7ff 100644 --- a/jdluc/__tests__/emit_test.py +++ b/jdluc/__tests__/emit_test.py @@ -445,15 +445,15 @@ def test_get_dset_for_output_appends_each_band_its_units() -> None: def test_get_linear_discounted() -> None: - # A boundary year takes the weight of the span it closes, and a pixel with no year takes none - # even where its value is missing + # 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]] ), darray=get_darray_for_data(data=[[8, 8, 8, 8, 8, numpy.nan, numpy.nan]]), ) - numpy.testing.assert_allclose(result.data, [[0.1, 0.1, 0.3, 0.5, 0.7, 0, 0]]) + numpy.testing.assert_allclose(result.data, [[0.02, 0.18, 0.22, 0.42, 0.78, 0, 0]]) def test_get_hectares_per_pixel() -> None: diff --git a/jdluc/__tests__/jurisdictional_direct_proxy_test.py b/jdluc/__tests__/jurisdictional_direct_proxy_test.py index 827af07..45b0524 100644 --- a/jdluc/__tests__/jurisdictional_direct_proxy_test.py +++ b/jdluc/__tests__/jurisdictional_direct_proxy_test.py @@ -26,8 +26,8 @@ ) from jdluc.statistical import Livestock -# A conversion in 2012 takes the 2010-2015 weight of 0.0625, so 160 and 1600 discount to 10 and 100 -CONVERSION_YEAR = 2012 +# A conversion in 2013 takes a factor of 0.0625, so 160 and 1600 discount to 10 and 100 +CONVERSION_YEAR = 2013 VEGETATION = 160.0 SOIL = 1600.0 LAND_CLASS_COMPONENT_NAMES = [ diff --git a/jdluc/emit.py b/jdluc/emit.py index 2c083d0..df0f12a 100644 --- a/jdluc/emit.py +++ b/jdluc/emit.py @@ -233,16 +233,13 @@ def get_mineral_soil_emissions( range(ASSESSMENT_YEAR - LOOKBACK_YEARS, ASSESSMENT_YEAR + 1) ) -# The GHGP linear discount weight by conversion year: the mean of the LSR Guidance's annual factors -# (Table 7.11) over the five-year span holding that year -YEAR_TO_LINEAR_DISCOUNT_WEIGHT: dict[int, float] = ( - dict.fromkeys(range(2001, 2006), 0.0125) - | dict.fromkeys(range(2006, 2011), 0.0375) - | dict.fromkeys(range(2011, 2016), 0.0625) - | dict.fromkeys(range(2016, 2021), 0.0875) -) -# NB: nothing can depart in the lookback's first year, which has no year before it to compare with -assert tuple(YEAR_TO_LINEAR_DISCOUNT_WEIGHT) == LOOKBACK_YEARS_RANGE[1:] +# 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) From 5bb2b14679fd0ca1fb6ac8132b1361fbe13b38a9 Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Mon, 9 Nov 2026 11:47:30 -0800 Subject: [PATCH 6/7] update docs 2 --- docs/executive_summary.md | 2 +- docs/methodology.md | 2 +- docs/validation.md | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/executive_summary.md b/docs/executive_summary.md index 32592ff..81f0aa2 100644 --- a/docs/executive_summary.md +++ b/docs/executive_summary.md @@ -9,7 +9,7 @@ This document is aimed at readers assessing how the methodology maps to the LSRS | LSRS requirement | How this methodology addresses it | Status | | ----------------------------------------------------------- | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | -------------------------------------------------- | | **Land-use-change emissions in scope** | Detects conversions out of forest, rangeland and pasture per ~30 m pixel from annual layers spanning 2000–2020, and quantifies the carbon lost. | Met | -| **20-year assessment period, linearly discounted (§7.2.1)** | All fluxes are allocated over a 20-year linearly-discounted lookback; more recent conversions are weighted more heavily. A conversion is dated to the year its source class ended, and each year carries the mean discount weight of the five-year span holding it. | Met | +| **20-year assessment period, linearly discounted (§7.2.1)** | All fluxes are allocated over a 20-year linearly-discounted lookback; more recent conversions are weighted more heavily. A conversion is dated to the year its source class ended and takes that year's linear amortization factor. | Met | | **Carbon pools (vegetation + soil)** | Includes key carbon pools described on the LSRS as follows: Above-ground biomass (Harris 2021), below-ground/root biomass (Huang 2021), dead organic matter (CDM AR-TOOL-12, forest only), mineral soil organic carbon (IPCC 2019 Tier 1 stock-change, 0–30 cm) and complements with peatland emissions (newly drained + long term drained). | Met | | **Jurisdictional direct vs. statistical land use change** | Two legs share one emissions core: jdLUC via spatial intersection with high-resolution crop maps (USDA CDL), and sLUC via each crop's share of local cropland expansion (IFPRI MapSPAM). | Met | | **Allocation to commodities** | Product expansion approach is used for sLUC calculation; pasture's share is divided among livestock commodities by grazing pressure. | Met | diff --git a/docs/methodology.md b/docs/methodology.md index 81b9f96..a4a9b29 100644 --- a/docs/methodology.md +++ b/docs/methodology.md @@ -174,7 +174,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 each year carries the mean of the GHGP's per-year weights over the five-year span holding it (`YEAR_TO_LINEAR_DISCOUNT_WEIGHT` in `emit.py`): 1.25% for 2001–2005, 3.75% for 2006–2010, 6.25% for 2011–2015 and 8.75% for 2016–2020, 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; `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. 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. diff --git a/docs/validation.md b/docs/validation.md index c8bf266..afa7632 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -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). `emit.YEAR_TO_LINEAR_DISCOUNT_WEIGHT` gives each loss year the mean of those factors over its five-year block, so its block sums equal WRI's exactly and only the slope within a block differs. +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. **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. From d5f54a2b0d2b34be7107b0622f2366d971a78a1f Mon Sep 17 00:00:00 2001 From: Jeffrey Seifried Date: Tue, 10 Nov 2026 09:35:41 -0800 Subject: [PATCH 7/7] Make get_linear_discounted more robust --- jdluc/__tests__/emit_test.py | 32 +++++++++++++++++++++++++++++ jdluc/emit.py | 39 +++++++++++++++++++++++++++--------- 2 files changed, 61 insertions(+), 10 deletions(-) diff --git a/jdluc/__tests__/emit_test.py b/jdluc/__tests__/emit_test.py index bc0a7ff..b22548b 100644 --- a/jdluc/__tests__/emit_test.py +++ b/jdluc/__tests__/emit_test.py @@ -44,6 +44,7 @@ get_hectares_per_pixel, get_last_departure_year, get_linear_discounted, + get_linear_discounted_values, get_mineral_soil_emissions, get_peatland_occupation_emissions, ) @@ -456,6 +457,37 @@ def test_get_linear_discounted() -> None: numpy.testing.assert_allclose(result.data, [[0.02, 0.18, 0.22, 0.42, 0.78, 0, 0]]) +@pytest.mark.parametrize( + "conversion_year", + ( + pytest.param(2021, id="a-year-after-the-assessment-year"), + pytest.param(2000, id="the-lookback-year-nothing-can-depart-in"), + pytest.param(0, id="no-year-at-all"), + pytest.param(-1, id="a-negative-year-that-would-index-from-the-end"), + pytest.param(2012.5, id="a-fractional-year-that-would-truncate"), + ), +) +def test_get_linear_discounted_values_rejects_emissions_it_cannot_date( + conversion_year: float, +) -> None: + with pytest.raises(AssertionError): + get_linear_discounted_values( + conversion_years=numpy.array([2012, conversion_year]), + values=numpy.array([1.0, 5.0]), + ) + + +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( + conversion_year=conversion_year, + darray=get_darray_for_data(data=[[1, 1]]).assign_coords( + x=conversion_year.x + 1 + ), + ) + + def test_get_hectares_per_pixel() -> None: result = get_hectares_per_pixel(darray=get_darray_for_data(data=[[0] * 2] * 2)) assert result.name == "ha" diff --git a/jdluc/emit.py b/jdluc/emit.py index df0f12a..ac75789 100644 --- a/jdluc/emit.py +++ b/jdluc/emit.py @@ -554,21 +554,40 @@ def get_peatland_occupation_emissions( ).rename("tco2e-per-ha") +def get_linear_discounted_values( + conversion_years: numpy.ndarray, values: numpy.ndarray +) -> 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(): + lookup[year] = weight + weights = lookup[indices] + # NB: emissions without a weight -- no year, or one before the lookback -- would be zeroed + # below, discarding their carbon; only a missing value may go without one + is_undated = (weights == 0) & numpy.isfinite(values) & (values != 0) + assert not is_undated.any(), "Emissions have no conversion year in the lookback" + # NB: NaN times a zero weight is still NaN, so a value with no year is zeroed, not weighed + return numpy.where(weights > 0, values * weights, 0) + + def get_linear_discounted( conversion_year: xarray.DataArray, darray: xarray.DataArray ) -> xarray.DataArray: """`darray` times the discount weight of each pixel's conversion year, or zero without one.""" - lookup = numpy.zeros(ASSESSMENT_YEAR + 1, dtype=numpy.float32) - for year, weight in YEAR_TO_LINEAR_DISCOUNT_WEIGHT.items(): - lookup[year] = weight - weights = xarray.apply_ufunc( - lookup.__getitem__, - conversion_year.fillna(0).astype(numpy.int64), + return xarray.apply_ufunc( + get_linear_discounted_values, + conversion_year.fillna(0), + darray, dask="parallelized", - output_dtypes=[numpy.float32], - ) - # NB: NaN times a zero weight is still NaN, so a pixel with no year is zeroed, not weighed - return (darray * weights).where(weights > 0, other=0).rename(darray.name) + # NB: arithmetic would align on an inner join, silently dropping pixels the grids don't share + join="exact", + output_dtypes=[numpy.result_type(darray.dtype, numpy.float32)], + ).rename(darray.name) def get_hectares_per_pixel(darray: xarray.DataArray) -> xarray.DataArray: