大致完成了仪器效应改正的代码, 修改了一些csstdata的小bug, deepcr部分还没添加

csst/core/data.py

@@ -2,6 +2,7 @@ from collections import OrderedDict
 
 import astropy.io.fits as fits
 from astropy.io.fits import HDUList, PrimaryHDU
+import numpy as np
 
 from csst.core.exception import CsstException
 
@@ -104,9 +105,13 @@ class CsstData:
         print("save L1 image to a fits file with name " + filename)
         try:
             self._l1hdr_global.set('TYPE', imgtype, 'Type of Level 1 data')
-            pri_hdu = PrimaryHDU(header=self._l1hdr_global)
-            hdulist = HDUList([pri_hdu, self._l1data[imgtype]])
-            hdulist.writeto(filename)
+            hdulist = fits.HDUList(
+                [
+                    fits.PrimaryHDU(header=self._l1hdr_global),
+                    fits.ImageHDU(header=self._l1data[imgtype].header, data=self._l1data[imgtype].data),
+                ]
+            )
+            hdulist.writeto(filename, overwrite=True)
         except Exception as e:
             print(e)
 

csst/msc/data.py

@@ -3,6 +3,7 @@ import astropy.io.fits as fits
 from astropy.io.fits import HDUList, PrimaryHDU, ImageHDU
 from astropy.io.fits.header import Header
 from ..core.data import CsstData, INSTRUMENT_LIST
+import numpy as np
 
 
 __all__ = ["CsstMscData", "CsstMscImgData"]
@@ -12,10 +13,8 @@ class CsstMscData(CsstData):
     _l1img_types = {'sci': True, 'weight': True, 'flag': True}
 
     def __init__(self, priHDU, imgHDU, **kwargs):
-        super(CsstData, self).__init__(priHDU, imgHDU, **kwargs)
+        super(CsstMscData, self).__init__(priHDU, imgHDU, **kwargs)
         self._l1hdr_global = priHDU.header.copy()
-        #         self._l1hdr_global['SIMPLE']  =  'T' #/ conforms to FITS standard
-        #         self._l1hdr_global['NAXIS']  =  0kkjk
         self._l1data['sci'] = ImageHDU()
         self._l1data['weight'] = ImageHDU()
         self._l1data['flag'] = ImageHDU()
@@ -74,8 +73,18 @@ class CsstMscData(CsstData):
     def set_l1data(self, imgtype, img):
         """ set L1 data """
         try:
-            if self._l1img_types[imgtype]:
-                self._l1data[imgtype].data = img.copy()
+            if imgtype == 'sci':
+                self._l1data[imgtype].header['EXTNAME'] = 'img'
+                self._l1data[imgtype].header['BUNIT'] = 'e/s'
+                self._l1data[imgtype].data = img.astype(np.float32) / self._l1hdr_global['exptime']
+            elif imgtype == 'weight':
+                self._l1data[imgtype].header['EXTNAME'] = 'wht'
+                self._l1data[imgtype].data = img.astype(np.float32)
+            elif imgtype == 'flag':
+                self._l1data[imgtype].header['EXTNAME'] = 'flg'
+                self._l1data[imgtype].data = img.astype(np.uint16)
+            else:
+                raise TypeError('unknow type image')
         except Exception as e:
             print(e)
         print('save image data to l1data')
@@ -93,7 +102,7 @@ class CsstMscData(CsstData):
 class CsstMscImgData(CsstMscData):
     def __init__(self, priHDU, imgHDU, **kwargs):
         # print('create CsstMscImgData')
-        super(CsstMscData, self).__init__(priHDU, imgHDU, **kwargs)
+        super(CsstMscImgData, self).__init__(priHDU, imgHDU, **kwargs)
 
     def __repr__(self):
         return "<CsstMscImgData: {} {}>".format(self.instrument, self.detector)
