Merge pull request #672 from cal-adapt/bug/wl-start-year-fix

Fixing start year WL slicing for LOCA data

Author: neilSchroeder

climakitae/core/constants.py

@@ -68,3 +68,10 @@
 ]
 
 CALISO_AREA_THRESHOLD = 100
+
+# Start and end years for LOCA and WRF data
+LOCA_START_YEAR = 1950
+LOCA_END_YEAR = 2100
+
+WRF_START_YEAR = 1981
+WRF_END_YEAR = 2099

climakitae/core/data_load.py

@@ -1297,7 +1297,7 @@ def read_catalog_from_select(selections: "DataParameters") -> xr.DataArray:
     """
 
     if selections.approach == "Warming Level":
-        selections.time_slice = (1980, 2100)  # Retrieve entire time period
+        selections.time_slice = (1950, 2100)  # Retrieve entire time period
 
     # Raise appropriate errors for time-based retrieval
     if selections.approach == "Time":

climakitae/new_core/processors/warming_level.py

@@ -20,7 +20,9 @@
     DataProcessor,
     register_processor,
 )
-from climakitae.util.utils import read_csv_file
+
+# from climakitae.new_core.processors.processor_utils import _determine_is_complete_wl
+from climakitae.util.utils import _determine_is_complete_wl, read_csv_file
 
 
 @register_processor("warming_level", priority=10)
@@ -50,9 +52,9 @@ class WarmingLevel(DataProcessor):
 
     Notes
     -----
-    The input data must span from 1980-2100 and include historical climate data
-    for proper warming level calculations. Data should have simulation and scenario
-    dimensions or be properly configured for stacking.
+    The input data must span from 1950-2100 or 1981-2100, depending on the downscaling method,
+    and include historical climate data for proper warming level calculations. Data should
+    have simulation and scenario dimensions or be properly configured for stacking.
     """
 
     def __init__(self, value: Dict[str, Any]):
