rdrgen.py

import ctypes
import os
import tempfile
from collections import OrderedDict
from contextlib import contextmanager
from datetime import datetime

import attr
import h5py
import numpy as np

import edosl0util
from edosl0util import compat
from edosl0util.jpssrdr import (
    StaticHeader, Apid as ApidListItem, PacketTracker, decode_rdr_blob)
from edosl0util.stream import jpss_packet_stream
from edosl0util.timecode import cds_to_iet, iet_to_dt


def packets_to_rdrs(sat, l0_files, **kwargs):
    def iter_pkts(l0_files):
        for l0_file in l0_files:
            with open(l0_file, 'rb') as l0_file_obj:
                for pkt in jpss_packet_stream(l0_file_obj):
                    yield pkt

    build_rdr(sat, iter_pkts(l0_files), **kwargs)


def build_rdr(sat, pkt_iter, output_dir='.', aggr_type='idps', aggr_level=None,
              diary_cushion=10000000, attr_overrides={}):
    """Construct RDR file(s) from L0 packets

    Default aggregation behavior uses file boundaries computed in the same way
    as for IDPS assuming the default aggregation level depending on the instrument
    (e.g. 85 sec for VIIRS, 8 min for CrIS). A different aggregation level can
    be given as an integer number of granule (e.g. 4 to get roughly 5.5 min files
    for VIIRS). Finally, aggr_type can be set to 'full' to request that a single
    file gets produced (thus aggregating all granules containing the given
    packets).
    """

    # divy packets up into temp files organized by granule
    file_mgr = BinnedTemporaryFileManager()
    get_jpss_packet_time = GetJpssPacketTime()
    gran_infos = set()  # (rdr_type, gran_iet) pairs
    for pkt in pkt_iter:
        rdr_type = get_rdr_type(pkt.apid)
        pkt_iet = get_jpss_packet_time(pkt)
        gran_iet = get_granule_start(sat, rdr_type.gran_len, pkt_iet)
        gran_info = (rdr_type, gran_iet)
        gran_infos.add(gran_info)
        file_mgr.add_data(gran_info, pkt.bytes())

    # determine what RDR files we'll be producing based on the packets we've seen
    rdr_types = set(rdr_type for rdr_type, gran_iet in gran_infos)
    primary_type, packaged_type = process_rdr_types(rdr_types, force_packaging=False)
    rdr_types = sorted(rdr_types, key=(lambda t: 1 if t is primary_type else 2))
    if aggr_type == 'idps':
        aggr_level = aggr_level or primary_type.default_aggregation
        primary_aggr_iets = sorted(set(
            get_aggregate_start(sat, primary_type.gran_len, aggr_level, gran_iet)
            for (rdr_type, gran_iet) in gran_infos if rdr_type is primary_type))
    elif aggr_type == 'full':
        # produce a single output file, ignoring IDPS-style aggregation boundaries
        assert aggr_level is None
        first_gran_iet = min(gran_iet for (rdr_type, gran_iet) in gran_infos
                             if rdr_type is primary_type)
        last_gran_iet = max(gran_iet for (rdr_type, gran_iet) in gran_infos
                             if rdr_type is primary_type)
        aggr_level = (last_gran_iet - first_gran_iet) // primary_type.gran_len + 1
        primary_aggr_iets = [first_gran_iet]
    else:
        raise ValueError('aggr_type must be idps or input')

    # now generate the RDRs
    rdr_files = []
    for aggr_iet in primary_aggr_iets:
        rdr_writer = RdrWriter(sat, rdr_types, aggr_iet, aggr_level, output_dir,
                               **attr_overrides)
        rdr_writer.write_aggregate(primary_type, aggr_iet, aggr_level)
        gran_iets = [aggr_iet + i * primary_type.gran_len for i in range(aggr_level)]
        for gran_iet in gran_iets:
            with file_mgr.process_file((primary_type, gran_iet)) as pkt_file:
                pkts = list(jpss_packet_stream(pkt_file))
            blob = build_rdr_blob(sat, pkts, primary_type, gran_iet)
            rdr_writer.write_granule(primary_type, gran_iet, blob)
        if packaged_type:
            packaged_gran_iets = get_overlapping_granules(
                sat, packaged_type.gran_len, aggr_iet - diary_cushion,
                                             aggr_iet + aggr_level * primary_type.gran_len + diary_cushion + 1)
            rdr_writer.write_aggregate(
                packaged_type, packaged_gran_iets[0], len(packaged_gran_iets))
            for gran_iet in packaged_gran_iets:
                with file_mgr.process_file((packaged_type, gran_iet)) as pkt_file:
                    pkts = list(jpss_packet_stream(pkt_file))
                blob = build_rdr_blob(sat, pkts, packaged_type, gran_iet)
                rdr_writer.write_granule(packaged_type, gran_iet, blob)
        rdr_writer.close()
        rdr_files.append(rdr_writer.file_name)
    file_mgr.clean_up()

    return rdr_files


