diff --git a/edosl0util/rdrgen.py b/edosl0util/rdrgen.py index 49eb7424265a9982d59c007de6d39a1914322c11..b16fdcfad47a212c22d5e11fabc3900c0e08985f 100644 --- a/edosl0util/rdrgen.py +++ b/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):