Merge branch 'main' into deprecate-sdba

CHANGELOG.rst

@@ -4,7 +4,7 @@ Changelog
 
 v0.55.0 (unreleased)
 --------------------
-Contributors to this version: Juliette Lavoie (:user:`juliettelavoie`), Trevor James Smith (:user:`Zeitsperre`), Sascha Hofmann (:user:`saschahofmann`), Pascal Bourgault (:user:`aulemahal`).
+Contributors to this version: Juliette Lavoie (:user:`juliettelavoie`), Trevor James Smith (:user:`Zeitsperre`), Sascha Hofmann (:user:`saschahofmann`), Pascal Bourgault (:user:`aulemahal`), Éric Dupuis (:user:`coxipi`).
 
 Breaking changes
 ^^^^^^^^^^^^^^^^
@@ -20,15 +20,17 @@ 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`).
+* Added ``xclim.indices.vapor_pressure_deficit`` to compute the vapor pressure deficit from temperature and relative humidity. (:issue:`1917`, :pull:`2072`).
 
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-* New function ``ensemble.partition.general_partition`` (:pull:`2035`)
+* 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`).
 * `xclim` now tracks energy usage and carbon emissions ("last run", "average", and "total") during CI workflows using the `eco-ci-energy-estimation` GitHub Action. (:pull:`2046`).
 * ``xclim.testing.helpers.test_timeseries`` now accepts a `calendar` argument that is forwarded to ``xr.cftime_range``. (:pull:`2019`).
 * New ``xclim.indices.fao_allen98``, exporting the FAO-56 Penman-Monteith equation for potential evapotranspiration (:issue:`2004`, :pull:`2067`).
 * Missing values method "pct" and "at_least_n" now accept a new "subfreq" option that allows to compute the missing mask in two steps. When given, the algorithm is applied at this "subfreq" resampling frequency first and then the result is resampled at the target indicator frequency. In the output, a period is invalid if any of its subgroup where flagged as invalid by the chosen method. (:pull:`2058`, :issue:`1820`).
+* `rv_continuous` functions can now be given directly as the `dist` argument in ``standardized_precipitation_index`` and ``standardized_precipitation_evapotranspiration_index`` indicators. This includes `lmoments3` distribution when specifying `method="PWM"`. (:issue:`2043`, :pull:`2045`).
 
 Internal changes
 ^^^^^^^^^^^^^^^^
@@ -39,6 +41,7 @@ Internal changes
     - Objects are initialized with their options and then called with the data, input frequency, target frequency and indexer.
     - Subclasses receive non-resampled DataArray in their ``is_missing`` methods.
     - ``MissingWMO`` now uses ``xclim.indices.helpers.resample_map`` which should greatly improve performance in a dask context.
+* There is now a warning stating that `fitkwargs` are not employed when using the `lmoments3` distribution. One exception is the use of `'floc'` which is allowed with the gamma distribution. `'floc'` is used to shift the distribution before computing fitting parameters with the `lmoments3` distribution since ``loc=0`` is always assumed in the library. (:issue:`2043`, :pull:`2045`).
 * `xclim.sdba` now emits a `FutureWarning` on import to inform users that the submodule is being deprecated. (:pull:`2073`).
 
 Bug fixes

docs/notebooks/customize.ipynb

