Add support for method argument in ensemble_percentile (#1732)

<!--Please ensure the PR fulfills the following requirements! -->
<!-- If this is your first PR, make sure to add your details to the
AUTHORS.rst! -->
### Pull Request Checklist:
- [x] This PR addresses an already opened issue (for bug fixes /
features)
    - This PR fixes #1694 
- [x] Tests for the changes have been added (for bug fixes / features)
- [ ] (If applicable) Documentation has been added / updated (for bug
fixes / features)
- [x] CHANGES.rst has been updated (with summary of main changes)
- [x] Link to issue (:issue:`number`) and pull request (:pull:`number`)
has been added

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

* Add method argument to ensemble_percentile. I chose to expose `method`
rather than `alpha`, `beta` as it felt more user friendly, since this is
what numpy and xarray use.

### Does this PR introduce a breaking change?
No

### Other information:

Raise an error if used with weights. 

Note that `test_calc_perc` takes 26 sec. on my laptop to run. I think
the test could be sped up by restricting computations to the things that
are ultimately checked.

CHANGES.rst

@@ -13,6 +13,9 @@ Announcements
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 * Indicator ``xclim.atmos.potential_evapotranspiration`` and indice ``xclim.indices.potential_evapotranspiration`` now accept a new value (`DA02`) for argument `method` implementing potential evapotranspiration based on Droogers and Allen (2002). (:issue:`1710`, :pull:`1723`).
+* `ensemble_percentiles` now takes a `method` argument, accepting
+  'interpolated_inverted_cdf', 'hazen', 'weibull', 'linear'  (default),
+  'median_unbiased' and 'normal_unbiased'. (:issue:`1694`, :pull:`1732`).
 
 New indicators
 ^^^^^^^^^^^^^^

tests/test_ensembles.py

@@ -176,6 +176,20 @@ def test_calc_perc(self, transpose, ensemble_dataset_objects, open_dataset):
             ),
             out1["tg_mean_p90"].isel(time=0, lon=5, lat=5),
         )
+
+        # Specify method
+        np.testing.assert_array_almost_equal(
+            mquantiles(
+                ens["tg_mean"].isel(time=0, lon=5, lat=5),
+                alphap=0.5,
+                betap=0.5,
+                prob=0.90,
+            ),
+            ensembles.ensemble_percentiles(
+                ens.isel(time=0, lon=5, lat=5), values=[90], method="hazen"
+            ).tg_mean_p90,
+        )
+
         assert np.all(out1["tg_mean_p90"] > out1["tg_mean_p50"])
         assert np.all(out1["tg_mean_p50"] > out1["tg_mean_p10"])
 

xclim/ensembles/_base.py

@@ -8,7 +8,7 @@
 from collections.abc import Sequence
 from glob import glob
 from pathlib import Path
-from typing import Any
+from typing import Any, Literal
 
 import numpy as np
 import xarray as xr
@@ -17,6 +17,16 @@
 from xclim.core.formatting import update_history
 from xclim.core.utils import calc_perc
 
+# The alpha and beta parameters for the quantile function
+_quantile_params = {
+    "interpolated_inverted_cdf": (0, 1),
+    "hazen": (0.5, 0.5),
+    "weibull": (0, 0),
+    "linear": (1, 1),
+    "median_unbiased": (1 / 3, 1 / 3),
+    "normal_unbiased": (3 / 8, 3 / 8),
+}
+
 
 def create_ensemble(
     datasets: Any,
@@ -207,6 +217,14 @@ def ensemble_percentiles(
     min_members: int | None = 1,
     weights: xr.DataArray | None = None,
     split: bool = True,
+    method: Literal[
+        "linear",
+        "interpolated_inverted_cdf",
+        "hazen",
+        "weibull",
+        "median_unbiased",
+        "normal_unbiased",
+    ] = "linear",
 ) -> xr.DataArray | xr.Dataset:
     """Calculate ensemble statistics between a results from an ensemble of climate simulations.
 
@@ -230,10 +248,13 @@ def ensemble_percentiles(
         The default (1) essentially skips this check.
     weights : xr.DataArray, optional
         Weights to apply along the 'realization' dimension. This array cannot contain missing values.
-        When given, the function uses xarray's quantile method which is slower than xclim's NaN-optimized algorithm.
+        When given, the function uses xarray's quantile method which is slower than xclim's NaN-optimized algorithm,
+        and does not support `method` values other than `linear`.
     split : bool
         Whether to split each percentile into a new variable
         or concatenate the output along a new "percentiles" dimension.
+    method : {"linear", "interpolated_inverted_cdf", "hazen", "weibull", "median_unbiased", "normal_unbiased"}
+        Method to use for estimating the percentile, see the `numpy.percentile` documentation for more information.
 
     Returns
     -------
@@ -274,6 +295,7 @@ def ensemble_percentiles(
                     split=split,
                     min_members=min_members,
                     weights=weights,
+                    method=method,
                 )
                 for da in ens.data_vars.values()
                 if "realization" in da.dims
@@ -310,20 +332,25 @@ def ensemble_percentiles(
             ens = ens.chunk({"realization": -1})
 
     if weights is None:
+        alpha, beta = _quantile_params[method]
+
         out = xr.apply_ufunc(
             calc_perc,
             ens,
             input_core_dims=[["realization"]],
             output_core_dims=[["percentiles"]],
             keep_attrs=True,
-            kwargs=dict(
-                percentiles=values,
-            ),
+            kwargs=dict(percentiles=values, alpha=alpha, beta=beta),
             dask="parallelized",
             output_dtypes=[ens.dtype],
             dask_gufunc_kwargs=dict(output_sizes={"percentiles": len(values)}),
         )
     else:
+        if method != "linear":
+            raise ValueError(
+                "Only the 'linear' method is supported when using weights."
+            )
+
         with xr.set_options(keep_attrs=True):
             # xclim's calc_perc does not support weighted arrays, so xarray's native function is used instead.
             qt = np.array(values) / 100  # xarray requires values between 0 and 1
@@ -350,7 +377,7 @@ def ensemble_percentiles(
             out = out.rename(name_dict={p: f"{ens.name}_p{int(p):02d}"})
 
     out.attrs["history"] = update_history(
-        f"Computation of the percentiles on {ens.realization.size} ensemble members.",
+        f"Computation of the percentiles on {ens.realization.size} ensemble members using method {method}.",
         ens,
     )