Find leaps seconds using the tz database. Convert UTC to TAI, GPS, TT(TAI)

steveha.py

#!/usr/bin/env python3
"""From
http://stackoverflow.com/questions/19332902/extract-historic-leap-seconds-from-tzdata
"""
import struct
import time

TZFILE_MAGIC = 'TZif'.encode('US-ASCII')

def leap_seconds(f):
    """
    Return a list of tuples of this format: (timestamp, number_of_seconds)
        timestamp: a 32-bit timestamp, seconds since the UNIX epoch
        number_of_seconds: how many leap-seconds occur at timestamp

    """
    fmt = ">4s c 15x 6l"
    size = struct.calcsize(fmt)
    (tzfile_magic, tzfile_format, ttisgmtcnt, ttisstdcnt, leapcnt, timecnt,
        typecnt, charcnt) =  struct.unpack(fmt, f.read(size))

    # Make sure it is a tzfile(5) file
    assert tzfile_magic == TZFILE_MAGIC, (
            "Not a tzfile; file magic was: '{}'".format(tzfile_magic))

    # comments below show struct codes such as "l" for 32-bit long integer
    offset = (timecnt*4  # transition times, each "l"
        + timecnt*1  # indices tying transition time to ttinfo values, each "B"
        + typecnt*6  # ttinfo structs, each stored as "lBB"
        + charcnt*1)  # timezone abbreviation chars, each "c"

    f.seek(offset, 1) # seek offset bytes from current position

    fmt = '>{}l'.format(leapcnt*2)
    size = struct.calcsize(fmt)
    data = struct.unpack(fmt, f.read(size))
    lst = [(data[i], data[i+1]) for i in range(0, len(data), 2)]
    assert all(lst[i][0] < lst[i+1][0] for i in range(len(lst)-1))
    assert all(lst[i][1] == lst[i+1][1]-1 for i in range(len(lst)-1))
    return lst

# From
# http://stackoverflow.com/questions/19332902/extract-historic-leap-seconds-from-tzdata
def print_leaps(leap_lst):
    print("         Leap second	Difference *after* the leap second")
    # leap_lst is tuples: (timestamp, num_leap_seconds)
    for tai_delta, (timestamp, leapcnt) in enumerate(leap_lst, start=11):
        assert tai_delta == leapcnt + 10
        time_tuple = time.gmtime(timestamp - leapcnt)
        assert time_tuple.tm_sec == 59
        print("%s\tTAI = UTC + %d" % (
            time.strftime("%Y-%m-%d %H:%M:60Z", time_tuple), tai_delta))


if __name__ == '__main__':
    # demo
    with open('/usr/share/zoneinfo/right/UTC', 'rb') as file:
        print_leaps(leap_seconds(file))
"""      Leap second	Difference *after* the leap second
1972-06-30 23:59:60Z	TAI = UTC + 11
1972-12-31 23:59:60Z	TAI = UTC + 12
1973-12-31 23:59:60Z	TAI = UTC + 13
1974-12-31 23:59:60Z	TAI = UTC + 14
1975-12-31 23:59:60Z	TAI = UTC + 15
1976-12-31 23:59:60Z	TAI = UTC + 16
1977-12-31 23:59:60Z	TAI = UTC + 17
1978-12-31 23:59:60Z	TAI = UTC + 18
1979-12-31 23:59:60Z	TAI = UTC + 19
1981-06-30 23:59:60Z	TAI = UTC + 20
1982-06-30 23:59:60Z	TAI = UTC + 21
1983-06-30 23:59:60Z	TAI = UTC + 22
1985-06-30 23:59:60Z	TAI = UTC + 23
1987-12-31 23:59:60Z	TAI = UTC + 24
1989-12-31 23:59:60Z	TAI = UTC + 25
1990-12-31 23:59:60Z	TAI = UTC + 26
1992-06-30 23:59:60Z	TAI = UTC + 27
1993-06-30 23:59:60Z	TAI = UTC + 28
1994-06-30 23:59:60Z	TAI = UTC + 29
1995-12-31 23:59:60Z	TAI = UTC + 30
1997-06-30 23:59:60Z	TAI = UTC + 31
1998-12-31 23:59:60Z	TAI = UTC + 32
2005-12-31 23:59:60Z	TAI = UTC + 33
2008-12-31 23:59:60Z	TAI = UTC + 34
2012-06-30 23:59:60Z	TAI = UTC + 35
"""

