From 60d3122bd5b0e54dae6d2eda38da50a229f7b540 Mon Sep 17 00:00:00 2001
From: Ryan Kingsbury <ryan.kingsbury@planet.com>
Date: Thu, 28 Apr 2016 15:26:59 +0000
Subject: [PATCH] Use OS tzdata for leap second calculation.

The tzdata database, also known as the Olsen, database seems to be the
definitive source for leap second data.  It is maintained regularly by folks
who follow this stuff.  The suggestion to use this database came from:
https://stackoverflow.com/questions/19332902/extract-historic-leap-seconds-from-tzdata

An alterate approach would be to include the tz repository as a submodule but
this would quickly go stale.  Therefore, I opted to pull the file from tzdata
which should be regularly updated by the OS.  The current location is suitable
for Ubuntu, not sure about other OSes.
---
 peregrine/gps_time.py | 90 +++++++++++++++++++++++++++++++++++++------
 1 file changed, 78 insertions(+), 12 deletions(-)

diff --git a/peregrine/gps_time.py b/peregrine/gps_time.py
index 9649fa0..7aac347 100644
--- a/peregrine/gps_time.py
+++ b/peregrine/gps_time.py
@@ -1,5 +1,78 @@
 import datetime
 
+# NTP timestamps are in units of seconds since the NTP epoch
+# See leap-seconds.list for further details.
+NTP_EPOCH = datetime.datetime(1900,1,1,0,0,0)
+
+def sec_since_NTP_epoch(dt):
+    """Return the number of seconds since the NTP epoch."""
+    return (dt-NTP_EPOCH).total_seconds()
+
+def load_leapsecond_table(f="/usr/share/zoneinfo/leap-seconds.list"):
+    """"
+    Loads leap second table from system.
+
+    The default file location is appropriate for Ubuntu and is provided by
+    the tzdata package.  Refer to the leap-seconds.list file for formatting
+    information.
+
+    Parameters
+    ----------
+    f : string
+      Path to a leap-seconds.list file
+
+    Returns: List of tuples in chronological order each containing an
+    epoch start time and number of leap seconds applicable for that epoch.
+    """
+
+    # Check the expiration date in the file
+    with open(f,'r') as fp:
+        for line in fp:
+            if line[0:2] != "#@":
+                continue
+
+            expiration_time = int(line.split('\t')[-1])
+            now = sec_since_NTP_epoch(datetime.datetime.now())
+
+            if expiration_time<now:
+                raise UserWarning("leap-seconds.list file has expired, update tzdata")
+
+    # Load the load the table
+    with open(f,'r') as fp:
+        table = []
+        for line in fp:
+            if line[0]=="#":
+                continue
+            raw = line.split('\t')
+            table.append( (int(raw[0]), int(raw[1])) )
+
+    # Add the expiration time to the end of the table
+    table.append( (expiration_time, None) )
+
+    return table
+
+def get_gpst_leap_seconds(dt):
+    """Returns the number of leap seconds at the time provided."""
+
+    # Load leap second data from system tzdata file
+    lstable = load_leapsecond_table()
+
+    t = sec_since_NTP_epoch(dt)
+
+    ls = None
+    for epoch_start,leapseconds in lstable:
+        if t>=epoch_start:
+            ls = leapseconds
+
+    # Raise warning if specified time is after expiry date of leap second table
+    if ls==None:
+        raise UserWarning("Specified datetime is after expiry time of leap second table.")
+
+    # GPS leap seconds relative to a Jan 1 1980, where TAI-UTC was 19 seconds.
+    gps_leap_seconds = ls - 19
+
+    return datetime.timedelta(seconds = gps_leap_seconds)
+
 def datetime_to_tow(t, mod1024 = True):
     """
     Convert a Python datetime object to GPS Week and Time Of Week.
@@ -39,12 +112,7 @@ def gpst_to_utc(t_gpst):
     -------
     t_utc : datetime
     """
-    t_utc = t_gpst - datetime.timedelta(seconds = 16)
-    if t_utc <= datetime.datetime(2012, 7, 1) or \
-       t_utc >= datetime.datetime(2015, 7, 1):
-        raise ValueError("Don't know how many leap seconds to use.  " +
-                         "Please implement something better in gpst_to_utc()")
-    return t_utc
+    return t_gpst - get_gpst_leap_seconds(t_gpst)
 
 def utc_to_gpst(t_utc):
     """
@@ -59,9 +127,7 @@ def utc_to_gpst(t_utc):
     -------
     t_gpst : datetime
     """
-    t_gpst = t_utc + datetime.timedelta(seconds = 16)
-    if t_utc <= datetime.datetime(2012, 7, 1) or \
-       t_utc >= datetime.datetime(2015, 7, 1):
-        raise ValueError("Don't know how many leap seconds to use.  " +
-                         "Please implement something better in utc_to_gpst()")
-    return t_gpst
+    # TODO: there may be an unhandled corner case here for conversions
+    # where the GPST-to-UTCT interval spans a leap second.
+
+    return t_utc + get_gpst_leap_seconds(t_utc)