[AstroPy] How to compute TAI-UTC using Astropy?

[Russell Owen]

> On Aug 25, 2019, at 3:42 AM, Aldcroft, Thomas <aldcroft at head.cfa.harvard.edu> wrote:
> 
> 
> 
> On Sat, Aug 24, 2019 at 7:24 PM Russell Owen <rowen at uw.edu <mailto:rowen at uw.edu>> wrote:
> Thank you for the suggestion, but I’m afraid I don’t see how it works.. The given date I have to work from is a date in UTC. By specifying “Time(date, scale=“tai”)” I would be lying to the system: saying the date was TAI when it is UTC. I would expect that to give the wrong answer very near a UTC leap second.
> 
> I think the code works as expected within a leap second, at the expense of an ERFA warning that you would need to catch:
> 
> In [2]: dts = [0, 0.5, 1.0, 1.5, 1.99, 2.0, 2.5, 3.0] * u.s
> In [3]: dates_utc = Time(Time('2015-06-30 23:59:59') + dts, scale='utc')
> In [4]: dates_tai = Time(dates_utc.iso, scale='tai')
> WARNING: ErfaWarning: ERFA function "dtf2d" yielded 3 of "time is after end of day (Note 5)" [astropy._erfa.core]
> In [5]: for date_utc, date_tai in zip(dates_utc, dates_tai):
>    ...:     print(f'{date_utc} {(date_utc - date_tai).sec:.2f}')
>    ...:     
> 2015-06-30 23:59:59.000 35.00
> 2015-06-30 23:59:59.500 35.00
> 2015-06-30 23:59:60.000 35.00
> 2015-06-30 23:59:60.500 35.00
> 2015-06-30 23:59:60.990 35.00
> 2015-07-01 00:00:00.000 36.00
> 2015-07-01 00:00:00.500 36.00
> 2015-07-01 00:00:01.000 36.00
>  
> 
> My code didn’t work at leap second either, but this code does:
> 
> def tai_from_utc(utc):
>     """Return TAI in unix seconds, given UTC in unix seconds.
>     """
>     astropy_utc = astropy.time.Time(utc, scale="utc", format="unix")
>     dt_utc = astropy_utc.utc.to_datetime()
>     dt_tai = astropy_utc.tai.to_datetime()
>     tai_minus_utc = (dt_tai - dt_utc).total_seconds()
>     return utc + tai_minus_utc
> 
> On this, I would caution against using the astropy unix format for calculations where the details matter.  From the TimeUnix API docs <https://docs.astropy.org/en/stable/api/astropy.time.TimeUnix.html#astropy.time.TimeUnix>:
> 
> NOTE: this quantity is not exactly unix time and differs from the strict POSIX definition by up to 1 second on days with a leap second. POSIX unix time actually jumps backward by 1 second at midnight on leap second days while this class value is monotonically increasing at 86400 seconds per UTC day.

I had missed that warning. Ouch!

I think I have been approaching this the wrong way. We are representing timestamps in TAI unix seconds. I have two cases:

1) Generate a TAI unix timestamp using the current time.

I really want to get this one right.

Originally I had hoped to use `Time.now().tai.unix`. Unfortunately that gives UTC, not TAI.

My fallback was `tai_from_utc(ttime.time()) but this is looking even worse, since it’s not completely clear what time.time() will return near the leap second (at least based on the wikipedia article on unix timestamps).

The wikipedia article on unix timestamps points out that the POSIX standard specifies one behavior for system time, but that systems that use NTP typically use a different standard (though one that includes useful flags). We’ll be using PTP or NTP. Unfortunately it may not be possible to use the time, datetime, and/or astropy.time packages as none of them appear to pay attention to those flags.

2) Convert a normal unix timestamp (as a float, any other information such as leap seconds flags is gone) to a TAI unix timestamp.

In this case we could probably live with a 1 second error at or near a leap second because the information is only used to measure delays in our messaging system (DDS). But if I can get it right, so much the better.

> I'm not sure about your test case, but the Python `datetime` does not support leap seconds.  I got this exception trying to run your code within a leap second, but maybe I did something differently:
> 
> ValueError: Time (array(2015, dtype=int32), array(6, dtype=int32), array(30, dtype=int32), array(23, dtype=int32), array(59, dtype=int32), array(60, dtype=int32), array(0, dtype=int32)) is within a leap second but datetime does not support leap seconds

That is very helpful. My test was missing “fraction of a second before seconds=60”. With that I get the same error. Of course I can subtract a second from the two times when that error occurs, but it doesn’t help with the fundamental issues.

Regards,

Russell

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/astropy/attachments/20190825/d0f7e910/attachment-0001.html>