Add Christmas Day indices, add bivariate_count_occurrences

.github/workflows/main.yml

@@ -19,7 +19,7 @@ on:
       - submitted
 
 env:
-  XCLIM_TESTDATA_BRANCH: v2024.8.23
+  XCLIM_TESTDATA_BRANCH: add-snw-prsn
 
 concurrency:
   # For a given workflow, if we push to the same branch, cancel all previous builds on that branch except on main.

CHANGELOG.rst

@@ -2,6 +2,16 @@
 Changelog
 =========
 
+
+v0.55.0 (unreleased)
+--------------------
+Contributors to this version: Trevor James Smith (:user:`Zeitsperre`).
+
+New indicators
+^^^^^^^^^^^^^^
+* Added ``xclim.indices.holiday_snow_days`` to compute the number of days with snow on the ground during holidays ("Christmas Days"). (:issue:`2029`, :pull:`2030`).
+* Added ``xclim.indices.holiday_snow_and_snowfall_days`` to compute the number of days with snow on the ground and measurable snowfall during holidays ("Perfect Christmas Days"). (:issue:`2029`, :pull:`2030`).
+
 v0.54.0 (2024-12-16)
 --------------------
 Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Éric Dupuis (:user:`coxipi`), Sascha Hofmann (:user:`saschahofmann`).

src/xclim/data/fr.json

@@ -1033,6 +1033,18 @@
     "title": "Jours avec neige",
     "abstract": "Nombre de jours où la neige est entre une borne inférieure et supérieure."
   },
