From c596819dffe33609254c3c55f99f01f509845196 Mon Sep 17 00:00:00 2001
From: Marko Mecina <marko.mecina@univie.ac.at>
Date: Tue, 10 Sep 2024 13:38:40 +0200
Subject: [PATCH] update L0 processing
---
.../dataprocessing/smile_L0b_converter.py | 17 ++-
.../dataprocessing/smile_raw_ce_converter.py | 124 ++++++++++++++----
2 files changed, 112 insertions(+), 29 deletions(-)
diff --git a/Ccs/tools/dataprocessing/smile_L0b_converter.py b/Ccs/tools/dataprocessing/smile_L0b_converter.py
index bae364b..6bd4789 100755
--- a/Ccs/tools/dataprocessing/smile_L0b_converter.py
+++ b/Ccs/tools/dataprocessing/smile_L0b_converter.py
@@ -56,6 +56,8 @@ PRODUCT_IDS = {0: 'SXI-SCI-ED',
SCI_PRODUCTS = {0: 'ED', 1: 'UNKNOWN', 2: 'FT', 3: 'UV', 4: 'FF'}
+UNKNOWN_PROD = 'UNKNOWN'
+
MODES = tuple(PRODUCT_IDS.values())
FT_NODES = ('FT_CCD_NODE_0', 'FT_CCD_NODE_1', 'FT_CCD_NODE_2', 'FT_CCD_NODE_3')
@@ -419,7 +421,7 @@ def get_ce_id(pkt):
obsid = int.from_bytes(pkt[SDU_DATA_OFF + 12: SDU_DATA_OFF + 16], 'big')
cecnt = int.from_bytes(pkt[SDU_DATA_OFF + 16: SDU_DATA_OFF + 18], 'big')
- product = PRODUCT_IDS[pkt[SDU_DATA_OFF + 28]]
+ product = PRODUCT_IDS.get(pkt[SDU_DATA_OFF + 28], UNKNOWN_PROD)
# product = pkt[SDU_DATA_OFF + 28]
return '{:010d}_{:05d}_{:09d}{:06d}_{:05d}_{}'.format(obsid, cecnt, coarse, fine, pktseqcnt, product)
@@ -535,7 +537,7 @@ def parse_pkts(fd):
elif pkt[SST_OFF] == 4:
if tx:
logging.warning('aborted downlink at {}'.format(get_pkt_time(pkt)))
- txpkts.append(extract_ce_data(pkt))
+ txpkts.append(extract_ce_data(pkt, check_seq=False))
bad_ces[ce_id] = b''.join(txpkts)
tx = False
else:
@@ -575,6 +577,7 @@ def extract(infile, outdir):
logging.warning('skipped {} bytes because of wrong CRCs'.format(trashcnt))
except Exception as err:
logging.exception(err)
+ good_ces, bad_ces, hks = [], [], []
logging.info('extracted {} files'.format(len(good_ces)))
@@ -1075,11 +1078,17 @@ def process_file(infile, outdir):
decompressed = {mode: [] for mode in MODES}
for ce in ces:
try:
+
+ if ce.endswith('_{}.ce'.format(UNKNOWN_PROD)):
+ logging.exception('Skipped decompression of {}'.format(ce))
+ continue
+
fitspath = decompress(ce, outdir)
if os.path.isfile(fitspath):
decompressed = sort_by_mode(decompressed, fitspath)
+
except Exception as err:
- # logging.error('Decompression failed for {}'.format(ce))
+ logging.error('Decompression failed for {}'.format(ce))
logging.exception(err)
# merged = merge_fits(decompressed, infile)
@@ -1122,7 +1131,7 @@ def setup_logging(output_dir):
if __name__ == '__main__':
# setup_logging('/home/marko/space/smile/cedata/proc')
- # process_file('/home/marko/space/smile/cedata/datapools/UL_flatsat_08072024_1156_rev_clk_dgen.bin', '/home/marko/space/smile/cedata/proc')
+ # process_file('/home/marko/ucloud/madrid/2024/X_BAND_SCOE/SXI_TMTC_RRC1_PS_X_20240729_merged.bin', '/home/marko/ucloud/madrid/2024/X_BAND_SCOE')
# sys.exit()
infile = sys.argv[1]
diff --git a/Ccs/tools/dataprocessing/smile_raw_ce_converter.py b/Ccs/tools/dataprocessing/smile_raw_ce_converter.py
index ed957aa..347358e 100755
--- a/Ccs/tools/dataprocessing/smile_raw_ce_converter.py
+++ b/Ccs/tools/dataprocessing/smile_raw_ce_converter.py
@@ -24,8 +24,9 @@ SIZE_FF = NROWS_FF * NCOLS_FF * 2
SIZE_FT = NROWS_FT * NCOLS_FT * 2 # 1 node
SIZE_UV = NROWS_UV * NCOLS_UV * 2 # 1 node
SIZE_ED = 64 # 1 event
+NODE_NUM = 4 # total number of nodes
-SCI_PRODUCTS = {0: 'ED', 1: 'UNKNOWN', 2: 'FT', 3: 'UV', 4: 'FF'}
+SCI_PRODUCTS = {0: 'ED', 1: 'UNKNOWN', 2: 'FT', 3: 'UV', 4: 'FF', 5: 'ED_WP'}
FILE_PREFIX = 'SMILE_SXI_L1'
@@ -49,11 +50,18 @@ def convert_ce(cefile, fitsfile=None, guess=False):
else:
print('Cannot determine product type for CE of length {}, aborting.'.format(len(cedata)))
sys.exit()
+
+ obsid = 0
+
else:
try:
ce = CompressionEntity(cedata)
prod = ce.header.items.product
+ obsid = ce.header.items.obsid
+
+ assert ce.meta_frame.shape[0] == ce.header.items.sdp_group_members
+
if SCI_PRODUCTS.get(prod) == 'FF':
mode, hdl = mk_ff(ce)
elif SCI_PRODUCTS.get(prod) == 'FT':
@@ -72,7 +80,7 @@ def convert_ce(cefile, fitsfile=None, guess=False):
if fitsfile is None:
outdir = os.path.dirname(os.path.abspath(cefile))
- fitsfile = os.path.join(outdir, '{}_{}_{}.fits'.format(FILE_PREFIX, mode, _mk_ts()))
+ fitsfile = os.path.join(outdir, '{}_{}_{:010d}_{}.fits'.format(FILE_PREFIX, mode, obsid, _mk_ts()))
hdl.writeto(fitsfile, overwrite=True)
@@ -106,17 +114,27 @@ def mk_ft(data):
arr = np.frombuffer(data.scidata, dtype='>H').reshape(-1, NROWS_FT, NCOLS_FT)
hdl = _mk_hdl('FT', data.header)
- for n in range(arr.shape[0]):
- node = fits.ImageHDU(data=arr[n, :, :], name='FT_CCD_NODE_{}'.format(n))
- node.add_checksum()
- hdl.append(node)
+ n_nodes = arr.shape[0]
+
+ if n_nodes % NODE_NUM:
+ print('Total number of nodes ({}) is not a multiple of 4!'.format(n_nodes))
+
+ for n in range(arr.shape[0]):
+ node = fits.ImageHDU(data=arr[n, :, :], name='FT_CCD_NODE_{}'.format(n))
+ node.add_checksum()
+ hdl.append(node)
+ else:
+ for n in range(NODE_NUM):
+ node = fits.ImageHDU(data=arr[n::NODE_NUM, :, :], name='FT_CCD_NODE_{}'.format(n))
+ node.add_checksum()
+ hdl.append(node)
# arrange all nodes to full CCD
- if arr.shape[0] == 4:
+ if n_nodes % 4 == 0:
nn = _assemble_ft_frames_to_fp_view(arr)
- hdl.append(fits.ImageHDU(data=nn, name='FULLCCD'))
+ hdl.append(fits.ImageHDU(data=nn, name='FULLARRAY'))
group_table = fits.BinTableHDU(data=data.meta_group, name='GROUP_HK')
frame_table = fits.BinTableHDU(data=data.meta_frame, name='FRAME_HK')
@@ -134,17 +152,25 @@ def mk_uv(data):
arr = np.frombuffer(data.scidata, dtype='>H').reshape(-1, NROWS_UV, NCOLS_UV)
hdl = _mk_hdl('UV', data.header)
- for n in range(arr.shape[0]):
- node = fits.ImageHDU(data=arr[n, :, :], name='UV_CCD_NODE_{}'.format(n))
- node.add_checksum()
- hdl.append(node)
+ n_nodes = arr.shape[0]
+
+ if n_nodes % NODE_NUM:
+ for n in range(arr.shape[0]):
+ node = fits.ImageHDU(data=arr[n, :, :], name='UV_CCD_NODE_{}'.format(n))
+ node.add_checksum()
+ hdl.append(node)
+ else:
+ for n in range(NODE_NUM):
+ node = fits.ImageHDU(data=arr[n::NODE_NUM, :, :], name='UV_CCD_NODE_{}'.format(n))
+ node.add_checksum()
+ hdl.append(node)
# arrange all nodes to full CCD
- if arr.shape[0] == 4:
+ if n_nodes % 4 == 0:
nn = _assemble_ft_frames_to_fp_view(arr)
- hdl.append(fits.ImageHDU(data=nn, name='FULLCCD'))
+ hdl.append(fits.ImageHDU(data=nn, name='FULLARRAY'))
group_table = fits.BinTableHDU(data=data.meta_group, name='GROUP_HK')
frame_table = fits.BinTableHDU(data=data.meta_frame, name='FRAME_HK')
@@ -162,10 +188,34 @@ def mk_ed(data):
# reshape into array of evt packets
arr = np.frombuffer(data.scidata, dtype='>H').reshape(-1, SIZE_ED // 2)
+ frmevtcnt = data.meta_frame['LastFrameEvtCount']
+
+ # handle wandering pixel mask if present
+ wpm = []
+ if frmevtcnt.sum() != arr.shape[0]:
+ print('Inconsistent LastFrameEvtCount {} (sum). {} in data block.'.format(frmevtcnt.sum(), arr.shape[0]))
+ if data.header.items.sdp_group_members == (arr.shape[0] - frmevtcnt.sum()):
+ print('Extracting {} wandering pixel(s).'.format(data.header.items.sdp_group_members))
+ wpm = arr[-data.header.items.sdp_group_members:]
+ arr = arr[:-data.header.items.sdp_group_members]
+
+ if data.meta_frame['LastFrameEvtCount'].sum() != data.meta_group['NOfEvtDet'][0]:
+ print('Inconsistent CE frame count: LastFrameEvtCount_sum = {}, NOfEvtDet = {}'.format(frmevtcnt.sum(), data.meta_group['NOfEvtDet'][0]))
+
hdl = _mk_hdl('ED', data.header)
# ts = int(hdl['PRIMARY'].header['OBS_ID'])
ts = data.meta_frame['sdpProductStarttimeCrs'] + data.meta_frame['sdpProductStarttimeFine'] / 1e6
- bindata = np.array([_mk_bin_entry(evt, ts) for evt in arr], dtype=ED_BIN_DTYPE)
+ ts_ext = ts.repeat(frmevtcnt)
+ bindata = np.array([_mk_bin_entry(evt, t, wpm_raise=True) for evt, t in zip(arr, ts_ext)], dtype=ED_BIN_DTYPE)
+
+ if len(wpm) != 0:
+ try:
+ wpmdata = np.array([_mk_bin_entry(evt, t) for evt, t in zip(arr, ts)], dtype=ED_BIN_DTYPE)
+ except Exception as err:
+ wpmdata = None
+ print('Failed extracting wandering pixel masks ({})'.format(err))
+ else:
+ wpmdata = None
# also add an HDU with event map
nodes = np.zeros((2, 2, NROWS_FT, NCOLS_FT))
@@ -176,7 +226,7 @@ def mk_ed(data):
except:
print(col, row, 'FAILED')
- nodes[nodes == 0] = np.nan
+ # nodes[nodes == 0] = np.nan
ed_img = _assemble_ft_frames_to_fp_view(nodes)
@@ -189,6 +239,11 @@ def mk_ed(data):
group_table = fits.BinTableHDU(data=data.meta_group, name='GROUP_HK')
frame_table = fits.BinTableHDU(data=data.meta_frame, name='FRAME_HK')
+ if wpmdata is not None:
+ wpm_table = fits.BinTableHDU(data=wpmdata, name='WNDRNG_PXL')
+ wpm_table.add_checksum()
+ hdl.append(wpm_table)
+
# checksums
group_table.add_checksum()
frame_table.add_checksum()
@@ -198,12 +253,16 @@ def mk_ed(data):
return 'ED', hdl
-def _mk_bin_entry(data, timestamp):
+def _mk_bin_entry(data, timestamp, wpm_raise=False):
spw, dlen, dtyp, fc, sc, col, row, *evts = data
ccdnr = (dtyp >> 4) & 1
node = (dtyp >> 5) & 0b11
+ if wpm_raise:
+ if (dtyp & 0b11) == 3:
+ raise ValueError('WPM detected')
+
return timestamp, fc, ccdnr, col, row, node, evts
@@ -214,23 +273,31 @@ def _assemble_ft_frames_to_fp_view(arrnd):
# FT
# according to MSSL-IF-122 the nodes arrive in the following order: E2, F2, E4, F4
if arrnd.ndim == 3:
- n00 = arrnd[1, ::-1, ::-1] # CCD2 F-side (upper right in FP view)
- n01 = arrnd[0, ::-1, :] # CCD2 E-side (upper left in FP view)
- n10 = arrnd[3, :, :] # CCD4 F-side (lower left in FP view)
- n11 = arrnd[2, :, ::-1] # CCD4 E-side (lower right in FP view)
+ n00 = arrnd[1::NODE_NUM, ::-1, ::-1] # CCD2 F-side (upper right in FP view)
+ n01 = arrnd[0::NODE_NUM, ::-1, :] # CCD2 E-side (upper left in FP view)
+ n10 = arrnd[3::NODE_NUM, :, :] # CCD4 F-side (lower left in FP view)
+ n11 = arrnd[2::NODE_NUM, :, ::-1] # CCD4 E-side (lower right in FP view)
+
# ED
elif arrnd.ndim == 4:
n00 = arrnd[0, 0, :, :][::-1, ::-1] # CCD2 F-side
n01 = arrnd[0, 1, :, :][::-1, :] # CCD2 E-side
n10 = arrnd[1, 0, :, :] # CCD4 F-side
n11 = arrnd[1, 1, :, :][:, ::-1] # CCD4 E-side
+
else:
return
- n0 = np.concatenate((n01, n00), axis=1) # CCD2
- n1 = np.concatenate((n10, n11), axis=1) # CCD4
+ # in case multiple frames are in the FT/UV CE
+ if n00.ndim == 3:
+ ax = 1
+ else:
+ ax = 0
- return np.concatenate((n1, n0), axis=0)
+ n0 = np.concatenate((n01, n00), axis=1 + ax) # CCD2
+ n1 = np.concatenate((n10, n11), axis=1 + ax) # CCD4
+
+ return np.concatenate((n1, n0), axis=0 + ax)
def _mk_hdl(dmode, dhead):
@@ -241,6 +308,8 @@ def _mk_hdl(dmode, dhead):
phdu.header['INSTRUME'] = 'SXI'
phdu.header['DATAMODE'] = dmode
phdu.header['OBS_ID'] = dhead.items.obsid
+ phdu.header['PRODTYP'] = 'single CE'
+ phdu.header['SDPGRPN'] = (dhead.items.sdp_group_members, 'number of SDP group members')
phdu.header['SOFTVER'] = dhead.items.version_number
phdu.header['BUILD'] = dhead.items.build_number
@@ -444,7 +513,7 @@ class CeHeader(ctypes.Union):
def __str__(self):
return '\n'.join(['{}: {}'.format(n, getattr(self.items, n)) for n, _ in self.items._fields_])
- def show(self):
+ def info(self):
print(self.__str__())
@@ -481,6 +550,11 @@ class CompressionEntity:
if __name__ == "__main__":
+ # cefile="/home/marko/space/smile/betatest/eqm_flatsat/newsci/LDT_207_0207139797_026819.ce"
+ # convert_ce(cefile, fitsfile=None)
+ #
+ # sys.exit()
+
if len(sys.argv) > 2:
cefile, fitsfile = sys.argv[1:3]
else:
--
GitLab