Exportable function for fao_allen98 (#2067)

### Pull Request Checklist: - [x] This PR addresses an already opened issue (for bug fixes / features) - This PR fixes #2004 - Also see discussion #1853 for more information - [x] Tests for the changes have been added (for bug fixes / features) - [x] (If applicable) Documentation has been added / updated (for bug fixes / features) - [x] CHANGELOG.rst has been updated (with summary of main changes) - [x] Link to issue (:issue:`number`) and pull request (:pull:`number`) has been added ### What kind of change does this PR introduce? * Adds a new function `fao_allen98` that implements the FAO-56 Penman-Monteith equation and is used inside `potential_evapotranspiration` for the `allen98` method. This allows providing some variables directly instead of using estimates e.g. tas mean or the different vapour pressures can be taken from datasets directly instead of being estimated. See discussion #1853 or issue #2004. * Fixes a bug where hurs provided in `%` lead to wrong results for PET when using `allen98` ### Does this PR introduce a breaking change? No. ### Questions: I would like to fix the units of the input for the new function exactly. Is there an easy way to do this? And how can I provide a scalar value as a default with units?
Ouranosinc · Feb 6, 2025 · 9845b7f · 9845b7f
2 parents f615ef7 + 34dded5
commit 9845b7f
Show file tree

Hide file tree

Showing 4 changed files with 84 additions and 13 deletions.
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -22,6 +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`).
+* New ``xclim.indices.fao_allen98``, exporting the FAO-56 Penman-Monteith equation for potential evapotranspiration (:issue:`2004`, :pull:`2067`).
 
 Internal changes
 ^^^^^^^^^^^^^^^^
@@ -32,6 +33,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
@@ -31,6 +31,7 @@
 
 __all__ = [
     "clausius_clapeyron_scaled_precipitation",
+    "fao_allen98",
     "heat_index",
     "humidex",
     "longwave_upwelling_radiation_from_net_downwelling",
@@ -1319,6 +1320,68 @@ def _get_D_from_M(time):  # noqa: N802
     )
 
 
+@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 MJ m-2 day-1"):
+    r"""
+    FAO-56 Penman-Monteith equation.
+
+    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`.
+
+    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 : float, 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].
+
+    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")
+    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")
+    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).assign_attrs(units="mm day-1")
+
+
 @declare_units(
     tasmin="[temperature]",
     tasmax="[temperature]",
@@ -1569,7 +1632,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.")
 
@@ -1586,25 +1649,17 @@ 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")
 
-            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]
 
-            # 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 = (
-                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, f"{gamma} kPa degC")
 
     else:
         raise NotImplementedError(f"'{method}' method is not implemented.")

diff --git 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
@@ -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)