diff --git a/edosl0util/__init__.py b/edosl0util/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..3d096a1e02125cd6a672db736d777dd9eeccd041 100644
--- a/edosl0util/__init__.py
+++ b/edosl0util/__init__.py
@@ -0,0 +1,3 @@
+from pkg_resources import get_distribution
+
+__version__ = get_distribution(__name__).version
\ No newline at end of file
diff --git a/edosl0util/jpssrdr.py b/edosl0util/jpssrdr.py
index beef1e2f2bd22c6106d990423725a4946cf52a92..4062c741a1b9048fa709acbb6bb8fac8b8148b36 100644
--- a/edosl0util/jpssrdr.py
+++ b/edosl0util/jpssrdr.py
@@ -192,7 +192,7 @@ def _rdrs_for_packet_dataset(group):
if group:
for name, buf in _generate_packet_datasets(group):
try:
- rdr = _decode_rdr_blob(buf)
+ rdr = decode_rdr_blob(buf)
except ValueError as e:
LOG.warn('{} ({})'.format(e.message, name))
continue
@@ -200,7 +200,7 @@ def _rdrs_for_packet_dataset(group):
return rdrs
-def _decode_rdr_blob(buf):
+def decode_rdr_blob(buf):
if buf.shape[0] < c.sizeof(StaticHeader):
raise ValueError('Not enough bytes for static header')
header = StaticHeader.from_buffer(buf)
diff --git a/edosl0util/rdrgen.py b/edosl0util/rdrgen.py
index 73ac856be35ec8018f59ee9ac574b85535337a17..420699039fe6834319bc42564de93f49d57a743e 100644
--- a/edosl0util/rdrgen.py
+++ b/edosl0util/rdrgen.py
@@ -1,27 +1,114 @@
import ctypes
+import itertools
from collections import OrderedDict
+from datetime import datetime
import attr
-from astropy.time import Time, TimeDelta
+import h5py
import numpy as np
+from astropy.time import Time, TimeDelta
-from edosl0util.jpssrdr import StaticHeader, Apid as ApidListItem, PacketTracker
+import edosl0util
+from edosl0util.jpssrdr import (
+ StaticHeader, Apid as ApidListItem, PacketTracker, decode_rdr_blob)
-def prototype():
- pds_file = 'P1571289CRISSCIENCEA6T17230135400001.PDS'
- satellite = 'snpp'
- rdr_type = 'CRIS-SCIENCE-RDR'
- granule_iet = 1881755831079000
+def write_rdr(sat, blob,
+ distributor=None, origin=None, domain=None, orbit_number=0,
+ creation_time=None, gran_creation_time=None, software_ver=None):
+ # TODO: write out the user block too??
+ distributor = distributor or origin or default_origin
+ origin = origin or default_origin
+ domain = domain or default_domain
+ creation_time = creation_time or datetime.now()
+ gran_creation_time = gran_creation_time or creation_time
+ software_ver = (software_ver
+ or edosl0util.__name__ + '-' + edosl0util.__version__)
+ blob_info = decode_rdr_blob(blob)
+ rdr_type = SpacecraftDiaryRdrType() # TODO: select via APID
+ gran_iet = blob_info.header.start_boundary
+ gran_end_iet = blob_info.header.end_boundary
+ gran_time = iet_to_datetime(gran_iet)
+ gran_end_time = iet_to_datetime(gran_end_iet)
+ gran_id = make_granule_id(sat, gran_iet)
+ gran_ver = 'A1'
+ format_date = lambda t: t.strftime('%Y%m%d')
+ format_time = lambda t: t.strftime('%H%M%S.%fZ')
+ with h5py.File('rdr.h5', 'w') as h5_file:
+ set_h5_attrs(h5_file, {
+ 'Distributor': distributor,
+ 'Mission_Name': 'S-NPP/JPSS',
+ 'N_Dataset_Source': origin,
+ 'N_HDF_Creation_Date': format_date(creation_time),
+ 'N_HDF_Creation_Time': format_time(creation_time),
+ 'Platform_Short_Name': platform_short_names[sat]})
+ all_grp = h5_file.create_group('All_Data')
+ raw_grp = all_grp.create_group(rdr_type.short_name + '_All')
+ raw_ds = raw_grp.create_dataset(
+ 'RawApplicationPackets_0', data=blob, maxshape=[None])
+ prod_grp = h5_file.create_group('Data_Products')
+ gran_grp = prod_grp.create_group(rdr_type.short_name)
+ 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': domain})
+ gran_ds = gran_grp.create_dataset(
+ rdr_type.short_name + '_Gran_0', [1], h5_regionref_dtype, maxshape=[None])
+ gran_ds[0] = raw_ds.regionref[:]
+ set_h5_attrs(gran_ds, {
+ 'Beginning_Date': format_date(gran_time),
+ 'Beginning_Time': format_time(gran_time),
+ 'Ending_Date': format_date(gran_end_time),
+ 'Ending_Time': format_time(gran_end_time),
+ 'N_Beginning_Orbit_Number': np.uint64(orbit_number), # FIXME: detect?
+ 'N_Beginning_Time_IET': np.uint64(gran_iet),
+ 'N_Creation_Date': format_date(gran_creation_time),
+ 'N_Creation_Time': format_time(gran_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': 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(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': software_ver})
+ aggr_ds = gran_grp.create_dataset(
+ rdr_type.short_name + '_Aggr', [1], h5_ref_dtype, maxshape=[None])
+ aggr_ds[0] = raw_grp.ref
+ set_h5_attrs(aggr_ds, {
+ 'AggregateBeginningDate': format_date(gran_time),
+ 'AggregateBeginningGranuleID': gran_id,
+ 'AggregateBeginningOrbitNumber': np.uint64(orbit_number),
+ 'AggregateBeginningTime': format_time(gran_time),
+ 'AggregateEndingDate': format_date(gran_end_time),
+ 'AggregateEndingGranuleID': gran_id,
+ 'AggregateEndingOrbitNumber': np.uint64(orbit_number),
+ 'AggregateEndingTime': format_time(gran_end_time),
+ 'AggregateNumberGranules': np.uint64(1)})
-class RdrBlobBuilder(object):
- def __init__(self, rdr_type, sat, granule_iet):
- self.sat = sat
+def set_h5_attrs(h5_obj, 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
-def build_rdr_blob(pkt_stream, rdr_type, sat, granule_iet):
- granule_iet_end = granule_iet + rdr_type.granularity
+def build_rdr_blob(sat, pkt_stream):
+ first_pkt = next(pkt_stream) # FIXME: what if there are no packets?
+ rdr_type = SpacecraftDiaryRdrType() # TODO: select based on first packet APID
+ granule_iet = calc_iet_granule(satellite_base_times[sat], rdr_type.gran_len,
+ timecode_to_iet(first_pkt.secondary_header.timecode))
+ granule_iet_end = granule_iet + rdr_type.gran_len
total_pkt_size = 0
apid_info = OrderedDict()
@@ -35,12 +122,13 @@ def build_rdr_blob(pkt_stream, rdr_type, sat, granule_iet):
'pkt_info': [{} for _ in range(apid.max_expected)]}
total_trackers += apid.max_expected
- for pkt in pkt_stream:
+ for pkt in itertools.chain([first_pkt], pkt_stream):
if pkt.apid not in apid_info:
- continue
+ raise ValueError(
+ 'APID {} not expected for {}'.format(pkt.apid, rdr_type.short_name))
pkt_iet = timecode_to_iet(pkt.secondary_header.timecode)
if not granule_iet <= pkt_iet < granule_iet_end:
- continue
+ 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
@@ -58,8 +146,8 @@ def build_rdr_blob(pkt_stream, rdr_type, sat, granule_iet):
buf = np.zeros([buf_size], np.uint8) # zeros needed to null-pad strings
header = StaticHeader.from_buffer(buf)
- header.satellite = sat # XXX: mapping?
- header.sensor = rdr_type.sensor
+ header.satellite = platform_short_names[sat]
+ header.sensor = instrument_short_names[rdr_type.sensor]
header.type_id = rdr_type.type_id
header.num_apids = len(apid_info)
header.apid_list_offset = apid_list_offset
@@ -140,14 +228,67 @@ class ApidSpec(object):
class CrisScienceRdrType(object):
short_name = 'CRIS-SCIENCE-RDR'
- granularity = 31997000
- apids = get_cris_science_apids()
- sensor = 'CrIS'
+ gran_len = 31997000
+ sensor = 'cris'
type_id = 'SCIENCE'
+ document = '474-00448-02-03_JPSS-DD-Vol-II-Part-3_0200B.pdf'
+ apids = get_cris_science_apids()
+
+
+class SpacecraftDiaryRdrType(object):
+ short_name = 'SPACECRAFT-DIARY-RDR'
+ gran_len = 20000000
+ sensor = None
+ type_id = 'DIARY'
+ document = '474-00448-02-08_JPSS-DD-Vol-II-Part-8_0200B.pdf'
+ apids = [ApidSpec(0, 'CRITICAL', max_expected=21),
+ ApidSpec(8, 'ADCS_HKH', max_expected=21),
+ ApidSpec(11, 'DIARY', max_expected=21)]
def timecode_to_iet(timecode):
# FIXME: move to timecode.py
- ccsds_epoch = iet_epoch = Time('1958-01-01', scale='tai')
+ 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)
+
+
+def iet_to_datetime(iet):
+ return (iet_epoch + TimeDelta(iet * 1e-6, format='sec')).utc.datetime
+
+
+def calc_iet_granule(base_time, gran_len, iet):
+ return (iet - base_time) // gran_len * gran_len + base_time
+
+
+def calc_iet_aggregate(base_time, gran_len, grans_per_aggr, iet):
+ # see "DDS Aggregation Methodology" in CDFCB vol I
+ Y = gran_len
+ N = grans_per_aggr
+ G_s = calc_iet_granule(base_time, Y, iet)
+ A_n = G_s // (N * Y)
+ O = G_s % Y
+ A_s = A_n * (Y * N) + O
+ return A_s
+
+
+def make_granule_id(sat, gran_iet):
+ tenths = (gran_iet - satellite_base_times[sat]) // 100000
+ return '{}{:012d}'.format(platform_short_names[sat], tenths)
+
+
+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
+
+
+iet_epoch = Time('1958-01-01', scale='tai')
+satellite_base_times = {'snpp': 1698019234000000}
+platform_short_names = {'snpp': 'NPP'}
+instrument_short_names = {'cris': 'CrIS', None: 'SPACECRAFT'}
+default_origin = 'ssec'
+default_domain = 'dev'
+
+h5_ref_dtype = h5py.special_dtype(ref=h5py.Reference)
+h5_regionref_dtype = h5py.special_dtype(ref=h5py.RegionReference)