Skip to content

Commit

Permalink
fix freq/offset/indexer attrs in SPI (#1538)
Browse files Browse the repository at this point in the history
<!--Please ensure the PR fulfills the following requirements! -->
<!-- If this is your first PR, make sure to add your details to the
AUTHORS.rst! -->
### Pull Request Checklist:
- [x] This PR addresses an already opened issue (for bug fixes /
features)
    - This PR fixes #xyz
- [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] CHANGES.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?

* Fixes wrong attributes in SPI/SPEI functions, allow to save files in
netcdf files.

### Does this PR introduce a breaking change?
No

### Other information:
  • Loading branch information
coxipi authored Dec 1, 2023
2 parents c748df5 + 68e3f88 commit d6f09d4
Show file tree
Hide file tree
Showing 4 changed files with 24 additions and 9 deletions.
7 changes: 4 additions & 3 deletions CHANGES.rst
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ Contributors to this version: Juliette Lavoie (:user:`juliettelavoie`), Pascal B

Announcements
^^^^^^^^^^^^^
* To circumvent issues stemming from changes to the frequency code convention in `pandas` v2.2, we have pinned `xarray` (<2023.11.0) and `pandas` (< 2.2) for this release. This change will be reverted in `xclim` v0.48.0 to support the newer versions (`xarray >= 2023.11.0` and `pandas >= 2.2`).
* To circumvent issues stemming from changes to the frequency code convention in `pandas` v2.2, we have pinned `xarray` (< 2023.11.0) and `pandas` (< 2.2) for this release. This change will be reverted in `xclim` v0.48.0 to support the newer versions (`xarray>= 2023.11.0` and `pandas>= 2.2`).
* `xclim` v0.47.0 will be the last release supporting Python3.8.

New features and enhancements
Expand All @@ -18,9 +18,10 @@ New features and enhancements

Bug fixes
^^^^^^^^^
* Fixed a bug with ``n_escore=-1`` in ``xclim.sdba.adjustment.NpdfTransform``. (:issue:`1515`, :pull:`1515`).
* In the documentation, fix the tooltips in the indicator search results. (:issue:`1524`, :pull:`1527`).
* Fixed a bug with ``n_escore=-1`` in ``xclim.sdba.adjustment.NpdfTransform``. (:issue:`1515`, :pull:`1516`).
* In the documentation, fixed the tooltips in the indicator search results. (:issue:`1524`, :pull:`1527`).
* If chunked inputs are passed to indicators ``mean_radiant_temperature`` and ``potential_evapotranspiration``, sub-calculations of the solar angle will also use the same chunks, instead of a single one of the same size as the data. (:issue:`1536`, :pull:`1542`).
* Fix wrong attributes in ``xclim.indices.standardized_precipitation_index``, ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`1537`, :pull:`1538`).

Internal changes
^^^^^^^^^^^^^^^^
Expand Down
2 changes: 2 additions & 0 deletions tests/test_indices.py
Original file line number Diff line number Diff line change
Expand Up @@ -653,6 +653,7 @@ def test_standardized_index_modularity(self, open_dataset):
dist=dist,
method=method,
offset="1 mm/d",
month=[2, 3],
)
spei1 = xci.standardized_precipitation_evapotranspiration_index(
wb.sel(time=slice("1998", "2000")), params=params
Expand All @@ -667,6 +668,7 @@ def test_standardized_index_modularity(self, open_dataset):
offset="1 mm/d",
cal_start="1950",
cal_end="1980",
month=[2, 3],
).sel(time=slice("1998", "2000"))

# In the previous computation, the first {window-1} values are NaN because the rolling is performed on the period [1998,2000].
Expand Down
13 changes: 10 additions & 3 deletions xclim/indices/_agro.py
Original file line number Diff line number Diff line change
Expand Up @@ -1174,6 +1174,7 @@ def standardized_precipitation_index(
... cal_end=cal_end,
... ) # Computing SPI-3 months using a gamma distribution for the fit
>>> # Fitting parameters can also be obtained ...
>>> from xclim.indices.stats import standardized_index_fit_params
>>> params = standardized_index_fit_params(
... pr.sel(time=slice(cal_start, cal_end)),
... freq="MS",
Expand All @@ -1189,7 +1190,12 @@ def standardized_precipitation_index(
:cite:cts:`mckee_relationship_1993`
"""
if params is not None and pr_cal is None:
freq, window = (params.attrs[s] for s in ["freq", "window"])
freq, window, indexer = (
params.attrs[s] for s in ["freq", "window", "time_indexer"]
)
# Unpack attrs to None and {} if needed
freq = None if freq == "" else freq
indexer = {} if indexer[0] == "" else {indexer[0]: indexer[1:]}
if cal_start or cal_end:
warnings.warn(
"Expected either `cal_{start|end}` or `params`, got both. The `params` input overrides other inputs."
Expand Down Expand Up @@ -1230,8 +1236,9 @@ def standardized_precipitation_index(
if paramsd != template.sizes:
params = params.broadcast_like(template)

spi = standardized_index(pr, params, **indexer)
spi = standardized_index(pr, params)
spi.attrs = params.attrs
spi.attrs["freq"] = freq or xarray.infer_freq(spi.time)
spi.attrs["units"] = ""
return spi

Expand Down Expand Up @@ -1324,7 +1331,7 @@ def standardized_precipitation_evapotranspiration_index(
"Proceeding with the value given in `params`."
)
offset = params_offset
offset = 0 if offset is None else convert_units_to(offset, wb, context="hydro")
offset = 0 if offset == "" else convert_units_to(offset, wb, context="hydro")
# Allowed distributions are constrained by the SPI function
if dist in ["gamma", "fisk"] and offset <= 0:
raise ValueError(
Expand Down
11 changes: 8 additions & 3 deletions xclim/indices/stats.py
Original file line number Diff line number Diff line change
Expand Up @@ -724,15 +724,20 @@ def standardized_index_fit_params(
)
params.attrs = {
"calibration_period": cal_range,
"freq": freq,
"freq": freq or "",
"window": window,
"scipy_dist": dist,
"method": method,
"group": group,
"time_indexer": indexer,
"units": "",
"offset": offset,
"offset": offset or "",
}
if indexer != {}:
method, args = indexer.popitem()
else:
method, args = "", []
params.attrs["time_indexer"] = (method, *args)

return params


Expand Down

0 comments on commit d6f09d4

Please sign in to comment.