From e7eee3b13972127940d2ce4a0fb6f5ffe2cf0d33 Mon Sep 17 00:00:00 2001
From: Marko Mecina <marko.mecina@univie.ac.at>
Date: Fri, 3 Nov 2023 17:54:50 +0100
Subject: [PATCH] implement pipeline to generate FITS files from S13 data
---
Ccs/ccs_function_lib.py | 19 +-
Ccs/decompression.py | 14 +-
Ccs/packet_config_SMILE.py | 1 +
.../dataprocessing/smile_raw_ce_converter.py | 162 ++++++++++++++++++
4 files changed, 189 insertions(+), 7 deletions(-)
create mode 100755 Ccs/tools/dataprocessing/smile_raw_ce_converter.py
diff --git a/Ccs/ccs_function_lib.py b/Ccs/ccs_function_lib.py
index 62f1d0f..43dd85a 100644
--- a/Ccs/ccs_function_lib.py
+++ b/Ccs/ccs_function_lib.py
@@ -5507,6 +5507,7 @@ def add_tst_import_paths():
sys.path.append(cfg.get('paths', 'tst') + '/sketch_desk')
sys.path.append(cfg.get('paths', 'tst') + '/test_specs')
sys.path.append(cfg.get('paths', 'tst') + '/testing_library')
+ sys.path.append(cfg.get('paths', 'tst') + '/testing_library/testlib')
# insert this to import the tst view.py, not the one in .local folder
sys.path.insert(0, cfg.get('paths', 'tst') + '/tst')
@@ -5688,10 +5689,15 @@ def collect_13(pool_name, starttime=None, endtime=None, startidx=None, endidx=No
continue
else:
pkts = [a.raw[S13_HEADER_LEN_TOTAL:-PEC_LEN] for a in tm_132.filter(DbTelemetry.idx > i.idx, DbTelemetry.idx < j.idx)]
+
+ # check for padding bytes in last packet
+ datalen = int.from_bytes(j.raw[S13_DATALEN_PAR_OFFSET:S13_DATALEN_PAR_OFFSET + S13_DATALEN_PAR_SIZE], 'big')
+ lastpktdata = j.raw[S13_HEADER_LEN_TOTAL:S13_HEADER_LEN_TOTAL + datalen]
+
if join:
- ces[float(i.timestamp[:-1])] = i.raw[S13_HEADER_LEN_TOTAL:-PEC_LEN] + b''.join(pkts) + j.raw[S13_HEADER_LEN_TOTAL:-PEC_LEN]
+ ces[float(i.timestamp[:-1])] = i.raw[S13_HEADER_LEN_TOTAL:-PEC_LEN] + b''.join(pkts) + lastpktdata
else:
- ces[float(i.timestamp[:-1])] = [i.raw[S13_HEADER_LEN_TOTAL:-PEC_LEN]] + [b''.join(pkts)] + [j.raw[S13_HEADER_LEN_TOTAL:-PEC_LEN]]
+ ces[float(i.timestamp[:-1])] = [i.raw[S13_HEADER_LEN_TOTAL:-PEC_LEN]] + pkts + [lastpktdata]
k += 2
cnt += 1
@@ -5738,6 +5744,11 @@ def dump_large_data(pool_name, starttime=0, endtime=None, outdir="", dump_all=Fa
try:
obsid, ctime, ftime, ctr = s13_unpack_data_header(ldt_dict[buf])
+ except NotImplementedError:
+ obsid = int(datetime.datetime.utcnow().strftime('%j'))
+ ctime = int(buf)
+ ftime = 0
+ ctr = crc(ldt_dict[buf])
except ValueError as err:
logger.error('Incompatible definition of S13 data header.')
raise err
@@ -6940,7 +6951,11 @@ try:
SDU_PAR_LENGTH = _s13_info[0][-1]
# length of PUS + source header in S13 packets (i.e. data to be removed when collecting S13)
S13_HEADER_LEN_TOTAL = TM_HEADER_LEN + sum([p[-1] for p in _s13_info])
+ S13_DATALEN_PAR_OFFSET, S13_DATALEN_PAR_SIZE = TM_HEADER_LEN + sum([x[1] for x in _s13_info[:-1]]), _s13_info[-1][1]
except (SQLOperationalError, NotImplementedError, IndexError):
logger.warning('Could not get S13 info from MIB, using default values')
SDU_PAR_LENGTH = 1
S13_HEADER_LEN_TOTAL = 21
+ S13_DATALEN_PAR_OFFSET = 19
+ S13_DATALEN_PAR_SIZE = 2
+
diff --git a/Ccs/decompression.py b/Ccs/decompression.py
index 7966865..c97bcce 100644
--- a/Ccs/decompression.py
+++ b/Ccs/decompression.py
@@ -76,9 +76,9 @@ def convert_fullframe_to_cheopssim(fname):
hdulist.writeto(fname[:-5] + '_CHEOPSSIM.fits')
-def ce_decompress(outdir, pool_name=None, sdu=None, starttime=None, endtime=None, startidx=None, endidx=None,
+def ce_decompress(pool_name, outdir, sdu=None, starttime=None, endtime=None, startidx=None, endidx=None,
ce_exec=None):
- decomp = CeDecompress(outdir, pool_name=pool_name, sdu=sdu, starttime=starttime, endtime=endtime, startidx=startidx,
+ decomp = CeDecompress(pool_name, outdir, sdu=sdu, starttime=starttime, endtime=endtime, startidx=startidx,
endidx=endidx, ce_exec=ce_exec)
decomp.start()
@@ -94,7 +94,7 @@ def ce_decompress_stop(name=None):
class CeDecompress:
- def __init__(self, outdir, pool_name=None, sdu=None, starttime=None, endtime=None, startidx=None, endidx=None,
+ def __init__(self, pool_name, outdir, sdu=None, starttime=None, endtime=None, startidx=None, endidx=None,
ce_exec=None):
self.outdir = outdir
self.pool_name = pool_name
@@ -156,7 +156,8 @@ class CeDecompress:
fitspath = cefile[:-2] + 'fits'
if os.path.isfile(fitspath):
subprocess.run(["rm", fitspath])
- subprocess.run([self.ce_exec, cefile, fitspath], stdout=open(cefile[:-2] + 'log', 'w'))
+ with open(cefile[:-2] + 'log', 'w') as logfd:
+ subprocess.run([self.ce_exec, cefile, fitspath], stdout=logfd, stderr=logfd)
# first, get all TM13s already complete in pool
try:
@@ -164,8 +165,11 @@ class CeDecompress:
outdir=self.outdir, dump_all=True, sdu=self.sdu, startidx=self.startidx,
endidx=self.endidx)
for ce in filedict:
- self.last_ce_time = ce
decompress(filedict[ce])
+ self.last_ce_time = ce
+
+ self.last_ce_time += self.ldt_minimum_ce_gap
+
except (ValueError, TypeError, AttributeError) as err:
ce_decompressors.pop(self.init_time)
raise err
diff --git a/Ccs/packet_config_SMILE.py b/Ccs/packet_config_SMILE.py
index ed9753f..802c4e4 100644
--- a/Ccs/packet_config_SMILE.py
+++ b/Ccs/packet_config_SMILE.py
@@ -793,4 +793,5 @@ _S13_HEADER_FMT = S13_FMT_OBSID + S13_FMT_TIME + S13_FMT_FTIME + S13_FMT_COUNTER
def s13_unpack_data_header(buf):
+ raise NotImplementedError
return struct.unpack('>' + _S13_HEADER_FMT, buf[:struct.calcsize(_S13_HEADER_FMT)])
diff --git a/Ccs/tools/dataprocessing/smile_raw_ce_converter.py b/Ccs/tools/dataprocessing/smile_raw_ce_converter.py
new file mode 100755
index 0000000..4930f74
--- /dev/null
+++ b/Ccs/tools/dataprocessing/smile_raw_ce_converter.py
@@ -0,0 +1,162 @@
+#!/usr/bin/env python
+
+"""
+Convert unprocessed CE raw data to FITS files. Product type is determined (guessed) based on CE size.
+"""
+
+import datetime
+import os
+import sys
+
+import numpy as np
+from astropy.io import fits
+
+# expected CE sizes in bytes
+SIZE_FF = 4511 * 4608 * 2
+SIZE_FT = 639 * 384 * 2 # 1 node
+SIZE_ED = 64 # 1 event
+
+FILE_PREFIX = 'SMILE_SXI_L1'
+
+ED_BIN_DTYPE = np.dtype(
+ [('TIME', '>f8'), ('CCDFRAME', '>i4'), ('CCDNR', 'u1'), ('RAWX', '>i2'), ('RAWY', '>i2'), ('AMP', 'u1'),
+ ('PHAS', '>i2', (25,))])
+
+
+def convert_ce(cefile, fitsfile=None):
+
+ cedata = open(cefile, 'rb').read()
+
+ # guess product based on CE size
+ if len(cedata) == SIZE_FF:
+ mode, hdl = mk_ff(cedata)
+ elif len(cedata) // SIZE_FT in [1, 2, 4]:
+ mode, hdl = mk_ft(cedata)
+ elif len(cedata) % SIZE_ED == 0:
+ mode, hdl = mk_ed(cedata)
+ else:
+ print('Cannot determine product type for CE of length {}, aborting.'.format(len(cedata)))
+ sys.exit()
+
+ if fitsfile is None:
+ outdir = os.path.dirname(os.path.abspath(cefile))
+ fitsfile = os.path.join(outdir, '{}_{}_{}.fits'.format(FILE_PREFIX, mode, _mk_ts()))
+ else:
+ fitsfile = fitsfile.replace('.fits', '_{}.fits'.format(mode))
+
+ hdl.writeto(fitsfile, overwrite=True)
+
+
+def mk_ff(data):
+ # create uint16 array from raw data and reshape
+ arr = np.frombuffer(data, dtype='>H').reshape(4511, -1)
+ fnode = arr[:, ::2]
+ enode = arr[:, 1::2][:, ::-1]
+ ff = np.concatenate((fnode, enode), axis=1)
+
+ # write array to FITS file
+ hdl = _mk_hdl('FF')
+ hdl.append(fits.ImageHDU(data=ff, name='FULLFRAME'))
+
+ return 'FF', hdl
+
+
+def mk_ft(data):
+ arr = np.frombuffer(data, dtype='>H').reshape(-1, 639, 384)
+
+ hdl = _mk_hdl('FT')
+ for n in range(arr.shape[0]):
+ hdl.append(fits.ImageHDU(data=arr[n, :, :], name='FT_CCD_NODE_{}'.format(n)))
+
+ # arrange all nodes to full CCD view
+ if arr.shape[0] == 4:
+ n00 = arr[0, :, :] # CCD2 F-side
+ n01 = arr[1, :, ::-1] # CCD2 E-side
+ n10 = arr[2, ::-1, :] # CCD4 F-side
+ n11 = arr[3, ::-1, ::-1] # CCD4 E-side
+
+ n0 = np.concatenate((n00, n01), 1)
+ n1 = np.concatenate((n10, n11), 1)
+
+ nn = np.concatenate((n0, n1), 0) # [::-1, :]
+
+ hdl.append(fits.ImageHDU(data=nn, name='FULLCCD'))
+
+ return 'FT', hdl
+
+
+def mk_ed(data):
+ # reshape into array of evt packets
+ arr = np.frombuffer(data, dtype='>H').reshape(-1, 32)
+
+ hdl = _mk_hdl('FT')
+ ts = int(hdl['PRIMARY'].header['OBS_ID'])
+ bindata = np.array([_mk_bin_entry(evt, ts) for evt in arr], dtype=ED_BIN_DTYPE)
+
+ # also add an HDU with event map
+ nodes = np.zeros((2, 2, 639, 384))
+ nodes[:] = np.nan
+
+ for _, _, ccd, col, row, node, fx in bindata:
+ try:
+ nodes[ccd, node, row - 2:row + 3, col - 2:col + 3] = fx.reshape(5, 5)
+ except:
+ print(col, row, 'FAILED')
+
+ n00 = nodes[0, 0, :, :] # CCD2 F-side
+ n01 = nodes[0, 1, :, :][:, ::-1] # CCD2 E-side
+ n10 = nodes[1, 0, :, :][::-1, :] # CCD4 F-side
+ n11 = nodes[1, 1, :, :][::-1, ::-1] # CCD4 E-side
+
+ n0 = np.concatenate((n00, n01), 1)
+ n1 = np.concatenate((n10, n11), 1)
+
+ ed_img = np.concatenate((n0, n1), 0) # [::-1, :]
+
+ hdl.append(fits.ImageHDU(data=ed_img, name='EVTMAP'))
+
+ hdl.append(fits.BinTableHDU(data=bindata, name='EVENTS'))
+
+ return 'ED', hdl
+
+
+def _mk_bin_entry(data, timestamp):
+ spw, dlen, dtyp, fc, sc, col, row, *evts = data
+
+ ccdnr = (dtyp >> 4) & 1
+ node = (dtyp >> 5) & 0b11
+
+ return timestamp, fc, ccdnr, col, row, node, evts
+
+
+def _mk_hdl(dmode):
+ hdl = fits.HDUList()
+ phdu = fits.PrimaryHDU()
+
+ phdu.header['TELESCOP'] = 'SMILE'
+ phdu.header['INSTRUME'] = 'SXI'
+ phdu.header['DATAMODE'] = dmode
+ phdu.header['OBS_ID'] = datetime.datetime.utcnow().strftime('%s')
+
+ hdl.append(phdu)
+
+ return hdl
+
+
+def _mk_ts(cefile=None):
+
+ if cefile is None:
+ return datetime.datetime.utcnow().strftime('%j_%H%M%S%f')[:-3]
+ else:
+ return cefile.split('_')[-1]
+
+
+if __name__ == "__main__":
+
+ if len(sys.argv) > 2:
+ cefile, fitsfile = sys.argv[1:3]
+ else:
+ cefile = sys.argv[1]
+ fitsfile = None
+
+ convert_ce(cefile, fitsfile=fitsfile)
--
GitLab