ENH: Lagrangian interpolation for UT1-UTC in IERS data #18036
Conversation
|
Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.
|
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
There was a problem hiding this comment.
Welcome to Astropy 👋 and congratulations on your first pull request! 🎉
A project member will respond to you as soon as possible; in the meantime, please have a look over the Checklist for Contributed Code and make sure you've addressed as many of the questions there as possible.
If you feel that this pull request has not been responded to in a timely manner, please send a message directly to the development mailing list. If the issue is urgent or sensitive in nature (e.g., a security vulnerability) please send an e-mail directly to the private e-mail [email protected].
|
CI is failing so I think this is too risky to go into v7.1 given how close we are to feature freeze. I have turned this into draft for and milestoned it for v7.2. Thank you for your interest and understanding. |
|
No worries, thanks for looking at this. |
mhvk
left a comment
mhvk
left a comment
There was a problem hiding this comment.
Thanks for the PR! Some smaller comments are in-line.
Bigger ones: As @pllim notes, first thing would be to make sure that the code passes tests. Note that the failures may well be because the answer actually changes a little. But we should not simply adjust test, but rather expand them for two choices, of either linear interpolation ignoring tides, or quadratic including them (see also below).
More generally: if we go this route, we will also want to interpolate other columns with higher order, so it may make more sense to try to have one interpolation routine rather than special-case UT1-UTC. The only reason UT1_UTC is different from X and Y is because of the leap seconds, so perhaps one can find some generic way to deal with that, then interpolate either linearly or quadratically, and then take into account the leap second again.
Another general note: I think that without the tides there is not so much point in doing the quadratic interpolation. Also, eventually, since the tides is a relatively heavy computation, we probably want a ScienceState configuration to decide how to do it.
Finally, just a note: we are very close to feature freeze, so I won't have time to help much the coming weeks. I do really like the idea of closing that very old issue, though!
astropy/utils/iers/iers.py
Outdated
| for j in range(len(x_vals)): | ||
| if i != j: | ||
| diff = x_vals[i] - x_vals[j] | ||
| diff = np.where(diff == 0, 1e-10, diff) |
There was a problem hiding this comment.
Within the IERS table, it is guaranteed that the x values are never the same, so this shouldn't be necessary.
astropy/utils/iers/iers.py
Outdated
|
|
||
| def _ut1_utc_interp(self, i0, i1, mjd): | ||
| indices = np.array([i0 - 1, i0, i1, i1 + 1]) | ||
| indices = np.clip(indices, 0, len(self) - 1) |
There was a problem hiding this comment.
I think you'll need to set up the clipping of i1 on l.462 (pre-change) in interpolate differently instead of redoing the clipping here. Note that if you look at https://hpiers.obspm.fr/iers/models/interp.f, it is not the indices that are clipped - rather it is made sure that there are always three values.
astropy/utils/iers/iers.py
Outdated
| vals[1:] += leap_diff | ||
|
|
||
| val = self._lagrange_interp(mjds, vals, mjd) | ||
| val -= leap_diff[0] |
There was a problem hiding this comment.
I think in the general case you cannot be sure index 0 is the right one to subtract. Which one it is depends on how the 3 surrounding points are selected, so this likely is an array.
Probably also need benchmarking then? |
Yes, but only if the tides are actually included (which is not yet the case). |
|
Thank you @mhvk ; these all make sense! I'm starting by generalizing the routine to allow (optionally) interpolation for any column / order, which should clear a lot of the issues, next I'll add the tides logic. |
…e leap second handling for UT1-UTC
|
Sorry for the slight delay and thanks again @mhvk, the pointers were really helpful! I have provisioned for allowing Lagrange interpolation over all columns, and the order can be changed from the default (cubic), but it's the same for all columns; I suppose we would want different orders for different columns (higher for X and Y perhaps), but what would they be (at least the defaults)? I think I've addressed the minor points--the leap second handling and index clipping is now more sophisticated, and works with the order. One can select to use linear interpolation and it is enabled by default, so the tests pass as usual. Thinking of adding new tests for the other interpolation methods but perhaps we can clarify the details and add the tides first. If this looks like the right direction I can go ahead and add the tides! |
|
The failing tests are for coverage (which will be added of course), and the changelog -- not sure if anything needs to be done there! |
|
I've made some improvements, including a fix to the leap second offset index, and added some parametrized coverage for interpolation to |
|
Hi @mhvk, I've added adjustment for the luni-solar tides! I'll double check again the constants as there were so many, and then maybe should add a boolean parameter for this? Would appreciate a look at this. Thank you! |
4 participants
Description
This PR implements cubic Lagrangian interpolation (instead of linear interpolation) for only UT1-UTC in IERS data as suggested in #1803. Currently does not correct for luni-solar tides (step 2 and 4).
I'm new to this project, so hope this aligns with the feature request! There was some confusion around verifying the procedure in the linked issue it seems, so open to discussing further--thank you for reviewing!
Fixes #1803.