Add potential evapotranspiration method based on Droogers and Allen (2002)

.zenodo.json

@@ -129,6 +129,11 @@
     {
       "name": "Castro, Dante",
       "affiliation": "Helmholtz-Zentrum Hereon, Geesthacht, Germany"
+    },
+    {
+      "name": "Diez-Sierra, Javier",
+      "affiliation": "Santander Meteorology Group, Instituto de Física de Cantabria (CSIC-UC), Santander, Spain",
+      "orcid": "0000-0001-9053-2542"
     }
   ],
   "keywords": [

AUTHORS.rst

@@ -42,3 +42,4 @@ Contributors
 * Christopher Whelan `@qwhelan <https://github.com/qwhelan>`_
 * Dante Castro <[email protected]> `@profesorpaiche <https://github.com/profesorpaiche>`_
 * Sascha Hofmann <[email protected]> `@saschahofmann <https://github.com/saschahofmann>`_
+* Javier Diez-Sierra <[email protected]> `@JavierDiezSierra <https://github.com/JavierDiezSierra>`_

CHANGES.rst

@@ -4,12 +4,16 @@ Changelog
 
 v0.49.0 (unreleased)
 --------------------
-Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user:`juliettelavoie`), David Huard (:user:`huard`), Gabriel Rondeau-Genesse (:user:`RondeauG`).
+Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user:`juliettelavoie`), David Huard (:user:`huard`), Gabriel Rondeau-Genesse (:user:`RondeauG`), Javier Diez-Sierra (:user:`JavierDiezSierra`).
 
 Announcements
 ^^^^^^^^^^^^^
 * `xclim` has migrated its development branch name from `master` to `main`. (:issue:`1667`, :pull:`1669`).
 
+New features and enhancements
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+* Indicator ``xclim.atmos.potential_evapotranspiration`` and indice ``xclim.indices.potential_evapotranspiration`` now accept a new value (`DA02`) for argument `method` implementing potential evapotranspiration based on Droogers and Allen (2002). (:issue:`1710`, :pull:`1723`).
+
 New indicators
 ^^^^^^^^^^^^^^
 * New ``snw_season_length`` and ``snd_season_length`` computing the duration between the start and the end of the snow season, both defined as the first day of a continuous period with snow above/under a threshold. Previous versions of these indicators were renamed ``snw_days_above`` and ``snd_days_above`` to better reflect what they computed : the number of days with snow above a given threshold (with no notion of continuity). (:issue:`1703`, :pull:`1708`).
@@ -19,6 +23,7 @@ Breaking changes
 ^^^^^^^^^^^^^^^^
 * The previously deprecated functions ``xclim.sdba.processing.construct_moving_yearly_window`` and ``xclim.sdba.processing.unpack_moving_yearly_window`` have been removed. These functions have been replaced by ``xclim.core.calendar.stack_periods`` and ``xclim.core.calendar.unstack_periods``. (:pull:`1717`).
 * Indicators ``snw_season_length`` and ``snd_season_length`` have been modified, see above.
+* The `hargeaves85`/`hg85` method for the ``potential_evapotranspiration`` indicator and indice has been modified for precision and consistency with recent academic literature. (:issue:`1710`, :pull:`1723`).
 
 Bug fixes
 ^^^^^^^^^

docs/references.bib