@@ -134,7 +136,7 @@ def execute(
         # and split the data by member_id
         ret = self.reformat_member_ids(result)
 
-        # extend the time domain of all ssp scenarios to 1980-2100
+        # extend the time domain of all ssp scenarios to cover historical period
         ret = self.extend_time_domain(ret)
 
         # first, extract the member IDs from the data
@@ -151,6 +153,9 @@ def execute(
         # get center years for each key for each warming level
         center_years = self.get_center_years(member_ids, ret.keys())
         retkeys = list(ret.keys())
+
+        dropped_slices_count = 0
+
         for key in retkeys:
             if key not in center_years:
                 del ret[key]
@@ -169,9 +174,18 @@ def execute(
                     )
                     continue
                 start_year = year - self.warming_level_window
-                start_year = max(start_year, 1981)
                 end_year = year + self.warming_level_window - 1
-                end_year = min(end_year, 2100)
+
+                is_full_wl = _determine_is_complete_wl(
+                    start_year, end_year, key, context["activity_id"], wl
+                )
+
+                if not is_full_wl:
+                    warnings.warn(
+                        f"\n\nIncomplete warming level for {key} at {wl}C. "
+                        "\nSkipping this warming level."
+                    )
+                    continue
 
                 da_slice = ret[key].sel(time=slice(f"{start_year}", f"{end_year}"))
 
@@ -211,6 +225,8 @@ def execute(
                 slices, dim="warming_level", join="outer", fill_value=np.nan
             )
         self.update_context(context)
+
+        print(f"Dropped slices count: {dropped_slices_count}")
         return ret
 
     def update_context(self, context: Dict[str, Any]) -> Dict[str, Any]:
@@ -287,7 +303,8 @@ def extend_time_domain(
         self, result: Dict[str, Union[xr.Dataset, xr.DataArray]]
     ) -> Dict[str, Union[xr.Dataset, xr.DataArray]]:
         """
-        Extend the time domain of the input data to cover 1980-2100.
+        Extend the time domain of the input data to cover the historical period,
+        either 1950-2100 for LOCA or 1981-2100 for WRF.
 
         This method ensures that all SSP scenarios have historical data
         included in the time series, allowing for proper warming level calculations.

climakitae/util/utils.py

@@ -1,6 +1,7 @@
 import calendar
 import datetime
 import os
+import warnings
 from typing import Any, Iterable, Union
 
 import geopandas as gpd
@@ -13,7 +14,14 @@
 from shapely.geometry import Point, mapping
 from timezonefinder import TimezoneFinder
 
-from climakitae.core.constants import SSPS, UNSET
+from climakitae.core.constants import (
+    LOCA_END_YEAR,
+    LOCA_START_YEAR,
+    SSPS,
+    UNSET,
+    WRF_END_YEAR,
+    WRF_START_YEAR,
+)
 
 # from climakitae.core.data_interface import DataParameters
 from climakitae.core.paths import DATA_CATALOG_URL, STATIONS_CSV_PATH
@@ -1370,3 +1378,51 @@ def clip_gpd_to_shapefile(
     print(f"Number of stations within area: {len(clipped)}")
 
     return clipped
+
+
+def _determine_is_complete_wl(
+    start_year: int,
+    end_year: int,
+    simulation_name: str,
+    downscaling_method: str,
+    level: float,
+) -> bool:
+    """
+    Determine if a complete warming level slice can be created for the given start and end years.
+    This checks if the years fall within valid ranges based on the downscaling method.
+
+    Parameters
+    ----------
+    start_year : int
+        The starting year of the warming level slice.
+    end_year : int
+        The ending year of the warming level slice.
+    simulation_name : str
+        The name of the simulation being evaluated.
+    downscaling_method : str
+        The downscaling method used ("Statistical" or "Dynamical").
+    level : float
+        The warming level being evaluated.
+
+    Returns
+    -------
+    bool
+        True if a complete warming level slice can be created, False otherwise.
+    """
+    valid_years = {
+        "LOCA2": (LOCA_START_YEAR, LOCA_END_YEAR),
+        "Statistical": (LOCA_START_YEAR, LOCA_END_YEAR),
+        "WRF": (WRF_START_YEAR, WRF_END_YEAR),
+        "Dynamical": (WRF_START_YEAR, WRF_END_YEAR),
+    }
+
+    min_year, max_year = valid_years.get(downscaling_method, (None, None))
+
+    if min_year and (start_year < min_year or end_year > max_year):
+        warnings.warn(
+            f"\n\nIncomplete warming level for {simulation_name} at {level}C. "
+            "\nSkipping this warming level."
+        )
+        return False
+
+    return True

climakitae/util/warming_levels.py

@@ -1,6 +1,7 @@
 """Helper functions related to applying a warming levels approach to a data object"""
 
 import calendar
+import warnings
 from typing import Union
 
 import intake
@@ -21,6 +22,7 @@
     SSP585_FILE,
 )
 from climakitae.util.utils import (
+    _determine_is_complete_wl,
     _get_cat_subset,
     read_csv_file,
     resolution_to_gridlabel,
@@ -118,16 +120,18 @@ def _get_sliced_data(
     # Dropping leap days before slicing time dimension because the window size can affect number of leap days per slice
     y = y.loc[~((y.time.dt.month == 2) & (y.time.dt.day == 29))]
 
+    # Getting start and end years for slicing if `center_time` is not NaN
     if not pd.isna(center_time):
-
-        # Find the centered year
         centered_year = pd.to_datetime(center_time).year
         start_year = centered_year - window
         end_year = centered_year + (window - 1)
 
-        if start_year < 1981:
-            start_year = 1981
+    # If the centered year is not NaN and a complete warming level slice can be created, proceed with slicing
+    if not pd.isna(center_time) and _determine_is_complete_wl(
+        start_year, end_year, y.simulation.item(), y.downscaling_method, level
+    ):
 
+        # Slicing data around the centered year
         sliced = y.sel(time=slice(str(start_year), str(end_year)))
 
         # Creating a mask for timestamps that are within the desired months

tests/conftest.py

@@ -164,3 +164,37 @@ def test_dataarray_dict():
             xr_dict["WRF.UCLA.EC-Earth3.ssp370.day.d03"] = ds
 
     return xr_dict
+
+
+@pytest.fixture
+def test_dataarray_dict_loca():
+    """Create test datasets using xarray for warming_level.py tests."""
+    xr_dict = {}
+    member_id = ["r1i1p1f1"]
+    hist_periods, ssp_periods = 65, 86
+    hist_time = pd.date_range(
+        "1950-01-01", periods=hist_periods, freq="YS"
+    )  # yearly start
+    ssp_time = pd.date_range(
+        "2015-01-01", periods=ssp_periods, freq="YS"
+    )  # yearly start
+    y = [0]
+    x = [0]
+
+    for pair in zip([hist_periods, ssp_periods], [hist_time, ssp_time]):
+        periods, timestamps = pair
+        ds = xr.Dataset(
+            {"t2": (("member_id", "y", "x", "time"), np.zeros((1, 1, 1, periods)))},
+            coords={
+                "member_id": member_id,
+                "y": y,
+                "x": x,
+                "time": timestamps,
+            },
+        )
+        if periods == hist_periods:
+            xr_dict["LOCA2.UCLA.ACCESS-CM2.historical.day.d03"] = ds
+        else:
+            xr_dict["LOCA2.UCLA.ACCESS-CM2.ssp585.day.d03"] = ds
+
+    return xr_dict

tests/new_core/processors/test_warming_level.py

@@ -78,6 +78,41 @@ def processor():
     yield wl_processor
 
 
+@pytest.fixture
+def edge_case_processor():
+    """Fixture for a WarmingLevel processor with edge case warming level time DataFrames."""
+    wl_processor = WarmingLevel(value={"warming_levels": [0.5, 4.5]})
+    wl_time_data = [
+        [
+            "GCM",
+            "run",
+            "scenario",
+            "0.5",
+            "4.5",
+        ],
+        [
+            "ACCESS-CM2",
+            "r1i1p1f1",
+            "ssp585",
+            "1955-01-01 00:00:00",
+            "2095-01-01 00:00:00",
+        ],
+    ]
+    wl_processor.warming_level_times = pd.DataFrame(
+        wl_time_data[1:], columns=wl_time_data[0]
+    ).set_index(["GCM", "run", "scenario"])
+
+    wl_time_idx_data = [
+        ["time", "ACCESS-CM2_r3i1p1f1_ssp585"],
+        ["2005-1", 1.0],
+    ]
+    wl_processor.warming_level_times_idx = pd.DataFrame(
+        wl_time_idx_data[1:], columns=wl_time_idx_data[0]
+    ).set_index("time")
+
+    yield wl_processor
+
+
 @pytest.fixture
 def full_processor():
     """Fixture for a full WarmingLevel processor with default warming level time DataFrames."""
@@ -343,7 +378,7 @@ def test_empty_warming_level_times_error(self, mock_read_csv_file, processor):
     def test_execute_updates_context(self, request, full_processor):
         """Test that execute updates context with warming level information."""
         test_result = request.getfixturevalue("test_dataarray_dict")
-        context = {}
+        context = {"activity_id": "WRF"}
         _ = full_processor.execute(result=test_result, context=context)
         assert full_processor.name in context[_NEW_ATTRS_KEY][full_processor.name]
         assert str(full_processor.value) in context[_NEW_ATTRS_KEY][full_processor.name]
@@ -378,15 +413,15 @@ def test_execute_skips_warming_level(self, request, full_processor):
                 match=f"No warming level data found",
             ),
         ):
-            ret = full_processor.execute(data, context={})
+            ret = full_processor.execute(data, context={"activity_id": "WRF"})
             for key in ret:
                 assert len(ret[key].warming_level) == 2
                 assert 5.8 not in ret[key].warming_level.values
 
     def test_execute_dims_correct(self, request, full_processor):
         """Test that execute returns a dict with expected keys and types."""
         test_result = request.getfixturevalue("test_dataarray_dict")
-        ret = full_processor.execute(result=test_result, context={})
+        ret = full_processor.execute(result=test_result, context={"activity_id": "WRF"})
         for key in ret:
             assert isinstance(ret[key], xr.Dataset)
             assert "warming_level" in ret[key].dims
@@ -397,7 +432,7 @@ def test_execute_years_correct(self, request, full_processor):
         """Test that execute manipulates the data to have correct dims and years."""
         test_result = request.getfixturevalue("test_dataarray_dict")
         test_key = "WRF.UCLA.EC-Earth3.ssp370.day.d03"
-        ret = full_processor.execute(result=test_result, context={})
+        ret = full_processor.execute(result=test_result, context={"activity_id": "WRF"})
         ret_key = "WRF.UCLA.EC-Earth3.ssp370.day.d03.r1i1p1f1"
 
         # Check that the warming_level coordinate matches the processor's warming_levels
@@ -417,3 +452,50 @@ def test_execute_years_correct(self, request, full_processor):
         assert isinstance(ret[ret_key].centered_year.item(), int)
         # Check that the centered_year is within expected range
         assert 1981 <= ret[ret_key].centered_year.item() <= 2100
+
+    def test_execute_loca_correct(self, request, full_processor):
+        """Test that execute works correctly for LOCA data."""
+        test_result = request.getfixturevalue("test_dataarray_dict_loca")
+        # The test fixture creates data for ACCESS-CM2 model
+        test_key = "LOCA2.UCLA.ACCESS-CM2.ssp585.day.d03"
+        ret = full_processor.execute(
+            result=test_result, context={"activity_id": "LOCA2"}
+        )
+        ret_key = "LOCA2.UCLA.ACCESS-CM2.ssp585.day.d03.r1i1p1f1"
+
+        # Check that the warming_level coordinate matches the processor's warming_levels
+        assert (
+            ret[ret_key].warming_level.values == full_processor.warming_levels
+        ).all()
+        # Check the length of the time_delta dimension
+        first_year = str(test_result[test_key].isel(time=0).time.dt.year.item())
+        # Find the number of elements in the first year of `ret[key]`
+        timesteps_per_year = (
+            test_result[test_key].sel(time=slice(first_year, first_year)).time.size
+        )
+        assert (
+            len(ret[ret_key].time_delta)
+            == timesteps_per_year * full_processor.warming_level_window * 2
+        )
+        assert isinstance(ret[ret_key].centered_year.item(), int)
+        # Check that the centered_year is within expected range
+        assert 1981 <= ret[ret_key].centered_year.item() <= 2100
+
+    def test_execute_edge_case_years(self, request, edge_case_processor):
+        """Test that execute handles edge case warming levels and years correctly."""
+        test_result = request.getfixturevalue("test_dataarray_dict_loca")
+        with pytest.warns(
+            UserWarning,
+        ) as record:
+            ret = edge_case_processor.execute(
+                result=test_result, context={"activity_id": "LOCA2"}
+            )
+            assert len(record) == 5
+            assert ret == {}
+            messages = [str(w.message) for w in record]
+            expected_substrings = [
+                "Incomplete warming level for LOCA2.UCLA.ACCESS-CM2.ssp585.day.d03.r1i1p1f1 at 0.5C",
+                "Incomplete warming level for LOCA2.UCLA.ACCESS-CM2.ssp585.day.d03.r1i1p1f1 at 4.5C",
+            ]
+            for substr in expected_substrings:
+                assert any(substr in msg for msg in messages), f"Missing: {substr}"