utc_to_tai.py

#!/usr/bin/env python3
from bisect import bisect_right
from calendar import timegm
from datetime import datetime, timedelta, timezone
from functools import lru_cache

# leap_seconds() from http://stackoverflow.com/a/19647860
from steveha import leap_seconds # leap epochs from tzdata

@lru_cache(maxsize=1)
def leap_second_epochs(tzfile='/usr/share/zoneinfo/right/UTC'):
    """posix_timestamp -> (TAI - UTC) delta"""
    with open(tzfile, 'rb') as file:
        L = [(timegm((1972, 1, 1, 0, 0, 0)), 10)] # TAI = UTC + 10
        L += [(ts-leapcnt+1, leapcnt + L[0][1])
              for ts, leapcnt in leap_seconds(file)]
    return tuple(zip(*L)) # transpose to return (timestamp, leapcnt)

def utc_datetime(posix_timestamp,
                 epoch=datetime(1970, 1, 1, tzinfo=timezone.utc)):
    """datetime.fromtimestamp(timestamp, timezine.utc) if gmtime(0) == 1970"""
    # it is undefined during a leap second (datetime can't represent it anyway)
    return epoch + timedelta(seconds=posix_timestamp)

def time_tai(utc_time_tuple):
    """UTC broken-down time -> TAI timestamp in seconds"""
    is_positive_leap_second = (utc_time_tuple[5] == 60)
    return utc_to_tai(timegm(utc_time_tuple)) - is_positive_leap_second

def utc_to_tai(posix_timestamp):
    """Return TAI timestamp.

    Accept dates no older than 1972-01-01

    Result for dates more than a half year into the future (starting
    with the last time the tz database was updated) may be incorrect (a
    leap second might occur)

    """
    timestamps, dATs = leap_second_epochs()
    if posix_timestamp < timestamps[0]:
        #TODO: add support for
# 1961 JAN  1 =JD 2437300.5  TAI-UTC=   1.4228180 S + (MJD - 37300.) X 0.001296 S
# 1961 AUG  1 =JD 2437512.5  TAI-UTC=   1.3728180 S + (MJD - 37300.) X 0.001296 S
# 1962 JAN  1 =JD 2437665.5  TAI-UTC=   1.8458580 S + (MJD - 37665.) X 0.0011232S
# 1963 NOV  1 =JD 2438334.5  TAI-UTC=   1.9458580 S + (MJD - 37665.) X 0.0011232S
# 1964 JAN  1 =JD 2438395.5  TAI-UTC=   3.2401300 S + (MJD - 38761.) X 0.001296 S
# 1964 APR  1 =JD 2438486.5  TAI-UTC=   3.3401300 S + (MJD - 38761.) X 0.001296 S
# 1964 SEP  1 =JD 2438639.5  TAI-UTC=   3.4401300 S + (MJD - 38761.) X 0.001296 S
# 1965 JAN  1 =JD 2438761.5  TAI-UTC=   3.5401300 S + (MJD - 38761.) X 0.001296 S
# 1965 MAR  1 =JD 2438820.5  TAI-UTC=   3.6401300 S + (MJD - 38761.) X 0.001296 S
# 1965 JUL  1 =JD 2438942.5  TAI-UTC=   3.7401300 S + (MJD - 38761.) X 0.001296 S
# 1965 SEP  1 =JD 2439004.5  TAI-UTC=   3.8401300 S + (MJD - 38761.) X 0.001296 S
# 1966 JAN  1 =JD 2439126.5  TAI-UTC=   4.3131700 S + (MJD - 39126.) X 0.002592 S
# 1968 FEB  1 =JD 2439887.5  TAI-UTC=   4.2131700 S + (MJD - 39126.) X 0.002592 S
# 1972 JAN  1 =JD 2441317.5  TAI-UTC=  10.0       S + (MJD - 41317.) X 0.0      S
        # e.g., 1970 is timedelta(seconds=8, microseconds=82)
        # see http://maia.usno.navy.mil/ser7/tai-utc.dat
        raise ValueError("Dates before %s are not supported, got %s" % (
            utc_datetime(timestamps[0]), utc_datetime(posix_timestamp)))
    return posix_timestamp + dATs[bisect_right(timestamps, posix_timestamp)-1]

