add holiday_snow_days and holiday_snow_and_snowfall_days

Signed-off-by: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com>

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::
 
@@ -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, "mm", context="hydro")
+    snow_depth_thresh = convert_units_to(snd_thresh, "mm", context="hydro")
+    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.all(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 = "10 mm",
+    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 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: 10 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, [bool]
+        Boolean array of years with Perfect Christmas Days.
+
+    References
+    ----------
+    https://www.canada.ca/en/environment-climate-change/services/weather-general-tools-resources/historical-christmas-snowfall-data.html
+    """
+    snowfall_rate = rate2amount(
+        convert_units_to(prsn, "mm/d", 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),
+    )
+
+    snow_depth = convert_units_to(snd, "mm", context="hydro")
+    snow_depth_thresh = convert_units_to(snd_thresh, "mm", context="hydro")
+    snow_depth_constrained = select_time(
+        snow_depth,
+        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.all(dim="time"))
+    )
+
+    perfect_xmas_days = perfect_xmas_days.assign_attrs({"units": "1"})
+    return perfect_xmas_days