Snow storm days indices (and change snw_season_XYZ default thresholds)

CHANGES.rst

@@ -6,6 +6,15 @@ v0.46.0 (unreleased)
 --------------------
 Contributors to this version: Éric Dupuis (:user:`coxipi`), Trevor James Smith (:user:`Zeitsperre`), David Huard (:user:`huard`), Pascal Bourgault (:user:`aulemahal`).
 
+New indicators
+^^^^^^^^^^^^^^
+* ``xclim.indices.snw_storm_days`` computes the number of days with snowfall amount accumulation above a threshold. (:pull:`1505`).
+
+Breaking changes
+^^^^^^^^^^^^^^^^
+* ``xclim.indices.winter_storm`` renamed to ``xclim.indices.snd_storm_days``.  (:pull:`1505`).
+* Default threshold in ``xclim.indices.snw_season_{start|length|end}`` changed form `20 kg m-2` to `4 kg m-2`. (:pull:`1505`).
+
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 * Add ``wind_power_potential`` to estimate the potential for wind power production given wind speed at the turbine hub height and turbine specifications, along with  ``wind_profile`` to estimate the wind speed at different heights based on wind speed at a reference height. (:issue:`1458`, :pull:`1471`)

tests/test_indices.py

@@ -2713,7 +2713,7 @@ def test_continous_snow_season_start(self, snd_series, snw_series):
 
         out = xci.snw_season_start(snw)
         assert len(out) == 2
-        np.testing.assert_array_equal(out, [snw.time.dt.dayofyear[0].data + 2, np.nan])
+        np.testing.assert_array_equal(out, [snw.time.dt.dayofyear[0].data + 1, np.nan])
         for attr in ["units", "is_dayofyear", "calendar"]:
             assert attr in out.attrs.keys()
         assert out.attrs["units"] == ""

xclim/data/fr.json

@@ -1069,7 +1069,20 @@
     "long_name": "Jours avec forte accumulation de neige",
     "description": "Nombre de jours où l'accumulation de neige est au-dessus ou égale à {thresh}.",
     "title": "Jours avec accumulation de neige supérieure à un seuil",
-    "abstract": "Nombre de jours où l'accumulation de neige est au-dessus ou égale à un seuil donné."
+    "abstract": "Nombre de jours où l'accumulation de neige est au-dessus ou égale à un seuil donné.",
+    "_deprecated_version": "0.46.0"
+  },
+  "SND_STORM_DAYS": {
+    "long_name": "Jours avec forte accumulation de l'épaisseur de neige",
+    "description": "Nombre de jours où l'accumulation de l'épaisseur de neige est au-dessus ou égale à {thresh}.",
+    "title": "Jours avec accumulation de l'épaisseur de neige supérieure à un seuil",
+    "abstract": "Nombre de jours où l'accumulation de l'épaisseur neige est au-dessus ou égale à un seuil donné."
+  },
+  "SNW_STORM_DAYS": {
+    "long_name": "Jours avec forte accumulation de la quantité de neige",
+    "description": "Nombre de jours où l'accumulation de la quantité de neige est au-dessus ou égale à {thresh}.",
+    "title": "Jours avec accumulation de la quantité de neige supérieure à un seuil",
+    "abstract": "Nombre de jours où l'accumulation de la quantité de neige est au-dessus ou égale à un seuil donné."
   },
   "HIGH_PRECIP_LOW_TEMP": {
     "long_name": "Jours avec des précipitations au dessus de {pr_thresh} et des températures sous {tas_thresh}",

xclim/indicators/land/_snow.py

@@ -10,6 +10,7 @@
     "snd_season_end",
     "snd_season_length",
     "snd_season_start",
+    "snd_storm_days",
     "snd_to_snw",
     "snow_depth",
     "snow_melt_we_max",
@@ -18,6 +19,7 @@
     "snw_season_end",
     "snw_season_length",
     "snw_season_start",
+    "snw_storm_days",
     "snw_to_snd",
     "winter_storm",
 ]
@@ -171,7 +173,28 @@ class SnowWithIndexing(ResamplingIndicatorWithIndexing):
     long_name="Days with snowfall at or above a given threshold",
     description="The {freq} number of days with snowfall accumulation above {thresh}.",
     units="days",