def utc_to_gps(posix_timestamp):
    """Return GPS time = TAI - 19 seconds"""
    return utc_to_tai(posix_timestamp) - 19

def utc_to_tt(posix_timestamp):
    """Return Terrestrial Time (TAI realization) -- TT(TAI) = TAI + 32.184s.

    http://en.wikipedia.org/wiki/Terrestrial_Time
    """
    return utc_to_tai(posix_timestamp) + 32.184

if __name__ == "__main__":
    import time

    print('%25s: %+10s -> %+10s' % (
        "UTC time", "POSIX", "TAI time"))
    for i, tt in enumerate([
            (2012, 6, 30, 23, 59, 60), # leap second
            (1970, 1, 1, 0, 0, 0),
            (1971, 12, 31, 23, 59, 59),
            (1972, 1, 1, 0, 0, 0),
            (1972, 7, 1, 0, 0, 0),
            (2006, 1, 1, 0, 0, 0),
            (2012, 6, 30, 23, 59, 59),
            (2012, 7, 1, 0, 0, 0),
            (2012, 7, 1, 0, 0, 1),
            (2012, 7, 1, 0, 0, 2),
    ]):
        ts = timegm(tt)
        s = str(utc_datetime(ts))
        leap_second = (i == 0)
        if leap_second:
            s = "2012-06-30 23:59:60+00:00"
        try:
            tai_ts = utc_to_tai(ts) - leap_second
            assert time_tai(tt) == tai_ts
            print('%s: %+10s -> %+10s' % (s, ts, tai_ts))
        except ValueError as e:
            print(e)
    t0 = timegm((1972, 1, 1, 0, 0, 0))
    assert ((utc_to_tai(ts) - utc_to_tai(t0) - (ts - t0))) == 25
    """
                     UTC time:      POSIX ->   TAI time
    2012-06-30 23:59:60+00:00: 1341100800 -> 1341100834
    Dates before 1972-01-01 00:00:00+00:00 are not supported, got 1970-01-01 00:00:00+00:00
    Dates before 1972-01-01 00:00:00+00:00 are not supported, got 1971-12-31 23:59:59+00:00
    1972-01-01 00:00:00+00:00:   63072000 ->   63072010
    1972-07-01 00:00:00+00:00:   78796800 ->   78796811
    2006-01-01 00:00:00+00:00: 1136073600 -> 1136073633
    2012-06-30 23:59:59+00:00: 1341100799 -> 1341100833
    2012-07-01 00:00:00+00:00: 1341100800 -> 1341100835
    2012-07-01 00:00:01+00:00: 1341100801 -> 1341100836
    2012-07-01 00:00:02+00:00: 1341100802 -> 1341100837
    """
    # example from http://cr.yp.to/proto/utctai.html
    old_ts = None
    for utc_s, expected in [
            ("1997-06-30 23:59:58 UTC", "1997-07-01 00:00:28 TAI"),
            ("1997-06-30 23:59:59 UTC", "1997-07-01 00:00:29 TAI"),
            ("1997-06-30 23:59:60 UTC", "1997-07-01 00:00:30 TAI"), # leap
            ("1997-07-01 00:00:00 UTC", "1997-07-01 00:00:31 TAI"),
    ]:
        tt = time.strptime(utc_s, "%Y-%m-%d %H:%M:%S UTC")
        # seconds since Epoch along can't represent a leap second
        is_leap_second = ("23:59:60" in utc_s)
        tai_ts = utc_to_tai(timegm(tt)) - is_leap_second
        assert tai_ts == time_tai(tt)
        if old_ts is not None:
            assert (tai_ts - old_ts) == 1, (old_ts, tai_ts)
        old_ts = tai_ts
        tai_time = time.gmtime(tai_ts)
        tai_s = time.strftime("%Y-%m-%d %H:%M:%S TAI", tai_time)
        assert tai_s == expected, (tai_s, expected)