Add ViirsGroupedPacketTimeTracker

edosl0util/rdrgen.py

@@ -16,6 +16,79 @@ from edosl0util.jpssrdr import (
 from edosl0util.stream import jpss_packet_stream
 
 
+class ViirsGroupedPacketTimeTracker(object):
+    grouped_apids = list(range(800, 824)) + [825]
+
+    @classmethod
+    def handles(cls, apid):
+        return apid in cls.grouped_apids
+
+    def __init_(self):
+        self.db = {}
+
+    def get_iet(self, pkt):
+        if not self.handles(pkt.apid):
+            raise ValueError('APID {} not a VIIRS grouped packet type'.format(pkt.apid))
+        viirs_iet = self.get_viirs_iet(pkt)
+        if pkt.is_first():
+            obs_iet = get_packet_iet(pkt)
+            self.db[pkt.apid] = (obs_iet, pkt.seq_id)
+            return obs_iet
+        else:
+            last_obs_iet, last_seq = self.db[pkt.apid]
+            group_size = self.get_group_size(pkt.apid)
+            if not self.check_sequence_number(pkt.seq_id, last_seq, group_size):
+                raise OrphanedViirsPacket(pkt)
+            if not self.check_packet_iet(self.get_viirs_iet(pkt), last_obs_iet):
+                raise OrphanedViirsPacket(pkt)
+            return last_obs_iet
+
+    @staticmethod
+    def get_viirs_iet(pkt):
+        idx = 20 if pkt.is_first() else 10
+        arr = np.frombuffer(pkt.bytes()[idx:idx+8], 'B')
+        days = arr[0:2].view('>u2')[0]
+        ms = arr[2:6].view('>u4')[0]
+        us = arr[6:8].view('>u2')[0]
+        return timecode_parts_to_iet(days, ms, us)
+
+    @staticmethod
+    def get_group_size(apid):
+        i_band_apids = [813, 817, 818, 819, 820]
+        cal_apid = 825
+        if apid in i_band_apids:
+            return 32
+        elif apid == cal_apid:
+            return 24
+        else:
+            return 16
+
+    @staticmethod
+    def check_sequence_number(nonfirst_seq_num, first_seq_num, group_size):
+        seq_limit = 2**14
+        group_end = first_seq_num + group_size
+        # the 2nd check below is needed to handle wrap-around
+        return (first_seq_num < nonfirst_seq_num < group_end
+                or first_seq_num < nonfirst_seq_num + seq_limit < group_end)
+
+    @staticmethod
+    def check_packet_iet(pkt_iet, obs_iet):
+        # this can be a pretty loose check since it is only needed to prevent
+        # the very rare case where the sequence number check yields a false
+        # positive
+        permitted_lag = 5  # seconds
+        lag = (pkt_iet - obs_iet) * 1e-6
+        return 0 <= lag <= permitted_lag
+
+
+class OrphanedViirsPacket(Exception):
+    def __init__(self, pkt):
+        self.packet = pkt
+
+    def __str__(self):
+        return repr(self.packet)
+
+
 def packets_to_rdrs(sat, pkt_files):
     # FIXME: refactor!!!
     rdr_pkt_files = {}
@@ -31,13 +104,11 @@ def packets_to_rdrs(sat, pkt_files):
     for rdr_pkt_file in rdr_pkt_files.values():
         rdr_pkt_file.seek(0)
         pkts = jpss_packet_stream(rdr_pkt_file)
-        pkts = sorted(pkts, key=(lambda p: (timecode_to_iet(p.secondary_header.timecode),
-                                            p.apid)))
+        pkts = sorted(pkts, key=(lambda p: (get_packet_iet(p), p.apid)))
         blob = build_rdr_blob(sat, pkts)
         write_rdr(sat, blob)
 
 
-
 def write_rdr(sat, blob, dir_path='.',
               distributor=None, origin=None, domain=None, compressed=False, orbit_num=0,
               creation_time=None, gran_creation_time=None, software_ver=None):
@@ -154,7 +225,7 @@ def build_rdr_blob(sat, pkt_stream):
         if pkt.apid not in apid_info:
             raise ValueError(
                 'APID {} not expected for {}'.format(pkt.apid, rdr_type.short_name))
-        pkt_iet = timecode_to_iet(pkt.secondary_header.timecode)
+        pkt_iet = get_packet_iet(pkt)
         if not granule_iet <= pkt_iet < granule_iet_end:
             raise ValueError('packet stream crosses granule boundary')
         info = apid_info[pkt.apid]
@@ -314,11 +385,16 @@ def make_rdr_filename(prod_id, sat, gran_time, gran_end_time, orbit_num, creatio
         o=format_origin(origin), m=domain)
 
 
-def timecode_to_iet(timecode):
+def get_packet_iet(pkt):
+    tc = pkt.secondary_header.timecode
+    return timecode_parts_to_iet(tc.days, tc.milliseconds, tc.microseconds)
+
+
+def timecode_parts_to_iet(days, ms, us):
     # FIXME: move to timecode.py
     ccsds_epoch = Time('1958-01-01', scale='tai')
-    day = Time(ccsds_epoch.jd + timecode.days, scale='utc', format='jd') - iet_epoch
-    return int(day.sec * 1e6 + timecode.milliseconds * 1e3 + timecode.microseconds)
+    day = Time(ccsds_epoch.jd + days, scale='utc', format='jd') - iet_epoch
+    return int(day.sec * 1e6 + ms * 1e3 + us)
 
 
 def iet_to_datetime(iet):
@@ -327,7 +403,7 @@ def iet_to_datetime(iet):
 
 def calc_rdr_granule(sat, rdr_type, pkt):
     return calc_iet_granule(satellite_base_times[sat], rdr_type.gran_len,
-                            timecode_to_iet(pkt.secondary_header.timecode))
+                            get_packet_iet(pkt))
 
 
 def calc_iet_granule(base_time, gran_len, iet):