From b64edb58a857db38dbbdf8cc4662035ed8c0d985 Mon Sep 17 00:00:00 2001 From: JavierDiezSierra Date: Mon, 22 Apr 2024 11:59:21 +0200 Subject: [PATCH 01/13] New method to calculate potential evapotranspiration based on Droogers and Allen 2002 --- docs/references.bib | 16 +++++++++++ xclim/indices/_conversion.py | 52 +++++++++++++++++++++++++++++++++--- 2 files changed, 65 insertions(+), 3 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index 579a3ef56..c74623835 100644 --- a/docs/references.bib +++ b/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} +} diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py index 0c6798cb6..7eda990ba 100644 --- a/xclim/indices/_conversion.py +++ b/xclim/indices/_conversion.py @@ -1281,6 +1281,7 @@ def clausius_clapeyron_scaled_precipitation( rlds="[radiation]", rlus="[radiation]", sfcWind="[speed]", + pr="[precipitation]", ) def potential_evapotranspiration( tasmin: xr.DataArray | None = None, @@ -1293,6 +1294,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 +1326,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 +1356,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,12 +1370,12 @@ 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` """ if lat is None: lat = _gather_lat(tasmin if tas is None else tas) @@ -1408,6 +1414,46 @@ 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") + if tas is None: + tas = (tasmin + tasmax) / 2 + else: + tas = convert_units_to(tas, "degC") + + # Monthly accumulated radiation + if xr.infer_freq(tasmin['time']) == 'D': + ra = extraterrestrial_solar_radiation( + tasmin.time, lat, chunks=tasmin.chunksizes + ) + tasmin = tasmin.resample(time="MS").mean(dim="time") + tasmax = tasmax.resample(time="MS").mean(dim="time") + tas = tas.resample(time="MS").mean(dim="time") + pr = pr.resample(time="MS").mean(dim="time") + + elif xr.infer_freq(tasmin['time']) == 'MS': + tasmin_day = tasmin.resample(time='D').ffill() + ra = extraterrestrial_solar_radiation( + tasmin_day.time, lat, chunks=tasmin_day.chunksizes + ) + + ra = convert_units_to(ra, "MJ m-2 d-1") + ra = ra.resample(time="MS").sum(dim="time") + ra = ra * 0.408 # Is used to convert the radiation to evaporation equivalents in mm (kg/MJ) + + 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") From d9c0093d9a33f360dce07c99264a1175e9d363fd Mon Sep 17 00:00:00 2001 From: JavierDiezSierra Date: Mon, 22 Apr 2024 12:11:20 +0200 Subject: [PATCH 02/13] New contributor added to the AUTHORS.rst file --- AUTHORS.rst | 1 + 1 file changed, 1 insertion(+) 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 `_ From 1797e9e140b3cdcb42e1b0efd7f9e16943cfa268 Mon Sep 17 00:00:00 2001 From: JavierDiezSierra Date: Mon, 22 Apr 2024 12:43:38 +0200 Subject: [PATCH 03/13] Changes added to the CHANGES.rst file --- CHANGES.rst | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index 1680d3857..da38e3119 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.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`). From 842999b04159942fac95a7206eb670445430097d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 22 Apr 2024 10:46:12 +0000 Subject: [PATCH 04/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xclim/indices/_conversion.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py index 7eda990ba..6e055a13f 100644 --- a/xclim/indices/_conversion.py +++ b/xclim/indices/_conversion.py @@ -1422,9 +1422,9 @@ def potential_evapotranspiration( tas = (tasmin + tasmax) / 2 else: tas = convert_units_to(tas, "degC") - + # Monthly accumulated radiation - if xr.infer_freq(tasmin['time']) == 'D': + if xr.infer_freq(tasmin["time"]) == "D": ra = extraterrestrial_solar_radiation( tasmin.time, lat, chunks=tasmin.chunksizes ) @@ -1433,26 +1433,26 @@ def potential_evapotranspiration( tas = tas.resample(time="MS").mean(dim="time") pr = pr.resample(time="MS").mean(dim="time") - elif xr.infer_freq(tasmin['time']) == 'MS': - tasmin_day = tasmin.resample(time='D').ffill() + elif xr.infer_freq(tasmin["time"]) == "MS": + tasmin_day = tasmin.resample(time="D").ffill() ra = extraterrestrial_solar_radiation( tasmin_day.time, lat, chunks=tasmin_day.chunksizes ) ra = convert_units_to(ra, "MJ m-2 d-1") ra = ra.resample(time="MS").sum(dim="time") - ra = ra * 0.408 # Is used to convert the radiation to evaporation equivalents in mm (kg/MJ) + ra = ( + ra * 0.408 + ) # Is used to convert the radiation to evaporation equivalents in mm (kg/MJ) tr = tasmax - tasmin - tr = tr.where(tr>0, 0) + 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 = 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 + out = out.clip(0) # mm/month elif method in ["mcguinnessbordne05", "MB05"]: if tas is None: From b8f4f65a8ff2f636910227865a1d548aff7a77db Mon Sep 17 00:00:00 2001 From: JavierDiezSierra Date: Mon, 22 Apr 2024 14:04:42 +0200 Subject: [PATCH 05/13] Add author to zenodo.json --- .zenodo.json | 5 +++++ 1 file changed, 5 insertions(+) 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": [ From df746842294c940dae068df4f2bd1d83aac454a9 Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Tue, 23 Apr 2024 09:36:04 -0400 Subject: [PATCH 06/13] Remove need for exact freq matching --- xclim/indices/_conversion.py | 77 +++++++++++++++--------------------- 1 file changed, 32 insertions(+), 45 deletions(-) diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py index 6e055a13f..d3375f182 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]", @@ -1376,7 +1392,8 @@ def potential_evapotranspiration( References ---------- :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) @@ -1417,33 +1434,24 @@ def potential_evapotranspiration( elif method in ["droogersallen02", "DA02"]: tasmin = convert_units_to(tasmin, "degC") tasmax = convert_units_to(tasmax, "degC") - pr = convert_units_to(pr, "mm/month") + pr = convert_units_to(pr, "mm/month", context="hydro") if tas is None: tas = (tasmin + tasmax) / 2 else: tas = convert_units_to(tas, "degC") - # Monthly accumulated radiation - if xr.infer_freq(tasmin["time"]) == "D": - ra = extraterrestrial_solar_radiation( - tasmin.time, lat, chunks=tasmin.chunksizes - ) - tasmin = tasmin.resample(time="MS").mean(dim="time") - tasmax = tasmax.resample(time="MS").mean(dim="time") - tas = tas.resample(time="MS").mean(dim="time") - pr = pr.resample(time="MS").mean(dim="time") - - elif xr.infer_freq(tasmin["time"]) == "MS": - tasmin_day = tasmin.resample(time="D").ffill() - ra = extraterrestrial_solar_radiation( - tasmin_day.time, lat, chunks=tasmin_day.chunksizes - ) + 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(dim="time") - ra = ( - ra * 0.408 - ) # Is used to convert the radiation to evaporation equivalents in mm (kg/MJ) + 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) @@ -1487,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 From 6f4a29fa89e3681881963329f46390f09aaca416 Mon Sep 17 00:00:00 2001 From: JavierDiezSierra Date: Wed, 24 Apr 2024 14:14:28 +0200 Subject: [PATCH 07/13] test for da02 method --- tests/test_indices.py | 35 +++++++++++++++++++++++++++++++++++ xclim/indices/_conversion.py | 9 +++++---- 2 files changed, 40 insertions(+), 4 deletions(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index 9fb3896cf..e86366227 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -3178,6 +3178,41 @@ def test_hargreaves(self, tasmin_series, tasmax_series, tas_series, lat_series): out.isel(lat=0, time=2), [3.962589 / 86400], rtol=1e-2 ) + 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-03-01", freq="MS", calendar="standard" + ) + tn = xr.DataArray( + np.array([0, 5, 10]), + dims=("time", "lat"), + coords={"time": time_std, "lat": lat}, + attrs={"units": "degC"}, + ) + tx = xr.DataArray( + np.array([10, 15, 20]), + dims=("time", "lat"), + coords={"time": time_std, "lat": lat}, + attrs={"units": "degC"}, + ) + tm = xr.DataArray( + np.array([5, 10, 15]), + dims=("time", "lat"), + coords={"time": time_std, "lat": lat}, + attrs={"units": "degC"}, + ) + pr = xr.DataArray( + np.array([30, 0, 60]), + dims=("time", "lat"), + coords={"time": time_std, "lat": lat}, + attrs={"units": "mm/month"}, + ) + + out = xci.potential_evapotranspiration(tn, tx, tm, pr=pr, lat=lat, method="DA02") + 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]) time_std = date_range( diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py index d3375f182..09c6a230d 100644 --- a/xclim/indices/_conversion.py +++ b/xclim/indices/_conversion.py @@ -1420,15 +1420,16 @@ 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"]: @@ -1447,7 +1448,7 @@ def potential_evapotranspiration( # Monthly accumulated radiation time_d = _get_D_from_M(tasmin.time) - ra = extraterrestrial_solar_radiation(time_d, lat, chunks=tasmin.chunksize) + 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) From 1f9bf22b94c472952bf89ecc8ad7ca888bbfa7c9 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 24 Apr 2024 12:15:38 +0000 Subject: [PATCH 08/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/test_indices.py | 8 ++++++-- xclim/indices/_conversion.py | 2 +- xclim/sdba/adjustment.py | 2 +- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index e86366227..3e53aa587 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -3178,7 +3178,9 @@ def test_hargreaves(self, tasmin_series, tasmax_series, tas_series, lat_series): out.isel(lat=0, time=2), [3.962589 / 86400], rtol=1e-2 ) - def test_droogersallen02(self, tasmin_series, tasmax_series, tas_series, pr_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-03-01", freq="MS", calendar="standard" @@ -3208,7 +3210,9 @@ def test_droogersallen02(self, tasmin_series, tasmax_series, tas_series, pr_seri attrs={"units": "mm/month"}, ) - out = xci.potential_evapotranspiration(tn, tx, tm, pr=pr, lat=lat, method="DA02") + out = xci.potential_evapotranspiration( + tn, tx, tm, pr=pr, lat=lat, method="DA02" + ) np.testing.assert_allclose( out.isel(lat=0, time=2), [2.32659206 / 86400], rtol=1e-2 ) diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py index 09c6a230d..50cc2ad32 100644 --- a/xclim/indices/_conversion.py +++ b/xclim/indices/_conversion.py @@ -1429,7 +1429,7 @@ def potential_evapotranspiration( ra = ra * 0.408 # Hargreaves and Samani (1985) formula - out = (0.0023 * ra * (tas + 17.8) * (tasmax - tasmin) ** 0.5) + out = 0.0023 * ra * (tas + 17.8) * (tasmax - tasmin) ** 0.5 out = out.clip(0) elif method in ["droogersallen02", "DA02"]: diff --git a/xclim/sdba/adjustment.py b/xclim/sdba/adjustment.py index 17d0d2068..253fcd457 100644 --- a/xclim/sdba/adjustment.py +++ b/xclim/sdba/adjustment.py @@ -109,7 +109,7 @@ def _check_inputs(cls, *inputs, group): ): raise ValueError( "Inputs have different multivariate coordinates " - f"({set(mv.name for mv in mvcrds)})." + f"({{mv.name for mv in mvcrds}})." ) if group.prop == "dayofyear" and ( From e8487565d11c5dfa944ec93c6cf763f7a4d49435 Mon Sep 17 00:00:00 2001 From: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> Date: Wed, 24 Apr 2024 16:42:20 -0400 Subject: [PATCH 09/13] Update CHANGES.rst --- CHANGES.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index 50ab33963..65521d8c0 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -12,7 +12,7 @@ Announcements 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`). +* 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 ^^^^^^^^^^^^^^ From aeb93a30866f91f5e6c3126920f1f822c8265581 Mon Sep 17 00:00:00 2001 From: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> Date: Wed, 24 Apr 2024 16:42:46 -0400 Subject: [PATCH 10/13] Update docs/references.bib --- docs/references.bib | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index c74623835..d8e0e8082 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -2125,8 +2125,5 @@ @article{droogers2002 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} + type = {Article} } From 474837f6a20ef3cf6d9d6915b2b8ec3e1ac69b20 Mon Sep 17 00:00:00 2001 From: Zeitsperre <10819524+Zeitsperre@users.noreply.github.com> Date: Wed, 24 Apr 2024 17:05:25 -0400 Subject: [PATCH 11/13] fix formatting --- xclim/sdba/adjustment.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/xclim/sdba/adjustment.py b/xclim/sdba/adjustment.py index 253fcd457..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 @@ -107,17 +107,17 @@ def _check_inputs(cls, *inputs, group): not all(mvcrds[0].equals(mv) for mv in mvcrds[1:]) or len(mvcrds) != len(inputs) ): + coords = {mv.name for mv in mvcrds} raise ValueError( - "Inputs have different multivariate coordinates " - f"({{mv.name for mv in mvcrds}})." + f"Inputs have different multivariate coordinates: {', '.join(coords)}." ) if group.prop == "dayofyear" and ( "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 From abb740ae9fa2d283316954c2cadc4f6e1f51dbd6 Mon Sep 17 00:00:00 2001 From: Zeitsperre <10819524+Zeitsperre@users.noreply.github.com> Date: Wed, 24 Apr 2024 17:13:45 -0400 Subject: [PATCH 12/13] add breaking change to CHANGES.rst --- CHANGES.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGES.rst b/CHANGES.rst index 65521d8c0..f6a0a7f32 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -23,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 ^^^^^^^^^ From 59f5d59528da0817d3a2b360499a4dbef18896a9 Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Mon, 29 Apr 2024 13:30:31 -0400 Subject: [PATCH 13/13] Fix tests and uniformize style --- tests/test_indices.py | 96 +++++++++++++++++++------------------------ 1 file changed, 42 insertions(+), 54 deletions(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index 3e53aa587..decd4b5b0 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -3175,43 +3175,28 @@ 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_droogersallen02( self, tasmin_series, tasmax_series, tas_series, pr_series, lat_series ): lat = lat_series([45]) - time_std = date_range( - "1990-01-01", "1990-03-01", freq="MS", calendar="standard" - ) - tn = xr.DataArray( - np.array([0, 5, 10]), - dims=("time", "lat"), - coords={"time": time_std, "lat": lat}, - attrs={"units": "degC"}, - ) - tx = xr.DataArray( - np.array([10, 15, 20]), - dims=("time", "lat"), - coords={"time": time_std, "lat": lat}, - attrs={"units": "degC"}, - ) - tm = xr.DataArray( - np.array([5, 10, 15]), - dims=("time", "lat"), - coords={"time": time_std, "lat": lat}, - attrs={"units": "degC"}, - ) - pr = xr.DataArray( - np.array([30, 0, 60]), - dims=("time", "lat"), - coords={"time": time_std, "lat": lat}, - attrs={"units": "mm/month"}, - ) + 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( - tn, tx, tm, pr=pr, lat=lat, method="DA02" + tasmin=tn, tasmax=tx, tas=tg, pr=pr, lat=lat, method="DA02" ) np.testing.assert_allclose( out.isel(lat=0, time=2), [2.32659206 / 86400], rtol=1e-2 @@ -3219,14 +3204,10 @@ def test_droogersallen02( def test_thornthwaite(self, tas_series, lat_series): lat = lat_series([45]) - 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"}, + 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 @@ -3286,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