🔧 Fix climatology season and weekofyear (#703)

* add failing test

* working dirty version

* refactor

* fix climatology weekofyear

CHANGELOG.rst

@@ -9,7 +9,11 @@ climpred unreleased (202x-xx-xx)
 
 Bug Fixes
 ---------
-- Fix when creating ``valid_time`` from ``lead.attrs["units"]`` in ``["seasons", "years"]`` with multi-month stride in ``init``. (:issue:`698`, :pr:`700`) `Aaron Spring`_.
+- Fix when creating ``valid_time`` from ``lead.attrs["units"]`` in
+  ``["seasons", "years"]`` with multi-month stride in ``init``.
+  (:issue:`698`, :pr:`700`) `Aaron Spring`_.
+- Fix ``seasonality="season"`` in ``reference="climatology"``.
+  (:issue:`641`, :pr:`703`) `Aaron Spring`_.
 
 New Features
 ------------
@@ -39,8 +43,8 @@ New Features
 - Implement Logarithmic Ensemble Skill Score :py:func:`~climpred.metrics._less`.
   (:issue:`239`, :pr:`687`) `Aaron Spring`_.
 - :py:meth:`~climpred.classes.HindcastEnsemble.remove_seasonality` and
-  :py:meth:`~climpred.classes.PerfectModelEnsemble.remove_seasonality` remove the seasonality
-  of all ``climpred`` datasets. (:issue:`530`, :pr:`688`) `Aaron Spring`_.
+  :py:meth:`~climpred.classes.PerfectModelEnsemble.remove_seasonality` remove the
+  seasonality of all ``climpred`` datasets. (:issue:`530`, :pr:`688`) `Aaron Spring`_.
 - Add keyword ``groupby`` in :py:meth:`~climpred.classes.HindcastEnsemble.verify`,
   :py:meth:`~climpred.classes.PerfectModelEnsemble.verify`,
   :py:meth:`~climpred.classes.HindcastEnsemble.bootstrap` and
@@ -73,7 +77,11 @@ New Features
   aligned based on the `alignment <alignment.html>`_ keyword. This may help
   understanding which dates are matched for the different ``alignment`` approaches.
   (:issue:`701`, :pr:`702`) `Aaron Spring`_.
-- Add ``attrs`` to new ``coordinates`` created by ``climpred``. (:issue:`695`, :pr:`697`) `Aaron Spring`_.
+- Add ``attrs`` to new ``coordinates`` created by ``climpred``.
+  (:issue:`695`, :pr:`697`) `Aaron Spring`_.
+- Add ``seasonality="weekofyear"`` in ``reference="climatology"``.
+  (:pr:`703`) `Aaron Spring`_.
+
 
 Internals/Minor Fixes
 ---------------------

climpred/reference.py

@@ -1,5 +1,6 @@
 import inspect
 
+import pandas as pd
 import xarray as xr
 
 from .alignment import return_inits_and_verif_dates
@@ -29,32 +30,81 @@
 )
 
 
+def _maybe_seasons_to_int(ds):
+    """set season str values or coords to int"""
+    seasonal = False
+    for season in ["DJF", "MAM", "JJA", "SON"]:
+        if season in ds:
+            seasonal = True
+    if seasonal:
+        ds = (
+            ds.str.replace("DJF", "1")
+            .str.replace("MAM", "2")
+            .str.replace("JJA", "3")
+            .str.replace("SON", "4")
+            .astype("int")
+        )
+    elif "season" in ds.coords:  # set season coords to int
+        seasonal = False
+        for season in ["DJF", "MAM", "JJA", "SON"]:
+            if season in ds.coords["season"]:
+                seasonal = True
+        if seasonal:
+            ds.coords["season"] = (
+                ds.coords.get("season")
+                .str.replace("DJF", "1")
+                .str.replace("MAM", "2")
+                .str.replace("JJA", "3")
+                .str.replace("SON", "4")
+                .astype("int")
+            )
+    return ds
+
+
 def persistence(verif, inits, verif_dates, lead):
     lforecast = verif.where(verif.time.isin(inits[lead]), drop=True)
     lverif = verif.sel(time=verif_dates[lead])
     return lforecast, lverif
 
 
 def climatology(verif, inits, verif_dates, lead):
+    init_lead = inits[lead].copy()
     seasonality_str = OPTIONS["seasonality"]
     if seasonality_str == "weekofyear":
-        raise NotImplementedError
+        # convert to datetime for weekofyear operations
+        from .utils import convert_cftime_to_datetime_coords
+
+        verif = convert_cftime_to_datetime_coords(verif, "time")
+        init_lead["time"] = init_lead["time"].to_index().to_datetimeindex()
+        init_lead = init_lead["time"]
     climatology_day = verif.groupby(f"time.{seasonality_str}").mean()
     # enlarge times to get climatology_forecast times
     # this prevents errors if verification.time and hindcast.init are too much apart
     verif_hind_union = xr.DataArray(
-        verif.time.to_index().union(inits[lead].time.to_index()), dims="time"
+        verif.time.to_index().union(init_lead.time.to_index()), dims="time"
     )
 
-    climatology_forecast = climatology_day.sel(
-        {seasonality_str: getattr(verif_hind_union.time.dt, seasonality_str)},
-        method="nearest",
-    ).drop(seasonality_str)
-
+    climatology_forecast = (
+        _maybe_seasons_to_int(climatology_day)
+        .sel(
+            {
+                seasonality_str: _maybe_seasons_to_int(
+                    getattr(verif_hind_union.time.dt, seasonality_str)
+                )
+            },
+            method="nearest",  # nearest may be a bit incorrect but doesnt error
+        )
+        .drop(seasonality_str)
+    )
     lforecast = climatology_forecast.where(
-        climatology_forecast.time.isin(inits[lead]), drop=True
+        climatology_forecast.time.isin(init_lead), drop=True
     )
     lverif = verif.sel(time=verif_dates[lead])
+    # convert back to CFTimeIndex if needed
+    if isinstance(lforecast["time"].to_index(), pd.DatetimeIndex):
+        lforecast = convert_time_index(lforecast, "time")
+    if isinstance(lverif["time"].to_index(), pd.DatetimeIndex):
+        lverif = convert_time_index(lverif, "time")
     return lforecast, lverif
 
 

climpred/tests/test_reference.py

@@ -0,0 +1,18 @@
+import pytest
+
+from climpred import set_options
+
+
+@pytest.mark.parametrize("seasonality", ["month", "season", "dayofyear", "weekofyear"])
+@pytest.mark.parametrize("reference", ["persistence", "climatology", "uninitialized"])
+def test_HindcastEnsemble_verify_reference(
+    hindcast_hist_obs_1d, seasonality, reference
+):
+    with set_options(seasonality=seasonality):
+        hindcast_hist_obs_1d.verify(
+            metric="mse",
+            comparison="e2o",
+            dim="init",
+            alignment="same_verifs",
+            reference=reference,
+        )