From b385dba4f103de70277e9a173c63e71a6d7cfae6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Mon, 13 Jan 2025 16:13:25 -0500 Subject: [PATCH 01/28] warn fitkwargs & lmoments, shift gamma with floc --- CHANGELOG.rst | 8 ++++++++ src/xclim/indices/stats.py | 15 +++++++++++++++ 2 files changed, 23 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index a4d922613..61da8f36e 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -2,6 +2,14 @@ Changelog ========= +v0.55.0 (unpublished) +-------------------- +Contributors to this version: Éric Dupuis (:user:`coxipi`). + +Internal changes +^^^^^^^^^^^^^^^^ +* There is now a warning stating that `fitkwargs` are not employed when using `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. + v0.54.0 (2024-12-16) -------------------- Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Éric Dupuis (:user:`coxipi`), Sascha Hofmann (:user:`saschahofmann`). diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index 76a018baa..e42e6d770 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -59,6 +59,21 @@ def _fitfunc_1d(arr, *, dist, nparams, method, **fitkwargs): # lmoments3 will raise an error if only dist.numargs + 2 values are provided if len(x) <= dist.numargs + 2: return np.asarray([np.nan] * nparams) + if (type(dist).__name__ != "GammaGen" and len(fitkwargs.keys()) != 0) or ( + type(dist).__name__ == "GammaGen" + and set(fitkwargs.keys()) - {"floc"} != set() + ): + warnings.warn( + "Lmoments3 does not use `fitkwargs` arguments, except for `floc` with the Gamma distribution." + ) + if "floc" in fitkwargs and type(dist).__name__ == "GammaGen": + # lmoments3 assumes `loc` is 0, so `x` may need to be shifted + # note that `floc` must already be in appropriate units for `x` + params = dist.lmom_fit(x - fitkwargs["floc"]) + params["loc"] = fitkwargs["floc"] + params = list(params.values()) + else: + params = list(dist.lmom_fit(x).values()) params = list(dist.lmom_fit(x).values()) elif method == "APP": args, kwargs = _fit_start(x, dist.name, **fitkwargs) From 6209134b70838215d21d408fe294f18873fe2420 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Mon, 13 Jan 2025 16:16:57 -0500 Subject: [PATCH 02/28] update CHANGELOG --- CHANGELOG.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 61da8f36e..b4157122a 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -8,7 +8,7 @@ Contributors to this version: Éric Dupuis (:user:`coxipi`). Internal changes ^^^^^^^^^^^^^^^^ -* There is now a warning stating that `fitkwargs` are not employed when using `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. +* There is now a warning stating that `fitkwargs` are not employed when using `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`). v0.54.0 (2024-12-16) -------------------- From 1fb25bd42f00ae9d987824f407733add4db6ea20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= <71575674+coxipi@users.noreply.github.com> Date: Tue, 14 Jan 2025 09:37:34 -0500 Subject: [PATCH 03/28] remove old line Co-authored-by: Pascal Bourgault --- src/xclim/indices/stats.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index e42e6d770..897fa091b 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -74,7 +74,6 @@ def _fitfunc_1d(arr, *, dist, nparams, method, **fitkwargs): params = list(params.values()) else: params = list(dist.lmom_fit(x).values()) - params = list(dist.lmom_fit(x).values()) elif method == "APP": args, kwargs = _fit_start(x, dist.name, **fitkwargs) kwargs.setdefault("loc", 0) From 8eb2b51db385d8adc0e6aaa639586756ea9a4fd3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 14 Jan 2025 09:46:43 -0500 Subject: [PATCH 04/28] fitkwargs&PWM=error, allow PWM in SPI/SPEI --- src/xclim/indices/_agro.py | 6 ++++-- src/xclim/indices/stats.py | 5 +++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index 52726c3f9..baa08662a 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -1150,6 +1150,7 @@ def standardized_precipitation_index( uses a deterministic function that does not involve any optimization. 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. @@ -1281,11 +1282,12 @@ def standardized_precipitation_evapotranspiration_index( 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'} + 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. 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. @@ -1311,7 +1313,7 @@ def standardized_precipitation_evapotranspiration_index( """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} + dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP", "PWM"]} if dist in dist_methods: if method not in dist_methods[dist]: raise NotImplementedError( diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index 897fa091b..b021857e7 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -63,7 +63,7 @@ def _fitfunc_1d(arr, *, dist, nparams, method, **fitkwargs): type(dist).__name__ == "GammaGen" and set(fitkwargs.keys()) - {"floc"} != set() ): - warnings.warn( + raise ValueError( "Lmoments3 does not use `fitkwargs` arguments, except for `floc` with the Gamma distribution." ) if "floc" in fitkwargs and type(dist).__name__ == "GammaGen": @@ -884,8 +884,9 @@ def standardized_index( The approximate method uses a deterministic function that doesn't involve any optimization. zero_inflated : bool If True, the zeroes of `da` are treated separately. - fitkwargs : dict + 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. From 2ef071c03d07ffc4670d427326a816c6340e4638 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 14 Jan 2025 15:58:21 -0500 Subject: [PATCH 05/28] add tests --- src/xclim/indices/_agro.py | 6 +-- src/xclim/indices/stats.py | 2 +- tests/test_indices.py | 75 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 79 insertions(+), 4 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index baa08662a..fd9355959 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -1145,7 +1145,7 @@ def standardized_precipitation_index( 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"} + 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. fitkwargs : dict, optional @@ -1219,7 +1219,7 @@ def standardized_precipitation_index( >>> spi_3_fitted = standardized_precipitation_index(pr, params=params) """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} + dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} if dist in dist_methods: if method not in dist_methods[dist]: raise NotImplementedError( @@ -1313,7 +1313,7 @@ def standardized_precipitation_evapotranspiration_index( """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP", "PWM"]} + dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} if dist in dist_methods: if method not in dist_methods[dist]: raise NotImplementedError( diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index b021857e7..477f8d7cf 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -808,7 +808,7 @@ def standardized_index_fit_params( "Pass a value for `floc` in `fitkwargs`." ) - dist_and_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} + dist_and_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} dist = get_dist(dist) if dist.name not in dist_and_methods: raise NotImplementedError(f"The distribution `{dist.name}` is not supported.") diff --git a/tests/test_indices.py b/tests/test_indices.py index 22978b49d..c8444fec4 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -660,6 +660,38 @@ class TestStandardizedIndices: [-1.08627775, -0.46491398, -0.77806462, 0.31759127, 0.03794528], 2e-2, ), + ( + "D", + 1, + "gamma", + "PWM", + [-0.13002, 1.346689, 0.965731, 0.245408, -0.427896], + 2e-2, + ), + ( + "D", + 12, + "gamma", + "PWM", + [-0.209411, -0.086357, 0.636851, 1.022608, 0.634409], + 2e-2, + ), + ( + "MS", + 1, + "gamma", + "PWM", + [1.364243, 1.478565, 1.915559, -3.055828, 0.905304], + 2e-2, + ), + ( + "MS", + 12, + "gamma", + "PWM", + [0.577214, 1.522867, 1.634222, 0.967847, 0.689001], + 2e-2, + ), ], ) def test_standardized_precipitation_index( @@ -671,6 +703,13 @@ def test_standardized_precipitation_index( and Version(__numpy_version__) < Version("2.0.0") ): pytest.skip("Skipping SPI/ML/D on older numpy") + + # change `dist` to a lmoments3 object if needed + if method == "PWM": + lmom = pytest.importorskip("lmoments3.distr") + scipy2lmom = {"gamma": "gam"} + dist = getattr(lmom, scipy2lmom[dist]) + ds = open_dataset("sdba/CanESM2_1950-2100.nc").isel(location=1) if freq == "D": # to compare with ``climate_indices`` @@ -922,6 +961,42 @@ def test_zero_inflated(self, open_dataset): np.all(np.not_equal(spid[False].values, spid[True].values)), True ) + def test_PWM_and_fitkwargs(self, open_dataset): + pr = ( + open_dataset("sdba/CanESM2_1950-2100.nc") + .isel(location=1) + .sel(time=slice("1950", "1980")) + ).pr + + lmom = pytest.importorskip("lmoments3.distr") + # for now, only one function used + scipy2lmom = {"gamma": "gam"} + dist = getattr(lmom, scipy2lmom["gamma"]) + fitkwargs = {"floc": 0} + input_params = dict( + freq=None, + window=1, + method="PWM", + dist=dist, + fitkwargs=fitkwargs, + # doy_bounds=(180, 180), + ) + # this should not cause a problem + params_d0 = xci.stats.standardized_index_fit_params(pr, **input_params).isel( + dayofyear=0 + ) + np.testing.assert_allclose( + params_d0, np.array([5.63e-01, 0, 3.37e-05]), rtol=0, atol=2e-2 + ) + # this should cause a problem + fitkwargs["fscale"] = 1 + input_params["fitkwargs"] = fitkwargs + with pytest.raises( + ValueError, + match="Lmoments3 does not use `fitkwargs` arguments, except for `floc` with the Gamma distribution.", + ): + xci.stats.standardized_index_fit_params(pr, **input_params) + class TestDailyFreezeThawCycles: @pytest.mark.parametrize( From dbce9669227f78a032d53e09781197d9c7ebb091 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 14 Jan 2025 16:05:20 -0500 Subject: [PATCH 06/28] update CHANGELOG --- CHANGELOG.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index e01f7af36..ea7a24d95 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -13,6 +13,7 @@ New features and enhancements Internal changes ^^^^^^^^^^^^^^^^ * There is now a warning stating that `fitkwargs` are not employed when using `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`). +* `PWM` method (`lmoments3`) is now available to be used with the `gamma` distribution in ``xclim.indices.standardized_precipitation_index`` and ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`2043`, :pull:`2045`). v0.54.0 (2024-12-16) -------------------- From c5d02593ec46de3633d9b02cab97dd4c0dbdc955 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Wed, 15 Jan 2025 13:18:00 -0500 Subject: [PATCH 07/28] remove conditions on PWM --- src/xclim/indices/_agro.py | 14 ++++++++------ src/xclim/indices/stats.py | 17 ++++++++++------- 2 files changed, 18 insertions(+), 13 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index fd9355959..9637a1d50 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -1143,8 +1143,8 @@ 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`). + dist : {'gamma', 'fisk'} + Name of the univariate distribution. (see :py:mod:`scipy.stats`). All possible distributions are allowed with 'PWM'. 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. @@ -1225,7 +1225,8 @@ def standardized_precipitation_index( raise NotImplementedError( f"{method} method is not implemented for {dist} distribution." ) - else: + # Constraints on distributions except for PWM + elif method != "PWM": raise NotImplementedError(f"{dist} distribution is not yet implemented.") # Precipitation is expected to be zero-inflated @@ -1281,7 +1282,7 @@ def standardized_precipitation_evapotranspiration_index( 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`). + Name of the univariate distribution. (see :py:mod:`scipy.stats`). All possible distributions are allowed with 'PWM'. 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. @@ -1313,13 +1314,14 @@ def standardized_precipitation_evapotranspiration_index( """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + 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: + # Constraints on distributions except for PWM + elif method != "PWM": raise NotImplementedError(f"{dist} distribution is not yet implemented.") # Water budget is not expected to be zero-inflated diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index 477f8d7cf..e2c222ad5 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -808,14 +808,17 @@ def standardized_index_fit_params( "Pass a value for `floc` in `fitkwargs`." ) - dist_and_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + dist_and_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} dist = get_dist(dist) - if dist.name not in dist_and_methods: - raise NotImplementedError(f"The distribution `{dist.name}` is not supported.") - if method not in dist_and_methods[dist.name]: - raise NotImplementedError( - f"The method `{method}` is not supported for distribution `{dist.name}`." - ) + if method != "PWM": + if dist.name not in dist_and_methods: + raise NotImplementedError( + f"The distribution `{dist.name}` is not supported." + ) + if method not in dist_and_methods[dist.name]: + raise NotImplementedError( + f"The method `{method}` is not supported for distribution `{dist.name}`." + ) da, group = preprocess_standardized_index(da, freq, window, **indexer) if zero_inflated: prob_of_zero = da.groupby(group).map( From 6052fe553c25254144545eef14a127d889cc55b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Mon, 20 Jan 2025 10:57:13 -0500 Subject: [PATCH 08/28] no PWM constraint for SPI either --- src/xclim/indices/_agro.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index 9637a1d50..0601c3f0c 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -1219,11 +1219,12 @@ def standardized_precipitation_index( >>> spi_3_fitted = standardized_precipitation_index(pr, params=params) """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + + 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." + f"{method} method is not implemented for {dist} distribution" ) # Constraints on distributions except for PWM elif method != "PWM": From 3323a9b9f0f7e23bc55bf2796418278bed5dcd00 Mon Sep 17 00:00:00 2001 From: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> Date: Tue, 11 Feb 2025 15:08:58 -0500 Subject: [PATCH 09/28] add vapour pressure deficit indice --- src/xclim/indices/_conversion.py | 45 +++++++++++++++++++++++++++++++- tests/test_indices.py | 28 +++++++++++++++++--- 2 files changed, 69 insertions(+), 4 deletions(-) diff --git a/src/xclim/indices/_conversion.py b/src/xclim/indices/_conversion.py index 329171041..f06086eaf 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -52,6 +52,7 @@ "tas", "uas_vas_2_sfcwind", "universal_thermal_climate_index", + "vapor_pressure_deficit", "wind_chill_index", "wind_power_potential", "wind_profile", @@ -356,7 +357,7 @@ def sfcwind_2_uas_vas( def saturation_vapor_pressure( tas: xr.DataArray, ice_thresh: Quantified | None = None, - method: str = "sonntag90", # noqa + method: str = "sonntag90", ) -> xr.DataArray: """ Saturation vapour pressure from temperature. @@ -492,6 +493,48 @@ def saturation_vapor_pressure( return e_sat +@declare_units(tas="[temperature]", hurs="[]", ice_thresh="[temperature]") +def vapor_pressure_deficit( + tas: xr.DataArray, + hurs: xr.DataArray, + ice_thresh: Quantified | None = None, + method: str = "sonntag90", +) -> xr.DataArray: + """ + Vapour pressure deficit. + + The measure of the moisture deficit of the air. + + Parameters + ---------- + tas : xarray.DataArray + Mean daily temperature. + hurs : xarray.DataArray + Relative humidity. + ice_thresh : Quantified, optional + Threshold temperature under which to switch to equations in reference to ice instead of water. + If None (default) everything is computed with reference to water. + method : {"goffgratch46", "sonntag90", "tetens30", "wmo08", "its90"} + Which method to use, see notes of :py:func:`saturation_vapor_pressure`. + Default is "sonntag90". + + Returns + ------- + xarray.DataArray, [Pa] + Vapour pressure deficit. + + See Also + -------- + :py:func:`xclim.indices.saturation_vapor_pressure` + """ + svp = saturation_vapor_pressure(tas, ice_thresh=ice_thresh, method=method) + + vpd = cast(xr.DataArray, (1 - (hurs / 100)) * svp) + + vpd = vpd.assign_attrs(units=svp.attrs["units"]) + return vpd + + @declare_units( tas="[temperature]", tdps="[temperature]", diff --git a/tests/test_indices.py b/tests/test_indices.py index aa9276bd3..64a4f95d9 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -2837,10 +2837,11 @@ def test_specific_humidity_from_dewpoint(tas_series, ps_series): @pytest.mark.parametrize( "ice_thresh,exp0", [(None, [125, 286, 568]), ("0 degC", [103, 260, 563])] ) -@pytest.mark.parametrize("units", ["degC", "degK"]) -def test_saturation_vapor_pressure(tas_series, method, ice_thresh, exp0, units): +@pytest.mark.parametrize("temp_units", ["degC", "degK"]) +def test_saturation_vapor_pressure(tas_series, method, ice_thresh, exp0, temp_units): tas = tas_series(np.array([-20, -10, -1, 10, 20, 25, 30, 40, 60]) + K2C) - tas = convert_units_to(tas, units) + tas = convert_units_to(tas, temp_units) + # Expected values obtained with the Sonntag90 method e_sat_exp = exp0 + [1228, 2339, 3169, 4247, 7385, 19947] @@ -2852,6 +2853,26 @@ def test_saturation_vapor_pressure(tas_series, method, ice_thresh, exp0, units): np.testing.assert_allclose(e_sat, e_sat_exp, atol=0.5, rtol=0.005) +@pytest.mark.parametrize( + "method", ["tetens30", "sonntag90", "goffgratch46", "wmo08", "its90"] +) +def test_vapor_pressure_deficit( + tas_series, hurs_series, method +): + tas = tas_series(np.array([-1, 10, 20, 25, 30, 40, 60]) + K2C) + hurs = hurs_series(np.array([0, 0.5, 0.8, 0.9, 0.95, 0.99, 1])) + + # Expected values obtained with the GroffGratch46 method + svp_exp = [567, 1220, 2317, 3136, 4200, 7300, 19717] + + vpd = xci.vapor_pressure_deficit( + tas=tas, + hurs=hurs, + method=method, + ) + np.testing.assert_allclose(vpd, svp_exp, atol=0.5, rtol=0.005) + + @pytest.mark.parametrize("method", ["tetens30", "sonntag90", "goffgratch46", "wmo08"]) @pytest.mark.parametrize( "invalid_values,exp0", [("clip", 100), ("mask", np.nan), (None, 188)] @@ -2860,6 +2881,7 @@ def test_relative_humidity( tas_series, hurs_series, huss_series, ps_series, method, invalid_values, exp0 ): tas = tas_series(np.array([-10, -10, 10, 20, 35, 50, 75, 95]) + K2C) + # Expected values obtained with the Sonntag90 method hurs_exp = hurs_series([exp0, 63.0, 66.0, 34.0, 14.0, 6.0, 1.0, 0.0]) ps = ps_series([101325] * 8) From 59d1fa4cde62098a444ca42ce0c6a5f4219aef9c Mon Sep 17 00:00:00 2001 From: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> Date: Tue, 11 Feb 2025 15:28:19 -0500 Subject: [PATCH 10/28] add indicators --- src/xclim/data/fr.json | 6 ++++++ src/xclim/indicators/atmos/_conversion.py | 20 ++++++++++++++++++++ 2 files changed, 26 insertions(+) diff --git a/src/xclim/data/fr.json b/src/xclim/data/fr.json index 548eb3cfb..3410c4695 100644 --- a/src/xclim/data/fr.json +++ b/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." }, + "E_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.", diff --git a/src/xclim/indicators/atmos/_conversion.py b/src/xclim/indicators/atmos/_conversion.py index 4d7ba5f9f..cb4929349 100644 --- a/src/xclim/indicators/atmos/_conversion.py +++ b/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,25 @@ def cfcheck(self, **das) -> None: compute=indices.specific_humidity_from_dewpoint, ) + +vapor_pressure_deficit = Converter( + title="Water vapour pressure deficit", + identifier="e_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", From 97dd3df2446bd3462a291e868e8475026a8203bf Mon Sep 17 00:00:00 2001 From: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> Date: Tue, 11 Feb 2025 15:31:54 -0500 Subject: [PATCH 11/28] update CHANGELOG.rst --- CHANGELOG.rst | 1 + tests/test_indices.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 5a371d480..c3381c81e 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -15,6 +15,7 @@ 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. New features and enhancements ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/tests/test_indices.py b/tests/test_indices.py index 64a4f95d9..d5a4496e4 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -2862,7 +2862,7 @@ def test_vapor_pressure_deficit( tas = tas_series(np.array([-1, 10, 20, 25, 30, 40, 60]) + K2C) hurs = hurs_series(np.array([0, 0.5, 0.8, 0.9, 0.95, 0.99, 1])) - # Expected values obtained with the GroffGratch46 method + # Expected values obtained with the GoffGratch46 method svp_exp = [567, 1220, 2317, 3136, 4200, 7300, 19717] vpd = xci.vapor_pressure_deficit( From 42e072856374fbb3613b82ae0402828fe988445d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 11 Feb 2025 20:36:56 +0000 Subject: [PATCH 12/28] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/xclim/data/fr.json | 4 ++-- src/xclim/indicators/atmos/_conversion.py | 3 ++- tests/test_indices.py | 4 +--- 3 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/xclim/data/fr.json b/src/xclim/data/fr.json index 3410c4695..73548d2ee 100644 --- a/src/xclim/data/fr.json +++ b/src/xclim/data/fr.json @@ -931,12 +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." }, - "E_DEFICIT": { + "E_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.", diff --git a/src/xclim/indicators/atmos/_conversion.py b/src/xclim/indicators/atmos/_conversion.py index cb4929349..83c510183 100644 --- a/src/xclim/indicators/atmos/_conversion.py +++ b/src/xclim/indicators/atmos/_conversion.py @@ -281,7 +281,8 @@ def cfcheck(self, **das) -> None: 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 "" diff --git a/tests/test_indices.py b/tests/test_indices.py index d5a4496e4..31f19cda3 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -2856,9 +2856,7 @@ def test_saturation_vapor_pressure(tas_series, method, ice_thresh, exp0, temp_un @pytest.mark.parametrize( "method", ["tetens30", "sonntag90", "goffgratch46", "wmo08", "its90"] ) -def test_vapor_pressure_deficit( - tas_series, hurs_series, method -): +def test_vapor_pressure_deficit(tas_series, hurs_series, method): tas = tas_series(np.array([-1, 10, 20, 25, 30, 40, 60]) + K2C) hurs = hurs_series(np.array([0, 0.5, 0.8, 0.9, 0.95, 0.99, 1])) From ac876d8879690761aa9fa34459532315d9e666f9 Mon Sep 17 00:00:00 2001 From: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> Date: Tue, 11 Feb 2025 15:43:28 -0500 Subject: [PATCH 13/28] fix see also --- src/xclim/indices/_conversion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/xclim/indices/_conversion.py b/src/xclim/indices/_conversion.py index f06086eaf..15cb9f0c0 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -525,7 +525,7 @@ def vapor_pressure_deficit( See Also -------- - :py:func:`xclim.indices.saturation_vapor_pressure` + saturation_vapor_pressure : Vapour pressure at saturation. """ svp = saturation_vapor_pressure(tas, ice_thresh=ice_thresh, method=method) From aa415da6fba4ae161f4b16c3c31a42030fbf97d9 Mon Sep 17 00:00:00 2001 From: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> Date: Tue, 11 Feb 2025 15:55:19 -0500 Subject: [PATCH 14/28] update CHANGELOG.rst --- CHANGELOG.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index f7a397257..b96da0c14 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -19,7 +19,7 @@ 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. +* 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 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ From 7133fcc1c4a42a626e75a421e48da6f648f4a6e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 11 Feb 2025 16:14:18 -0500 Subject: [PATCH 15/28] no constraints on any rv_continuous function --- src/xclim/indices/_agro.py | 51 +++++++++++++++++++------------------- 1 file changed, 26 insertions(+), 25 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index 0601c3f0c..c08d138e9 100644 --- a/src/xclim/indices/_agro.py +++ b/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,11 +1144,11 @@ 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`). All possible distributions are allowed with 'PWM'. + 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. @@ -1221,14 +1222,14 @@ def standardized_precipitation_index( 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" - ) - # Constraints on distributions except for PWM - elif method != "PWM": - 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 @@ -1257,7 +1258,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, @@ -1282,11 +1283,11 @@ 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`). All possible distributions are allowed with 'PWM'. - method : {'APP', 'ML', 'PWM'} + 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. @@ -1316,14 +1317,14 @@ def standardized_precipitation_evapotranspiration_index( 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" - ) - # Constraints on distributions except for PWM - elif method != "PWM": - 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 From 3f0f24b8f3a93c4e8164590ebb2c01ffcc3d76b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 11 Feb 2025 16:37:12 -0500 Subject: [PATCH 16/28] update CHANGELOG --- CHANGELOG.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 6a1df537e..7c66248ce 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -23,6 +23,8 @@ New features and enhancements * `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`). +* 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`). +* `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 ^^^^^^^^^^^^^^^^ @@ -30,7 +32,6 @@ Internal changes * Adjusted the ``TestOfficialYaml`` test to use a dynamic method for finding the installed location of `xclim`. (:pull:`2028`). * Adjusted two tests for better handling when running in Windows environments. (:pull:`2057`). * 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`). -* The `PWM` method (from `lmoments3`) is now available to be used with the `gamma` distribution in ``xclim.indices.standardized_precipitation_index`` and ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`2043`, :pull:`2045`). Bug fixes ^^^^^^^^^ From daf85192df42cfc600a2245d1d9a7ef790f50cf5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 11 Feb 2025 16:41:19 -0500 Subject: [PATCH 17/28] test rv_continuous --- tests/test_indices.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/tests/test_indices.py b/tests/test_indices.py index b6c0f692e..6de6fc626 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -22,6 +22,7 @@ from numpy import __version__ as __numpy_version__ from packaging.version import Version from pint import __version__ as __pint_version__ +from scipy import stats from xclim import indices as xci from xclim.core import ValidationError @@ -749,6 +750,34 @@ def test_standardized_precipitation_index( np.testing.assert_allclose(spi.values, values, rtol=0, atol=diff_tol) + @pytest.mark.slow + @pytest.mark.parametrize("dist", ["gamma", "fisk"]) + def test_str_vs_rv_continuous(self, open_dataset, dist): + ds = open_dataset("sdba/CanESM2_1950-2100.nc").isel(location=1) + window = 1 + method = "ML" + freq = "MS" + + pr = ds.pr.sel(time=slice("1998", "2000")) + pr_cal = ds.pr.sel(time=slice("1950", "1980")) + fitkwargs = {} + + out = [] + for dist0 in [dist, getattr(stats, dist)]: + params = xci.stats.standardized_index_fit_params( + pr_cal, + freq=freq, + window=window, + dist=dist0, + method=method, + fitkwargs=fitkwargs, + zero_inflated=True, + ) + spi = xci.standardized_precipitation_index(pr, params=params) + # Only a few moments before year 2000 are tested + out.append(spi.isel(time=slice(-11, -1, 2))) + assert all(out[0] == out[1]) + # See SPI version @pytest.mark.slow @pytest.mark.parametrize( From 4920932cb00399366054f45d7325d77d94d97d41 Mon Sep 17 00:00:00 2001 From: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> Date: Tue, 11 Feb 2025 17:24:05 -0500 Subject: [PATCH 18/28] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Éric Dupuis <71575674+coxipi@users.noreply.github.com> --- src/xclim/data/fr.json | 2 +- src/xclim/indicators/atmos/_conversion.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/xclim/data/fr.json b/src/xclim/data/fr.json index 73548d2ee..8efa0c85b 100644 --- a/src/xclim/data/fr.json +++ b/src/xclim/data/fr.json @@ -931,7 +931,7 @@ "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." }, - "E_DEFICIT": { + "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", diff --git a/src/xclim/indicators/atmos/_conversion.py b/src/xclim/indicators/atmos/_conversion.py index 83c510183..66f83cd8e 100644 --- a/src/xclim/indicators/atmos/_conversion.py +++ b/src/xclim/indicators/atmos/_conversion.py @@ -274,7 +274,7 @@ def cfcheck(self, **das) -> None: vapor_pressure_deficit = Converter( title="Water vapour pressure deficit", - identifier="e_deficit", + identifier="vapor_pressure_deficit", units="Pa", long_name='Vapour pressure deficit ("{method}" method)', standard_name="water_vapor_saturation_deficit_in_air", From 162f42a0e80dee252ab67f1917e1cf5833317bcd Mon Sep 17 00:00:00 2001 From: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> Date: Tue, 11 Feb 2025 17:26:05 -0500 Subject: [PATCH 19/28] update docstring --- src/xclim/indices/_conversion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/xclim/indices/_conversion.py b/src/xclim/indices/_conversion.py index 15cb9f0c0..98a4eb0ee 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -515,7 +515,7 @@ def vapor_pressure_deficit( Threshold temperature under which to switch to equations in reference to ice instead of water. If None (default) everything is computed with reference to water. method : {"goffgratch46", "sonntag90", "tetens30", "wmo08", "its90"} - Which method to use, see notes of :py:func:`saturation_vapor_pressure`. + Method used to calculate saturation vapour pressure, see notes of :py:func:`saturation_vapor_pressure`. Default is "sonntag90". Returns From 6a8a5a12a6b0d3821ec04a1a716026246d9180fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Wed, 12 Feb 2025 13:48:59 -0500 Subject: [PATCH 20/28] better doc --- src/xclim/indices/_agro.py | 16 +++++++++++----- src/xclim/indices/stats.py | 17 +++++++++++++++-- 2 files changed, 26 insertions(+), 7 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index cfa49faa1..e155f0c9b 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -1171,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 ---------- @@ -1313,6 +1317,8 @@ 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 {} diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index 085ff7e46..ca2b2874a 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -788,9 +788,11 @@ def standardized_index_fit_params( Notes ----- Supported combinations of `dist` and `method` are: - * Gamma ("gamma") : "ML", "APP", "PWM" + * Gamma ("gamma") : "ML", "APP" * Log-logistic ("fisk") : "ML", "APP" * "APP" method only supports two-parameter distributions. Parameter `loc` will be set to 0 (setting `floc=0` in `fitkwargs`). + * Otherwise, generic `rv_continuous` methods can be used. This includes distributions from `lmoments3` which should be used with + `method="PWM"`. When using the zero inflated option, : A probability density function :math:`\texttt{pdf}_0(X)` is fitted for :math:`X \neq 0` and a supplementary parameter :math:`\pi` takes into account the probability of :math:`X = 0`. The full probability density @@ -909,11 +911,22 @@ def standardized_index( xarray.DataArray, [unitless] Standardized Precipitation Index. + See Also + -------- + standardized_index_fit_params : Standardized Index Fit Params. + Notes ----- * 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. - * ``window``, ``dist``, ``method``, ``zero_inflated`` are only optional if ``params`` is given. + * ``window``, ``dist``, ``method``, ``zero_inflated`` are only optional if ``params`` is given. If `params` is given as input, + it overrides the `cal_start`, `cal_end`, `freq` and `window`, `dist` and `method` options. + * Supported combinations of `dist` and `method` are: + * Gamma ("gamma") : "ML", "APP" + * Log-logistic ("fisk") : "ML", "APP" + * "APP" method only supports two-parameter distributions. Parameter `loc` will be set to 0 (setting `floc=0` in `fitkwargs`). + * Otherwise, generic `rv_continuous` methods can be used. This includes distributions from `lmoments3` which should be used with + `method="PWM"`. References ---------- From a50929a7cd79872c1a1aa1a9219e310656640a5e Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Wed, 12 Feb 2025 15:12:06 -0500 Subject: [PATCH 21/28] Fix doc for custom missing methods --- docs/notebooks/customize.ipynb | 35 ++++++++++++++++++++++++++-------- 1 file changed, 27 insertions(+), 8 deletions(-) diff --git a/docs/notebooks/customize.ipynb b/docs/notebooks/customize.ipynb index 3d28e5f21..5a60e24e3 100644 --- a/docs/notebooks/customize.ipynb +++ b/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", + "
\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 preserverd, 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", + "
\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 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 the following arguments:\n", + "\n", + "- `null`, a `DataArray` of the mask of invalid values in the input data array.\n", + "- `count`, the number of days in each resampled periods\n", + "- `freq`, the resampling frequency.\n", + "\n", + "The `is_missing` method should return a boolean mask, 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` staticmethod, which validates those. 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,12 @@ "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, max_n=5):\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 null.resample(time=freq).map(longest_run, dim=\"time\") >= max_n\n", "\n", " @staticmethod\n", " def validate(max_n):\n", @@ -204,6 +216,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 +236,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.12.7" } }, "nbformat": 4, From bef557ce3418085e071e2fc741cda6712efd6e3d Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Wed, 12 Feb 2025 15:16:51 -0500 Subject: [PATCH 22/28] fix words --- docs/notebooks/customize.ipynb | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/notebooks/customize.ipynb b/docs/notebooks/customize.ipynb index 5a60e24e3..139b8521b 100644 --- a/docs/notebooks/customize.ipynb +++ b/docs/notebooks/customize.ipynb @@ -151,19 +151,19 @@ "\n", "
\n", "\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 preserverd, but custom subclasses might break with future changes.\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", "
\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 the following arguments:\n", + "Finally, it is possible for advanced users to register their own method. 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.\n", - "- `count`, the number of days in each resampled periods\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, 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", + "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` staticmethod, which validates those. See example below and the docstrings in the module.\n", + "To add additional arguments, one should override the `__init__` (receiving those arguments) and the `validate` static method, which validates them. 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. " ] From 9479f51d3a7116846b34ce9d479a61a889b3b37d Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Wed, 12 Feb 2025 15:31:00 -0500 Subject: [PATCH 23/28] Correct way to pass extra args --- docs/notebooks/customize.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/notebooks/customize.ipynb b/docs/notebooks/customize.ipynb index 139b8521b..90e5007e5 100644 --- a/docs/notebooks/customize.ipynb +++ b/docs/notebooks/customize.ipynb @@ -163,7 +163,7 @@ "\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. See example below and the docstrings in the module.\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. " ] @@ -185,9 +185,9 @@ " def __init__(self, max_n: int = 5):\n", " super().__init__(max_n=max_n)\n", "\n", - " def is_missing(self, null, count, freq, max_n=5):\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.resample(time=freq).map(longest_run, dim=\"time\") >= max_n\n", + " return null.resample(time=freq).map(longest_run, dim=\"time\") >= self.options['max_n']\n", "\n", " @staticmethod\n", " def validate(max_n):\n", From 364aa81147d511ad5fcd5ff209902f1af1d4f7cd Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Wed, 12 Feb 2025 15:31:44 -0500 Subject: [PATCH 24/28] Apply suggestions from code review Co-authored-by: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> --- docs/notebooks/customize.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/notebooks/customize.ipynb b/docs/notebooks/customize.ipynb index 90e5007e5..5ca84e5aa 100644 --- a/docs/notebooks/customize.ipynb +++ b/docs/notebooks/customize.ipynb @@ -155,13 +155,13 @@ "\n", "\n", "\n", - "Finally, it is possible for advanced users to register their own method. 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", + "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", + "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", From 78f8e76198f3ad89a21c668c7ab11013109486d1 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 12 Feb 2025 20:31:50 +0000 Subject: [PATCH 25/28] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- docs/notebooks/customize.ipynb | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/docs/notebooks/customize.ipynb b/docs/notebooks/customize.ipynb index 5ca84e5aa..cdb4d93bc 100644 --- a/docs/notebooks/customize.ipynb +++ b/docs/notebooks/customize.ipynb @@ -187,7 +187,10 @@ "\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.resample(time=freq).map(longest_run, dim=\"time\") >= self.options['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", From ec7722133892cdd192fa288686cb3ada9bded9c9 Mon Sep 17 00:00:00 2001 From: Ouranos Helper Bot Date: Wed, 12 Feb 2025 21:42:19 +0000 Subject: [PATCH 26/28] =?UTF-8?q?Bump=20version:=200.54.1-dev.11=20?= =?UTF-8?q?=E2=86=92=200.54.1-dev.12?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Ouranos Helper Bot --- pyproject.toml | 2 +- src/xclim/__init__.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 42f942171..782144f84 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -137,7 +137,7 @@ target-version = [ ] [tool.bumpversion] -current_version = "0.54.1-dev.11" +current_version = "0.54.1-dev.12" commit = true commit_args = "--no-verify --signoff" tag = false diff --git a/src/xclim/__init__.py b/src/xclim/__init__.py index a635b16ee..c73c39df8 100644 --- a/src/xclim/__init__.py +++ b/src/xclim/__init__.py @@ -13,7 +13,7 @@ __author__ = """Travis Logan""" __email__ = "logan.travis@ouranos.ca" -__version__ = "0.54.1-dev.11" +__version__ = "0.54.1-dev.12" with _resources.as_file(_resources.files("xclim.data")) as _module_data: From 9c608fff25ecd1261aed4addfc5078eb86f9d97b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Thu, 13 Feb 2025 00:25:10 -0500 Subject: [PATCH 27/28] remove comment --- tests/test_indices.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index 6de6fc626..835582407 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -1009,7 +1009,6 @@ def test_PWM_and_fitkwargs(self, open_dataset): method="PWM", dist=dist, fitkwargs=fitkwargs, - # doy_bounds=(180, 180), ) # this should not cause a problem params_d0 = xci.stats.standardized_index_fit_params(pr, **input_params).isel( From 00de70377b9f15c46d82adda2bff6c39da783388 Mon Sep 17 00:00:00 2001 From: Ouranos Helper Bot Date: Thu, 13 Feb 2025 05:37:57 +0000 Subject: [PATCH 28/28] =?UTF-8?q?Bump=20version:=200.54.1-dev.12=20?= =?UTF-8?q?=E2=86=92=200.54.1-dev.13?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Ouranos Helper Bot --- pyproject.toml | 2 +- src/xclim/__init__.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 782144f84..b344ba32d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -137,7 +137,7 @@ target-version = [ ] [tool.bumpversion] -current_version = "0.54.1-dev.12" +current_version = "0.54.1-dev.13" commit = true commit_args = "--no-verify --signoff" tag = false diff --git a/src/xclim/__init__.py b/src/xclim/__init__.py index c73c39df8..9becf7526 100644 --- a/src/xclim/__init__.py +++ b/src/xclim/__init__.py @@ -13,7 +13,7 @@ __author__ = """Travis Logan""" __email__ = "logan.travis@ouranos.ca" -__version__ = "0.54.1-dev.12" +__version__ = "0.54.1-dev.13" with _resources.as_file(_resources.files("xclim.data")) as _module_data: