Improve ERFA warnings about dubious years to give more user-friendly information

[astrofrog]

As discussed in some other issues (#9587 (comment); #9593 (comment)) it would be good if ERFA warnings could be more user-friendly. For example, this tells me nothing much as a user:

WARNING: ErfaWarning: ERFA function "utctai" yielded 3 of "dubious year (Note 3)" [astropy._erfa.core]

A better warning would be:

WARNING: conversion between the utc and tai scales is may be inaccurate as dates in the utc scale outside of the range ... to ... are not well defined

In general whenever I get ERFA warnings I basically never really understand them and also don't know how to fix them. Clearer messages such as the above would really help in that regard.

[mhvk]

Currently, the warning messages are auto-generated from the *.c files; clarifying some would not be too difficult (but obviously a bit of a hassle) by special-casing some functions.

As an example, the warning discussed most so far comes from dat.c; it is extracted from:

**  Returned (function value):
**            int      status (Note 5):
**                       1 = dubious year (Note 1)
**                       0 = OK
**                      -1 = bad year
**                      -2 = bad month
**                      -3 = bad day (Note 3)
**                      -4 = bad fraction (Note 4)
**                      -5 = internal error (Note 5)

where note 1 in dat.c is

**  1) UTC began at 1960 January 1.0 (JD 2436934.5) and it is improper
**     to call the function with an earlier date.  If this is attempted,
**     zero is returned together with a warning status.
**
**     Because leap seconds cannot, in principle, be predicted in
**     advance, a reliable check for dates beyond the valid range is
**     impossible.  To guard against gross errors, a year five or more
**     after the release year of the present function (see the constant
**     IYV) is considered dubious.  In this case a warning status is
**     returned but the result is computed in the normal way.
**
**     For both too-early and too-late years, the warning status is +1.
**     This is distinct from the error status -1, which signifies a year
**     so early that JD could not be computed.

[astrofrog]

One of the issues is that the notes don't mean anything to a user - i.e. we cannot realistically expect them to go and check the source code. I wonder whether we can still somehow automate the warnings but include the actual text of the note too, or whether that would be confusing.

[mhvk]

Agreed, My sense is that the easiest is simply to replace some of the more common warnings with custom ones. For instance, one could override specific entries in STATUS_CODES dict in _erfa/core.py.templ at the very end of the file.

[timj]

My biggest problem with the current warning situation is that it's impossible for me to tell which conversion is actually causing the problem. It would help a lot if the warning could come from user code instead of deep inside erfa core. In general I think that warnings should probably always set the stacklevel argument to warnings.warn if you want the warning to be actionable. Using stacklevel=3 pulls it out of ERFA but you still get something confusing to the user:

astropy/time/core.py:800: ErfaWarning: ERFA function "taiutc" yielded 1 of "dubious year (Note 4)"
    jd1, jd2 = conv_func(*args)

Making this change to core.py in ERFA makes the warning useful:

        stacklevel = 1
        import traceback
        stack = traceback.extract_stack()
        for i, s in enumerate(reversed(stack)):
            if "astropy" not in s.filename:
                stacklevel = i + 1
                break

        warnings.warn('ERFA function "' + func_name + '" yielded ' + wmsg, ErfaWarning, stacklevel=stacklevel)

since it now tells me which line of code outside of astropy was causing the warning. If all astropy warnings tried to report user code they become much more useful.

I'm not sure what people think about this sort of change though.

[mhvk]

That's such a tricky one! Numpy has been adjusting stack levels based on similar rationale, but for those the net result seems to be that often it comes out of another library that is called numpy (for me, usually astropy), which is not all that much more useful.

Also, the stacklevel is actually tricky: the erfa routines are now in pyerfa so from their perspective I don't see how one could do anything but stacklevel=1 - they do not know whether they are called from astropy or anything else...

Overall, I think most important would be to improve the actual warnings. Right now, these are simply automatically extracted from the erfa C code.

[timj]

With the code above I was finally able to track down where the elusive warning was coming from. It had be driving me bonkers for a while.

You are right that with pyerfa as a first class entity this is seemingly impossible to fix in the erfa code without astropy intercepting the warning and then reissuing it with a clever stacklevel itself. At Rubin we get this erfa warning a lot (sometimes simulation data is used for tests) and it invariably takes a lot of time to work out which particular time conversion was triggering it in our code. It would help enormously if it could at least report the first stack entry outside of astropy.time specifically. Yes, that may well not be useful if it's astropy.coordinates that is triggering it but it will help the people who are using astropy.time directly and seems like it's hard to trigger a false positive for that specific case.

[mhvk]

I tend to just turn warnings into errors, so that I get a proper trace and can enter the debugger. Would that work here?

I ask because the cleanest solution would seem to catch warnings in astropy and then adjust stacklevel or the warning message, but I've been getting a bit worried generally about the balance between trying to help users and code clarity - the latter obviously important if we want to get new developers.

[timj]

Turning all warnings into errors can't quite work because at the moment we have some warnings from third parties that happen first. I wish that python had an option to show stack trace for warnings (maybe only in developer mode -X dev).

[mhvk]

But one can just turn ErfaWarning into an error, no?

[timj]

I admit I have not looked deeply into all the options that I can use with -W.

[astrojuanlu]

I don't think this can be done with custom warning classes but I'd love to be proven wrong here...

[mhvk]

@timj, @astrojuanlu - Maybe I'm simply missing what is being tried/needed, but the following certainly error with a full traceback:

python3 -W "error:Erfa function" -c 'from astropy.time import Time; Time("J1950").ut1'
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/usr/lib/python3/dist-packages/astropy/time/core.py", line 1708, in __getattr__
    tm._set_scale(attr)
  File "/usr/lib/python3/dist-packages/astropy/time/core.py", line 800, in _set_scale
    jd1, jd2 = conv_func(*args)
  File "/usr/lib/python3/dist-packages/astropy/_erfa/core.py", line 15911, in taiutc
    check_errwarn(c_retval, 'taiutc')
  File "/usr/lib/python3/dist-packages/astropy/_erfa/core.py", line 104, in check_errwarn
    warnings.warn('ERFA function "' + func_name + '" yielded ' + wmsg, ErfaWarning)
astropy.utils.exceptions.ErfaWarning: ERFA function "taiutc" yielded 1 of "dubious year (Note 4)"

Or simply working interactively (on astropy 4.0, so astropy._erfa is relevant module):

import warnings
from astropy.time import Time


warnings.filterwarnings('error', module='astropy._erfa')

Time('J1950').ut1

[taldcroft]

So maybe this comes down the usual (and relatively easy!) documentation fix of including an example of tracking down these warnings.

@astrojuanlu - I have the same recollection about not being able to filter warnings for custom user warning classes from the command line. In some cases it will be possible to import the ErfaWarning class and use that in a call to filterwarnings. But for this particular case the solution from @mhvk looks good.

[mzechmeister]

Is there a possibility to update a leap second table? I couldn't figure out from where erfa retrieves the leap seconds. I got stucked at https://github.com/liberfa/erfa/blob/5233e79cd97c2d3fdae4b7a8f802d46be46867a5/src/dtf2d.c.
That would probably be a better approach, instead of suppressing the warning.

[mhvk]

@mzechmeister - the leapsecond table should be updated automatically in astropy - is it not for you?

[mzechmeister]

@mhvk How can I check, if the leap second table it up-to-date? I'm still running on 2.0.3, sorry. astropy.utils.iers.conf.auto_download returns True.
Here is what I tried. I'm not sure that table contains leap second information.

>>> from astropy.time import Time
>>> Time.now()
WARNING: ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "d2dtf" yielded 1 of "dubious year (Note 5)" [astropy._erfa.core]
<Time object: scale='utc' format='datetime' value=2022-01-17 08:05:23.061601>
>>>
>>> from astropy.utils import iers
>>> iers_table = iers.IERS_Auto.open()
WARNING: failed to download http://maia.usno.navy.mil/ser7/finals2000A.all, using local IERS-B: <urlopen error timed out> [astropy.utils.iers.iers]
>>> iers_table
<IERS_B length=19861>
 year month  day    MJD      PM_x     PM_y    UT1_UTC      LOD    dX_2000A dY_2000A  e_PM_x  e_PM_y e_UT1_UTC  e_LOD   e_dX_2000A e_dY_2000A
                     d      arcsec   arcsec      s          s      arcsec   arcsec   arcsec  arcsec     s        s       arcsec     arcsec  
int64 int64 int64 float64  float64  float64   float64    float64  float64  float64  float64 float64  float64  float64   float64    float64  
----- ----- ----- ------- --------- -------- ---------- --------- -------- -------- ------- ------- --------- -------- ---------- ----------
 1962     1     1 37665.0   -0.0127    0.213  0.0326338  0.001723      0.0      0.0    0.03    0.03     0.002   0.0014   0.004774      0.002
 1962     1     2 37666.0   -0.0159   0.2141  0.0320547  0.001669      0.0      0.0    0.03    0.03     0.002   0.0014   0.004774      0.002
 1962     1     3 37667.0    -0.019   0.2152  0.0315526  0.001582      0.0      0.0    0.03    0.03     0.002   0.0014   0.004774      0.002
 1962     1     4 37668.0 -0.021999 0.216301  0.0311435  0.001496      0.0      0.0    0.03    0.03     0.002   0.0014   0.004774      0.002
  ...   ...   ...     ...       ...      ...        ...       ...      ...      ...     ...     ...       ...      ...        ...        ...
 2016     5    13 57521.0  0.056405 0.487315 -0.1562674 0.0014815 -0.00012 0.000139   4e-05 3.7e-05   8.3e-06  1.9e-05    3.7e-05      5e-05
 2016     5    14 57522.0  0.058448  0.48823 -0.1578078 0.0015995 -9.3e-05 0.000112   4e-05 3.7e-05   8.3e-06 1.88e-05    3.7e-05    4.9e-05
 2016     5    15 57523.0   0.06049 0.488841 -0.1594424 0.0016863 -7.3e-05  9.2e-05   4e-05 3.7e-05   8.2e-06 1.86e-05    3.6e-05    4.7e-05
 2016     5    16 57524.0  0.062539 0.489126 -0.1611736 0.0017572 -5.4e-05  7.2e-05   4e-05 3.7e-05   8.1e-06 1.85e-05    3.5e-05    4.5e-05
 2016     5    17 57525.0  0.064598 0.489628 -0.1629401 0.0017687 -3.4e-05  5.2e-05 4.1e-05 3.7e-05   7.9e-06 1.83e-05    3.5e-05    4.4e-05

Using a more recent iers_auto_url does not help.

>>> from astropy.utils import iers
>>> iers.conf.iers_auto_url = 'https://datacenter.iers.org/data/9/finals2000A.all'
>>> iers_table = iers.IERS_Auto.open()
>>> iers_table
<IERS_Auto length=18254>
 year month  day    MJD   PolPMFlag_A  PM_x_A  e_PM_x_A  PM_y_A  ... dX_2000A_B dY_2000A_B  UT1_UTC   UT1Flag    PM_x     PM_y   PolPMFlag
                     d                 arcsec   arcsec   arcsec  ...  marcsec    marcsec       s                arcsec   arcsec           
int64 int64 int64 float64     str1    float64  float64  float64  ...  float64    float64    float64   unicode1 float64  float64   unicode1
----- ----- ----- ------- ----------- -------- -------- -------- ... ---------- ---------- ---------- -------- -------- -------- ---------
   73     1     2 41684.0           I 0.120733 0.009786 0.136966 ...    -18.637     -3.667  0.8078584        B   0.1235    0.123         B
   73     1     3 41685.0           I  0.11898 0.011039 0.135656 ...    -18.636     -3.571  0.8051525        B   0.1223    0.121         B
   73     1     4 41686.0           I 0.117227 0.011039 0.134348 ...    -18.669     -3.621  0.8024015        B 0.120901 0.119001         B
   73     1     5 41687.0           I 0.115473 0.009743 0.133044 ...    -18.751     -3.769  0.7995416        B 0.119601 0.117101         B
  ...   ...   ...     ...         ...      ...      ...      ... ...        ...        ...        ...      ...      ...      ...       ...
   22    12    20 59933.0           P 0.117546 0.017537 0.246806 ...      1e+20      1e+20 -0.0262065        P 0.117546 0.246806         P
   22    12    21 59934.0           P 0.115571 0.017564 0.247231 ...      1e+20      1e+20 -0.0261637        P 0.115571 0.247231         P
   22    12    22 59935.0           P 0.113605 0.017591 0.247687 ...      1e+20      1e+20 -0.0258937        P 0.113605 0.247687         P
   22    12    23 59936.0           P 0.111649 0.017618 0.248174 ...      1e+20      1e+20 -0.0254386        P 0.111649 0.248174         P
   22    12    24 59937.0           P 0.109703 0.017645 0.248692 ...      1e+20      1e+20   -0.02489        P 0.109703 0.248692         P
>>>
>>> from astropy.time import Time
>>> Time.now()
WARNING: ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "d2dtf" yielded 1 of "dubious year (Note 5)" [astropy._erfa.core]
<Time object: scale='utc' format='datetime' value=2022-01-17 08:08:37.522942>

[mhvk]

@mzechmeister - the error you get is because the leap-second table is out of date, not the IERS tables. Unfortunately, in astropy 2.0.3, there was not yet a way to update the leap seconds. I fear you'll just have to update astropy (or at least the erfa library, but if that is not a system library, then for 2.0.3 the erfa library is built in to astropy, so likely is more work).

[stroncod]

@mhvk solution to hide the error does not work on astropy '5.1.1'.

>>> import warnings
>>> from astropy.time import Time
>>> warnings.filterwarnings('error', module='astropy._erfa')
>>> Time('J1950').ut1
/Users/stroncoso/opt/miniconda3/envs/scheduler-env/lib/python3.10/site-packages/erfa/core.py:154: ErfaWarning: ERFA function "taiutc" yielded 1 of "dubious year (Note 4)"
  warnings.warn('ERFA function "{}" yielded {}'.format(func_name, wmsg)|

[mhvk]

In astropy 5.1.1, erfa is no longer bundled in, so I think one should filter on module erfa not astropy._erfa.