def process_rdr_types(given_rdr_types, force_packaging):
    """Determine the RDR type we'll be making based on the packets we've been given

    Return both a primary_type and packaged_type. The primary type could indicate
    either a science or spacecraft RDR. The packaged type is used when spacecraft
    packets are to be packaged along with a science RDR; it will be None if that's
    not the case. forceA_packaging is a boolean that if set will force spacecraft
    RDR structures to be written out even if no spacecraft packets are present.
    """
    rdr_types = set(given_rdr_types)
    if not rdr_types:
        raise RdrgenError('No RDR types to generate!')
    sci_type = {t for t in rdr_types if t.type_id == 'SCIENCE'}
    rdr_types -= sci_type
    if len(sci_type) > 1:
        raise RdrgenError(
            'Cannot process more than 1 science RDR type at once '
            + '(have {})'.format(', '.join(t.short_name for t in sci_type)))
    diary_type = {t for t in rdr_types if t is SpacecraftDiaryRdrType}
    rdr_types -= diary_type
    if rdr_types:
        raise RdrgenError(
            'Unsupported RDR type(s): ' + ', '.join(t.short_name for t in rdr_types))
    if sci_type:
        primary_type, = sci_type
        packaged_type = (
            SpacecraftDiaryRdrType if (diary_type or force_packaging) else None)
    else:
        primary_type = SpacecraftDiaryRdrType
        packaged_type = None
    return primary_type, packaged_type


class BinnedTemporaryFileManager(object):
    """Manage a set of append-mode temporary files each labeled with a 'bin key'

    Limits the number of file handles kept open at a time.
    """

    def __init__(self, parent_dir='.', max_open_files=32):
        self.max_open_files = max_open_files
        self.dir = tempfile.mkdtemp(dir=parent_dir)
        self._file_paths = {}
        self._file_objs = OrderedDict()

    def add_data(self, bin_key, data):
        """Append some data to the temp file associated with bin_key

        Creates the file if needed.
        """
        file_obj = self._file_objs.pop(bin_key, None)
        if not file_obj:
            file_path = self._file_paths.get(bin_key)
            if file_path:
                file_obj = open(file_path, 'a+b')
            else:
                file_obj = tempfile.NamedTemporaryFile(dir=self.dir, delete=False, mode='wb')
                file_path = file_obj.name
                self._file_paths[bin_key] = file_path
            if len(self._file_objs) == self.max_open_files:
                _, old_file_obj = self._file_objs.popitem(last=False)
                old_file_obj.close()
        self._file_objs[bin_key] = file_obj
        file_obj.write(data)

    @contextmanager
    def process_file(self, bin_key):
        """'Check out' a file for processing

        Returns a context manager that provides a read-mode file handle and
        removes the file after processing. Can be called with a bin_key that
        has no data yet which just results in an empty file.
        """
        file_obj = self._file_objs.pop(bin_key, None)
        if file_obj:
            file_obj.close()
        file_path = self._file_paths.pop(bin_key, '/dev/null')
        file_obj = open(file_path, 'rb')
        try:
            yield file_obj
        finally:
            file_obj.close()
            if file_path != '/dev/null':
                os.remove(file_path)

    def clean_up(self):
        """Call after all files are processed to avoid leaving a dir sitting around"""
        os.rmdir(self.dir)


