1#!/usr/bin/env python3
2
3"""Utilities for handling CCSDS packets and datatypes found in CCSDS data."""
4
5
6
7import mmap
8import logging
9from datetime import datetime, timedelta
10
11import numpy as np
12
13# an info structure would be better, listing APID, packet name, instrument, size
14# and frequency
15APID_TM = 1
16APID_REDUCED_NAV = 2
17APID_ADMIN = 6
18APID_AMSUA1 = 39
19APID_AMSUA2 = 40
20APID_AVHRR_HR_3A = 103
21APID_AVHRR_HR_3B = 104
22APID_GRAS = 480
23APID_SEM = 37
24APID_AVHRR_3A = 103
25APID_AVHRR_3B = 104
26
27ccsds_time_desc = np.dtype([('days', '>u2'), ('ms', '>u4'), ('us', '>u2')])
28ccsds_end_epoch = np.zeros((1,), dtype=ccsds_time_desc)
29ccsds_end_epoch['days'] = 65535
30ccsds_start_epoch = np.zeros((1,), dtype=ccsds_time_desc)
31
32np_uint16 = np.dtype('>u2') # big endian
33np_uint32 = np.dtype('>u4') # big endian
34
35DAYS_PER_YEAR = 365
36UTC_TO_JULIAN_YEARS = 2000 - 1970
37LEAP_YEARS_1970_2000 = 7
38SECONDS_PER_DAY = 86400
39UTC_TO_JULIAN_SECONDS = SECONDS_PER_DAY * (DAYS_PER_YEAR * UTC_TO_JULIAN_YEARS +
40 LEAP_YEARS_1970_2000)
41CCSDS_REFTIME = datetime(2000, 1, 1)
42
43
44def extract_grh_sensing_start(grh):
45 """Convert a numpy GRH structure to a datetime.
46 """
47
48 return CCSDS_REFTIME + timedelta(days=int(grh['RECORD_START_TIME']['DAYS']),
49 milliseconds=int(grh['RECORD_START_TIME']['MILLISECONDS']))
50
51
52def ccsds_to_datetime(ccsds):
53 """Extracting onboard UTC time from CCSDS packet
54 """
55
56 #day = ccsds['day']
57 #ms = ccsds['ms']
58 #us = ccsds['us']
59 day = ccsds[0]
60 ms = ccsds[1]
61 us = ccsds[2]
62 offt = timedelta(days=int(day), milliseconds=int(ms), microseconds=int(us))
63 return CCSDS_REFTIME + offt
64
65
66def gentime_to_ccsds(gentime):
67 """Convert an EPS gentime ('20100912120000') to a datetime.
68 """
69
70 year = int(gentime[0:4])
71 month = int(gentime[4:6])
72 day = int(gentime[6:8])
73 hour = int(gentime[8:10])
74 minute = int(gentime[10:12])
75 second = int(gentime[12:14])
76 dt = datetime(year, month, day, hour, minute, second)
77 diff = dt - CCSDS_REFTIME
78 return (diff.days, diff.seconds * 1000, 0)
79
80
81# def datetime_to_ccsds(dt):
82# """Convert a datetime to a CCSDS time triple of (days since 2001-01-01, seconds of day,
83# milliseconds).
84# Obsolete function - does not appear to be used?
85# """
86
87# diff = dt - CCSDS_REFTIME
88# #print 'diff ', diff.seconds, diff.microseconds
89# return (diff.days, diff.seconds * 1000 + diff.microseconds / 1000, diff.microseconds % 1000)
90
91def datetime_to_ccsds_time(dt):
92 """Convert a datetime into a binary CCSDS packet UTC time (days, ms, us)
93 accurate to microseconds and binary compatible with METOP/MSG timestamps
94 """
95 ccsds_time = np.zeros((1,), dtype=ccsds_time_desc)
96 diff = dt - CCSDS_REFTIME
97 ccsds_time['days'] = diff.days
98 ccsds_time['ms'] = diff.seconds * 1000 + diff.microseconds // 1000
99 ccsds_time['us'] = diff.microseconds % 1000
100 return ccsds_time
101
102
103def get_apid(ccsds):
104 """Extract APID from a CCSDS packet.
105 """
106
107 #print mdr_payload.size
108 pi = ccsds[0:2]
109 return (pi[0] & 7) << 8 | pi[1]
110
111
112def gen_gras_subpackets(timestamp, ccsds, subtypes):
113 """Generate a list of GRAS subpackets from `ccsds`, filtering for subpacket
114 IDs in `subtypes`.
115 """
116
117 #ccsds_obt = ccsds[14:20]
118 ccsds_obt = ((int(ccsds[14]) << 32) +
119 (int(ccsds[15]) << 24) +
120 (int(ccsds[16]) << 16) +
121 (int(ccsds[17]) << 8) +
122 int(ccsds[18]))
123 # logging.debug('obt '+str(ccsds_obt) + ' subtypes ' + str(subtypes))
124 i = 24
125 while i < ccsds.size - 2:
126 len_typ = ccsds[i:i + 2]
127 len_typ.dtype = np_uint16
128 sub_type = (len_typ[0]) & 0xf
129 sub_len = (len_typ[0]) >> 4
130
131 if sub_len == 0:
132 return
133
134 # logging.debug('subpacket type '+str(sub_type)+' len '+str(sub_len))
135 if sub_type in subtypes:
136 if sub_type == 3: # temperature and voltage
137 obt = ccsds[i + 2:i + 6]
138 obt.dtype = np_uint32
139 delta = ccsds_obt - (obt[0] * 256)
140 delta_td = timedelta(microseconds=delta * 1000000.0 / 65536.0)
141 #logging.debug('subpacket obt '+str(obt[0]*256) + ' delta '+\
142 #str(delta)+' delta td '+str(delta_td))
143 yield timestamp - delta_td, ccsds[i:i + sub_len]
144 else:
145 raise Exception('Subpacket ' + str(sub_type) + ' not handled')
146
147 # logging.debug('subpacket done')
148 i += sub_len
149
150 # logging.debug('packet done')
151
152
153class CCSDSReader:
154 """Read a stream of CCSDS packets.
155 """
156
157 # these must match the names in distiller.distillers.
158 instruments = {
159 'AMSA': (39, 40),
160 'AVHR': (103, 104, 66, 67, 68, 69),
161 'GOME': (384, ),
162 'GRAS': (448, 480),
163 'HIRS': (38, ),
164 'HKTM': (1, 2, 3, 6),
165 'MHSx': (34, ),
166 'SEMx': (37, ),
167 }
168
169 def __init__(self, filename, multiple_instrument_check='throw'):
170 """
171 `packets` should be the filename to be read.
172 `multiple_instrument_check` can be set to:
173 - 'throw' :: throw an exception if data from multiple instruments is found/
174 - 'ignore' :: the first packet of the file will lock the reader to that instrument.
175 Data from other instruments will be discarded.
176 - 'disable' :: allow data from multiple instruments to be read.
177 """
178
179 self.filename = filename
180 self.multiple_instrument_check = multiple_instrument_check
181
182 fd = filename.open('rb')
183 self.buf = mmap.mmap(fd.fileno(), 0, access=mmap.ACCESS_READ)
184
185 self.instrument = None
186 self.iter = self._gen_ccsds()
187 self.hold = next(self.iter)
188
189 def get_scid(self):
190 """We define this function because some distillers will call it.
191 In practise they don't actually need a correct SCID though."""
192
193 return None
194
195 scid = property(get_scid)
196
197 def get_sensing_start(self):
198 """Sensing time of first snack."""
199 return self.hold[0]
200
201 sensing_start = property(get_sensing_start)
202
203 def gen_ccsds(self):
204 """Yield all snacks."""
205 yield self.hold
206 for x in self.iter:
207 yield x
208
209 def _gen_ccsds(self):
210 """Yield a stream of tuples of (timestamp, payload),
211 where `timestamp` is the UTC time of the packet.
212 """
213
214 np_uint16_ne = np.dtype('>u2')
215
216 pos = 0
217 while pos < len(self.buf):
218 utc_time = ccsds_to_datetime(
219 np.frombuffer(
220 buffer=self.buf,
221 dtype=ccsds_time_desc,
222 count=1,
223 offset=pos + 6)[0])
224
225 payload = None
226
227 length = np.frombuffer(
228 buffer=self.buf,
229 dtype=np_uint16_ne,
230 count=1,
231 offset=pos + 4)[0]
232
233 payload = np.frombuffer(buffer=self.buf, dtype=np.uint8, count=length + 7, offset=pos)
234
235 apid = get_apid(payload)
236
237 found = False
238 for k, v in CCSDSReader.instruments.items():
239 if apid in v:
240 if self.instrument is None:
241 self.instrument = k
242 logging.info('Detected instrument ' + k)
243
244 elif self.instrument != k:
245 logging.error('Found packet of instrument ' + k + ' when previous data '
246 'had instrument ' + self.instrument)
247
248 found = True
249
250 if self.instrument is not None and not found:
251 logging.error('Found data of unknown instrument (APID {apid}) when previous data '
252 'has instrument {inst}'.format(apid=apid, inst=self.instrument))
253
254 yield utc_time, payload
255
256 pos += length + 7
257
258
259def main():
260 """Command line entry point."""
261
262 from chart.common.args import ArgumentParser
263
264 parser = ArgumentParser()
265 parser.add_argument('files', nargs='*')
266 args = parser.parse_args()
267
268 for f in args.files:
269 prod = CCSDSReader(f)
270 print('Instrument ', prod.instrument)
271 print('SCID ', prod.scid)
272 for timestamp, payload in prod.gen_ccsds():
273 print(timestamp, get_apid(payload), len(payload))
274
275if __name__ == '__main__':
276 main()