Add Christmas Day indices, add bivariate_count_occurrences (#2030)

### What kind of change does this PR introduce? * Adds two new indices and indicators: * `holiday_snow_days` → units = "days" * `holiday_snow_and_snowfall_days` → units = "days" * Added `bivariate_count_occurrences` to generic indices. ### Does this PR introduce a breaking change? No. ### Other information: `prsn` and `snd` data has been added to the `xclim-testdata` repository: Ouranosinc/xclim-testdata#30
Ouranosinc · Jan 14, 2025 · 4cc8b2e · 4cc8b2e
2 parents eda9f51 + 93b5839
commit 4cc8b2e
Show file tree

Hide file tree

Showing 12 changed files with 451 additions and 43 deletions.
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
@@ -19,7 +19,7 @@ on:
       - submitted
 
 env:
-  XCLIM_TESTDATA_BRANCH: v2024.8.23
+  XCLIM_TESTDATA_BRANCH: v2025.1.8
 
 concurrency:
   # For a given workflow, if we push to the same branch, cancel all previous builds on that branch except on main.

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -4,11 +4,21 @@ Changelog
 
 v0.55.0 (unreleased)
 --------------------
-Contributors to this version: Juliette Lavoie (:user:`juliettelavoie`)
+Contributors to this version: Juliette Lavoie (:user:`juliettelavoie`), 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`).
 
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 * New function ``ensemble.partition.general_partition`` (:pull:`2035`)
+* Added a new ``xclim.indices.generic.bivariate_count_occurrences`` function to count instances where operations and performed and validated for two variables. (:pull:`2030`).
+
+Internal changes
+^^^^^^^^^^^^^^^^
+* `sphinx-codeautolink` and `pygments` have been temporarily pinned due to breaking API changes. (:pull:`2030`).
 
 v0.54.0 (2024-12-16)
 --------------------

diff --git a/environment.yml b/environment.yml
@@ -29,14 +29,14 @@ dependencies:
   - lmoments3 >=1.0.7  # Required for some Jupyter notebooks
   - pot >=0.9.4
   # Testing and development dependencies
-  - black ==24.10.0
-  - blackdoc ==0.3.9
+  - black =24.10.0
+  - blackdoc =0.3.9
   - bump-my-version >=0.28.1
   - cairosvg >=2.6.0
-  - codespell ==2.3.0
+  - codespell =2.3.0
   - coverage >=7.5.0
   - coveralls >=4.0.1  # Note: coveralls is not yet compatible with Python 3.13
-  - deptry ==0.20.0
+  - deptry =0.21.2
   - distributed >=2.0
   - flake8 >=7.1.1
   - flake8-rst-docstrings >=0.3.0
@@ -45,10 +45,10 @@ dependencies:
   - h5netcdf >=1.3.0
   - ipykernel
   - ipython >=8.5.0
-  - isort ==5.13.2
+  - isort =5.13.2
   - matplotlib >=3.6.0
   - mypy >=1.10.0
-  - nbconvert <7.14  # Pinned due to directive errors in sphinx. See: https://github.com/jupyter/nbconvert/issues/2092
+  - nbconvert >=7.16.4
   - nbqa >=1.8.2
   - nbsphinx >=0.9.5
   - nbval >=0.11.0
@@ -60,6 +60,7 @@ dependencies:
   - pooch >=1.8.0
   - pre-commit >=3.7
   - pybtex >=0.24.0
+  - pygments <2.19  # FIXME: temporary fix for sphinx-codeautolink
   - pylint >=3.3.1
   - pytest >=8.0.0
   - pytest-cov >=5.0.0
@@ -69,15 +70,15 @@ dependencies:
   - sphinx >=7.0.0
   - sphinx-autobuild >=2024.4.16
   - sphinx-autodoc-typehints
-  - sphinx-codeautolink
+  - sphinx-codeautolink >=0.15.2,!=0.16.0  # FIXME: temporary fix for sphinx-codeautolink
   - sphinx-copybutton
   - sphinx-mdinclude
   - sphinxcontrib-bibtex
   - sphinxcontrib-svg2pdfconverter
   - tokenize-rt >=5.2.0
   - tox >=4.21.2
   - tox-gh >=1.4.4
-  - vulture ==2.13
+  - vulture =2.14
   - xdoctest >=1.1.5
   - yamllint >=1.35.1
   - pip >=24.2.0

diff --git a/pyproject.toml b/pyproject.toml
@@ -72,7 +72,7 @@ dev = [
   "ipython >=8.5.0",
   "isort ==5.13.2",
   "mypy >=1.10.0",
-  "nbconvert <7.14", # Pinned due to directive errors in sphinx. See: https://github.com/jupyter/nbconvert/issues/2092
+  "nbconvert >=7.16.4",
   "nbqa >=1.8.2",
   "nbval >=0.11.0",
   "numpydoc >=1.8.0",
@@ -104,10 +104,11 @@ docs = [
   "nc-time-axis >=1.4.1",
   "pooch >=1.8.0",
   "pybtex >=0.24.0",
+  "pygments <2.19", # FIXME: temporary fix for sphinx-codeautolink
   "sphinx >=7.0.0",
   "sphinx-autobuild >=2024.4.16",
   "sphinx-autodoc-typehints",
-  "sphinx-codeautolink",
+  "sphinx-codeautolink >=0.15.2,!=0.16.0", # FIXME: temporary fix for sphinx-codeautolink
   "sphinx-copybutton",
   "sphinx-mdinclude",
   "sphinxcontrib-bibtex",

diff --git a/src/xclim/data/fr.json b/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 de Noël",
+    "description": "Nombre de jours de neige durant les jours de Noël.",
+    "title": "Jours de neige durant les jours de Noël",
+    "abstract": "Nombre de jours de neige durant les jours de Noël."
+  },
+  "HOLIDAY_SNOW_AND_SNOWFALL_DAYS": {
+    "long_name": "Nombre de jours de neige et de jours de chute de neige durant les jours de Noël",
+    "description": "Nombre de jours de neige et de jours de chute de neige durant les jours de Noël.",
+    "title": "Jours de neige et de chute de neige durant les jours de Noël",
+    "abstract": "Nombre de jours de neige et de jours de chute de neige durant les jours de Noël."
+  },
   "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.",

diff --git a/src/xclim/indicators/land/_snow.py b/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="days",
+    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="days",
+    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,
+)
diff --git a/src/xclim/indices/_threshold.py b/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,
@@ -22,7 +22,9 @@
 )
 from xclim.indices import run_length as rl
 from xclim.indices.generic import (
+    bivariate_count_occurrences,
     compare,
+    count_occurrences,
     cumulative_difference,
     domain_count,
     first_day_threshold_reached,
@@ -67,6 +69,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 +380,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 +2771,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 +2826,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 +3167,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::
 
@@ -3454,7 +3459,7 @@ def wet_spell_frequency(
         Resampling frequency.
     resample_before_rl : bool
         Determines if the resampling should take place before or after the run length encoding (or a similar algorithm) is applied to runs.
-    op : {"sum","min", "max", "mean"}
+    op : {"sum", "min", "max", "mean"}
         Operation to perform on the window.
         Default is "sum", which checks that the sum of accumulated precipitation over the whole window is more than the
         threshold.
@@ -3633,3 +3638,145 @@ 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",
+    op: str = ">=",
+    date_start: str = "12-25",
+    date_end: str | None = None,
+    freq: str = "YS",
+) -> xarray.DataArray:  # numpydoc ignore=SS05
+    """
+    Christmas Days.
+
+    Whether there is a significant amount of snow on the ground on December 25th (or a given date range).
+
+    Parameters
+    ----------
+    snd : xarray.DataArray
+        Surface snow depth.
+    snd_thresh : Quantified
+        Threshold snow amount. Default: 20 mm.
+    op : {">", "gt", ">=", "ge"}
+        Comparison operation. Default: ">=".
+    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
+    """
+    snd_constrained = select_time(
+        snd,
+        date_bounds=(date_start, date_start if date_end is None else date_end),
+    )
+
+    xmas_days = count_occurrences(
+        snd_constrained, snd_thresh, freq, op, constrain=[">=", ">"]
+    )
+
+    xmas_days = to_agg_units(xmas_days, snd, "count")
+    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 mm",
+    snd_op: str = ">=",
+    prsn_op: str = ">=",
+    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 daily snowfall liquid-water equivalent thickness. Default: 1 mm.
+    snd_op : {">", "gt", ">=", "ge"}
+        Comparison operation for snow depth. Default: ">=".
+    prsn_op : {">", "gt", ">=", "ge"}
+        Comparison operation for snowfall flux. Default: ">=".
+    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
+    """
+    snd_constrained = select_time(
+        snd,
+        date_bounds=(date_start, date_start if date_end is None else date_end),
+    )
+
+    prsn_mm = rate2amount(
+        convert_units_to(prsn, "mm day-1", context="hydro"), out_units="mm"
+    )
+    prsn_mm_constrained = select_time(
+        prsn_mm,
+        date_bounds=(date_start, date_start if date_end is None else date_end),
+    )
+
+    perfect_xmas_days = bivariate_count_occurrences(
+        data_var1=snd_constrained,
+        data_var2=prsn_mm_constrained,
+        threshold_var1=snd_thresh,
+        threshold_var2=prsn_thresh,
+        op_var1=snd_op,
+        op_var2=prsn_op,
+        freq=freq,
+        var_reducer="all",
+        constrain_var1=[">=", ">"],
+        constrain_var2=[">=", ">"],
+    )
+
+    perfect_xmas_days = to_agg_units(perfect_xmas_days, snd, "count")
+    return perfect_xmas_days