-    compute=xci.winter_storm,
+    compute=xci.snd_storm_days,
+    _version_deprecated="0.46.0",
+)
+
+snd_storm_days = SnowWithIndexing(
+    title="Winter storm days",
+    identifier="snd_storm_days",
+    var_name="{freq}_snd_storm_days",
+    long_name="Days with snowfall depth at or above a given threshold",
+    description="The {freq} number of days with snowfall depth accumulation above {thresh}.",
+    units="days",
+    compute=xci.snd_storm_days,
+)
+
+snw_storm_days = SnowWithIndexing(
+    title="Winter storm days",
+    identifier="snw_storm_days",
+    var_name="{freq}_snw_storm_days",
+    long_name="Days with snowfall amount at or above a given threshold",
+    description="The {freq} number of days with snowfall amount accumulation above {thresh}.",
+    units="days",
+    compute=xci.snw_storm_days,
 )
 
 

xclim/indices/_threshold.py

@@ -80,11 +80,13 @@
     "snd_season_end",
     "snd_season_length",
     "snd_season_start",
+    "snd_storm_days",
     "snowfall_frequency",
     "snowfall_intensity",
     "snw_season_end",
     "snw_season_length",
     "snw_season_start",
+    "snw_storm_days",
     "tg_days_above",
     "tg_days_below",
     "tn_days_above",
