Chill Units from Utah Model, Chill portions from Dynamic Model

[saschahofmann]

Addressing a Problem?

I am currently implementing an indicator that computes utah model chilling units based of Richardson EA, Seeley SD, Walker DR (1974) A model for estimating the completion of rest for Redhaven and Elberta peach trees. HortScience 9(4), 331-332.

The computation is quite straightforward but requires hourly temperature data. Since we don't have this I am implementing some utility function to compute hourly temperature from daily maximum and minimum temperature.

I was wondering whether you would be interested in me integrating this into xclim. My code would need some more adaptations to work with different calendars etc. but I think it could be done.

Additional context

Computing the hourly temperature requires computing time of sunrise and sunset, which is probably the biggest headache. My computation is based on this https://en.wikipedia.org/wiki/Sunrise_equation. This paper is suggesting how to compute the hourly temperature, it assumes that minimum temperature is reached just before sunrise, that maximum temperature is reached two hours before sunset and that during daylight the relation ship is sinoid and during nighttime logarithmic.

Most of the logic is needed because in UTC today's sunrise could be yesterday or tomorrow for locations close to -180 or 180 respectively. Additionally, the sun equation requires julian dates that I so far have found difficult to compute without several calendar conversions (we are using a 360day calendar from cftime and I couldn't figure out how to use cftime.date2num to compute the julian day number).

Contribution

• [X] I would be willing/able to open a Pull Request to contribute this feature.

Code of Conduct

• [X] I agree to follow this project's Code of Conduct

[coxipi]

This might help you: xclim.indices.helpers.cosine_of_solar_zenith_angle ? The case with sunlit=True?

[saschahofmann]

Do I understand this right that this would be equivalent to this ?

[coxipi]

# cosine_of_solar_zenith_angle
    if sunlit:
        # hour angle of sunset (eq. 2.15), with NaNs inside the polar day/night
        tantan = -np.tan(lat) * np.tan(declination)
        h_ss = np.arccos(tantan.where(abs(tantan) <= 1))

I would say no, it looks like it implements the first formula on the wiki page:

The code in stat == "instant" uses the formula: see https://en.wikipedia.org/wiki/Solar_zenith_angle which is equivalent to your formula, with the altitude not fixed, but to be computed. The output cosine_of_solar_zenith_angle is the cosine of solar zenith angle = sin of solar altitude.

I'm a bit confused, is this a fixed or a variable quantity? @aulemahal , any thoughts?

[aulemahal]

@coxipi I'm not sure I understand your question. Which quantity are you referring to ?

@saschahofmann Indeed cosine_of_solar_zenith_angle might help you, but maybe not as is, more like an inspiration.

(All functions referred to in this comment are in the xclim.indices.helpers module)

In the function, for the finer-than-daily case, we compute the local hour angle of the sun, of the time coordinate as:

S_IN_D = 24 * 60 * 60   # Number of seconds in a day
if time.dtype == "O":  # cftime
    # time_as_s : seconds elapsed since Jan 1st 1970 at 00:00.
    time_as_s = time.copy(data=xr.CFTimeIndex(time.values).asi8 / 1e6)
else:  # numpy
    time_as_s = time.copy(data=time.astype(float) / 1e9)
h_s_utc = (
    (
        (time_as_s % S_IN_D)   # Number of seconds since last midnight
        / S_IN_D               # Fraction of a day since last midnight
    ) * 2 * np.pi              # as radians (1 day = 2 pi)
   - np.pi                     # Midnight is -pi, noon is 0, next midnight is pi (see note)
).assign_attrs( units="rad")
h_s = h_s_utc + lon            # Assuming time as UTC we get local hour angle by simply adding longitude in radians

note : h_s is in radians, so it is defined modulo 2 pi. At the last step of h_s_utc, the current function does + np.pi, it's the same once we ensure radians are wrapped. See _wrap_radians.

Thus we have the local hour angle.

The (local) sunrise and sunset hour angles are further below, as shown by Eric:

if sunlit:
    # hour angle of sunset (eq. 2.15), with NaNs inside the polar day/night
    tantan = -np.tan(lat) * np.tan(declination)
    h_ss = np.arccos(tantan.where(abs(tantan) <= 1))

Where, declination was computed with solar_declination, two methods are implemented.

There is also a time_correction_for_solar_angle function that was taken from this article : https://link.springer.com/article/10.1007/s00484-020-01900-5, which cites this dead website : https://web.archive.org/web/20130217133755/http://www2.ncdc.noaa.gov/docs/gviug/html/l/app-l.htm. I have not dug deep enough to fully understand what that does, but it was needed to get similar results the first article, when implementing the cosine_of_solar_zenith_angle function for MRT and UTCI. The correction is used when we want the instant cosine, which means you might want it if you are comparing "instant" hour angles.

The cosine_of_solar_zenith_angle actually does an integral when stat is not "instant", which is why (?) the correction is not always required and why some of the trigonometric formulations are different : they are already integrated.

[saschahofmann]

I actually figured that for my needs its enough to have the daily temperature distribution. The true times don't really matter that much. So I used xclim's day_lengths function and computed the hourly temperature but put the sunrise globally at midnight and sunrise would then be after daylength hours. Since I am using this to compute daily aggregated values like chilling units the true timestamps don't really matter. Would you be happy with an implementation like that ?

[saschahofmann]

I now also have a chill portions/dynamic model implementation. If you want I could port them to xclim. If so, where should I put them (I can see that there is a module called _agro.py).

And about the interface: Both functions require hourly data. Should the functions perform a check on the resolution of the data? Usually they will be computed for the dormancy period. The functions could either simply compute each indicator for the whole period provided or accept a range of days and compute it for only these ranged within the whole dataset.

Both functions could also be understood as accumulation over these periods or should they return the indicator value per hour?

[Zeitsperre]

@saschahofmann I think _agro.py is likely the best place to put this.

For the hourly time step, we don't really support that level of granularity in our Indicators, though I imagine that most indices would technically work for hourly data. @aulemahal correct me if I'm wrong, but supporting some proper hourly Indicators would likely be as simple as inheriting an Indicator class and modifying the methods responsible for checking the time, right?

[aulemahal]

Exact. The src_freq attribute of the Indicator class can simply be modified to ["H"] to (only) accept hourly data.

[saschahofmann]

Alright great. I would also add a functionality that estimates hourly temperatures from daily min max. But my suggestion is that this should be called outside the chill portion function