diff --git a/.zenodo.json b/.zenodo.json index 9242df632..bc29d900e 100644 --- a/.zenodo.json +++ b/.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": [ diff --git a/AUTHORS.rst b/AUTHORS.rst index b8fd617c7..17582a6b6 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -42,3 +42,4 @@ Contributors * Christopher Whelan `@qwhelan `_ * Dante Castro `@profesorpaiche `_ * Sascha Hofmann `@saschahofmann `_ +* Javier Diez-Sierra `@JavierDiezSierra `_ diff --git a/CHANGES.rst b/CHANGES.rst index bf9636021..f6a0a7f32 100644 --- a/CHANGES.rst +++ b/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 ^^^^^^^^^ diff --git a/docs/references.bib b/docs/references.bib index 579a3ef56..d8e0e8082 100644 --- a/docs/references.bib +++ b/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} +} diff --git a/tests/test_indices.py b/tests/test_indices.py index 9fb3896cf..decd4b5b0 100644 --- a/tests/test_indices.py +++ b/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 diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py index 0c6798cb6..50cc2ad32 100644 --- a/xclim/indices/_conversion.py +++ b/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 diff --git a/xclim/sdba/adjustment.py b/xclim/sdba/adjustment.py index 7a2a5f3d7..8436d278b 100644 --- a/xclim/sdba/adjustment.py +++ b/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