@@ -399,13 +401,14 @@ def snd_season_end(
 @declare_units(snw="[mass]/[area]", thresh="[mass]/[area]")
 def snw_season_end(
     snw: xarray.DataArray,
-    thresh: Quantified = "20 kg m-2",
+    thresh: Quantified = "UNSET",
     window: int = 14,
     freq: str = "AS-JUL",
 ) -> xarray.DataArray:
     r"""End date of continuous snow water cover.
 
-    First day after the start of the continuous snow water cover when snow water is below a threshold (default: 2 cm)
+    First day after the start of the continuous snow water cover
+    when snow water is below a threshold (Current default:  20 kg m-2. xclim >=0.47.0 default: 4 kg m-2)
     for at least `N` (default: 14) consecutive days.
 
     Warnings
@@ -434,6 +437,11 @@ def snw_season_end(
     ----------
     :cite:cts:`chaumont_elaboration_2017`
     """
+    if thresh == "UNSET":
+        warnings.warn(
+            "The default value for this threshold will change in xclim>=0.47.0,  from `20 kg m-2` to `4 kg m-2`. Using `20 kg m-2` for now."
+        )
+        thresh = "20 kg m-2"
     valid = at_least_n_valid(snw.where(snw > 0), n=1, freq=freq)
 
     thresh = convert_units_to(thresh, snw)
@@ -507,13 +515,13 @@ def snd_season_start(
 @declare_units(snw="[mass]/[area]", thresh="[mass]/[area]")
 def snw_season_start(
     snw: xarray.DataArray,
-    thresh: Quantified = "20 kg m-2",
+    thresh: Quantified = "UNSET",
     window: int = 14,
     freq: str = "AS-JUL",
 ) -> xarray.DataArray:
     r"""Start date of continuous snow water cover.
 
-    Day of year when snow water is above or equal to a threshold (default: 2 cm)
+    Day of year when snow water is above or equal to a threshold  (Current default:  20 kg m-2. xclim >=0.47.0 default: 4 kg m-2)
     for at least `N` (default: 14) consecutive days.
 
     Warnings
@@ -541,6 +549,11 @@ def snw_season_start(
     ----------
     :cite:cts:`chaumont_elaboration_2017`
     """
+    if thresh == "UNSET":
+        warnings.warn(
+            "The default value for this threshold will change in xclim>=0.47.0,  from `20 kg m-2` to `4 kg m-2`. Using `20 kg m-2` for now."
+        )
+        thresh = "20 kg m-2"
     valid = at_least_n_valid(snw.where(snw > 0), n=1, freq=freq)
 
     thresh = convert_units_to(thresh, snw)
@@ -560,6 +573,90 @@ def snw_season_start(
     return out.where(~valid)
 
 
+@declare_units(snd="[length]", thresh="[length]")
+def snd_storm_days(
+    snd: xarray.DataArray, thresh: Quantified = "25 cm", freq: str = "AS-JUL"
+) -> xarray.DataArray:
+    """Days with snowfall over threshold.
+
+    Number of days with snowfall depth accumulation greater or equal to threshold (default: 25 cm).
+
+    Warnings
+    --------
+    The default `freq` is valid for the northern hemisphere.
+
+    Parameters
+    ----------
+    snd : xarray.DataArray
+        Surface snow depth.
+    thresh : Quantified
+        Threshold on snowfall depth accumulation require to label an event a `snd storm`.
+    freq : str
+        Resampling frequency.
+
+    Returns
+    -------
+    xarray.DataArray
+        Number of days per period identified as winter storms.
+
+    Notes
+    -----
+    Snowfall accumulation is estimated by the change in snow depth.
+    """
+    thresh = convert_units_to(thresh, snd)
+
+    # Compute daily accumulation
+    acc = snd.diff(dim="time")
+
+    # Winter storm condition
+    out = threshold_count(acc, ">=", thresh, freq)
+
+    out.attrs["units"] = to_agg_units(out, snd, "count")
+    return out
+
+
+@declare_units(snw="[mass]/[area]", thresh="[mass]/[area]")
+def snw_storm_days(
+    snw: xarray.DataArray, thresh: Quantified = "10 kg m-2", freq: str = "AS-JUL"
+) -> xarray.DataArray:
+    """Days with snowfall over threshold.
+
+    Number of days with snowfall amount accumulation greater or equal to threshold (default: 10 kg m-2).
+
+    Warnings
+    --------
+    The default `freq` is valid for the northern hemisphere.
+
+    Parameters
+    ----------
+    snw : xarray.DataArray
+        Surface snow amount.
+    thresh : Quantified
+        Threshold on snowfall amount accumulation require to label an event a `snw storm`.
+    freq : str
+        Resampling frequency.
+
+    Returns
+    -------
+    xarray.DataArray
+        Number of days per period identified as winter storms.
+
+    Notes
+    -----
+    Snowfall accumulation is estimated by the change in snow amount.
+    """
+    thresh = convert_units_to(thresh, snw)
+
+    # Compute daily accumulation
+    acc = snw.diff(dim="time")
+
+    # Winter storm condition
+    out = threshold_count(acc, ">=", thresh, freq)
+
+    out.attrs["units"] = to_agg_units(out, snw, "count")
+    return out
+
+
 @declare_units(pr="[precipitation]", thresh="[precipitation]")
 def daily_pr_intensity(
     pr: xarray.DataArray,
@@ -2078,13 +2175,13 @@ def snd_season_length(
 @declare_units(snw="[mass]/[area]", thresh="[mass]/[area]")
 def snw_season_length(
     snw: xarray.DataArray,
-    thresh: Quantified = "20 kg m-2",
+    thresh: Quantified = "UNSET",
     freq: str = "AS-JUL",
     op: str = ">=",
 ) -> xarray.DataArray:
     """The number of days with snow water above a threshold.
 
-    Number of days where surface snow water is greater or equal to given threshold (default: 2 cm).
+    Number of days where surface snow water is greater or equal to given threshold (Current default:  20 kg m-2. xclim >=0.47.0 default: 4 kg m-2).
 
     Warnings
     --------
@@ -2106,6 +2203,11 @@ def snw_season_length(
     xarray.DataArray, [time]
         Number of days where snow water is greater than or equal to threshold.
     """
+    if thresh == "UNSET":
+        warnings.warn(
+            "The default value for this threshold will change in xclim>=0.47.0,  from `20 kg m-2` to `4 kg m-2`. Using `20 kg m-2` for now."
+        )
+        thresh = "20 kg m-2"
     valid = at_least_n_valid(snw, n=1, freq=freq)
     thresh = convert_units_to(thresh, snw)
     out = threshold_count(snw, op, thresh, freq)
@@ -2986,6 +3088,8 @@ def winter_storm(
     Warnings
     --------
     The default `freq` is valid for the northern hemisphere.
+    The `winter_storm` indice is being deprecated in favour of `snd_storm_days`. This indice will
+    be removed in `xclim>=0.47.0`.
 
     Parameters
     ----------
@@ -3005,6 +3109,12 @@ def winter_storm(
     -----
     Snowfall accumulation is estimated by the change in snow depth.
     """
+    warnings.warn(
+        "The `winter_storm` indice is being deprecated in favour of `snd_storm_days`"
+        "This indice will be removed in `xclim>=0.47.0`. Please update your scripts accordingly.",
+        DeprecationWarning,
+        stacklevel=3,
+    )
     thresh = convert_units_to(thresh, snd)
 
     # Compute daily accumulation