From 23ff158404c059924b68369a583daadbc36ea8c4 Mon Sep 17 00:00:00 2001 From: Ivan Ruiz Manuel <72193617+irm-codebase@users.noreply.github.com> Date: Fri, 13 Mar 2026 22:31:10 +0100 Subject: [PATCH 1/2] Add optional CH4 effects --- .../dynamic_characterization.py | 42 +++- dynamic_characterization/prospective/agwp.py | 185 +++++++++++---- .../prospective/radiative_forcing.py | 57 ++--- tests/test_prospective.py | 224 ++++++------------ 4 files changed, 267 insertions(+), 241 deletions(-) diff --git a/dynamic_characterization/dynamic_characterization.py b/dynamic_characterization/dynamic_characterization.py index ad690d1..29ec5f8 100644 --- a/dynamic_characterization/dynamic_characterization.py +++ b/dynamic_characterization/dynamic_characterization.py @@ -1,6 +1,5 @@ import json import os -import warnings from collections.abc import Collection from datetime import datetime from typing import Callable, Dict, Tuple @@ -40,6 +39,7 @@ def characterize( characterization_function_co2: Callable = None, time_varying_re: bool = False, fallback_to_ipcc: bool = True, + indirect_effects: str = "all", ) -> pd.DataFrame: """ Characterizes the dynamic inventory, formatted as a Dataframe, by evaluating each emission (row in DataFrame) using given dynamic characterization functions. @@ -88,6 +88,11 @@ def characterize( fallback_to_ipcc: bool, optional Only for prospective metrics. If True (default), use IPCC AR6 characterization functions for GHGs not available in the Watanabe module (e.g., CO and other GHGs). If False, only GHGs available in Watanabe (CO2, CH4, N2O) are characterized. + indirect_effects : str + Indirect effect treatment. One of: + - all: all indirect effects. Matches Watanabe et al. Default. + - no_carbon_cyle: all indirect effects except carbon cycle feedbacks. + - none: only direct effects of the gasae. Returns ------- @@ -163,6 +168,7 @@ def characterize( original_time_horizon=time_horizon, dynamic_time_horizon=dynamic_time_horizon, time_varying_re=time_varying_re, + indirect_effects=indirect_effects, ) ) @@ -184,6 +190,7 @@ def characterize( row=row, time_horizon=dynamic_time_horizon, time_varying_re=time_varying_re, + indirect_effects=indirect_effects, ) ) @@ -421,6 +428,7 @@ def _characterize_pgwp( original_time_horizon, dynamic_time_horizon, time_varying_re: bool = False, + indirect_effects: str = "all", ) -> CharacterizedRow: """ Calculate prospective GWP using Watanabe et al. (2026) characterization. @@ -428,16 +436,23 @@ def _characterize_pgwp( Uses AGWP_gas / AGWP_CO2 to calculate kg CO2 equivalent. Emission year is extracted from the row's date. """ - # Get emission year from the row date - emission_date = row.date - emission_year = int(str(emission_date.to_numpy())[:4]) # Calculate AGWP for the gas using its characterization function - radiative_forcing_ghg = characterization_functions[row.flow]( - row, - dynamic_time_horizon, - time_varying_re=time_varying_re, - ) + char_func = characterization_functions[row.flow] + if char_func == prospective_characterize_ch4: + # Only CH4 supports indirect_effects + radiative_forcing_ghg = char_func( + row, + dynamic_time_horizon, + time_varying_re=time_varying_re, + indirect_effects=indirect_effects, + ) + else: + radiative_forcing_ghg = char_func( + row, + dynamic_time_horizon, + time_varying_re=time_varying_re, + ) # Calculate reference AGWP for 1 kg of CO2 radiative_forcing_co2 = prospective_characterize_co2( @@ -533,6 +548,7 @@ def _characterize_prospective_radiative_forcing( row, time_horizon, time_varying_re: bool = False, + indirect_effects: str = "all", ) -> CharacterizedRow: """ Calculate prospective radiative forcing using Watanabe et al. (2026) characterization. @@ -550,6 +566,14 @@ def _characterize_prospective_radiative_forcing( prospective_characterize_n2o, ) + if char_func == prospective_characterize_ch4: + # CH4 supports multiple treatments of indirect effects + return char_func( + row, + time_horizon, + time_varying_re=time_varying_re, + indirect_effects=indirect_effects, + ) if char_func in prospective_functions: # Watanabe functions support time_varying_re return char_func(row, time_horizon, time_varying_re=time_varying_re) diff --git a/dynamic_characterization/prospective/agwp.py b/dynamic_characterization/prospective/agwp.py index 78f3e3c..da16e79 100644 --- a/dynamic_characterization/prospective/agwp.py +++ b/dynamic_characterization/prospective/agwp.py @@ -13,8 +13,8 @@ import numpy as np -from .config import get_scenario -from .data_loader import ( +from dynamic_characterization.prospective.config import get_scenario +from dynamic_characterization.prospective.data_loader import ( load_irf_ch4, load_irf_co2, load_irf_n2o, @@ -46,6 +46,134 @@ CONST_CH4 = 1.0 / PPB_TO_KG_CH4 # ppb/kg CONST_N2O = 1.0 / PPB_TO_KG_N2O # ppb/kg +CH4_INDIRECT_RE = 0.00014 + 0.00004 + + +def _check_indirect_effects_mode(mode: str) -> str: + """Validate the selected mode for indirect effects.""" + valid_modes = {"none", "no_carbon_cycle", "all"} + if mode not in valid_modes: + raise ValueError( + f"Invalid indirect effect mode '{mode}'. Valid modes: {sorted(valid_modes)}" + ) + return mode + + +def _get_re_value( + re_series: np.ndarray, + year_idx: int, + step: int, + time_varying_re: bool, +) -> float: + """Configuration-aware radiate efficiency getter.""" + if time_varying_re: + # RE evolves: use RE at emission_year + t + re_idx = min(year_idx + step, len(re_series) - 1) + re_t = re_series[re_idx] + else: + re_t = re_series[year_idx] + return re_t + + +def _agwp_ch4_cumulative( + emission_year: int, + time_horizon: int = 100, + time_varying_re: bool = False, + indirect_effects: str = "all", +) -> np.ndarray: + """ + Return cumulative CH4 AGWP values for each year of the selected horizon. + + Based on Watanabe's Code #4 CH4 AGWP function. + + Parameters + ---------- + emission_year : int + Year of emission (2030-2100, clamped if outside) + time_horizon : int + Integration period in years + time_varying_re : bool + If True, use RE that evolves over the decay period. + indirect_effects : str + Indirect effects to include. Defaults to "all". + + Returns + ------- + np.ndarray + Cumulative AGWP in W*yr/m^2/kg + """ + indirect_effects = _check_indirect_effects_mode(indirect_effects) + + scenario = get_scenario() + iam, ssp, rcp = scenario["iam"], scenario["ssp"], scenario["rcp"] + + # Load data + re_data_ch4 = load_re_ch4(iam) + re_data_co2 = load_re_co2(iam) + irf_series = load_irf_ch4() + + years = re_data_ch4["_years"] + re_series_ch4 = re_data_ch4[ssp][rcp] + re_series_co2 = re_data_co2[ssp][rcp] + + year_idx = _get_year_index(emission_year, years) + max_years = min(time_horizon, len(irf_series)) + + if max_years <= 0: + # Year 0 -> no effect + return np.array([0], dtype="float64") + + total_rf = np.zeros(max_years, dtype="float64") + + temp_s = np.zeros(max_years, dtype="float64") + temp_d = np.zeros(max_years, dtype="float64") + c1 = np.zeros(max_years, dtype="float64") + c2 = np.zeros(max_years, dtype="float64") + c3 = np.zeros(max_years, dtype="float64") + + for t in range(max_years): + re_ch4 = _get_re_value(re_series_ch4, year_idx, t, time_varying_re) + re_co2 = _get_re_value(re_series_co2, year_idx, t, time_varying_re) + + # Apply indirect effects if requested + re_ch4 = ( + re_ch4 + CH4_INDIRECT_RE + if indirect_effects in {"no_carbon_cycle", "all"} + else re_ch4 + ) + factor1 = re_ch4 * irf_series[t] * CONST_CH4 + factor1_model = factor1 * 1e12 + + factor2 = 0.0 + if indirect_effects == "all" and t > 0: + temp_s_prev = temp_s[t - 1] if t > 0 else 0 + temp_d_prev = temp_d[t - 1] if t > 0 else 0 + c1_prev = c1[t - 1] if t > 0 else 0 + + temp_s[t] = ( + temp_s_prev + + factor1_model / 7.7 + - 1.31 * temp_s_prev / 7.7 + - 0.88 * (temp_s_prev - temp_d_prev) / 7.7 + ) + temp_d[t] = temp_d_prev + (0.88 / 1.03) * (temp_s_prev - temp_d_prev) / 147 + + c1[t] = temp_s[t] * 11.06 * 0.6368 + c1_prev * np.exp(-1 / 2.376) + c2[t] = temp_s[t] * 11.06 * 0.3322 + c1_prev * np.exp(-1 / 30.14) + c3[t] = temp_s[t] * 11.06 * 0.031 + c1_prev * np.exp(-1 / 490.1) + + factor2 = (c1[t] + c2[t] + c3[t]) * re_co2 * CONST_CO2 + + total_rf[t] = factor1 + factor2 + + cumulative_agwp = np.zeros(max_years, dtype="float64") + for t in range(1, max_years): + cumulative_agwp[t] = ( + cumulative_agwp[t - 1] + (total_rf[t - 1] + total_rf[t]) / 2 + ) + + return cumulative_agwp + def _get_year_index(emission_year: int, years: np.ndarray) -> int: """ @@ -119,14 +247,7 @@ def agwp_co2( for t in range(max_years): irf_t = irf_series[t] - if time_varying_re: - # RE evolves: use RE at emission_year + t - re_idx = min(year_idx + t, len(re_series) - 1) - re_t = re_series[re_idx] - else: - # Fixed RE from emission year - re_t = re_series[year_idx] - + re_t = _get_re_value(re_series, year_idx, t, time_varying_re) # AGWP contribution: RE * IRF * conversion * dt (dt=1 year) agwp += re_t * irf_t * CONST_CO2 @@ -137,6 +258,7 @@ def agwp_ch4( emission_year: int, time_horizon: int = 100, time_varying_re: bool = False, + indirect_effects: str = "all", ) -> float: """ Calculate AGWP for 1 kg CH4. @@ -149,39 +271,21 @@ def agwp_ch4( Integration period in years time_varying_re : bool If True, use RE that evolves over the decay period. + indirect_effects : str + Indirect effects to include. Defaults to "all". Returns ------- float AGWP in W*yr/m^2/kg """ - scenario = get_scenario() - iam, ssp, rcp = scenario["iam"], scenario["ssp"], scenario["rcp"] - - # Load data - re_data = load_re_ch4(iam) - irf_series = load_irf_ch4() - - years = re_data["_years"] - re_series = re_data[ssp][rcp] # W/m^2/ppb - - year_idx = _get_year_index(emission_year, years) - - max_years = min(time_horizon, len(irf_series)) - - agwp = 0.0 - for t in range(max_years): - irf_t = irf_series[t] - - if time_varying_re: - re_idx = min(year_idx + t, len(re_series) - 1) - re_t = re_series[re_idx] - else: - re_t = re_series[year_idx] - - agwp += re_t * irf_t * CONST_CH4 - - return agwp + cumulative_agwp = _agwp_ch4_cumulative( + emission_year=emission_year, + time_horizon=time_horizon, + time_varying_re=time_varying_re, + indirect_effects=indirect_effects, + ) + return float(cumulative_agwp[-1]) def agwp_n2o( @@ -224,12 +328,7 @@ def agwp_n2o( for t in range(max_years): irf_t = irf_series[t] - if time_varying_re: - re_idx = min(year_idx + t, len(re_series) - 1) - re_t = re_series[re_idx] - else: - re_t = re_series[year_idx] - + re_t = _get_re_value(re_series, year_idx, t, time_varying_re) agwp += re_t * irf_t * CONST_N2O return agwp diff --git a/dynamic_characterization/prospective/radiative_forcing.py b/dynamic_characterization/prospective/radiative_forcing.py index ffc8efc..77fd7a2 100644 --- a/dynamic_characterization/prospective/radiative_forcing.py +++ b/dynamic_characterization/prospective/radiative_forcing.py @@ -13,8 +13,9 @@ from dynamic_characterization.classes import CharacterizedRow -from .config import get_scenario -from .data_loader import ( +from dynamic_characterization.prospective.agwp import _agwp_ch4_cumulative +from dynamic_characterization.prospective.config import get_scenario +from dynamic_characterization.prospective.data_loader import ( load_irf_ch4, load_irf_co2, load_irf_n2o, @@ -189,6 +190,7 @@ def characterize_ch4( period: int = 100, cumulative: bool = False, time_varying_re: bool = False, + indirect_effects: str = "all", ) -> CharacterizedRow: """ Calculate radiative forcing time series for 1 kg CH4 emission. @@ -196,8 +198,8 @@ def characterize_ch4( Uses scenario-based radiative efficiencies from Watanabe et al. (2026). Scenario must be set via prospective.set_scenario() before calling. - This includes only direct effects. For indirect effects (ozone, water vapor), - multiply the result by the appropriate factor (~1.43 for GWP100). + This can include either only direct effects, indirect effects w/o carbon feedback or all effects. + Only including direct effects will reduce values by ~1.43 for GWP100. Parameters ---------- @@ -211,6 +213,11 @@ def characterize_ch4( time_varying_re : bool If True, use RE that evolves over the decay period. If False, use fixed RE from emission year (IPCC standard, default). + indirect_effects : str + Treatment of indirect CH4 effects. One of: + - all: all indirect effects. Matches Watanabe et al. (2026). Default. + - no_carbon_cyle: all indirect effects except carbon cycle feedbacks. + - none: only direct effects of CH4. Returns ------- @@ -218,48 +225,20 @@ def characterize_ch4( namedtuple with date, amount, flow, activity arrays. Amount is in W*yr/m^2/kg CH4. """ - scenario = get_scenario() - iam, ssp, rcp = scenario["iam"], scenario["ssp"], scenario["rcp"] - - # Load data - re_data = load_re_ch4(iam) - irf_series = load_irf_ch4() - - years = re_data["_years"] - re_series = re_data[ssp][rcp] # W/m^2/ppb - - # Get emission year from series date date_beginning = series.date.to_numpy() emission_year = int(str(date_beginning)[:4]) - year_idx = _get_year_index(emission_year, years) - # Limit time horizon to available IRF data - max_years = min(period, len(irf_series)) + forcing = _agwp_ch4_cumulative( + emission_year=emission_year, + time_horizon=period, + time_varying_re=time_varying_re, + indirect_effects=indirect_effects, + ) - # Create date array dates_characterized = date_beginning + np.arange( - start=0, stop=max_years, dtype="timedelta64[Y]" + start=0, stop=len(forcing), dtype="timedelta64[Y]" ).astype("timedelta64[s]") - # Calculate cumulative radiative forcing at each time step - # Start from t=1 so forcing[0]=0 (no time elapsed = no forcing yet) - forcing = np.zeros(max_years, dtype="float64") - cumulative_forcing = 0.0 - - for t in range(1, max_years): - irf_t = irf_series[t] - - if time_varying_re: - re_idx = min(year_idx + t, len(re_series) - 1) - re_t = re_series[re_idx] - else: - re_t = re_series[year_idx] - - delta_forcing = re_t * irf_t * CONST_CH4 - cumulative_forcing += delta_forcing - forcing[t] = cumulative_forcing - - # Scale by emission amount forcing = forcing * series.amount if not cumulative: diff --git a/tests/test_prospective.py b/tests/test_prospective.py index 4a3eef1..4c44714 100644 --- a/tests/test_prospective.py +++ b/tests/test_prospective.py @@ -6,6 +6,7 @@ import os import sys import types +import pandas as pd # Set up fake parent package structure to enable relative imports in prospective modules _prospective_dir = os.path.join( @@ -204,31 +205,15 @@ def test_agwp_n2o_basic(): # --- pGWP100 Validation Tests Against SI Reference Tables --- -# + # IMPORTANT NOTE ON INDIRECT EFFECTS: # The Watanabe SI tables (es5c12391_si_003.xlsx for CH4, es5c12391_si_004.xlsx for N2O) # report pGWP100 values that INCLUDE indirect effects for CH4: # - Tropospheric ozone formation from CH4 oxidation # - Stratospheric water vapor from CH4 oxidation -# -# Per IPCC AR6 Chapter 7 Table 7.15: -# - Direct CH4 GWP100: ~19.3 -# - Including indirect effects: ~27.9 (factor of ~1.45) -# -# Our AGWP implementation calculates DIRECT radiative forcing only. -# Therefore, our pGWP100 = AGWP_CH4 / AGWP_CO2 will be ~30% lower than SI table values. -# -# The indirect effects factor can be estimated by comparing: -# SI_pGWP100 / calculated_direct_pGWP100 ~ 1.4-1.5 -# -# For N2O, indirect effects are minimal, so the match should be closer. - -import pandas as pd - -# Indirect effects multiplier for CH4 (based on IPCC AR6) -# The factor varies slightly by scenario/year, but ~1.43 is typical -CH4_INDIRECT_EFFECTS_FACTOR = 1.43 - +# Our AGWP implementation of CH4 can calculate both DIRECT and INDIRECT radiative forcing. +# Our AGWP implementation of NO2 calculates DIRECT radiative forcing only. +# For N2O, indirect effects are minimal, so the match should be close. def _find_si_row(df: pd.DataFrame, iam: str, ssp: str, rcp: str) -> pd.Series: """ @@ -250,7 +235,6 @@ def _find_si_row(df: pd.DataFrame, iam: str, ssp: str, rcp: str) -> pd.Series: pd.Series Matching row, or empty Series if not found """ - # Match based on components in the scenario string mask = ( df["Scenario"].str.contains(iam, case=False, na=False) & df["Scenario"].str.contains(ssp, case=False, na=False) @@ -272,106 +256,93 @@ def _find_si_row(df: pd.DataFrame, iam: str, ssp: str, rcp: str) -> pd.Series: ) -def test_pgwp100_ch4_direct_image_ssp1_26(): - """ - Direct pGWP100 for CH4 (without indirect effects) for IMAGE-SSP1-2.6. - Our AGWP implementation calculates direct radiative forcing only. - The SI table values INCLUDE indirect effects (~1.43x multiplier). - This test verifies: - 1. Direct pGWP is reasonable (between IPCC direct value ~19-22) - 2. When adjusted for indirect effects, matches SI table within 10% - """ +def _reference_agwp_co2(emission_year: int, time_horizon: int = 100) -> float: + """Reference CO2 AGWP using the trapezoidal integration from Watanabe SI Code #3.""" + scenario = config.get_scenario() + iam, ssp, rcp = scenario["iam"], scenario["ssp"], scenario["rcp"] + + re_data = data_loader.load_re_co2(iam) + years = re_data["_years"] + year_idx = np.searchsorted(years, emission_year) + re_value = re_data[ssp][rcp][year_idx] + + # Watanabe's SI code selects the CO2 IRF column by scenario position, not by RCP label. + # Reproducing the SI tables requires mirroring that behavior in the test reference. + scenario_index = list(re_data[ssp].keys()).index(rcp) + irf_df = pd.read_excel(os.path.join(_DATA_DIR, "es5c12391_si_010.xlsx")) + irf_values = irf_df.iloc[:time_horizon, 1 + scenario_index].to_numpy() + + rf_values = re_value * irf_values * agwp.CONST_CO2 + agwp_values = [(rf_values[0] + rf_values[1]) / 2] + for t in range(1, len(rf_values) - 1): + agwp_values.append(agwp_values[-1] + (rf_values[t] + rf_values[t + 1]) / 2) + agwp_values.append(agwp_values[-1] + rf_values[-1]) + return float(agwp_values[-1]) + +def _assert_pgwp100_ch4_matches_si(iam: str, ssp: str, rcp: str, year: int): + """Reusable helper to check pGWP values.""" ref_df = pd.read_excel(os.path.join(_DATA_DIR, "es5c12391_si_003.xlsx")) + row = _find_si_row(ref_df, iam, ssp, rcp) + assert not row.empty, f"{iam}-{ssp}-{rcp} scenario not found in SI table" - config.set_scenario(iam="IMAGE", ssp="SSP1", rcp="2.6") + config.set_scenario(iam=iam, ssp=ssp, rcp=rcp) - # Test year 2030 - agwp_ch4_val = agwp.agwp_ch4(emission_year=2030, time_horizon=100) - agwp_co2_val = agwp.agwp_co2(emission_year=2030, time_horizon=100) - calculated_direct_pgwp = agwp_ch4_val / agwp_co2_val + calculated_pgwp = agwp.agwp_ch4( + emission_year=year, time_horizon=100, indirect_effects="all" + ) / _reference_agwp_co2(emission_year=year, time_horizon=100) - row = _find_si_row(ref_df, "IMAGE", "SSP1", "2.6") - assert not row.empty, "IMAGE-SSP1-2.6 scenario not found in SI table" + assert calculated_pgwp == pytest.approx(row[year], abs=0.2) - si_value = row[2030] - # Verify direct pGWP is in expected range (~19-22) - assert 18 < calculated_direct_pgwp < 25, ( - f"Direct pGWP100 CH4 outside expected range: {calculated_direct_pgwp:.1f}" - ) +def test_agwp_ch4_indirect_effects_are_monotonic(): + """The three CH4 modes should monotonically add warming contributions.""" + config.set_scenario(iam="IMAGE", ssp="SSP1", rcp="2.6") - # Verify that applying indirect effects factor brings us close to SI value - adjusted_pgwp = calculated_direct_pgwp * CH4_INDIRECT_EFFECTS_FACTOR - assert adjusted_pgwp == pytest.approx(si_value, rel=0.10), ( - f"pGWP100 CH4 with indirect adjustment mismatch for IMAGE-SSP1-2.6 at 2030: " - f"got {adjusted_pgwp:.1f} (direct: {calculated_direct_pgwp:.1f} x {CH4_INDIRECT_EFFECTS_FACTOR}), " - f"expected {si_value:.1f}" + direct = agwp.agwp_ch4(emission_year=2030, time_horizon=100, indirect_effects="none") + no_carbon_cycle = agwp.agwp_ch4( + emission_year=2030, time_horizon=100, indirect_effects="no_carbon_cycle" ) + all_indirect = agwp.agwp_ch4(emission_year=2030, time_horizon=100, indirect_effects="all") + assert direct < no_carbon_cycle < all_indirect -def test_pgwp100_ch4_direct_aim_ssp3_60(): - """ - Direct pGWP100 for CH4 (without indirect effects) for AIM-SSP3-6.0. - """ - ref_df = pd.read_excel(os.path.join(_DATA_DIR, "es5c12391_si_003.xlsx")) - - config.set_scenario(iam="AIM", ssp="SSP3", rcp="6.0") - - agwp_ch4_val = agwp.agwp_ch4(emission_year=2030, time_horizon=100) - agwp_co2_val = agwp.agwp_co2(emission_year=2030, time_horizon=100) - calculated_direct_pgwp = agwp_ch4_val / agwp_co2_val - - row = _find_si_row(ref_df, "AIM", "SSP3", "6.0") - assert not row.empty, "AIM-SSP3-6.0 scenario not found in SI table" - - si_value = row[2030] - # Verify direct pGWP is in expected range - assert 18 < calculated_direct_pgwp < 25 +def test_pgwp100_ch4_image_ssp1_26_matches_si(): + """CH4 pGWP100 should match Watanabe Table S.2 for IMAGE-SSP1-2.6.""" + _assert_pgwp100_ch4_matches_si("IMAGE", "SSP1", "2.6", 2030) - # Verify that applying indirect effects factor brings us close to SI value - adjusted_pgwp = calculated_direct_pgwp * CH4_INDIRECT_EFFECTS_FACTOR - assert adjusted_pgwp == pytest.approx(si_value, rel=0.10), ( - f"pGWP100 CH4 with indirect adjustment mismatch for AIM-SSP3-6.0 at 2030: " - f"got {adjusted_pgwp:.1f}, expected {si_value:.1f}" - ) +def test_pgwp100_ch4_aim_ssp3_60_matches_si(): + """CH4 pGWP100 should match Watanabe Table S.2 for AIM-SSP3-6.0.""" + _assert_pgwp100_ch4_matches_si("AIM", "SSP3", "6.0", 2030) -def test_pgwp100_ch4_direct_multiple_years(): - """ - Direct pGWP100 for CH4 across multiple years (2030, 2040, 2050). - - Tests that our direct AGWP calculation is consistent across years. - Note: The indirect effects factor may vary slightly by year. - """ - ref_df = pd.read_excel(os.path.join(_DATA_DIR, "es5c12391_si_003.xlsx")) - - config.set_scenario(iam="IMAGE", ssp="SSP1", rcp="2.6") - - row = _find_si_row(ref_df, "IMAGE", "SSP1", "2.6") - assert not row.empty +def test_pgwp100_ch4_multiple_years_match_si(): + """CH4 pGWP100 should match Watanabe Table S.2 across pulse years.""" for year in [2030, 2040, 2050]: - agwp_ch4_val = agwp.agwp_ch4(emission_year=year, time_horizon=100) - agwp_co2_val = agwp.agwp_co2(emission_year=year, time_horizon=100) - calculated_direct_pgwp = agwp_ch4_val / agwp_co2_val - - si_value = row[year] + _assert_pgwp100_ch4_matches_si("IMAGE", "SSP1", "2.6", year) - # Direct pGWP should be in expected range - assert 18 < calculated_direct_pgwp < 28, ( - f"Direct pGWP100 CH4 at {year} outside expected range: {calculated_direct_pgwp:.1f}" - ) - # With indirect effects adjustment, should be within 15% of SI - # (allow more tolerance since indirect factor varies by year) - adjusted_pgwp = calculated_direct_pgwp * CH4_INDIRECT_EFFECTS_FACTOR - assert adjusted_pgwp == pytest.approx(si_value, rel=0.15), ( - f"pGWP100 CH4 with indirect adjustment at {year}: " - f"got {adjusted_pgwp:.1f}, expected {si_value:.1f}" - ) +@pytest.mark.parametrize( + "iam,ssp,rcp", + [ + ("IMAGE", "SSP1", "2.6"), + ("IMAGE", "SSP1", "4.5"), + ("AIM", "SSP3", "4.5"), + ("AIM", "SSP3", "6.0"), + ("GCAM4", "SSP4", "2.6"), + ("GCAM4", "SSP4", "4.5"), + ("MESSAGE", "SSP2", "4.5"), + ("MESSAGE", "SSP2", "6.0"), + ("REMIND", "SSP5", "4.5"), + ("REMIND", "SSP5", "8.5"), + ], +) +def test_pgwp100_ch4_all_scenarios_match_si(iam, ssp, rcp): + """CH4 pGWP100 should match Watanabe Table S.2 across scenarios.""" + _assert_pgwp100_ch4_matches_si(iam, ssp, rcp, 2030) def test_pgwp100_n2o_image_ssp1_26(): @@ -459,53 +430,6 @@ def test_pgwp100_n2o_multiple_years(): ) -@pytest.mark.parametrize( - "iam,ssp,rcp", - [ - ("IMAGE", "SSP1", "2.6"), - ("IMAGE", "SSP1", "4.5"), - ("AIM", "SSP3", "4.5"), - ("AIM", "SSP3", "6.0"), - ("GCAM4", "SSP4", "2.6"), - ("GCAM4", "SSP4", "4.5"), - ("MESSAGE", "SSP2", "4.5"), - ("MESSAGE", "SSP2", "6.0"), - ("REMIND", "SSP5", "4.5"), - ("REMIND", "SSP5", "8.5"), - ], -) -def test_pgwp100_ch4_direct_all_scenarios(iam, ssp, rcp): - """ - Direct pGWP100 for CH4 across all scenarios at year 2030. - - Tests that our direct AGWP implementation, when adjusted for indirect effects, - matches the SI table values within 15% tolerance. - """ - ref_df = pd.read_excel(os.path.join(_DATA_DIR, "es5c12391_si_003.xlsx")) - - config.set_scenario(iam=iam, ssp=ssp, rcp=rcp) - - agwp_ch4_val = agwp.agwp_ch4(emission_year=2030, time_horizon=100) - agwp_co2_val = agwp.agwp_co2(emission_year=2030, time_horizon=100) - calculated_direct_pgwp = agwp_ch4_val / agwp_co2_val - - row = _find_si_row(ref_df, iam, ssp, rcp) - assert not row.empty, f"{iam}-{ssp}-{rcp} scenario not found in SI table" - - si_value = row[2030] - - # Direct pGWP should be in expected range (~18-25) - assert 15 < calculated_direct_pgwp < 30, ( - f"Direct pGWP100 CH4 for {iam}-{ssp}-{rcp} outside expected range: {calculated_direct_pgwp:.1f}" - ) - - # With indirect effects adjustment, should be within 15% of SI - adjusted_pgwp = calculated_direct_pgwp * CH4_INDIRECT_EFFECTS_FACTOR - assert adjusted_pgwp == pytest.approx(si_value, rel=0.15), ( - f"pGWP100 CH4 with indirect adjustment for {iam}-{ssp}-{rcp} at 2030: " - f"got {adjusted_pgwp:.1f} (direct: {calculated_direct_pgwp:.1f}), expected {si_value:.1f}" - ) - @pytest.mark.parametrize( "iam,ssp,rcp", @@ -1167,8 +1091,8 @@ def test_pgwp_characterization(): # Allow slightly larger tolerance due to numerical differences in integration methods assert calculated_pgwp == pytest.approx(expected_pgwp, rel=0.1) - # pGWP for CH4 should be in reasonable range (direct effect ~21-25) - assert 15 < expected_pgwp < 40 + # pGWP for CH4 should be in reasonable range for the Watanabe default mode + assert 25 < expected_pgwp < 40 def test_pgtp_characterization(): From 611ef5d6efb75767c1cc2b664ceb48fe4565a0a9 Mon Sep 17 00:00:00 2001 From: TimoDiepers Date: Tue, 28 Apr 2026 10:21:29 +0200 Subject: [PATCH 2/2] require bw2data only when characterizing from method --- dynamic_characterization/dynamic_characterization.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/dynamic_characterization/dynamic_characterization.py b/dynamic_characterization/dynamic_characterization.py index 29ec5f8..71eb667 100644 --- a/dynamic_characterization/dynamic_characterization.py +++ b/dynamic_characterization/dynamic_characterization.py @@ -5,10 +5,8 @@ from typing import Callable, Dict, Tuple from loguru import logger -import bw2data as bd import numpy as np import pandas as pd -from bw2data.utils import UnknownObject from dynamic_characterization.classes import CharacterizedRow from dynamic_characterization.ipcc_ar6.radiative_forcing import ( @@ -251,6 +249,16 @@ def create_characterization_functions_from_method( characterization_functions = dict() + try: + import bw2data as bd + from bw2data.utils import UnknownObject + except ModuleNotFoundError as exc: + raise ModuleNotFoundError( + "bw2data is required to build default characterization functions from an LCIA method. " + "Install bw2data with its dependencies, or provide custom characterization_functions " + "when calling characterize()." + ) from exc + filepath = os.path.join( os.path.dirname(os.path.abspath(__file__)), "ipcc_ar6",