Add chill portion and chill units

[saschahofmann]

Closes #1753

TODO:

• Tests for the changes have been added (for bug fixes / features)
  • (If applicable) Documentation has been added / updated (for bug fixes / features)
  • Test make_hourly_temperature
  • Test chill_portions
  • Test chill_units
• CHANGELOG.rst has been updated (with summary of main changes)
  • Link to issue (:issue:number) and pull request (:pull:number) has been added
• Check make_hourly_temperature is working with other calendars than 360_day

What kind of change does this PR introduce?

• Adds a chill_portion indicator based on the Dynamic model
• Adds a chill_unit indicator based on the Utah model
• Adds a helper function to estimate hourly temperatures from daily tasmin and tasmax

[Zeitsperre]

This looks great, thanks so much! I'll leave another reviewer to give the OK, though.

[saschahofmann]

I am still trying to figure out what the functions should return. Right now, the chill_portions returns an aggregated value per year but chill_units a value per hour that can be aggregated at will. Ideally, they would work the same but chill_portions are significantly harder to aggregate since an hours value depends on the previous hour's value. This also has implications on computing chill portions for a season. For chill units, you can simply select the range of days or hours afterwards and aggregate but for chill portions you need to select the period before and I could imagine this messing with the ResamplingIndicator class that fills missing data (@aulemahal ?) so I guess I need to do it inside the function?

As you can see in the current implementation of chill_portions, I am doing a yearly grouping to get the values per year but this would fail if you're interested in a season crossing new years (e.g. Nov-Feb). We do a trick where we shift the day so the season is within a year but that's only easy if you assume a 360_day calendar.

Update: chill_units now resamples to a year as well and sums the values but the question above remains: what do you think would be the right mechanism to operate on only a selected time range?

[saschahofmann]

I am also not sure what the IndexingIndicator class exactly is doing?

injects "indexer" kwargs to subset the inputs before computation.

sounds a little like what I am suggesting to do but I couldn't find how exactly it works?

In general, I was expecting these classes to rename and add the units specified but it seems they working more as a check? e.g. I made cu the varname for the chill units index but after the operation it's still tas. I guess I need to manually rename it?

And finally about testing: I am currently adding tests to the TestAgroclimaticIndices class in test_indices.py but I can see that there are also tests in test_temperature.py. What's the use case for the two?

[aulemahal]

By "over years", you mean periods like "YS-JUL" that contain the change in calendar year ? I'm not I see how this is problematic for chill_portions as you don't explicitly make use of the datetime information ?

[saschahofmann]

Haha. Its only problematic if you had no idea that something like YS-JUL existed ... Alright, that solves a lot. Will implement the indexer solution for chill_portions and test it.

[saschahofmann]

also linking @eikeluedeling who formalised the method for chill portions.

[aulemahal]

This looks good to go now! Is there a reason the PR was closed ?

[saschahofmann]

I must have misclicked! Reopening

[aulemahal]

I'll try to contribute the french translations soon.

[aulemahal]

Sorry for the delay. I think these translations are good enough. I have asked the opinion of a colleague but we could merge away when the tests pass.

[aulemahal]

For the record, I have found only one french language document discussing the chill units and portions, and they didn't translate the names : https://archipel.uqam.ca/7003/1/M13555.pdf
My colleague nonetheless thinks we should translate them, which I did in the last commit. Especially in the context of xclim, the translated attributes never appear by themselves, so there's no ambiguity.
But I'm guessing you don't really care about all that, and I can understand. This requirement of french translations makes less and less sense as the number of non-Ouranosinc contributors grows...

Hopefully this last commit also triggers the correct set of CI jobs and we can merge.

[coveralls]

coverage: 89.435% (+0.06%) from 89.374%
when pulling 427a313 on saschahofmann:chilling
into 7502953 on Ouranosinc:main.

[aulemahal]

Thanks @saschahofmann !

[saschahofmann]

Awesome thank you so much !

[baptistehamon]

Hi, I've been using the Utah chill unit model for my research but my code is a mess and time-consuming so I'm happy to see it implemented in xclim (thanks @saschahofmann for that)! I did some tests with my data today and wanted to inform you about two things I've noticed.

Memory issue

First, it's more a warning than an issue. Because chill_unit uses hourly data, some memory errors can occur using an important dataset while the computation of other indices or indicators with the same dataset (but daily resolution) is fine. This can be handled easily by iterating over the years and it's what I've done. However, I think it could be relevant to add a warning message when using make_hourly_temperature. What do you think about this?

Output bug

The second is make_hourly_temperature adds one more day to compute the hourly temperature but does not remove it after the computation, thus returning an array with a shape different to the inputs. See below for data corresponding to a complete year with a shape (time: 365, lat: 56, lon: 54) as input:

>>> data_yr['tasmin'].shape, data_yr['tasmax'].shape
((365, 56, 54), (365, 56, 54))
>>> make_hourly_temperature(data_yr['tasmin'], data_yr['tasmax']).resample(time='D').mean().shape
(366, 56, 54)

This is particularly annoying if the input data ends on the last day of the year and a resampling YS is done, creating a year at the end considering only one day.

>>> hourly_temp = make_hourly_temperature(data_yr['tasmin'], data_yr['tasmax'])
>>> cu = chill_units(hourly_temp, freq='YS')
>>> print(cu)
<xarray.DataArray (time: 2, lat: 56, lon: 54)> Size: 48kB
dask.array<stack, shape=(2, 56, 54), dtype=float64, chunksize=(1, 56, 54), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 448B -47.38 -47.12 -46.88 ... -34.12 -33.88 -33.62
  * lon      (lon) float64 432B 166.1 166.4 166.6 166.9 ... 178.9 179.1 179.4
  * time     (time) datetime64[ns] 16B 2100-01-01 2101-01-01
Attributes:
    units:    1

I hope my explanations are clear enough! Thanks !