ENH: Lagrangian interpolation for UT1-UTC in IERS data - #18036

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.

• Do the proposed changes actually accomplish desired goals?
• Do the proposed changes follow the Astropy coding guidelines?
• Are tests added/updated as required? If so, do they follow the Astropy testing guidelines?
• Are docs added/updated as required? If so, do they follow the Astropy documentation guidelines?
• Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for rebase and squash.
• Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the bot.
• Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
• Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
• At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) before merge.

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.

tides is a relatively heavy computation

Probably also need benchmarking then?

tides is a relatively heavy computation

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.

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 test_iers.py. Seems the tests pass now!

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!

Hi humans 👋 - this pull request hasn't had any new commits for approximately 5 months. I plan to close this in a month if the pull request doesn't have any new commits by then.

In lieu of a stalled pull request, please consider closing this and open an issue instead if a reminder is needed to revisit in the future. Maintainers may also choose to add keep-open label to keep this PR open but it is discouraged unless absolutely necessary.

If this PR still needs to be reviewed, as an author, you can rebase it to reset the clock.

If you believe I commented on this pull request incorrectly, please report this here.

@MaanasArora - I'm very sorry that this completely slipped through the cracks... I only got reminded because the bot is trying to close this PR (I've removed the close? label for now).

I think the code itself is good, but I wonder about how best to expose it to the user. Generally, someone does not use this code directly, but it gets called when transforming a time or coordinate. Passing up arguments seems a non-starter; instead, I think these belong as configuration items in iers.conf -- see the Conf class. I think it is not super helpful to make this incredibly flexible. E.g., I could imagine having an iers_interpolation_scheme configuration option that can take values of

• "linear" or "astropy03": to get the old scheme
• "tides" or "astropy81" ("astropy90"?): new lagrange + tides

I guess in principle one could also have "lagrange" which omits the tides, but it is not obvious to me that is useful: is it actually closer to the truth than linear? It would need testing. I do not think we should expose things like order

