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.indices.potential_evapotranspiration`` now accepts 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`).

docs/references.bib

@@ -2114,3 +2114,19 @@ @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},
+	publication_stage = {Final},
+	source = {Scopus},
+	note = {Cited by: 747}
+}

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)
 
@@ -1408,6 +1431,37 @@ def potential_evapotranspiration(
         out = (0.0023 * ra * (tas + 17.8) * (tasmax - tasmin) ** 0.5) / lv
         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, chunks=tasmin.chunksize)
+        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 +1495,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