class RdrWriter(object):
    def __init__(self, sat, rdr_types, aggr_iet, aggr_level, dir_path='.',
                 distributor=None, origin=None, domain=None, compressed=False,
                 orbit_num=0, creation_time=None, software_ver=None):
        self._sat = sat
        origin = origin or self.default_origin
        distributor = distributor or origin
        self._domain = domain or self.default_domain
        self._orbit_num = orbit_num
        self._creation_time = creation_time or datetime.now()
        self._software_ver = (
                software_ver or edosl0util.__name__ + '-' + edosl0util.__version__)

        aggr_end_iet = aggr_iet + aggr_level * rdr_types[0].gran_len
        self.file_name = make_rdr_filename(
            rdr_types, self._sat, aggr_iet, aggr_end_iet, self._orbit_num,
            self._creation_time, origin, self._domain, compressed)
        self._h5_file = h5py.File(os.path.join(dir_path, self.file_name), 'w')
        self._write_skeleton(rdr_types, distributor, origin)
        self._aggr_starts = {}

    def _write_skeleton(self, rdr_types, distributor, origin):
        self._set_h5_attrs(self._h5_file, {
            'Distributor': distributor,
            'Mission_Name': 'S-NPP/JPSS',  # FIXME: what will this be for J1?
            'N_Dataset_Source': origin,
            'N_HDF_Creation_Date': self._format_date_attr(self._creation_time),
            'N_HDF_Creation_Time': self._format_time_attr(self._creation_time),
            'Platform_Short_Name': platform_short_names[self._sat]})
        all_grp = self._h5_file.create_group('All_Data')
        prod_grp = self._h5_file.create_group('Data_Products')
        for rdr_type in rdr_types:
            all_grp.create_group(rdr_type.short_name + '_All')
            gran_grp = prod_grp.create_group(rdr_type.short_name)
            self._set_h5_attrs(gran_grp, {
                'Instrument_Short_Name': instrument_short_names[rdr_type.sensor],
                'N_Collection_Short_Name': rdr_type.short_name,
                'N_Dataset_Type_Tag': 'RDR',
                'N_Processing_Domain': self._domain})

    def write_aggregate(self, rdr_type, aggr_iet, num_grans):
        self._aggr_starts[rdr_type] = aggr_iet
        grp = self._h5_file['Data_Products'][rdr_type.short_name]
        ds = grp.create_dataset(
            rdr_type.short_name + '_Aggr', [1], self._h5_ref_dtype, maxshape=[None])
        ds[0] = self._h5_file['All_Data/{}_All'.format(rdr_type.short_name)].ref
        aggr_end_iet = aggr_iet + num_grans * rdr_type.gran_len
        last_gran_iet = aggr_end_iet - rdr_type.gran_len
        self._set_h5_attrs(ds, {
            'AggregateBeginningDate': self._format_date_attr(aggr_iet),
            'AggregateBeginningGranuleID': make_granule_id(self._sat, aggr_iet),
            'AggregateBeginningOrbitNumber': np.uint64(self._orbit_num),
            'AggregateBeginningTime': self._format_time_attr(aggr_iet),
            'AggregateEndingDate': self._format_date_attr(aggr_end_iet),
            'AggregateEndingGranuleID': make_granule_id(self._sat, last_gran_iet),
            'AggregateEndingOrbitNumber': np.uint64(self._orbit_num),
            'AggregateEndingTime': self._format_time_attr(aggr_end_iet),
            'AggregateNumberGranules': np.uint64(num_grans)})

    def write_granule(self, rdr_type, gran_iet, blob, creation_time=None):
        raw_grp = self._h5_file['All_Data/{}_All'.format(rdr_type.short_name)]
        gran_idx = (gran_iet - self._aggr_starts[rdr_type]) // rdr_type.gran_len
        raw_ds = raw_grp.create_dataset(
            'RawApplicationPackets_{}'.format(gran_idx), data=blob, maxshape=[None])
        gran_grp = self._h5_file['Data_Products'][rdr_type.short_name]
        gran_ds = gran_grp.create_dataset(
            rdr_type.short_name + '_Gran_{}'.format(gran_idx),
            [1], self._h5_regionref_dtype, maxshape=[None])
        gran_ds[0] = raw_ds.regionref[:]
        gran_end_iet = gran_iet + rdr_type.gran_len
        creation_time = creation_time or self._creation_time or datetime.now()
        gran_id = make_granule_id(self._sat, gran_iet)
        gran_ver = 'A1'
        blob_info = decode_rdr_blob(blob)
        self._set_h5_attrs(gran_ds, {
            'Beginning_Date': self._format_date_attr(gran_iet),
            'Beginning_Time': self._format_time_attr(gran_iet),
            'Ending_Date': self._format_date_attr(gran_end_iet),
            'Ending_Time': self._format_time_attr(gran_end_iet),
            'N_Beginning_Orbit_Number': np.uint64(self._orbit_num),
            'N_Beginning_Time_IET': np.uint64(gran_iet),
            'N_Creation_Date': self._format_date_attr(creation_time),
            'N_Creation_Time': self._format_time_attr(creation_time),
            'N_Ending_Time_IET': np.uint64(gran_end_iet),
            'N_Granule_ID': gran_id,
            'N_Granule_Status': 'N/A',
            'N_Granule_Version': gran_ver,
            'N_IDPS_Mode': self._domain,
            'N_JPSS_Document_Ref': rdr_type.document,
            'N_LEOA_Flag': 'Off',
            'N_Packet_Type': [a.name for a in blob_info.apids],
            'N_Packet_Type_Count': [np.uint64(a.pkts_received) for a in blob_info.apids],
            'N_Percent_Missing_Data': np.float32(self._calc_percent_missing(blob_info)),
            'N_Primary_Label': 'Primary',  # TODO: find out what this is
            'N_Reference_ID': ':'.join([rdr_type.short_name, gran_id, gran_ver]),
            'N_Software_Version': self._software_ver})

    def close(self):
        self._h5_file.close()

    @staticmethod
    def _set_h5_attrs(h5_obj, attrs):
        attrs = _encodeall(attrs)
        # for some reason all the H5 attributes are in 2-D arrays in IDPS files
        for name, value in attrs.items():
            if isinstance(value, list):
                value = [[x] for x in value]
            else:
                value = [[value]]
            h5_obj.attrs[name] = value

    @staticmethod
    def _format_date_attr(t):
        return iet_to_datetime(t).strftime('%Y%m%d')

    @staticmethod
    def _format_time_attr(t):
        return iet_to_datetime(t).strftime('%H%M%S.%fZ')

    @staticmethod
    def _calc_percent_missing(common_rdr):
        total_received = sum(a.pkts_received for a in common_rdr.apids)
        total_reserved = sum(a.pkts_reserved for a in common_rdr.apids)
        return 100. * (total_reserved - total_received) / total_reserved

    default_origin = 'ssec'
    default_domain = 'dev'
    _h5_ref_dtype = h5py.special_dtype(ref=h5py.Reference)
    _h5_regionref_dtype = h5py.special_dtype(ref=h5py.RegionReference)