@@ -149,14 +149,23 @@
    "source": [
     "This method checks that there are less than `nm=11` invalid values in a month and that there are no consecutive runs of `nc>=5` invalid values. Thus, every month is now valid.\n",
     "\n",
-    "Finally, it is possible for advanced users to register their own method. Xclim's missing methods are in fact based on class instances. Thus, to create a custom missing class, one should implement a subclass based on `xclim.core.checks.MissingBase` and overriding at least the `is_missing` method. The method should take a `null` argument and  a `count` argument.\n",
+    "<div class=\"alert alert-warning\">\n",
     "\n",
-    "- `null` is a `DataArrayResample` instance of the resampled mask of invalid values in the input data array.\n",
-    "- `count` is the number of days in each resampled periods and any number of other keyword arguments.\n",
+    "The content that follows is based on an experimental part of xclim, the way missing methods are implemented might change in the near future. Access of missing methods as functions of ``xclim.core.missing`` and usage of these algorithms in the indicators will be preserved, but custom subclasses might break with future changes.\n",
     "\n",
-    "The `is_missing` method should return a boolean mask, at the same frequency as the indicator output (same as `count`), where True values are for elements that are considered missing and masked on the output.\n",
+    "</div>\n",
     "\n",
-    "When registering the class with the `xclim.core.checks.register_missing_method` decorator, the keyword arguments will be registered as options for the missing method. One can also implement a `validate` static method that receives only those options and returns whether they should be considered valid or not."
+    "Finally, it is possible for advanced users to register their own methods. Xclim's missing methods are in fact class-based. To create a custom missing class, one should implement a subclass of `xclim.core.checks.MissingBase` and override at least the `is_missing` method. This method should take the following arguments:\n",
+    "\n",
+    "- `null`, a `DataArray` of the mask of invalid values in the input data array (with the same time coordinate as the raw data).\n",
+    "- `count`, `DataArray` of the number of days in each resampled periods\n",
+    "- `freq`, the resampling frequency.\n",
+    "\n",
+    "The `is_missing` method should return a boolean mask, resampled at the `freq` frequency, the same as the indicator output (same as `count`), where `True` values are for elements that are considered missing and masked on the output.\n",
+    "\n",
+    "To add additional arguments, one should override the `__init__` (receiving those arguments) and the `validate` static method, which validates them. The options are then stored in the `options` property of the instance. See example below and the docstrings in the module.\n",
+    "\n",
+    "When registering the class with the `xclim.core.checks.register_missing_method` decorator, the keyword arguments will be registered as options for the missing method. "
    ]
   },
   {
@@ -173,9 +182,15 @@
     "class MissingConsecutive(MissingBase):\n",
     "    \"\"\"Any period with more than max_n consecutive missing values is considered invalid\"\"\"\n",
     "\n",
-    "    def is_missing(self, null, count, max_n=5):\n",
+    "    def __init__(self, max_n: int = 5):\n",
+    "        super().__init__(max_n=max_n)\n",
+    "\n",
+    "    def is_missing(self, null, count, freq):\n",
     "        \"\"\"Return a boolean mask where True values are for elements that are considered missing and masked on the output.\"\"\"\n",
-    "        return null.map(longest_run, dim=\"time\") >= max_n\n",
+    "        return (\n",
+    "            null.resample(time=freq).map(longest_run, dim=\"time\")\n",
+    "            >= self.options[\"max_n\"]\n",
+    "        )\n",
     "\n",
     "    @staticmethod\n",
     "    def validate(max_n):\n",
@@ -204,6 +219,13 @@
     "    )  # compute monthly max tasmax\n",
     "tx_mean.sel(time=\"2013\", lat=75, lon=200)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -217,7 +239,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.9"
+   "version": "3.12.7"
   }
  },
  "nbformat": 4,

pyproject.toml

@@ -137,7 +137,7 @@ target-version = [
 ]
 
 [tool.bumpversion]
-current_version = "0.54.1-dev.11"
+current_version = "0.54.1-dev.13"
 commit = true
 commit_args = "--no-verify --signoff"
 tag = false

src/xclim/__init__.py

@@ -13,7 +13,7 @@
 
 __author__ = """Travis Logan"""
 __email__ = "[email protected]"
-__version__ = "0.54.1-dev.11"
+__version__ = "0.54.1-dev.13"
 
 
 with _resources.as_file(_resources.files("xclim.data")) as _module_data:

src/xclim/data/fr.json

@@ -931,6 +931,12 @@
     "title": "Humidité spécifique calculée à partir de la température du point de rosée et de la pression",
     "abstract": "Humidité spécifique calculée à partir de la température du point de rosée et de la pression à l'aide de la pression de vapeur saturante."
   },
