Add standard deviation to properties, update xclim.sdba.properties.tr…

…end documentation (#1445)

### What kind of change does this PR introduce?

* Add new function ``xclim.sdba.properties.std``to calculate the
standard deviation of a variable over all years at a given time
resolution.
* Add a test for the above.
* Amended the documentation of ``xclim.sdba.properties.trend`` to
document already existing functionality of calculating the return values
of scipy.stats.linregress.
* Modified the test for the above to test all 6 possible return values
of scipy.stats.linregress.

### Does this PR introduce a breaking change?

No.

.zenodo.json

@@ -120,6 +120,11 @@
     {
       "name": "Whelan, Chistopher",
       "affiliation": "Independent Researcher, United States"
+    },
+    {
+      "name": "Braun, Marco",
+      "affiliation": "Ouranos Consortium, Montréal, Québec, Canada",
+      "orcid": "0000-0001-5061-3217"
     }
   ],
   "keywords": [

AUTHORS.rst

@@ -25,6 +25,7 @@ Contributors
 * Clair Barnes <clair.barnes.16@ucl.ac.uk> `@clairbarnes <https://github.com/clairbarnes>`_
 * Alexis Beaupré-Laperrière `@Beauprel <https://github.com/Beauprel>`_
 * Sébastien Biner `@sbiner <https://github.com/sbiner>`_
+* Marco Braun <Braun.Marco@ouranos.ca> `@vindelico <https://github.com/vindelico>`_
 * David Caron `@davidcaron <https://github.com/davidcaron>`_
 * Carsten Ehbrecht <ehbrecht@dkrz.de> `@cehbrecht <https://github.com/cehbrecht>`_
 * Jeremy Fyke `@jeremyfyke <https://github.com/jeremyfyke>`_

CHANGES.rst

@@ -4,12 +4,14 @@ Changelog
 
 v0.45.0 (unreleased)
 --------------------
-Contributors to this version: David Huard (:user:`huard`), Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user: `juliettelavoie`), Gabriel Rondeau-Genesse (:user:`RondeauG`).
+Contributors to this version: David Huard (:user:`huard`), Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user: `juliettelavoie`), Gabriel Rondeau-Genesse (:user:`RondeauG`), Marco Braun (:user:`vindelico`).
 
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 * Added ``ensembles.hawkins_sutton`` method to partition the uncertainty sources in a climate projection ensemble. (:issue:`771`, :pull:`1262`).
 * New function ``xclim.core.calendar.convert_doy`` to transform day-of-year data between calendars. Also accessible from ``convert_calendar`` with ``doy=True``. (:issue:`1283`, :pull:`1406`).
+* Added new function ``xclim.sdba.properties.std`` to calculate the standard deviation of a variable over all years at a given time resolution. (:pull:`1445`).
+* Amended the documentation of ``xclim.sdba.properties.trend`` to document already existing functionality of calculating the return values of scipy.stats.linregress. (:pull:`1445`).
 * Add support for setting optional variables through the `ds` argument. (:issue:`1432`, :pull:`1435`).
 
 Bug fixes

tests/test_sdba/test_properties.py

@@ -44,6 +44,24 @@ def test_var(self, open_dataset):
         assert out_season.long_name.startswith("Variance")
         assert out_season.units == "kg^2 m-4 s-2"
 
+    def test_std(self, open_dataset):
+        sim = (
+            open_dataset("sdba/CanESM2_1950-2100.nc")
+            .sel(time=slice("1950", "1980"), location="Vancouver")
+            .pr
+        ).load()
+
+        out_year = sdba.properties.std(sim)
+        np.testing.assert_array_almost_equal(out_year.values, [5.08770208398345e-05])
+
+        out_season = sdba.properties.std(sim, group="time.season")
+        np.testing.assert_array_almost_equal(
+            out_season.values,
+            [6.2666411e-05, 3.5410259e-05, 4.3654352e-05, 5.3643853e-05],
+        )
+        assert out_season.long_name.startswith("Standard deviation")
+        assert out_season.units == "kg m-2 s-1"
+
     def test_skewness(self, open_dataset):
         sim = (
             open_dataset("sdba/CanESM2_1950-2100.nc")
@@ -352,19 +370,63 @@ def test_trend(self, open_dataset):
         ).load()
 
         slope = sdba.properties.trend(simt).values
+        intercept = sdba.properties.trend(simt, output="intercept").values
+        rvalue = sdba.properties.trend(simt, output="rvalue").values
         pvalue = sdba.properties.trend(simt, output="pvalue").values