def _encodeall(v):
    """
    Recurse v encoding stringy values. dict keys are not encoded.
    """
    if isinstance(v, dict):
        dtmp = {}
        for k, v in v.items():
            dtmp[k] = _encodeall(v)
        return dtmp
    if isinstance(v, list):
        return [_encodeall(x) for x in v]
    if isinstance(v, compat.str_types):
        return v.encode()
    return v


def build_rdr_blob(sat, pkt_stream, rdr_type, granule_iet):
    get_jpss_packet_time = GetJpssPacketTime()
    granule_iet_end = granule_iet + rdr_type.gran_len

    total_pkt_size = 0
    apid_info = OrderedDict()
    total_trackers = 0
    all_pkts = []
    for apid in rdr_type.apids:
        apid_info[apid.num] = {
            'name': apid.name.encode(),
            'pkts_reserved': apid.max_expected,
            'pkts_received': 0,
            'first_tracker_index': total_trackers,
            'pkt_info': [{} for _ in range(apid.max_expected)]}
        total_trackers += apid.max_expected

    for pkt in pkt_stream:
        if pkt.apid not in apid_info:
            raise ValueError(
                'APID {} not expected for {}'.format(pkt.apid, rdr_type.short_name))
        pkt_iet = get_jpss_packet_time(pkt)
        if not granule_iet <= pkt_iet < granule_iet_end:
            raise ValueError('packet stream crosses granule boundary')
        info = apid_info[pkt.apid]
        pkt_info = info['pkt_info'][info['pkts_received']]
        pkt_info['obs_time'] = pkt_iet
        pkt_info['seq_num'] = pkt.seqid
        pkt_info['size'] = pkt.size
        pkt_info['offset'] = total_pkt_size
        info['pkts_received'] += 1
        total_pkt_size += pkt.size
        all_pkts.append(pkt)

    apid_list_offset = ctypes.sizeof(StaticHeader)
    pkt_tracker_offset = apid_list_offset + len(apid_info) * ctypes.sizeof(ApidListItem)
    ap_storage_offset = pkt_tracker_offset + total_trackers * ctypes.sizeof(PacketTracker)
    buf_size = ap_storage_offset + total_pkt_size
    buf = np.zeros([buf_size], np.uint8)  # zeros needed to null-pad strings

    header = StaticHeader.from_buffer(buf)
    header.satellite = platform_short_names[sat].encode()
    header.sensor = instrument_short_names[rdr_type.sensor].encode()
    header.type_id = rdr_type.type_id.encode()
    header.num_apids = len(apid_info)
    header.apid_list_offset = apid_list_offset
    header.pkt_tracker_offset = pkt_tracker_offset
    header.ap_storage_offset = ap_storage_offset
    header.next_pkt_pos = total_pkt_size
    header.start_boundary = granule_iet
    header.end_boundary = granule_iet_end

    for i, (apid, info) in enumerate(apid_info.items()):
        offset = header.apid_list_offset + i * ctypes.sizeof(ApidListItem)
        item = ApidListItem.from_buffer(buf, offset)
        item.name = info['name']
        item.value = apid
        item.pkt_tracker_start_idx = info['first_tracker_index']
        item.pkts_reserved = info['pkts_reserved']
        item.pkts_received = info['pkts_received']

        for j, pkt_info in enumerate(info['pkt_info']):
            offset = (header.pkt_tracker_offset
                      + (info['first_tracker_index'] + j) * ctypes.sizeof(PacketTracker))
            tracker = PacketTracker.from_buffer(buf, offset)
            if pkt_info:
                tracker.obs_time = pkt_info['obs_time']
                tracker.sequence_number = pkt_info['seq_num']
                tracker.size = pkt_info['size']
                tracker.offset = pkt_info['offset']
                tracker.fill_percent = 0
            else:
                tracker.offset = -1

    buf[ap_storage_offset:] = bytearray().join(pkt.bytes() for pkt in all_pkts)

    return buf