+  "VAPOR_PRESSURE_DEFICIT": {
+    "long_name": "Déficit de pression de vapeur (méthode \"{method}\")",
+    "description": "Déficit de pression de vapeur calculé à partir de la température et de l'humidité relative à l'aide de la pression de vapeur saturante, laquelle fut calculée en suivant la méthode {method}.",
+    "title": "Déficit de pression de vapeur calculé à partir de la température et de l'humidité relative",
+    "abstract": "Déficit de pression de vapeur calculé à partir de la température et de l'humidité relative à l'aide de la pression de vapeur saturante."
+  },
   "FIRST_DAY_TG_BELOW": {
     "long_name": "Premier jour de l'année avec une température moyenne quotidienne sous {thresh} durant au moins {window} jours",
     "description": "Premier jour de l'année avec une température moyenne quotidienne sous {thresh} durant au moins {window} jours.",

src/xclim/indicators/atmos/_conversion.py

@@ -24,6 +24,7 @@
     "specific_humidity_from_dewpoint",
     "tg",
     "universal_thermal_climate_index",
+    "vapor_pressure_deficit",
     "water_budget",
     "water_budget_from_tas",
     "wind_chill_index",
@@ -270,6 +271,26 @@ def cfcheck(self, **das) -> None:
     compute=indices.specific_humidity_from_dewpoint,
 )
 