@@ -2114,3 +2114,16 @@ @article{brode_utci_2012
     doi = {10.1007/s00484-011-0454-1},
     url = {https://doi.org/10.1007/s00484-011-0454-1}
 }
+
+@article{droogers2002,
+	author = {Droogers, Peter and Allen, Richard G.},
+	title = {Estimating reference evapotranspiration under inaccurate data conditions},
+	year = {2002},
+	journal = {Irrigation and Drainage Systems},
+	volume = {16},
+	number = {1},
+	pages = {33 – 45},
+	doi = {10.1023/A:1015508322413},
+	url = {https://www.scopus.com/inward/record.uri?eid=2-s2.0-0036464359&doi=10.1023%2fA%3a1015508322413&partnerID=40&md5=7322aaa4c6874878f5b1dab3c73c1718},
+	type = {Article}
+}

tests/test_indices.py

@@ -3175,19 +3175,39 @@ def test_hargreaves(self, tasmin_series, tasmax_series, tas_series, lat_series):
 
         out = xci.potential_evapotranspiration(tn, tx, tm, lat=lat, method="HG85")
         np.testing.assert_allclose(
-            out.isel(lat=0, time=2), [3.962589 / 86400], rtol=1e-2
+            out.isel(lat=0, time=2), [4.030339 / 86400], rtol=1e-2
         )
 
-    def test_thornthwaite(self, tas_series, lat_series):
+    def test_droogersallen02(
+        self, tasmin_series, tasmax_series, tas_series, pr_series, lat_series
+    ):
         lat = lat_series([45])
-        time_std = date_range(
-            "1990-01-01", "1990-12-01", freq="MS", calendar="standard"
+        tn = tasmin_series(
+            np.array([0, 5, 10]), start="1990-01-01", freq="MS", units="degC"
+        ).expand_dims(lat=lat)
+        tx = tasmax_series(
+            np.array([10, 15, 20]), start="1990-01-01", freq="MS", units="degC"
+        ).expand_dims(lat=lat)
+        tg = tas_series(
+            np.array([5, 10, 15]), start="1990-01-01", freq="MS", units="degC"
+        ).expand_dims(lat=lat)
+        pr = pr_series(
+            np.array([30, 0, 60]), start="1990-01-01", freq="MS", units="mm/month"
+        ).expand_dims(lat=lat)
+
+        out = xci.potential_evapotranspiration(
+            tasmin=tn, tasmax=tx, tas=tg, pr=pr, lat=lat, method="DA02"
         )
-        tm = xr.DataArray(
-            np.ones((time_std.size, 1)),
-            dims=("time", "lat"),
-            coords={"time": time_std, "lat": lat},
-            attrs={"units": "degC"},
+        np.testing.assert_allclose(
+            out.isel(lat=0, time=2), [2.32659206 / 86400], rtol=1e-2
+        )
+
+    def test_thornthwaite(self, tas_series, lat_series):
+        lat = lat_series([45])
+        tm = (
+            tas_series(np.ones(12), start="1990-01-01", freq="MS", units="degC")
+            .expand_dims(lat=lat)
+            .assign_coords(lat=lat)
         )
 
         # find lat implicitly
@@ -3247,31 +3267,38 @@ def test_allen(
         )
 
 
-def test_water_budget_from_tas(pr_series, tasmin_series, tasmax_series, lat_series):
+def test_water_budget_from_tas(
+    pr_series, tasmin_series, tasmax_series, tas_series, lat_series
+):
     lat = lat_series([45])
     pr = pr_series(np.array([10, 10, 10])).expand_dims(lat=lat)
     pr.attrs["units"] = "mm/day"
-    tn = tasmin_series(np.array([0, 5, 10]) + K2C).expand_dims(lat=lat)
-    tx = tasmax_series(np.array([10, 15, 20]) + K2C).expand_dims(lat=lat)
+    tn = (
+        tasmin_series(np.array([0, 5, 10]) + K2C)
+        .expand_dims(lat=lat)
+        .assign_coords(lat=lat)
+    )
+    tx = (
+        tasmax_series(np.array([10, 15, 20]) + K2C)
+        .expand_dims(lat=lat)
+        .assign_coords(lat=lat)
+    )
 
     out = xci.water_budget(pr, tasmin=tn, tasmax=tx, lat=lat, method="BR65")
     np.testing.assert_allclose(out[0, 2], 6.138921 / 86400, rtol=2e-3)
 
     out = xci.water_budget(pr, tasmin=tn, tasmax=tx, lat=lat, method="HG85")
-    np.testing.assert_allclose(out[0, 2], 6.037411 / 86400, rtol=2e-3)
-
-    time_std = date_range("1990-01-01", "1990-12-01", freq="MS", calendar="standard")
-    tm = xr.DataArray(
-        np.ones((time_std.size, 1)),
-        dims=("time", "lat"),
-        coords={"time": time_std, "lat": lat},
-        attrs={"units": "degC"},
-    )
-    prm = xr.DataArray(
-        np.ones((time_std.size, 1)) * 10,
-        dims=("time", "lat"),
-        coords={"time": time_std, "lat": lat},
-        attrs={"units": "mm/day"},
+    np.testing.assert_allclose(out[0, 2], 5.969661 / 86400, rtol=2e-3)
+
+    tm = (
+        tas_series(np.ones(12), start="1990-01-01", freq="MS", units="degC")
+        .expand_dims(lat=lat)
+        .assign_coords(lat=lat)
+    )
+    prm = (
+        pr_series(np.ones(12) * 10, start="1990-01-01", freq="MS", units="mm/day")
+        .expand_dims(lat=lat)
+        .assign_coords(lat=lat)
     )
 
     # find lat implicitly

xclim/indices/_conversion.py

@@ -5,7 +5,6 @@
 import xarray as xr
 from numba import float32, float64, vectorize  # noqa
 
-from xclim.core.calendar import date_range
 from xclim.core.units import (
     amount2rate,
     convert_units_to,
@@ -1270,6 +1269,23 @@ def clausius_clapeyron_scaled_precipitation(
     return pr_out
 
 
+def _get_D_from_M(time):
+    start = time[0].dt.strftime("%Y-%m-01").item()
+    yrmn = time[-1].dt.strftime("%Y-%m").item()
+    end = f"{yrmn}-{time[-1].dt.daysinmonth.item()}"
+    return xr.DataArray(
+        xr.date_range(
+            start,
+            end,
+            freq="D",
+            calendar=time.dt.calendar,
+            use_cftime=(time.dtype == "O"),
+        ),
+        dims="time",
+        name="time",
+    )
+
+
 @declare_units(
     tasmin="[temperature]",
     tasmax="[temperature]",
@@ -1281,6 +1297,7 @@ def clausius_clapeyron_scaled_precipitation(
     rlds="[radiation]",
     rlus="[radiation]",
     sfcWind="[speed]",
+    pr="[precipitation]",
 )
 def potential_evapotranspiration(
     tasmin: xr.DataArray | None = None,
@@ -1293,6 +1310,7 @@ def potential_evapotranspiration(
     rlds: xr.DataArray | None = None,
     rlus: xr.DataArray | None = None,
     sfcWind: xr.DataArray | None = None,
+    pr: xr.DataArray | None = None,
     method: str = "BR65",
     peta: float = 0.00516409319477,
     petb: float = 0.0874972822289,
@@ -1324,7 +1342,9 @@ def potential_evapotranspiration(
         Surface Upwelling Longwave Radiation
     sfcWind : xarray.DataArray, optional
         Surface wind velocity (at 10 m)
-    method : {"baierrobertson65", "BR65", "hargreaves85", "HG85", "thornthwaite48", "TW48", "mcguinnessbordne05", "MB05", "allen98", "FAO_PM98"}
+    pr : xarray.DataArray
+        Mean daily precipitation flux.
+    method : {"baierrobertson65", "BR65", "hargreaves85", "HG85", "thornthwaite48", "TW48", "mcguinnessbordne05", "MB05", "allen98", "FAO_PM98", "droogersallen02", "DA02"}
         Which method to use, see notes.
     peta : float
         Used only with method MB05 as :math:`a` for calculation of PET, see Notes section.
@@ -1352,6 +1372,8 @@ def potential_evapotranspiration(
       (optional: tas can be given instead of tasmin and tasmax).
     - "allen98" or "FAO_PM98", based on :cite:t:`allen_crop_1998`. Modification of Penman-Monteith method.
       Requires tasmin and tasmax, relative humidity, radiation flux and wind speed (10 m wind will be converted to 2 m).
+    - "droogersallen02" or "DA02", based on :cite:t:`droogers2002`.
+      Requires tasmin, tasmax and precipitation, monthly [MS] or daily [D] freq. (optional: tas can be given in addition of tasmin and tasmax).
 
     The McGuinness-Bordne :cite:p:`mcguinness_comparison_1972` equation is:
 
@@ -1364,13 +1386,14 @@ def potential_evapotranspiration(
     with :math:`a=0.0147` and :math:`b=0.07353`. The default parameters used here are calibrated for the UK,
     using the method described in :cite:t:`tanguy_historical_2018`.
 
-    Methods "BR65", "HG85" and "MB05" use an approximation of the extraterrestrial radiation.
+    Methods "BR65", "HG85", "MB05" and "DA02" use an approximation of the extraterrestrial radiation.
     See :py:func:`~xclim.indices._helpers.extraterrestrial_solar_radiation`.
 
     References
     ----------
-    :cite:cts:`baier_estimation_1965,george_h_hargreaves_reference_1985,tanguy_historical_2018,thornthwaite_approach_1948,mcguinness_comparison_1972,allen_crop_1998`
-    """
+    :cite:cts:`baier_estimation_1965,george_h_hargreaves_reference_1985,tanguy_historical_2018,thornthwaite_approach_1948,mcguinness_comparison_1972,allen_crop_1998,droogers2002`
+    """  # noqa: E501
+    # ^ Ignoring "line too long" as it comes from un-splittable constructs
     if lat is None:
         lat = _gather_lat(tasmin if tas is None else tas)
 
@@ -1397,17 +1420,49 @@ def potential_evapotranspiration(
         else:
             tas = convert_units_to(tas, "degC")
 
-        lv = 2.5  # MJ/kg
-
         ra = extraterrestrial_solar_radiation(
             tasmin.time, lat, chunks=tasmin.chunksizes
         )
         ra = convert_units_to(ra, "MJ m-2 d-1")
 
+        # Is used to convert the radiation to evaporation equivalents in mm (kg/MJ)
+        ra = ra * 0.408
+
         # Hargreaves and Samani (1985) formula
-        out = (0.0023 * ra * (tas + 17.8) * (tasmax - tasmin) ** 0.5) / lv
+        out = 0.0023 * ra * (tas + 17.8) * (tasmax - tasmin) ** 0.5
         out = out.clip(0)
 
+    elif method in ["droogersallen02", "DA02"]:
+        tasmin = convert_units_to(tasmin, "degC")
+        tasmax = convert_units_to(tasmax, "degC")
+        pr = convert_units_to(pr, "mm/month", context="hydro")
+        if tas is None:
+            tas = (tasmin + tasmax) / 2
+        else:
+            tas = convert_units_to(tas, "degC")
+
+        tasmin = tasmin.resample(time="MS").mean()
+        tasmax = tasmax.resample(time="MS").mean()
+        tas = tas.resample(time="MS").mean()
+        pr = pr.resample(time="MS").mean()
+
+        # Monthly accumulated radiation
+        time_d = _get_D_from_M(tasmin.time)
+        ra = extraterrestrial_solar_radiation(time_d, lat)
+        ra = convert_units_to(ra, "MJ m-2 d-1")
+        ra = ra.resample(time="MS").sum()
+        # Is used to convert the radiation to evaporation equivalents in mm (kg/MJ)
+        ra = ra * 0.408
+
+        tr = tasmax - tasmin
+        tr = tr.where(tr > 0, 0)
+
+        # Droogers and Allen (2002) formula
+        ab = tr - 0.0123 * pr
+        out = 0.0013 * ra * (tas + 17.0) * ab**0.76
+        out = xr.where(np.isnan(ab**0.76), 0, out)
+        out = out.clip(0)  # mm/month
+
     elif method in ["mcguinnessbordne05", "MB05"]:
         if tas is None:
             tasmin = convert_units_to(tasmin, "degC")
@@ -1441,30 +1496,9 @@ def potential_evapotranspiration(
         tas = tas.clip(0)
         tas = tas.resample(time="MS").mean(dim="time")
 
-        start = "-".join(
-            [
-                str(tas.time[0].dt.year.values),
-                f"{tas.time[0].dt.month.values:02d}",
-                "01",
-            ]
-        )
-
-        end = "-".join(
-            [
-                str(tas.time[-1].dt.year.values),
-                f"{tas.time[-1].dt.month.values:02d}",
-                str(tas.time[-1].dt.daysinmonth.values),
-            ]
-        )
-
-        time_v = xr.DataArray(
-            date_range(start, end, freq="D", calendar="standard"),
-            dims="time",
-            name="time",
-        )
-
         # Thornthwaite measures half-days
-        dl = day_lengths(time_v, lat) / 12
+        time_d = _get_D_from_M(tas.time)
+        dl = day_lengths(time_d, lat) / 12
         dl_m = dl.resample(time="MS").mean(dim="time")
 
         # annual heat index

xclim/sdba/adjustment.py

@@ -74,7 +74,7 @@ def __init__(self, *args, _trained=False, **kwargs):
             super().__init__(*args, **kwargs)
         else:
             raise ValueError(
-                "As of xclim 0.29, Adjustment object should be initialized through their `train` or  `adjust` methods."
+                "As of xclim 0.29, Adjustment object should be initialized through their `train` or `adjust` methods."
             )
 
     @classmethod
@@ -116,8 +116,8 @@ def _check_inputs(cls, *inputs, group):
             "default" in calendars or "standard" in calendars
         ):
             warn(
-                "Strange results could be returned when using dayofyear grouping "
-                "on data defined in the proleptic_gregorian calendar "
+                "Strange results could be returned when using `dayofyear` grouping "
+                "on data defined in the 'proleptic_gregorian' calendar."
             )
 
     @classmethod