From f6e05d18e54b539e991820c5616b0206b5584617 Mon Sep 17 00:00:00 2001 From: Aaron Spring Date: Mon, 6 Dec 2021 17:31:28 +0100 Subject: [PATCH] =?UTF-8?q?=F0=9F=94=A7=20Fix=20climatology=20`season`=20a?= =?UTF-8?q?nd=20`weekofyear`=20(#703)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * add failing test * working dirty version * refactor * fix climatology weekofyear --- CHANGELOG.rst | 16 ++++++-- climpred/reference.py | 66 ++++++++++++++++++++++++++++---- climpred/tests/test_reference.py | 18 +++++++++ 3 files changed, 88 insertions(+), 12 deletions(-) create mode 100644 climpred/tests/test_reference.py diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 5fdc2f662..8b6edb091 100644 --- a/CHANGELOG.rst +++ b/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 `_ 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 --------------------- diff --git a/climpred/reference.py b/climpred/reference.py index fe0133f71..7a528450e 100644 --- a/climpred/reference.py +++ b/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,6 +30,37 @@ ) +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]) @@ -36,25 +68,43 @@ def persistence(verif, inits, verif_dates, lead): 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 diff --git a/climpred/tests/test_reference.py b/climpred/tests/test_reference.py new file mode 100644 index 000000000..33112a8f2 --- /dev/null +++ b/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, + )