On "the truth", we do need something to compare with! I guess easiest would be to run interp.f. (Indeed, I now wonder if we cannot just use numpy's f2py to translate that to a python module, even if just for testing.)

@MaanasArora - turned out interp.f was actually fairly easily compilable with f2py (see below). We can use this to verify results (I even wondered about just using it, but do not know we can start relying on a fortran compiler being present...)

Interpolation of IERS tables

The file interp.f was downloaded from
https://hpiers.obspm.fr/iers/models/interp.f, untabified in emacs, and
the following patches applied:

29a30
>       use, intrinsic :: ieee_arithmetic, only: IEEE_Value, IEEE_QUIET_NAN
31c32
<       integer i,n
---
>       integer n
35a37,42
>       IF (rjd_int .lt. RJD(1) .or. rjd_int .gt. RJD(n)) THEN
>          x_int = IEEE_VALUE(x_int, IEEE_QUIET_NAN)
>          y_int = x_int
>          ut1_int = x_int
>          RETURN
>       END IF
86a94
>       k = 0
90a99
>       if ( k .eq. 0 ) return

A draft signature file was made by

f2py interp.f -m interp_mod -h interp.pyf

and then edited to make clear what variables are input and output
ones, leading to:

!    -*- f90 -*-
! Note: the context of this file is case sensitive.

python module interp_mod ! in
    interface  ! in :interp
        subroutine interp(rjd,x,y,ut1,n,rjd_int,x_int,y_int,ut1_int) ! in :interp:interp.f
            double precision dimension(n), intent(in) :: rjd
            double precision dimension(n), intent(in) :: x
            double precision dimension(n), intent(in) :: y
            double precision dimension(n), intent(in) :: ut1
            integer, intent(hide), depend(rjd) :: n=shape(rjd, 0)
            double precision :: rjd_int
            double precision, intent(out) :: x_int
            double precision, intent(out) :: y_int
            double precision, intent(out) :: ut1_int
        end subroutine interp
        subroutine lagint(x,y,n,xint,yout) ! in :interp:interp.f
            real*8 dimension(n), intent(in) :: x
            real*8 dimension(n), intent(in) :: y
            integer, intent(hide), depend(x) :: n=shape(x, 0)
            real*8, intent(in) :: xint
            real*8, intent(out) :: yout
        end subroutine lagint
        subroutine pmut1_oceans(rjd,cor_x,cor_y,cor_ut1,cor_lod) ! in :interp:interp.f
            real*8, intent(in) :: rjd
            real*8, intent(out) :: cor_x
            real*8, intent(out) :: cor_y
            real*8, intent(out) :: cor_ut1
            real*8, intent(out) :: cor_lod
        end subroutine pmut1_oceans
        subroutine pm_gravi(rjd,cor_x,cor_y) ! in :interp:interp.f
            double precision, intent(in) :: rjd
            double precision, intent(out) :: cor_x
            double precision, intent(out) :: cor_y
        end subroutine pm_gravi
    end interface
end python module interp

! This file was auto-generated with f2py (version:2.5.1).
! and then edited by MHvK
! See:
! https://web.archive.org/web/20140822061353/http://cens.ioc.ee/projects/f2py2e

Compilation to a module is then with,

f2py -c --f77flags="-ffixed-line-length-0" interp.f interp.pyf

Then, as a trial,

from interp_mod import interp
from astropy.time import Time
from astropy.utils.iers import IERS_Auto
iers = IERS_Auto.open()

t = Time("2025-04-08T09:10:11")
iers.pm_xy(t)
# [<Quantity 0.06410133 arcsec>, <Quantity 0.38834011 arcsec>]
iers.ut1_utc(t)
# <Quantity 0.03676799 s>

interp(iers["MJD"], iers["PM_x"], iers["PM_y"], iers["UT1_UTC"], t.mjd)
(0.06414316493735002, 0.38862301914940645, 0.03677844300346548)

@mhvk - no worries, thank you for looking at this, and the helpful pointers! (Sorry, this comment ended up long.)

About exposing the feature: yes, agreed, providing the interpolation method in configuration seems much better than passing arguments. I added a configuration parameter to the IERS configuration class. It accepts linear (old scheme) and tides (with lagrange interpolation and tidal corrections).

I was able to compile interp.f with f2py as well using your steps and get results (more below). But while debugging an issue during compilation, not sure how, I noticed my implementation was actually very garbled! It seemed to mainly just have the arg such as chi, omega right, but the algorithm following was incorrectly mixed between routines and had coefficient mistakes. So, I had to go back and rewrite the whole thing before testing:

1. As part of this I converted the correct coefficient matrix (the relevant one was very large) into a CSV file, mainly by replacing the special characters (&, extra whitespace, and slashes), and adding a header. Then I had a small Python script write it as NumPy code, below for future reference. I'm wondering if having such a large array in the code is acceptable, though? I'm happy to use a CSV, npy, or another data file directly, I'm not sure what the preference is.

import pandas as pd
import numpy as np

csv = pd.read_csv("interp_pmut1_oceans.csv", index_col=False)
arg = csv.iloc[:, :6].to_numpy()
val = csv.iloc[:, 6:].to_numpy()

# format as Python code to create a numpy array
print("OCEANS_ARG = np.array([")
for row in arg:
    print("    [", end="")
    print(", ".join(f"{x}" for x in row), end="")
    print("],")
print("])")

print("OCEANS_VAL = np.array([")
for row in val:
    print("    [", end="")
    print(", ".join(f"{x}" for x in row), end="")
    print("],")
print("])")

2. Since it seems that the PM_GRAVI routine only impacts the x/y columns and not UT1-UTC, I have excluded that part of interp.f from this PR. This does mean that I only use the tides method for UT1-UTC currently, even if it is selected, but I think the code can be expanded quite easily as the foundation is there. We would only need to convert that routine and add it to the pipeline if the other columns are involved.

On testing, and whether just a lagrange option is worth it: I used this script to compare each algorithm with the f2py reference:

from interp_mod import interp
from astropy.time import Time
from astropy.utils import iers as iers_module
from astropy.utils.iers import IERS_Auto
import numpy as np


def get_expected_values(iers, ts):
    return np.array(
        [
            interp(iers["MJD"], iers["PM_x"], iers["PM_y"], iers["UT1_UTC"], t.mjd)
            for t in ts
        ]
    ).transpose()


def get_values_from_iers(iers, ts, method):
    iers_module.conf.ut1_utc_interpolation = method
    x, y = iers.pm_xy(ts)
    ut1_utc = iers.ut1_utc(ts)
    return np.array([x, y, ut1_utc])


def compare_interp_methods(iers, ts):
    ts_expected = get_expected_values(iers, ts)[2]
    ts_linear = get_values_from_iers(iers, ts, "linear")[2]
    ts_tides = get_values_from_iers(iers, ts, "tides")[2]

    error_linear = np.abs(ts_linear - ts_expected)
    error_tides = np.abs(ts_tides - ts_expected)

    return error_linear, error_tides


if __name__ == "__main__":
    iers = IERS_Auto.open()

    # sample 1000 points between 2024 and 2026
    start_mjd = Time("2024-01-01").mjd
    end_mjd = Time("2026-01-01").mjd
    ts = Time(np.linspace(start_mjd, end_mjd, 1000), format="mjd")

    error_linear, error_tides = compare_interp_methods(iers, ts)

    print(f"Mean absolute error (linear): {np.mean(error_linear):.3e} seconds")
    print(f"Mean absolute error (tides): {np.mean(error_tides):.3e} seconds")

    print(f"Max absolute error (linear): {np.max(error_linear):.3e} seconds")
    print(f"Max absolute error (tides): {np.max(error_tides):.3e} seconds")

With the current state of the branch, the results are very close to the reference:

Mean absolute error (linear): 2.510e-05 seconds
Mean absolute error (tides): 6.638e-13 seconds

Max absolute error (linear): 1.093e-04 seconds
Max absolute error (tides): 2.094e-12 seconds

The gains without tidal corrections seem much more marginal. So, I suppose it is not really worth it, so I'm not exposing a lagrange option, nor an order parameter:

Mean absolute error (linear): 2.510e-05 seconds
Mean absolute error (tides): 2.285e-05 seconds

Max absolute error (linear): 1.093e-04 seconds
Max absolute error (tides): 8.199e-05 seconds

Thanks, that is great! So, it looks like we just have two config options.

I agree with your hesitance about large parameter arrays inside the code. Since you already made the .csv, I think we should just include that and load it on demand. After all, we are already loading the (far larger) IERS tables on demand as well, so it is not like this will meaningfully affect performance.

However, I think we should use the Table machinery (or np.loadtxt) to load the .csv; we do not want to add a dependence on pandas. Alternatively, we can include the interp.f file and extract the arguments programmatically from that. A bit of playing around shows this works:

from astropy.table import Table
from astropy.utils.data import get_pkg_data_contents

interp = get_pkg_data_contents("src/interp.f")
pre, _, gravi = interp.rpartition("data(")
oceans = pre.partition("data(")[2]

def parse_data_block(data):
    # Get header into a single line
    hdr, _, post  = data.partition(")/\n")
    # remove bits that are irrelevant
    hdr = hdr.replace("     &", "").replace("),j=1,nlines", "").strip()
    hdr = hdr.replace("(j,", "").replace("(j", "").replace("),", ",").replace("\n", "")
    # Get value lines
    values = post.partition("/\n")[0].replace("     &", "").replace(",\n", "\n")
    return Table.read("\n".join([hdr, values]), format="ascii.csv")

ocean = parse_data_block(oceans)
gravi = parse_data_block(gravi)

p.s. Should have added: since getting it from interp.f works, to me it makes most sense to just include that file. That also means that if a fortran compiler is available, we can easily re-check numbers.

Thanks, yes I just included the interp.f file in a src/ directory in the module (I thought of data/, but this is moreso code). The snippet you shared worked nearly as-is. Yes, it is good to have the source for verification in the future, and this also means we can easily add gravi when need be! I'd mark this as ready for review now, except codecov doesn't seem to recognize testing coverage for the new config option.

Thanks @mhvk! np.newaxis used, makes sense. I also added the location procedure. The code was shared mostly so it didn't take up too many more lines. That said, I think I might have found a mistake in the ARG preparation, not entirely sure why ARG(3) (lp) is different in the oceans versus gravi procedure:

PMUT1_OCEANS:

      ARG(3) = -0.00001149d0*T**4 + 0.000136d0*T**3
     .  -  0.5532d0*T**2
     .  + 129596581.0481d0*T + 1287104.79305d0

PM_GRAVI:

      ARG(3) = -0.00001149d0*T**4 - 0.000136d0*T**3
     .  -  0.5532d0*T**2
     .  + 129596581.0481d0*T + 1287104.79305d0

The sign of the second term (0.000136d0*T**3) seems flipped. I've not inherited that, but happy to go either direction if you know what should be done here. To be fair even it's an error, might not be very harmful to be consistent if we're unsure. And yes, we can move to a better configuration in ScienceState whichever way you like!

That's a nice catch! Looking at https://hpiers.obspm.fr/eop-pc/models/models.html_nov08, I'm led to https://hpiers.obspm.fr/eop-pc/models/PM/PM_oceanic_hf_tab.html and via interp.f I infer that ARG(3) is lp, which presumably means l', for which the bottom of the web page lists,

Mean anomaly of the Sun  : 
   l'= 357°.529 109 18 +   129 596 581.0481" t - 0".553 2 t2 -  0".000 136 t3 - 0".000 011 49 t4

So, the term should be negative in OCEANS, while GRAVI is correct. Indeed, that can also be inferred from the derivative right below in OCEANS.

I think that more or less settles it, but let me send an e-mail to point out the inconsistency. That also gives me a nice opportunity to verify the procedure, and ask whether there are some reference numbers we can compare with.

Thanks @mhvk for the email, that is likely correct indeed, so I've updated the sign to be negative!