+        stderr = sdba.properties.trend(simt, output="stderr").values
+        intercept_stderr = sdba.properties.trend(simt, output="intercept_stderr").values
+
         np.testing.assert_array_almost_equal(
-            [slope, pvalue], [-1.33720000e-01, 0.154605951862107], 4
+            [slope, intercept, rvalue, pvalue, stderr, intercept_stderr],
+            [
+                -0.133711111111111,
+                288.762132222222222,
+                -0.9706433333333333,
+                0.1546344444444444,
+                0.033135555555555,
+                0.042776666666666,
+            ],
+            4,
         )
 
         slope = sdba.properties.trend(simt, group="time.month").sel(month=1)
+        intercept = (
+            sdba.properties.trend(simt, output="intercept", group="time.month")
+            .sel(month=1)
+            .values
+        )
+        rvalue = (
+            sdba.properties.trend(simt, output="rvalue", group="time.month")
+            .sel(month=1)
+            .values
+        )
         pvalue = (
             sdba.properties.trend(simt, output="pvalue", group="time.month")
             .sel(month=1)
             .values
         )
+        stderr = (
+            sdba.properties.trend(simt, output="stderr", group="time.month")
+            .sel(month=1)
+            .values
+        )
+        intercept_stderr = (
+            sdba.properties.trend(simt, output="intercept_stderr", group="time.month")
+            .sel(month=1)
+            .values
+        )
+
         np.testing.assert_array_almost_equal(
-            [slope.values, pvalue], [0.8254349999999988, 0.6085783558202086], 4
+            [slope.values, intercept, rvalue, pvalue, stderr, intercept_stderr],
+            [
+                0.8254511111111111,
+                281.76353222222222,
+                0.576843333333333,
+                0.6085644444444444,
+                1.1689105555555555,
+                1.509056666666666,
+            ],
+            4,
         )
 
         assert slope.long_name.startswith("Slope of the interannual linear trend")

xclim/sdba/properties.py

@@ -176,6 +176,41 @@ def _var(da: xr.DataArray, *, group: str | Grouper = "time") -> xr.DataArray:
 )
 
 
+def _std(da: xr.DataArray, *, group: str | Grouper = "time") -> xr.DataArray:
+    """Standard Deviation.
+
+    Standard deviation of the variable over all years at the time resolution.
+
+    Parameters
+    ----------
+    da : xr.DataArray
+      Variable on which to calculate the diagnostic.
+    group : {'time', 'time.season', 'time.month'}
+      Grouping of the output.
+      E.g. If 'time.month', the standard deviation is performed separately for each month.
+
+    Returns
+    -------
+    xr.DataArray,
+      Standard deviation of the variable.
+    """
+    units = da.units
+    if group.prop != "group":
+        da = da.groupby(group.name)
+    out = da.std(dim=group.dim)
+    out.attrs["units"] = units
+    return out
+
+
+std = StatisticalProperty(
+    identifier="std",
+    aspect="marginal",
+    cell_methods="time: std",
+    compute=_std,
+    measure="xclim.sdba.measures.RATIO",
+)
+
+
 def _skewness(da: xr.DataArray, *, group: str | Grouper = "time") -> xr.DataArray:
     """Skewness.
 
@@ -823,11 +858,14 @@ def _trend(
     ----------
     da : xr.DataArray
         Variable on which to calculate the diagnostic.
-    output : {'slope', 'pvalue'}
-        Attributes of the linear regression to return.
+    output : {'slope', 'intercept', 'rvalue', 'pvalue', 'stderr', 'intercept_stderr'}
+        The attributes of the linear regression to return, as defined in scipy.stats.linregress:
         'slope' is the slope of the regression line.
-        'pvalue' is  for a hypothesis test whose null hypothesis is that the slope is zero,
-        using Wald Test with t-distribution of the test statistic.
+        'intercept' is the intercept of the regression line.
+        'rvalue' is the The Pearson correlation coefficient. The square of rvalue is equal to the coefficient of determination.
+        'pvalue' is the p-value for a hypothesis test whose null hypothesis is that the slope is zero, using Wald Test with t-distribution of the test statistic.
+        'stderr' is the standard error of the estimated slope (gradient), under the assumption of residual normality.
+        'intercept_stderr' is the standard error of the estimated intercept, under the assumption of residual normality.
     group : {'time', 'time.season', 'time.month'}
         Grouping on the output.