From d52d11742ccbe184cb153a711a721201b2a06343 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Wed, 5 Feb 2025 17:20:18 +0100 Subject: [PATCH 01/10] Exportable function for fao_allen98 --- src/xclim/indices/_conversion.py | 59 ++++++++++++++++++++++++++++---- tests/test_atmos.py | 14 ++++++++ tests/test_indices.py | 2 +- 3 files changed, 67 insertions(+), 8 deletions(-) diff --git a/src/xclim/indices/_conversion.py b/src/xclim/indices/_conversion.py index b8868baef..6d083907f 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -1319,6 +1319,55 @@ def _get_D_from_M(time): # noqa: N802 ) +# TODO: How can I check exact units for inputs +# @declare_units( +# net_radiation="[MJ m-2 day-1]", +# tas="[degC]", +# wind="[m s-1]", +# es="[kPa]", +# ea="[kPa]", +# delta_svp="[kPa degC-1]", +# gamma="[kPa degC]", +# G="[MJ m-2 day-1]", +# ) +def fao_allen98(net_radiation, tas, wind, es, ea, delta_svp, gamma, G=0): + r""" + FAO-56 Penman-Monteith equation. + Estimates reference evapotranspiration (ETo) from a hypothetical short grass reference surface. + Based on equation 6 in based on :cite:t:`allen_crop_1998`. + + Parameters + ---------- + net_radiation : xarray.DataArray + Net radiation at crop surface [MJ m-2 day-1]. + tas : xarray.DataArray + Air temperature at 2 m height [degC]. + wind : xarray.DataArray + Wind speed at 2 m height [m s-1]. + es : xarray.DataArray + Saturation vapour pressure [kPa]. + ea : xarray.DataArray + Actual vapour pressure [kPa]. + delta_svp : xarray.DataArray + Slope of saturation vapour pressure curve [kPa degC-1]. + gamma : xarray.DataArray + Psychrometric constant [kPa deg C]. + G : xarray.DataArray, optional + Soil heat flux (G) [MJ m-2 day-1] (For daily default to 0). + + Returns + ------- + xarray.DataArray + Potential Evapotranspiration from a hypothetical grass reference surface [mm day-1]. + """ + tasK = convert_units_to(tas, "K") + a1 = 0.408 * delta_svp * (net_radiation - G) + a2 = gamma * 900 / (tas + 273) * wind * (es - ea) + a3 = delta_svp + (gamma * (1 + 0.34 * wind)) + + return (a1 + a2) / a3 + + @declare_units( tasmin="[temperature]", tasmax="[temperature]", @@ -1563,7 +1612,7 @@ def potential_evapotranspiration( elif method in ["allen98", "FAO_PM98"]: _tasmax = convert_units_to(tasmax, "degC") _tasmin = convert_units_to(tasmin, "degC") - + _hurs = convert_units_to(hurs, "1") if sfcWind is None: raise ValueError("Wind speed is required for Allen98 method.") @@ -1580,14 +1629,13 @@ def potential_evapotranspiration( ) es = convert_units_to(es, "kPa") # mean actual vapour pressure [kPa] - ea = hurs * es + ea = _hurs * es # slope of saturation vapour pressure curve [kPa degC-1] delta = 4098 * es / (tas_m + 237.3) ** 2 # net radiation Rn = convert_units_to(rsds - rsus - (rlus - rlds), "MJ m-2 d-1") - G = 0 # Daily soil heat flux density [MJ m-2 d-1] P = 101.325 # Atmospheric pressure [kPa] gamma = 0.665e-03 * P # psychrometric const = C_p*P/(eps*lam) [kPa degC-1] @@ -1595,10 +1643,7 @@ def potential_evapotranspiration( # height = 0.12m, surface resistance = 70 s m-1, albedo = 0.23 # Surface resistance implies a ``moderately dry soil surface resulting from # about a weekly irrigation frequency'' - pet = ( - 0.408 * delta * (Rn - G) - + gamma * (900 / (tas_m + 273)) * wa2 * (es - ea) - ) / (delta + gamma * (1 + 0.34 * wa2)) + pet = fao_allen98(Rn, tas_m, wa2, es, ea, delta, gamma) else: raise NotImplementedError(f"'{method}' method is not implemented.") diff --git a/tests/test_atmos.py b/tests/test_atmos.py index 6f3a10848..6bcd289c0 100644 --- a/tests/test_atmos.py +++ b/tests/test_atmos.py @@ -343,6 +343,8 @@ def test_convert_units(self, atmosds): tnC.attrs["units"] = "degC" tmC = tm - K2C tmC.attrs["units"] = "degC" + hurs_pct = hurs * 100 + hurs_pct.attrs["units"] = "%" pet_br65 = atmos.potential_evapotranspiration(tn, tx, method="BR65") pet_br65C = atmos.potential_evapotranspiration(tnC, tx, method="BR65") @@ -374,12 +376,24 @@ def test_convert_units(self, atmosds): sfcWind=sfcWind, method="FAO_PM98", ) + pet_fao_pm98pct = atmos.potential_evapotranspiration( + tn, + tx, + hurs=hurs_pct, + rsds=rsds, + rsus=rsus, + rlds=rlds, + rlus=rlus, + sfcWind=sfcWind, + method="FAO_PM98", + ) np.testing.assert_allclose(pet_br65, pet_br65C, atol=1) np.testing.assert_allclose(pet_hg85, pet_hg85C, atol=1) np.testing.assert_allclose(pet_tw48, pet_tw48C, atol=1) np.testing.assert_allclose(pet_mb05, pet_mb05C, atol=1) np.testing.assert_allclose(pet_fao_pm98, pet_fao_pm98C, atol=1) + np.testing.assert_allclose(pet_fao_pm98, pet_fao_pm98pct, atol=1) def test_nan_values(self, atmosds): ds = atmosds diff --git a/tests/test_indices.py b/tests/test_indices.py index 6987959a3..aa9276bd3 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -3515,7 +3515,7 @@ def test_allen( tn = tasmin_series(np.array([0, 5, 10]) + 273.15).expand_dims(lat=lat) tx = tasmax_series(np.array([10, 15, 20]) + 273.15).expand_dims(lat=lat) tm = tas_series(np.array([5, 10, 15]) + 273.15).expand_dims(lat=lat) - hurs = hurs_series(np.array([0.8, 0.7, 0.73])).expand_dims(lat=lat) + hurs = hurs_series(np.array([80, 70, 73])).expand_dims(lat=lat) rsds = rsds_series(np.array([43.09, 43.57, 70.20])).expand_dims(lat=lat) rsus = rsus_series(np.array([12.51, 14.46, 20.36])).expand_dims(lat=lat) rlds = rlds_series(np.array([293.65, 228.96, 275.40])).expand_dims(lat=lat) From 2dd0374e06091f82436bae42f02744f29f0b6364 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Wed, 5 Feb 2025 17:40:30 +0100 Subject: [PATCH 02/10] Edit changelog --- CHANGELOG.rst | 2 ++ src/xclim/indices/_conversion.py | 10 ++++------ 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index f476f6d7e..0c27dfc68 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -17,6 +17,7 @@ New features and enhancements * 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`). +* ``xclim.indicators.atmos`` now exports the FAO-56 Penman-Monteith equation as `fao_allen98` (:issue:`2004`, :pull:`2067`). Internal changes ^^^^^^^^^^^^^^^^ @@ -27,6 +28,7 @@ Internal changes Bug fixes ^^^^^^^^^ * Fixed a bug in ``xclim.sdba.Grouper.get_index`` that didn't correctly interpolate seasonal values (:issue:`2014`, :pull:`2019`). +* Fixed a bug where ``xclim.indicators.atmos.potential_evapotranspiration`` would return wrong results when `hurs` was provided in `%` (:pull:`2067`). v0.54.0 (2024-12-16) -------------------- diff --git a/src/xclim/indices/_conversion.py b/src/xclim/indices/_conversion.py index 6d083907f..4f3b448e8 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -1333,7 +1333,9 @@ def _get_D_from_M(time): # noqa: N802 def fao_allen98(net_radiation, tas, wind, es, ea, delta_svp, gamma, G=0): r""" FAO-56 Penman-Monteith equation. - Estimates reference evapotranspiration (ETo) from a hypothetical short grass reference surface. + Estimates reference evapotranspiration (ETo) from a hypothetical short grass reference surface ( + height = 0.12m, surface resistance = 70 s m-1, albedo = 0.23 and a ``moderately dry soil surface resulting from + about a weekly irrigation frequency``). Based on equation 6 in based on :cite:t:`allen_crop_1998`. Parameters @@ -1352,7 +1354,7 @@ def fao_allen98(net_radiation, tas, wind, es, ea, delta_svp, gamma, G=0): Slope of saturation vapour pressure curve [kPa degC-1]. gamma : xarray.DataArray Psychrometric constant [kPa deg C]. - G : xarray.DataArray, optional + G : float, optional Soil heat flux (G) [MJ m-2 day-1] (For daily default to 0). Returns @@ -1639,10 +1641,6 @@ def potential_evapotranspiration( P = 101.325 # Atmospheric pressure [kPa] gamma = 0.665e-03 * P # psychrometric const = C_p*P/(eps*lam) [kPa degC-1] - # Penman-Monteith formula with reference grass: - # height = 0.12m, surface resistance = 70 s m-1, albedo = 0.23 - # Surface resistance implies a ``moderately dry soil surface resulting from - # about a weekly irrigation frequency'' pet = fao_allen98(Rn, tas_m, wa2, es, ea, delta, gamma) else: From 063f2f677298c03f5a0b1c59cf474ea86d1a6d89 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Wed, 5 Feb 2025 17:43:31 +0100 Subject: [PATCH 03/10] Fix using tasK in fao equation --- 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 4f3b448e8..c7b6b5a05 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -1364,7 +1364,7 @@ def fao_allen98(net_radiation, tas, wind, es, ea, delta_svp, gamma, G=0): """ tasK = convert_units_to(tas, "K") a1 = 0.408 * delta_svp * (net_radiation - G) - a2 = gamma * 900 / (tas + 273) * wind * (es - ea) + a2 = gamma * 900 / (tasK) * wind * (es - ea) a3 = delta_svp + (gamma * (1 + 0.34 * wind)) return (a1 + a2) / a3 From c199cf9cfbabe8883891eb7524ccffaccfd22612 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Wed, 5 Feb 2025 18:11:15 +0100 Subject: [PATCH 04/10] Fix linting --- src/xclim/indices/_conversion.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/xclim/indices/_conversion.py b/src/xclim/indices/_conversion.py index c7b6b5a05..d54dfca10 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -1333,6 +1333,7 @@ def _get_D_from_M(time): # noqa: N802 def fao_allen98(net_radiation, tas, wind, es, ea, delta_svp, gamma, G=0): r""" FAO-56 Penman-Monteith equation. + Estimates reference evapotranspiration (ETo) from a hypothetical short grass reference surface ( height = 0.12m, surface resistance = 70 s m-1, albedo = 0.23 and a ``moderately dry soil surface resulting from about a weekly irrigation frequency``). From 44580762ad95b6e337d2442f794518c5168fd1da Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Wed, 5 Feb 2025 22:18:26 +0100 Subject: [PATCH 05/10] Fix changelog and do unit conversion --- CHANGELOG.rst | 2 +- src/xclim/indices/_conversion.py | 38 +++++++++++++++++++------------- 2 files changed, 24 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 0c27dfc68..d5009fd4e 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -17,7 +17,7 @@ New features and enhancements * 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`). -* ``xclim.indicators.atmos`` now exports the FAO-56 Penman-Monteith equation as `fao_allen98` (:issue:`2004`, :pull:`2067`). +* ``xclim.indices.atmos`` now exports the FAO-56 Penman-Monteith equation as `fao_allen98` (:issue:`2004`, :pull:`2067`). Internal changes ^^^^^^^^^^^^^^^^ diff --git a/src/xclim/indices/_conversion.py b/src/xclim/indices/_conversion.py index d54dfca10..09458eb4c 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -1319,22 +1319,23 @@ def _get_D_from_M(time): # noqa: N802 ) -# TODO: How can I check exact units for inputs -# @declare_units( -# net_radiation="[MJ m-2 day-1]", -# tas="[degC]", -# wind="[m s-1]", -# es="[kPa]", -# ea="[kPa]", -# delta_svp="[kPa degC-1]", -# gamma="[kPa degC]", -# G="[MJ m-2 day-1]", -# ) -def fao_allen98(net_radiation, tas, wind, es, ea, delta_svp, gamma, G=0): +@declare_units( + net_radiation="[radiation]", + tas="[temperature]", + wind="[speed]", + es="[pressure]", + ea="[pressure]", + delta_svp="[pressure] / [temperature]", + # gamma="[pressure] * [temperature]", + # G="[radiation]", +) +def fao_allen98( + net_radiation, tas, wind, es, ea, delta_svp, gamma, G=0 # "0 MJ m-2 day-1" +): r""" FAO-56 Penman-Monteith equation. - Estimates reference evapotranspiration (ETo) from a hypothetical short grass reference surface ( + Estimates reference evapotranspiration from a hypothetical short grass reference surface ( height = 0.12m, surface resistance = 70 s m-1, albedo = 0.23 and a ``moderately dry soil surface resulting from about a weekly irrigation frequency``). Based on equation 6 in based on :cite:t:`allen_crop_1998`. @@ -1363,7 +1364,14 @@ def fao_allen98(net_radiation, tas, wind, es, ea, delta_svp, gamma, G=0): xarray.DataArray Potential Evapotranspiration from a hypothetical grass reference surface [mm day-1]. """ + net_radiation = convert_units_to(net_radiation, "MJ m-2 day-1") + wind = convert_units_to(wind, "m s-1") tasK = convert_units_to(tas, "K") + es = convert_units_to(es, "kPa") + ea = convert_units_to(ea, "kPa") + delta_svp = convert_units_to(delta_svp, "kPa degC-1") + # gamma = convert_units_to(gamma, "kPa degC-1") + # G = convert_units_to(G, "MJ m-2 day-1") a1 = 0.408 * delta_svp * (net_radiation - G) a2 = gamma * 900 / (tasK) * wind * (es - ea) a3 = delta_svp + (gamma * (1 + 0.34 * wind)) @@ -1632,10 +1640,10 @@ def potential_evapotranspiration( ) es = convert_units_to(es, "kPa") # mean actual vapour pressure [kPa] - ea = _hurs * es + ea = es * _hurs # slope of saturation vapour pressure curve [kPa degC-1] - delta = 4098 * es / (tas_m + 237.3) ** 2 + delta = (4098 * es / (tas_m + 237.3) ** 2).assign_attrs(units="kPa degC-1") # net radiation Rn = convert_units_to(rsds - rsus - (rlus - rlds), "MJ m-2 d-1") From 7e77690129e3ef21ee9589cf2f94f93ea3666ec5 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Thu, 6 Feb 2025 09:00:45 +0100 Subject: [PATCH 06/10] Fully unitted fao function --- src/xclim/indices/_conversion.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/xclim/indices/_conversion.py b/src/xclim/indices/_conversion.py index 09458eb4c..668b21d14 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -1326,12 +1326,10 @@ def _get_D_from_M(time): # noqa: N802 es="[pressure]", ea="[pressure]", delta_svp="[pressure] / [temperature]", - # gamma="[pressure] * [temperature]", - # G="[radiation]", + gamma="[pressure] [temperature]", + G="[radiation]", ) -def fao_allen98( - net_radiation, tas, wind, es, ea, delta_svp, gamma, G=0 # "0 MJ m-2 day-1" -): +def fao_allen98(net_radiation, tas, wind, es, ea, delta_svp, gamma, G="0 MJ m-2 day-1"): r""" FAO-56 Penman-Monteith equation. @@ -1370,13 +1368,13 @@ def fao_allen98( es = convert_units_to(es, "kPa") ea = convert_units_to(ea, "kPa") delta_svp = convert_units_to(delta_svp, "kPa degC-1") - # gamma = convert_units_to(gamma, "kPa degC-1") - # G = convert_units_to(G, "MJ m-2 day-1") + gamma = convert_units_to(gamma, "kPa degC") + G = convert_units_to(G, "MJ m-2 day-1") a1 = 0.408 * delta_svp * (net_radiation - G) a2 = gamma * 900 / (tasK) * wind * (es - ea) a3 = delta_svp + (gamma * (1 + 0.34 * wind)) - return (a1 + a2) / a3 + return ((a1 + a2) / a3).assign_attrs(units="mm day-1") @declare_units( @@ -1650,7 +1648,7 @@ def potential_evapotranspiration( P = 101.325 # Atmospheric pressure [kPa] gamma = 0.665e-03 * P # psychrometric const = C_p*P/(eps*lam) [kPa degC-1] - pet = fao_allen98(Rn, tas_m, wa2, es, ea, delta, gamma) + pet = fao_allen98(Rn, tas_m, wa2, es, ea, delta, f"{gamma} kPa degC") else: raise NotImplementedError(f"'{method}' method is not implemented.") From 68ac659e2bfc74e99c70d340d314c98894d2641e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= <71575674+coxipi@users.noreply.github.com> Date: Thu, 6 Feb 2025 09:43:07 -0500 Subject: [PATCH 07/10] add references section --- src/xclim/indices/_conversion.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/xclim/indices/_conversion.py b/src/xclim/indices/_conversion.py index c2bdf0b83..365d6846a 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -1361,6 +1361,10 @@ def fao_allen98(net_radiation, tas, wind, es, ea, delta_svp, gamma, G="0 MJ m-2 ------- xarray.DataArray Potential Evapotranspiration from a hypothetical grass reference surface [mm day-1]. + + References + ---------- + :cite:t:`allen_crop_1998` """ net_radiation = convert_units_to(net_radiation, "MJ m-2 day-1") wind = convert_units_to(wind, "m s-1") From 6c3f30a70b31630dffc39e98d8bbfc51e5b88cb3 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 6 Feb 2025 14:44:05 +0000 Subject: [PATCH 08/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- 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 365d6846a..e0d3dc3ce 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -1361,7 +1361,7 @@ def fao_allen98(net_radiation, tas, wind, es, ea, delta_svp, gamma, G="0 MJ m-2 ------- xarray.DataArray Potential Evapotranspiration from a hypothetical grass reference surface [mm day-1]. - + References ---------- :cite:t:`allen_crop_1998` From 09bc8ef6980b821bc48c080b163fbfb6ded5d497 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= <71575674+coxipi@users.noreply.github.com> Date: Thu, 6 Feb 2025 09:50:11 -0500 Subject: [PATCH 09/10] add fao_allen98 in _all_ --- src/xclim/indices/_conversion.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/xclim/indices/_conversion.py b/src/xclim/indices/_conversion.py index e0d3dc3ce..329171041 100644 --- a/src/xclim/indices/_conversion.py +++ b/src/xclim/indices/_conversion.py @@ -31,6 +31,7 @@ __all__ = [ "clausius_clapeyron_scaled_precipitation", + "fao_allen98", "heat_index", "humidex", "longwave_upwelling_radiation_from_net_downwelling", From 34dded5dc89401c2ba5a79dffc34ec235258cd83 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Thu, 6 Feb 2025 17:00:59 +0100 Subject: [PATCH 10/10] Update CHANGELOG.rst Co-authored-by: Pascal Bourgault --- CHANGELOG.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 3877145d3..5a371d480 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,7 +22,7 @@ New features and enhancements * 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`). -* ``xclim.indices.atmos`` now exports the FAO-56 Penman-Monteith equation as `fao_allen98` (:issue:`2004`, :pull:`2067`). +* New ``xclim.indices.fao_allen98``, exporting the FAO-56 Penman-Monteith equation for potential evapotranspiration (:issue:`2004`, :pull:`2067`). Internal changes ^^^^^^^^^^^^^^^^