+
+vapor_pressure_deficit = Converter(
+    title="Water vapour pressure deficit",
+    identifier="vapor_pressure_deficit",
+    units="Pa",
+    long_name='Vapour pressure deficit ("{method}" method)',
+    standard_name="water_vapor_saturation_deficit_in_air",
+    description=lambda **kws: (
+        "The difference between the saturation vapour pressure and the actual vapour pressure,"
+        "calculated from temperature and relative humidity according to the {method} method."
+    )
+    + (
+        " The computation was done in reference to ice for temperatures below {ice_thresh}."
+        if kws["ice_thresh"] is not None
+        else ""
+    ),
+    abstract="Difference between the saturation vapour pressure and the actual vapour pressure.",
+    compute=indices.vapor_pressure_deficit,
+)
+
 snowfall_approximation = Converter(
     title="Snowfall approximation",
     identifier="prsn",

src/xclim/indices/_agro.py

@@ -6,6 +6,7 @@
 
 import numpy as np
 import xarray
+from scipy.stats import rv_continuous
 
 import xclim.indices.run_length as rl
 from xclim.core import DateStr, DayOfYearStr, Quantified
@@ -1122,7 +1123,7 @@ def standardized_precipitation_index(
     pr: xarray.DataArray,
     freq: str | None = "MS",
     window: int = 1,
-    dist: str = "gamma",
+    dist: str | rv_continuous = "gamma",
     method: str = "ML",
     fitkwargs: dict | None = None,
     cal_start: DateStr | None = None,
@@ -1143,13 +1144,14 @@ def standardized_precipitation_index(
     window : int
         Averaging window length relative to the resampling frequency. For example, if `freq="MS"`,
         i.e. a monthly resampling, the window is an integer number of months.
-    dist : {"gamma", "fisk"}
-        Name of the univariate distribution. (see :py:mod:`scipy.stats`).
-    method : {"APP", "ML"}
+    dist : {'gamma', 'fisk'} or `rv_continuous` function
+        Name of the univariate distribution, or a callable `rv_continuous` (see :py:mod:`scipy.stats`).
+    method : {"APP", "ML", "PWM"}
         Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method
-        uses a deterministic function that does not involve any optimization.
+        uses a deterministic function that does not involve any optimization. `PWM` should be used with a `lmoments3` distribution.
     fitkwargs : dict, optional
         Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`).
+        If method is `PWM`, `fitkwargs` should be empty, except for `floc` with `dist`=`gamma` which is allowed.
     cal_start : DateStr, optional
         Start date of the calibration period. A `DateStr` is expected, that is a `str` in format `"YYYY-MM-DD"`.
         Default option `None` means that the calibration period begins at the start of the input dataset.
@@ -1169,18 +1171,22 @@ def standardized_precipitation_index(
     xarray.DataArray, [unitless]
         Standardized Precipitation Index.
 
+    See Also
+    --------
+    xclim.indices.stats.standardized_index : Standardized Index.
+    xclim.indices.stats.standardized_index_fit_params : Standardized Index Fit Params.
+
     Notes
     -----
     * N-month SPI / N-day SPI is determined by choosing the `window = N` and the appropriate frequency `freq`.
     * Supported statistical distributions are: ["gamma", "fisk"], where "fisk" is scipy's implementation of
-       a log-logistic distribution
+      a log-logistic distribution
     * Supported frequencies are daily ("D"), weekly ("W"), and monthly ("MS").
-      Weekly frequency will only work if the input array has a "standard" (non-cftime) calendar.
+    * Weekly frequency will only work if the input array has a "standard" (non-cftime) calendar.
     * If `params` is given as input, it overrides the `cal_start`, `cal_end`, `freq` and `window`, `dist` and `method` options.
     * "APP" method only supports two-parameter distributions. Parameter `loc` needs to be fixed to use method `APP`.
-    * The standardized index is bounded by ±8.21. 8.21 is the largest standardized index as constrained by the float64 precision in
-      the inversion to the normal distribution.
-    * The results from `climate_indices` library can be reproduced with `method = "APP"` and `fitwkargs = {"floc": 0}`
+    * The results from `climate_indices` library can be reproduced with `method = "APP"` and `fitwkargs = {"floc": 0}`, except for the maximum
+      and minimum values allowed which are greater in xclim ±8.21, . See `xclim.indices.stats.standardized_index`
 
     References
     ----------
@@ -1218,14 +1224,16 @@ def standardized_precipitation_index(
     >>> spi_3_fitted = standardized_precipitation_index(pr, params=params)
     """
     fitkwargs = fitkwargs or {}
+
     dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]}
-    if dist in dist_methods:
-        if method not in dist_methods[dist]:
-            raise NotImplementedError(
-                f"{method} method is not implemented for {dist} distribution."
-            )
-    else:
-        raise NotImplementedError(f"{dist} distribution is not yet implemented.")
+    if isinstance(dist, str):
+        if dist in dist_methods:
+            if method not in dist_methods[dist]:
+                raise NotImplementedError(
+                    f"{method} method is not implemented for {dist} distribution"
+                )
+        else:
+            raise NotImplementedError(f"{dist} distribution is not yet implemented.")
 
     # Precipitation is expected to be zero-inflated
     zero_inflated = True
@@ -1254,7 +1262,7 @@ def standardized_precipitation_evapotranspiration_index(
     wb: xarray.DataArray,
     freq: str | None = "MS",
     window: int = 1,
-    dist: str = "gamma",
+    dist: str | rv_continuous = "gamma",
     method: str = "ML",
     fitkwargs: dict | None = None,
     cal_start: DateStr | None = None,
@@ -1279,13 +1287,14 @@ def standardized_precipitation_evapotranspiration_index(
     window : int
         Averaging window length relative to the resampling frequency. For example, if `freq="MS"`, i.e. a monthly
         resampling, the window is an integer number of months.
-    dist : {'gamma', 'fisk'}
-        Name of the univariate distribution. (see :py:mod:`scipy.stats`).
-    method : {'APP', 'ML'}
+    dist : {'gamma', 'fisk'} or `rv_continuous` function
+        Name of the univariate distribution, or a callable `rv_continuous` (see :py:mod:`scipy.stats`).
+    method : {"APP", "ML", "PWM"}
         Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method
-        uses a deterministic function that doesn't involve any optimization.
+        uses a deterministic function that does not involve any optimization. `PWM` should be used with a `lmoments3` distribution.
     fitkwargs : dict, optional
         Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`).
+        If method is `PWM`, `fitkwargs` should be empty, except for `floc` with `dist`=`gamma` which is allowed.
     cal_start : DateStr, optional
         Start date of the calibration period. A `DateStr` is expected, that is a `str` in format `"YYYY-MM-DD"`.
         Default option `None` means that the calibration period begins at the start of the input dataset.
@@ -1308,17 +1317,20 @@ def standardized_precipitation_evapotranspiration_index(
     See Also
     --------
     standardized_precipitation_index : Standardized Precipitation Index.
+    xclim.indices.stats.standardized_index : Standardized Index.
+    xclim.indices.stats.standardized_index_fit_params : Standardized Index Fit Params.
     """
     fitkwargs = fitkwargs or {}
 
     dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]}
-    if dist in dist_methods:
-        if method not in dist_methods[dist]:
-            raise NotImplementedError(
-                f"{method} method is not implemented for {dist} distribution"
-            )
-    else:
-        raise NotImplementedError(f"{dist} distribution is not yet implemented.")
+    if isinstance(dist, str):
+        if dist in dist_methods:
+            if method not in dist_methods[dist]:
+                raise NotImplementedError(
+                    f"{method} method is not implemented for {dist} distribution"
+                )
+        else:
+            raise NotImplementedError(f"{dist} distribution is not yet implemented.")
 
     # Water budget is not expected to be zero-inflated
     zero_inflated = False