class ViirsScienceApidInfo(object):
    apids = list(x for x in range(800, 829) if x != 824)
    names = ['M04', 'M05', 'M03', 'M02', 'M01', 'M06', 'M07', 'M09', 'M10',
             'M08', 'M11', 'M13', 'M12', 'I04', 'M16', 'M15', 'M14', 'I05',
             'I01', 'I02', 'I03', 'DNB', 'DNB_MGS', 'DNB_LGS',
             'CAL', 'ENG', 'DNB_HGA', 'DNB_HGB']

    @classmethod
    def get_specs(cls):
        return [ApidSpec(apid, cls.get_name(apid), cls.get_max_expected(apid))
                for apid in cls.apids]

    @classmethod
    def get_name(cls, apid):
        return cls.names[cls.apids.index(apid)]

    @classmethod
    def get_max_expected(cls, apid):
        max_scans = 48
        return max_scans * cls.get_packets_per_scan(apid)

    @classmethod
    def get_packets_per_scan(cls, apid):
        name = cls.get_name(apid)
        if name == 'ENG':
            return 1
        elif name == 'CAL':
            return 24
        elif name.startswith('M') or name.startswith('DNB'):
            return 17
        else:
            return 33


class CrisScienceApidInfo(object):
    apids = [1289, 1290] + list(range(1315, 1396))

    @classmethod
    def get_specs(cls):
        return [ApidSpec(apid, cls.get_name(apid), cls.get_max_expected(apid))
                for apid in cls.apids]

    @classmethod
    def get_name(cls, apid):
        if apid == 1289:
            return 'EIGHT_S_SCI'
        elif apid == 1290:
            return 'ENG'
        else:
            offset = apid - 1315
            view_types = ['N', 'S', 'C']
            bands = ['LW', 'MW', 'SW']
            num_fovs = 9
            view_type = view_types[offset // (num_fovs * len(bands))]
            band = bands[offset // num_fovs % len(bands)]
            fov = str(offset % num_fovs + 1)
            return view_type + band + fov

    @classmethod
    def get_max_expected(cls, apid):
        name = cls.get_name(apid)
        if name == 'EIGHT_S_SCI':
            return 5
        elif name == 'ENG':
            return 1
        else:
            view_type = name[0]
            if view_type == 'N':
                return 121
            else:
                return 9


@attr.s
class ApidSpec(object):
    num = attr.ib()
    name = attr.ib()
    max_expected = attr.ib()


class RdrTypeManager(object):
    def __init__(self):
        self._types_by_apid = {}

    def register_type(self, cls):
        for apid_spec in cls.apids:
            if apid_spec.num in self._types_by_apid:
                raise ValueError('each APID can only be handled by one RDR type')
            self._types_by_apid[apid_spec.num] = cls
        return cls

    def get_type_for_apid(self, apid):
        try:
            return self._types_by_apid[apid]
        except KeyError:
            raise RdrgenError('Unsupported APID: {}'.format(apid))


rdr_type_mgr = RdrTypeManager()
rdr_type_spec = rdr_type_mgr.register_type
get_rdr_type = rdr_type_mgr.get_type_for_apid


@rdr_type_spec
class ViirsScienceRdrType(object):
    product_id = 'RVIRS'
    short_name = 'VIIRS-SCIENCE-RDR'
    gran_len = 85350000
    sensor = 'viirs'
    type_id = 'SCIENCE'
    document = '474-00448-02-06_JPSS-DD-Vol-II-Part-6_0200G.pdf'
    apids = ViirsScienceApidInfo.get_specs()
    default_aggregation = 1


@rdr_type_spec
class CrisScienceRdrType(object):
    product_id = 'RCRIS'
    short_name = 'CRIS-SCIENCE-RDR'
    gran_len = 31997000
    sensor = 'cris'
    type_id = 'SCIENCE'
    document = '474-00448-02-03_JPSS-DD-Vol-II-Part-3_0200B.pdf'
    apids = CrisScienceApidInfo.get_specs()
    default_aggregation = 15


@rdr_type_spec
class AtmsScienceRdrType(object):
    product_id = 'RATMS'
    short_name = 'ATMS-SCIENCE-RDR'
    gran_len = 31997000
    sensor = 'atms'
    type_id = 'SCIENCE'
    document = '474-00448-02-02_JPSS-DD-Vol-II-Part-2_0200B.pdf'
    apids = [ApidSpec(515, 'CAL', max_expected=5),
             ApidSpec(528, 'SCI', max_expected=1249),
             ApidSpec(530, 'ENG_TEMP', max_expected=13),
             ApidSpec(531, 'ENG_HS', max_expected=5)]
    default_aggregation = 15


@rdr_type_spec
class SpacecraftDiaryRdrType(object):
    product_id = 'RNSCA'
    short_name = 'SPACECRAFT-DIARY-RDR'
    gran_len = 20000000
    sensor = None
    type_id = 'DIARY'
    document = '474-00448-02-08_JPSS-DD-Vol-II-Part-8_0200F.pdf'
    apids = [ApidSpec(0, 'CRITICAL', max_expected=21),
             ApidSpec(8, 'ADCS_HKH', max_expected=21),
             ApidSpec(11, 'DIARY', max_expected=21)]
    default_aggregation = 303


class GetJpssPacketTime(object):
    def __init__(self):
        self._viirs_tracker = ViirsGroupedPacketTimeTracker()

    def __call__(self, pkt):
        if self._viirs_tracker.tracks_apid(pkt.apid):
            return self._viirs_tracker.get_iet(pkt)
        else:
            return get_packet_iet(pkt)


class ViirsGroupedPacketTimeTracker(object):
    grouped_apids = list(range(800, 824)) + [825]

    @classmethod
    def tracks_apid(cls, apid):
        return apid in cls.grouped_apids

    def __init__(self):
        self._db = {}

    def get_iet(self, pkt):
        if not self.tracks_apid(pkt.apid):
            raise ValueError('APID {} not a VIIRS grouped packet type'.format(pkt.apid))
        if pkt.is_first():
            obs_iet = get_packet_iet(pkt)
            self._db[pkt.apid] = (obs_iet, pkt.seqid)
            return obs_iet
        else:
            last_obs_iet, last_seq = self._db[pkt.apid]
            group_size = ViirsScienceApidInfo.get_packets_per_scan(pkt.apid)
            if not self.check_sequence_number(pkt.seqid, 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):
        if pkt.is_standalone():
            idx = 18
        elif pkt.is_first():
            idx = 20
        else:
            idx = 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 cds_to_iet(days, ms, us)

    @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 RdrgenError(Exception):
    pass


class OrphanedViirsPacket(RdrgenError):
    def __init__(self, pkt):
        self.packet = pkt

    def __str__(self):
        return repr(self.packet)


def make_rdr_filename(rdr_types, sat, aggr_begin, aggr_end, orbit_num, creation_time,
                      origin, domain, compressed):
    aggr_begin = iet_to_datetime(aggr_begin)
    aggr_end = iet_to_datetime(aggr_end)
    prod_ids = '-'.join(sorted(t.product_id for t in rdr_types))
    sat = {'snpp': 'npp', 'j01': 'j01'}[sat]
    if origin.endswith('-'):
        origin = origin[:-1] + ('c' if compressed else 'u')

    def format_time(t):
        return t.strftime('%H%M%S') + str(t.microsecond // 100000)

    return '{p}_{s}_d{d:%Y%m%d}_t{b}_e{e}_b{n:05d}_c{c:%Y%m%d%H%M%S%f}_{o}_{m}.h5'.format(
        p=prod_ids, s=sat, d=aggr_begin, b=format_time(aggr_begin),
        e=format_time(aggr_end), n=orbit_num, c=creation_time, o=origin, m=domain)


def make_granule_id(sat, gran_iet):
    tenths = (gran_iet - satellite_base_times[sat]) // 100000
    return '{}{:012d}'.format(platform_short_names[sat], tenths)


def get_packet_iet(pkt):
    tc = pkt.secondary_header.timecode
    return cds_to_iet(tc.days, tc.milliseconds, tc.microseconds)


def iet_to_datetime(iet):
    if isinstance(iet, datetime):
        return iet
    return iet_to_dt(iet)


def get_granule_start(sat, gran_len, iet):
    base_time = satellite_base_times[sat]
    return (iet - base_time) // gran_len * gran_len + base_time


def get_aggregate_start(sat, gran_len, grans_per_aggr, iet):
    # see "DDS Aggregation Methodology" in CDFCB vol I
    Y = gran_len
    N = grans_per_aggr
    G_s = get_granule_start(sat, Y, iet)
    A_n = G_s // (N * Y)
    O = G_s % Y
    A_s = A_n * (Y * N) + O
    return A_s


def get_aggregate_granule_times(sat, gran_len, aggr_level, aggr_iet):
    aggr_start = get_aggregate_start(sat, gran_len, aggr_level, aggr_iet)
    aggr_end = aggr_start + aggr_level * gran_len
    return get_overlapping_granules(sat, gran_len, aggr_start, aggr_end)


def get_overlapping_granules(sat, gran_len, start_iet, stop_iet):
    rv = []
    gran_iet = get_granule_start(sat, gran_len, start_iet)
    while gran_iet < stop_iet:
        rv.append(gran_iet)
        gran_iet += gran_len
    return rv


satellite_base_times = {
    'snpp': 1698019234000000,
    'j01': 1698019234000000,
}
platform_short_names = {
    'snpp': 'NPP',
    'j01': 'J01',
}
instrument_short_names = {'viirs': 'VIIRS', 'cris': 'CrIS', 'atms': 'ATMS', None: 'SPACECRAFT'}