@@ -124,14 +133,16 @@ class CsstMscImgData(CsstMscData):
         """
 
         try:
-            hl = fits.open(fp)
-            instrument = hl[0].header.get('INSTRUME')  # strip or not?
-            detector = hl[0].header.get('DETECTOR')  # strip or not?
-            print("@CsstMscImgData: reading data {} ...".format(fp))
-            assert instrument in INSTRUMENT_LIST
-            if instrument == 'MSC' and 6 <= int(detector[3:5]) <= 25:
-                # multi-band imaging
-                data = CsstMscImgData(hl[0], hl[1], instrument=instrument, detector=detector)
-                return data
+            with fits.open(fp) as hdulist:
+                instrument = hdulist[0].header.get('INSTRUME')  # strip or not?
+                detector = hdulist[0].header.get('DETECTOR')  # strip or not?
+                print("@CsstMscImgData: reading data {} ...".format(fp))
+                assert instrument in INSTRUMENT_LIST
+                if instrument == 'MSC' and 6 <= int(detector[3:5]) <= 25:
+                    # multi-band imaging
+                    hdu0 = hdulist[0].copy()
+                    hdu1 = hdulist[1].copy()
+                    data = CsstMscImgData(hdu0, hdu1, instrument=instrument, detector=detector)
+                    return data
         except Exception as e:
             print(e)

csst/msc/instrument.py

@@ -2,37 +2,107 @@ from collections import OrderedDict
 from abc import ABCMeta, abstractmethod
 from enum import Enum
 import numpy as np
+from deepCR import deepCR
 from ..core.processor import CsstProcStatus, CsstProcessor
 
 
 class CsstMscInstrumentProc(CsstProcessor):
     _status = CsstProcStatus.empty
-    _switches = {'crosstalk': False, 'nonlinear': False, 'deepcr': False, 'cti': False, 'brighterfatter': False}
+    _switches = {'crosstalk': False, 'nonlinear': False,
+                 'deepcr': False, 'cti': False, 'brighterfatter': False}
 
     def __init__(self):
         pass
 
-    def _do_crosstalk(self):
-        if self._switches['crosstalk']:
-            print('Crosstalk correction')
-
-    def _do_nonlinear(self):
-        if self._switches['nonlinear']:
-            print('Nonlinear effect correction')
-
-    def _do_deepcr(self):
-        if self._switches['deepcr']:
-            print('Deep CR operation')
-        else:
-            print('Laplace CR correction')
-
-    def _do_cti(self):
-        if self._switches['cti']:
-            print('CTI effect correction')
-
-    def _do_brighterfatter(self):
-        if self._switches['brighterfatter']:
-            print('Brighter-Fatter effect correction')
+    def _do_fix(self, raw, bias, dark, flat, exptime):
+        '''仪器效应改正
+
+        将raw扣除本底, 暗场, 平场. 并且避免了除0
+        Args:
+            raw: 科学图生图
+            bias: 本底
+            dark: 暗场
+            flat: 平场
+            exptime: 曝光时长
+        '''
+        self.__l1img = np.divide(
+            raw - bias - dark * exptime, flat,
+            out=np.zeros_like(raw, float),
+            where=(flat != 0),
+        )
+
+    def _do_badpix(self, flat):
+        '''坏像元标记
+
+        因为探测器本身原因造成不可用于有效科学研究的像元, 像元响应小于0.5中值或大于1.5中值
+
+        Args:
+            flat: 平场
+        '''
+        med = np.median(flat)
+        flg = (flat < 0.5 * med) | (1.5 * med < flat)
+        self.__flagimg = self.__flagimg | (flg * 1)
+
+    def _do_hot_and_warm_pix(self, dark, exptime, rdnoise):
+        '''热像元与暖像元标记
+
+        因为探测器本身原因造成的影响科学研究结果的像元. 像元在曝光时长积分时间内, 
+        热像元为: 暗流计数大于探测器平均读出噪声的平方.
+        暖像元为: 暗流计数大于0.5倍读出噪声的平方, 但小于1倍读出噪声的平方.
+
+        Args:
+            dark: 暗场
+            exptime: 曝光时长
+            rdnoise: 读出噪声
+        '''
+        tmp = dark * exptime
+        tmp[tmp < 0] = 0
+        flg = 1 * rdnoise ** 2 <= tmp  # 不确定是否包含 暂定包含
+        self.__flagimg = self.__flagimg | (flg * 2)
+        flg = (0.5 * rdnoise ** 2 < tmp) & (tmp < 1 * rdnoise ** 2)
+        self.__flagimg = self.__flagimg | (flg * 4)
+
+    def _do_over_saturation(self, raw):
+        '''饱和溢出像元标记
+
+        饱和像元及流量溢出污染的像元.
+
+        Args:
+            raw: 科学图生图
+        '''
+        flg = raw == 65535
+        self.__flagimg = self.__flagimg | (flg * 8)
+
+    def _do_cray(self):
+        '''宇宙线像元标记
+
+        宇宙线污染的像元, 用deepCR进行标记
+        '''
+
+
+        
+        flg, _ = deepCR(self.__l1img)
+
+
+
+        self.__flagimg = self.__flagimg | (flg * 16)
+
+    def _do_weight(self, bias, gain, rdnoise, exptime):
+        '''权重图
+
+        Args:
+            bias: 本底
+            gain: 增益
+            rdnoise: 读出噪声
+            exptime: 曝光时长
+        '''
+        data = self.__l1img.copy()
+        data[self.__l1img < 0] = 0
+        weight_raw = 1. / (gain * data + rdnoise ** 2)
+        bias_weight = np.std(bias)
+        weight = 1. / (1. / weight_raw + 1. / bias_weight) * exptime ** 2
+        weight[self.__flagimg > 0] = 0
+        self.__weightimg = weight
 
     def prepare(self, **kwargs):
         for name in kwargs:
@@ -40,19 +110,24 @@ class CsstMscInstrumentProc(CsstProcessor):
 
     def run(self, data):
         if type(data).__name__ == 'CsstMscImgData' or type(data).__name__ == 'CsstMscSlsData':
-            self.__l1img = data.get_l0data(copy=True)
-            self.__weightimg = np.random.uniform(0, 1, (9216, 9232))
-            self.__flagimg = np.random.uniform(0, 1, (9216, 9232))
-
+            raw = data.get_l0data()
+            self.__l1img = raw.copy()
+            self.__weightimg = np.zeros_like(raw)
+            self.__flagimg = np.zeros_like(raw, dtype=np.uint16)
+
+            exptime = data.get_l0keyword('pri', 'EXPTIME')
+            gain = data.get_l0keyword('img', 'GAIN1')
+            rdnoise = data.get_l0keyword('img' ,'RDNOISE1')
             flat = data.get_flat()
             bias = data.get_bias()
             dark = data.get_dark()
             print('Flat and bias correction')
-            self._do_crosstalk()
-            self._do_nonlinear()
-            self._do_cti()
-            self._do_deepcr()
-            self._do_brighterfatter()
+            self._do_fix(raw, bias, dark, flat, exptime)
+            self._do_badpix(flat)
+            self._do_hot_and_warm_pix(dark, exptime, rdnoise)
+            self._do_over_saturation(raw)
+            # self._do_cray()
+            self._do_weight(bias, gain, rdnoise, exptime)
 
             print('fake to finish the run and save the results back to CsstData')
 
@@ -61,7 +136,8 @@ class CsstMscInstrumentProc(CsstProcessor):
             data.set_l1data('flag', self.__flagimg)
 
             print('Update keywords')
-            data.set_l1keyword('SOMEKEY', 'some value', 'Test if I can append the header')
+            data.set_l1keyword('SOMEKEY', 'some value',
+                               'Test if I can append the header')
 
             self._status = CsstProcStatus.normal
         else:

local_test_instrument.py

+from csst.msc import CsstMscImgData
+from csst.msc.instrument import CsstMscInstrumentProc
+from astropy.io.fits import getdata
+
+
+fp = "/data/cali_20211012/L0/150s/MSC_MS_210525120000_100000000_13_raw.fits"
+bs = "/data/ref/MSC_CLB_210525190000_100000014_13_combine.fits"
+dk = "/data/ref/MSC_CLD_210525192000_100000014_13_combine.fits"
+ft = "/data/ref/MSC_CLF_210525191000_100000014_13_combine.fits"
+op = "/data/test/"
+data = CsstMscImgData.read(fp)
+# print(repr(data._l1data.header))
+data.set_bias(getdata(bs))
+data.set_dark(getdata(dk))
+data.set_flat(getdata(ft))
+
+proc = CsstMscInstrumentProc()
+proc.run(data)
+data.save_l1data('sci', op+'sci.fits')
+data.save_l1data('flag', op+'flag.fits')
+data.save_l1data('weight', op+'weight.fits')