+  "HOLIDAY_SNOW_DAYS": {
+    "long_name": "Nombre de jours de neige durant les jours fériés",
+    "description": "Nombre de jours de neige durant les jours fériés.",
+    "title": "Jours de neige durant les jours fériés",
+    "abstract": "Nombre de jours de neige durant les jours fériés."
+  },
+  "HOLIDAY_SNOW_AND_SNOWFALL_DAYS": {
+    "long_name": "Nombre de jours de neige et de jours de chute de neige durant les jours fériés",
+    "description": "Nombre de jours de neige et de jours de chute de neige durant les jours fériés.",
+    "title": "Jours de neige et de chute de neige durant les jours fériés",
+    "abstract": "Nombre de jours de neige et de jours de chute de neige durant les jours fériés."
+  },
   "SND_SEASON_LENGTH": {
     "long_name": "Durée de couvert de neige",
     "description": "La saison débute lorsque l'épaisseur de neige est au-dessus de {thresh} durant {window} jours et se termine lorsqu'elle redescend sous {thresh} durant {window} jours.",

src/xclim/indicators/land/_snow.py

@@ -8,6 +8,8 @@
 
 __all__ = [
     "blowing_snow",
+    "holiday_snow_and_snowfall_days",
+    "holiday_snow_days",
     "snd_days_above",
     "snd_max_doy",
     "snd_season_end",
@@ -254,3 +256,26 @@ class SnowWithIndexing(ResamplingIndicatorWithIndexing):
     abstract="Number of days when the snow amount is greater than or equal to a given threshold.",
     compute=xci.snw_days_above,
 )
+
+holiday_snow_days = Snow(
+    title="Christmas snow days",
+    identifier="holiday_snow_days",
+    units="1",
+    long_name="Number of holiday days with snow",
+    description="The total number of days where snow on the ground was greater than or equal to {snd_thresh} "
+    "occurring on {date_start} and ending on {date_end}.",
+    abstract="The total number of days where there is a significant amount of snow on the ground on December 25th.",
+    compute=xci.holiday_snow_days,
+)
+
+holiday_snow_and_snowfall_days = Snow(
+    title="Perfect Christmas snow days",
+    identifier="holiday_snow_and_snowfall_days",
+    units="1",
+    long_name="Number of holiday days with snow and snowfall",
+    description="The total number of days where snow on the ground was greater than or equal to {snd_thresh} "
+    "and snowfall was greater than or equal to {prsn_thresh} occurring on {date_start} and ending on {date_end}.",
+    abstract="The total number of days where there is a significant amount of snow on the ground "
+    "and a measurable snowfall occurring on December 25th.",
+    compute=xci.holiday_snow_and_snowfall_days,
+)

src/xclim/indices/_threshold.py

@@ -9,7 +9,7 @@
 import xarray
 
 from xclim.core import DayOfYearStr, Quantified
-from xclim.core.calendar import doy_from_string, get_calendar
+from xclim.core.calendar import doy_from_string, get_calendar, select_time
 from xclim.core.missing import at_least_n_valid
 from xclim.core.units import (
     convert_units_to,
@@ -67,6 +67,8 @@
     "growing_season_start",
     "heat_wave_index",
     "heating_degree_days",
+    "holiday_snow_and_snowfall_days",
+    "holiday_snow_days",
     "hot_spell_frequency",
     "hot_spell_max_length",
     "hot_spell_max_magnitude",
@@ -376,7 +378,8 @@ def snd_season_end(
     window : int
         Minimum number of days with snow depth below threshold.
     freq : str
-        Resampling frequency. The default value is chosen for the northern hemisphere.
+        Resampling frequency. Default: "YS-JUL".
+        The default value is chosen for the northern hemisphere.
 
     Returns
     -------
@@ -2766,7 +2769,7 @@ def maximum_consecutive_frost_days(
     Let :math:`\mathbf{t}=t_0, t_1, \ldots, t_n` be a minimum daily temperature series and :math:`thresh` the threshold
     below which a day is considered a frost day. Let :math:`\mathbf{s}` be the sorted vector of indices :math:`i`
     where :math:`[t_i < thresh] \neq [t_{i+1} < thresh]`, that is, the days where the temperature crosses the threshold.
-    Then the maximum number of consecutive frost days is given by
+    Then the maximum number of consecutive frost days is given by:
 
     .. math::
 
@@ -2821,7 +2824,7 @@ def maximum_consecutive_dry_days(
     Let :math:`\mathbf{p}=p_0, p_1, \ldots, p_n` be a daily precipitation series and :math:`thresh` the threshold
     under which a day is considered dry. Then let :math:`\mathbf{s}` be the sorted vector of indices :math:`i` where
     :math:`[p_i < thresh] \neq [p_{i+1} < thresh]`, that is, the days where the precipitation crosses the threshold.
-    Then the maximum number of consecutive dry days is given by
+    Then the maximum number of consecutive dry days is given by:
 
     .. math::
 
@@ -3162,7 +3165,7 @@ def degree_days_exceedance_date(
     -----
     Let :math:`TG_{ij}` be the daily mean temperature at day :math:`i` of period :math:`j`,
     :math:`T` is the reference threshold and :math:`ST` is the sum threshold. Then, starting
-    at day :math:i_0:, the degree days exceedance date is the first day :math:`k` such that
+    at day :math:i_0:, the degree days exceedance date is the first day :math:`k` such that:
 
     .. math::
 
@@ -3633,3 +3636,142 @@ def wet_spell_max_length(
         resample_before_rl=resample_before_rl,
         **indexer,
     )
+
+
+@declare_units(
+    snd="[length]",
+    prsn="[precipitation]",
+    snd_thresh="[length]",
+    prsn_thresh="[length]",
+)
+def holiday_snow_days(
+    snd: xarray.DataArray,
+    snd_thresh: Quantified = "20 mm",
+    date_start: str = "12-25",
+    date_end: str | None = None,
+    freq: str = "YS",
+) -> xarray.DataArray:  # numpydoc ignore=SS05
+    r"""
+    Christmas Days.
+
+    Whether there is a significant amount of snow on the ground on December 25th (or around that time).
+
+    Parameters
+    ----------
+    snd : xarray.DataArray
+        Surface snow depth.
+    snd_thresh : Quantified
+        Threshold snow amount. Default: 20 mm.
+    date_start : str
+        Beginning of analysis period. Default: "12-25" (December 25th).
+    date_end : str, optional
+        End of analysis period. If not provided, `date_start` is used.
+        Default: None.
+    freq : str
+        Resampling frequency. Default: "YS".
+        The default value is chosen for the northern hemisphere.
+
+    Returns
+    -------
+    xarray.DataArray, [bool]
+        Boolean array of years with Christmas Days.
+
+    References
+    ----------
+    https://www.canada.ca/en/environment-climate-change/services/weather-general-tools-resources/historical-christmas-snowfall-data.html
+    """
+    snow_depth = convert_units_to(snd, "m")
+    snow_depth_thresh = convert_units_to(snd_thresh, "m")
+    snow_depth_constrained = select_time(
+        snow_depth,
+        drop=True,
+        date_bounds=(date_start, date_start if date_end is None else date_end),
+    )
+
+    xmas_days = (
+        (snow_depth_constrained >= snow_depth_thresh)
+        .resample(time=freq)
+        .map(lambda x: x.sum(dim="time"))
+    )
+
+    xmas_days = xmas_days.assign_attrs({"units": "1"})
+    return xmas_days
+
+
+@declare_units(
+    snd="[length]",
+    prsn="[precipitation]",
+    snd_thresh="[length]",
+    prsn_thresh="[length]",
+)
+def holiday_snow_and_snowfall_days(
+    snd: xarray.DataArray,
+    prsn: xarray.DataArray | None = None,
+    snd_thresh: Quantified = "20 mm",
+    prsn_thresh: Quantified = "1 cm",
+    date_start: str = "12-25",
+    date_end: str | None = None,
+    freq: str = "YS-JUL",
+) -> xarray.DataArray:
+    r"""
+    Perfect Christmas Days.
+
+    Whether there is a significant amount of snow on the ground and measurable snowfall occurring on December 25th.
+
+    Parameters
+    ----------
+    snd : xarray.DataArray
+        Surface snow depth.
+    prsn : xarray.DataArray
+        Snowfall flux.
+    snd_thresh : Quantified
+        Threshold snow amount. Default: 20 mm.
+    prsn_thresh : Quantified
+        Threshold snowfall flux. Default: 1 mm.
+    date_start : str
+        Beginning of analysis period. Default: "12-25" (December 25th).
+    date_end : str, optional
+        End of analysis period. If not provided, `date_start` is used.
+        Default: None.
+    freq : str
+        Resampling frequency. Default: "YS-JUL".
+        The default value is chosen for the northern hemisphere.
+
+    Returns
+    -------
+    xarray.DataArray, [int]
+        The total number of days with snow and snowfall during the holiday.
+
+    References
+    ----------
+    https://www.canada.ca/en/environment-climate-change/services/weather-general-tools-resources/historical-christmas-snowfall-data.html
+    """
+    snow_depth = convert_units_to(snd, "m")
+    snow_depth_thresh = convert_units_to(snd_thresh, "m")
+    snow_depth_constrained = select_time(
+        snow_depth,
+        drop=True,
+        date_bounds=(date_start, date_start if date_end is None else date_end),
+    )
+
+    snowfall_rate = rate2amount(
+        convert_units_to(prsn, "mm day-1", context="hydro"), out_units="mm"
+    )
+    snowfall_thresh = convert_units_to(prsn_thresh, "mm", context="hydro")
+    snowfall_constrained = select_time(
+        snowfall_rate,
+        drop=True,
+        date_bounds=(date_start, date_start if date_end is None else date_end),
+    )
+
+    perfect_xmas_days = (
+        (
+            (snow_depth_constrained >= snow_depth_thresh)
+            & (snowfall_constrained >= snowfall_thresh)
+        )
+        .resample(time=freq)
+        .map(lambda x: x.sum(dim="time"))
+    )
+
+    perfect_xmas_days = perfect_xmas_days.assign_attrs({"units": "1"})
+    return perfect_xmas_days

tests/test_indices.py

@@ -1468,6 +1468,64 @@ def test_1d(
         np.testing.assert_allclose(hwml.values, expected)
 
 
+class TestHolidayIndices:
+
+    def test_xmas_days_simple(self, snd_series):
+        # 5ish years of data
+        snd = snd_series(np.zeros(365 * 5), units="cm")
+
+        # add snow on ground on December 25 for first 3 years
+        snd.loc["2000-12-25"] = 2
+        snd.loc["2001-12-25"] = 1.5  # not enough
+        snd.loc["2002-12-25"] = 2
+        snd.loc["2003-12-25"] = 0  # no snow
+        snd.loc["2004-12-25"] = 6
+
+        out = xci.holiday_snow_days(snd)
+        np.testing.assert_array_equal(out, [1, 0, 1, 0, 1])
+
+    def test_xmas_days_range(self, snd_series):
+        # 5ish years of data
+        snd = snd_series(np.zeros(365 * 5), units="cm")
+
+        # add snow on ground on December 25 for first 3 years
+        snd.loc["2000-12-25"] = 2
+        snd.loc["2001-12-25"] = 1.5  # not enough
+        snd.loc["2002-12-24"] = 10  # a réveillon miracle
+        snd.loc["2002-12-25"] = 2
+        snd.loc["2003-12-25"] = 0  # no snow
+        snd.loc["2004-12-25"] = 6
+
+        out = xci.holiday_snow_days(snd, date_start="12-24", date_end="12-25")
+        np.testing.assert_array_equal(out, [1, 0, 2, 0, 1])
+
+    def test_perfect_xmas_days_simple(self, snd_series, prsn_series):
+        # 5ish years of data
+        a = np.zeros(365 * 5)
+        snd = snd_series(a, units="mm")
+        prsn = prsn_series(a.copy(), units="cm day-1")
+
+        # add snow on ground on December 25
+        snd.loc["2000-12-25"] = 20
+        snd.loc["2001-12-25"] = 15  # not enough
+        snd.loc["2001-12-26"] = 30  # too bad it's Boxing Day
+        snd.loc["2002-12-25"] = 20
+        snd.loc["2003-12-25"] = 0  # no snow
+        snd.loc["2004-12-25"] = 60
+
+        # add snowfall on December 25
+        prsn.loc["2000-12-25"] = 5
+        prsn.loc["2001-12-25"] = 2
+        prsn.loc["2001-12-26"] = 30  # too bad it's Boxing Day
+        prsn.loc["2002-12-25"] = 0.5  # not quite enough
+        prsn.loc["2003-12-25"] = 0  # no snow
+        prsn.loc["2004-12-25"] = 10
+        prsn = convert_units_to(prsn, "kg m-2 s-1", context="hydro")
+
+        out = xci.holiday_snow_and_snowfall_days(snd, prsn)
+        np.testing.assert_array_equal(out, [1, 0, 0, 0, 1])
+
+
 class TestHotSpellFrequency:
     @pytest.mark.parametrize(
         "thresh,window,op,expected",