add L1 msc instrument pipeline modules

measurement_pipeline/L1_pipeline/csst_msc_instrument/__init__.py

+import os
+from .csst_msc_mbi_instrument import core_msc_l1_mbi_instrument
+__all__ = ["core_msc_l1_mbi_instrument"]
\ No newline at end of file

measurement_pipeline/L1_pipeline/csst_msc_instrument/format.py

+"""
+Identifier:     csst-l1/msc/csst_msc_instrument/csst_msc_instrument/format.py
+Name:           format.py
+Description:    Data format conversion.
+Author:         Li Shao (shaoli@nao.cas.cn)
+Created:        2023-11-16
+Modified-History:
+    2023-11-16, Li Shao, created
+    2023-11-22, Li Shao, split from common_detector
+    2023-11-23, Li Shao, change the way to get CCDArrayConfig
+    2023-12-21, Li Shao, add function convert_format_one_cube
+    2024-02-06, Li Shao, move to csst_msc_instrument
+    2024-03-05, Li Shao, merge all three functions into one
+"""
+
+import numpy as np
+
+from .config import get_array_config
+
+__all__ = ["convert_format", ]
+
+
+def convert_format(image: np.ndarray,
+                   output_format: str,
+                   array_config_name: str = 'CCDArrayConfig',
+                   ) -> np.ndarray:
+    """
+    Convert image format.
+
+    The function will automatically detect input format and convert to the output format.
+    The output format includes: "raw", "cube", "sky".
+
+    Parameters
+    ----------
+    image : numpy.ndarray
+        Input 2-D image.
+    output_format : str
+        The output format name: "raw", "cube" or "sky".
+    array_config_name : str, optional
+        The name of the pixel array configuration class. Default is "CCDArrayConfig".
+
+    Returns
+    -------
+    numpy.ndarray
+        Converted image. If convert from "sky" to other format, non-active region will be filled with 0.
+    """
+
+    config = get_array_config(array_config_name)
+    conf_raw = config('raw')
+    conf_cube = config('cube')
+    conf_sky = config('sky')
+    shape_input = image.shape
+    if np.array_equal(shape_input, conf_raw.get_shape(full=True)):
+        iconf = conf_raw
+    elif np.array_equal(shape_input, conf_cube.get_shape(full=True)):
+        iconf = conf_cube
+    elif np.array_equal(shape_input, conf_sky.get_shape(full=True)):
+        iconf = conf_sky
+    else:
+        raise ValueError('Input data shape is inconsistent with array configuration.')
+
+    if output_format == 'raw':
+        oconf = conf_raw
+    elif output_format == 'cube':
+        oconf = conf_cube
+    elif output_format == 'sky':
+        oconf = conf_sky
+    else:
+        raise ValueError('Invalid output type: "{}". Must be "raw", "cube" or "sky".'.format(output_format))
+
+    if iconf == oconf:
+        print('Warning: input and output is the same. No format conversion is performed.')
+        return image
+
+    if iconf.image_format == 'sky' or output_format == 'sky':
+        rtype = 'active'
+    else:
+        rtype = 'channel'
+    oshape = oconf.get_shape(full=True)
+    output = np.zeros(oshape, dtype=image.dtype)
+    for i in range(iconf.nchan):
+        cutout = iconf.get_cutout(i+1, image, rtype=rtype)
+        oconf.set_cutout(i+1, output, cutout, rtype=rtype)
+
+    return output

measurement_pipeline/L1_pipeline/csst_msc_instrument/image/__init__.py

+"""
+Identifier:     csst-l1/msc/csst_msc_instrument/csst_msc_instrument/image/__init__.py
+Name:           __init__.py
+Description:    Instrument effect module (basic operations).
+Author:         Li Shao (shaoli@nao.cas.cn)
+Created:        2023-11-27
+Modified-History:
+    2023-11-27, Li Shao, created
+    2024-02-06, Li Shao, move to csst_msc_instrument
+"""
+
+from .basic import subtract_bias, subtract_dark
+from .crosstalk import remove_crosstalk
+from .gain import apply_gain
+from .overscan import correct_overscan
+
+__all__ = ['correct_overscan', 'apply_gain', 'subtract_bias', 'remove_crosstalk', 'subtract_dark', ]
\ No newline at end of file

measurement_pipeline/L1_pipeline/csst_msc_instrument/image/basic.py

+"""
+Identifier:     csst-l1/msc/csst_msc_instrument/csst_msc_instrument/image/basic.py
+Name:           basic.py
+Description:    Basic instrument effect correction.
+Author:         Li Shao (shaoli@nao.cas.cn)
+Created:        2023-11-16
+Modified-History:
+    2023-11-16, Li Shao, created
+    2023-11-22, Li Shao, split from common_detector
+    2024-02-06, Li Shao, move to csst_msc_instrument
+"""
+
+import numpy as np
+
+__all__ = ["subtract_bias", "subtract_dark", ]
+
+
+def subtract_bias(image: np.ndarray,
+                  image_err: np.ndarray,
+                  bias: np.ndarray,
+                  bias_err: np.ndarray,
+                  ) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Subtract bias.
+
+    Subtract bias.
+
+    Parameters
+    ----------
+    image : numpy.ndarray
+        The input image to be corrected.
+    image_err : numpy.ndarray
+        The uncertainty of input image.
+    bias : numpy.ndarray
+        The input bias to be subtracted.
+    bias_err : numpy.ndarray
+        The uncertainty of input bias.
+
+    Returns
+    -------
+    output : np.ndarray
+        Output corrected image.
+    output_err : np.ndarray
+        Output uncertainty map.
+    """
+
+    output = image - bias
+    output_err = np.sqrt(image_err ** 2 + bias_err ** 2)
+
+    return output, output_err
+
+
+def subtract_dark(image: np.ndarray,
+                  image_err: np.ndarray,
+                  dark: np.ndarray,
+                  dark_err: np.ndarray,
+                  tdark_image: float | int = 1.0,
+                  ) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Subtract dark current.
+
+    Subtract dark current.
+
+    Parameters
+    ----------
+    image : numpy.ndarray
+        The input image to be corrected.
+    image_err : numpy.ndarray
+        The uncertainty of input image.
+    dark : numpy.ndarray
+        The input dark current image to be subtracted.
+    dark_err : numpy.ndarray
+        The uncertainty of input dark current.
+    tdark_image : float or int, optional
+        The effective dark current cumulation time of input image. Default value is 1.0.
+
+    Returns
+    -------
+    output : np.ndarray
+        Output corrected image.
+    output_err : np.ndarray
+        Output uncertainty map.
+    """
+
+    output = image - dark * tdark_image
+    output_err = np.sqrt(image_err ** 2 + (dark_err * tdark_image) ** 2)
+
+    return output, output_err

measurement_pipeline/L1_pipeline/csst_msc_instrument/image/crosstalk.py

+"""
+Identifier:     csst-l1/msc/csst_msc_instrument/csst_msc_instrument/image/crosstalk.py
+Name:           crosstalk.py
+Description:    Crosstalk correction.
+Author:         Tianmeng Zhang (zhangtm@nao.cas.cn)
+Created:        2023-11-16
+Modified-History:
+    2023-11-16, Tianmeng Zhang, created
+    2023-11-22, Li Shao, split from common_detector
+    2024-02-06, Li Shao, move to csst_msc_instrument
+"""
+
+import numpy as np
+
+__all__ = ["remove_crosstalk", ]
+
+
+def remove_crosstalk(data: np.ndarray,
+                     crosstalk_coe: np.ndarray,
+                     hdu_num: int = 16,
+                     switch: bool = False,
+                     ) -> np.ndarray:
+    """
+    Function to remove the crosstalk between different channel.
+
+    Function to remove the crosstalk between different channel.
+
+    Parameters
+    ----------
+    data : numpy.ndarray
+        The input data array which is needed to remove the crosstalk from each hdu, 3D-array.
+    crosstalk_coe : numpy.ndarray
+        The coefficients of crosstalk, in size [hdu_num, hdu_num].
+    hdu_num : int, optional
+        The number of extension for each detector.
+    switch : bool, optional
+        If True, do the crosstalk correction. If False, return original data. Default is False.
+
+    Returns
+    -------
+    numpy.ndarray
+        The data array after crosstalk correction, 3D-array.
+    """
+
+    data_cor = np.copy(data)  # backup original data array
+
+    if switch:
+        for i in range(hdu_num):
+            for j in range(hdu_num):
+                data_cor[i, :, :] = data[i, :, :] + data[j, :, :] * crosstalk_coe[i, j]
+    else:
+        data_cor = data
+        print('No crosstalk correction, return the original array')
+
+    return data_cor

measurement_pipeline/L1_pipeline/csst_msc_instrument/image/gain.py

+"""
+Identifier:     csst-l1/msc/csst_msc_instrument/csst_msc_instrument/image/gain.py
+Name:           gain.py
+Description:    Gain map related function.
+Author:         Li Shao (shaoli@nao.cas.cn)
+Created:        2023-12-22
+Modified-History:
+    2023-12-22, Li Shao, created
+    2024-02-06, Li Shao, move to csst_msc_instrument
+"""
+
+import numpy as np
+
+from ..config import get_array_config
+
+__all__ = ['apply_gain', 'make_gainmap_from_channel_value', ]
+
+
+def apply_gain(image: np.ndarray,
+               gainmap: np.ndarray,
+               ) -> np.ndarray:
+    """
+    Apply gain map.
+
+    The gain map is in unit of electrons/ADU. The output image will be in unit of electrons or electrons/s.
+
+    Parameters
+    ----------
+    image : numpy.ndarray
+        Input image, in unit of ADU or ADU/s.
+    gainmap : numpy.ndarray
+        Gain map, in unit of electrons/ADU.
+
+    Returns
+    -------
+    numpy.ndarray
+        Output image, in unit of electrons or electrons/s.
+    """
+
+    return image * gainmap
+
+
+def make_gainmap_from_channel_value(gain_values: list[float] | tuple[float] | np.ndarray,
+                                    array_config_name: str = 'CCDArrayConfig',
+                                    image_format: str = 'raw',
+                                    ) -> np.ndarray:
+    """
+    Make gain map from channel gain values.
+
+    Parameters
+    ----------
+    gain_values : list[float] or tuple[float] or numpy.ndarray
+        Gain value of each channel.
+    array_config_name : str, optional
+        The name of the pixel array configuration class. Default is "CCDArrayConfig".
+    image_format : str, optional
+        The output image format: "raw", "cube" or "sky". See config.CCDArrayConfig().
+
+    Returns
+    -------
+    numpy.ndarray
+        Output gain map.
+    """
+
+    aconf = get_array_config(array_config_name)(image_format)
+    if len(gain_values) != aconf.nchan:
+        raise ValueError('The length of input gain values does not match with number of channels.')
+
+    full_shape = aconf.get_shape(full=True)
+    output = np.zeros(full_shape, dtype=np.float32)
+    for i in range(aconf.nchan):
+        cutout = aconf.get_cutout(i+1, output, rtype='channel') + gain_values[i]
+        aconf.set_cutout(i+1, output, cutout, rtype='channel')
+
+    return output
\ No newline at end of file

measurement_pipeline/L1_pipeline/csst_msc_instrument/image/overscan.py

+"""
+Identifier:     csst-l1/msc/csst_msc_instrument/csst_msc_instrument/image/overscan.py
+Name:           overscan.py
+Description:    Overscan correction and data format conversion.
+Author:         Li Shao (shaoli@nao.cas.cn)
+Created:        2023-11-16
+Modified-History:
+    2023-11-16, Li Shao, created
+    2023-11-22, Li Shao, split from common_detector
+    2024-02-06, Li Shao, move to csst_msc_instrument
+"""
+
+from typing import Callable
+import numpy as np
+from scipy.ndimage import median_filter
+from scipy.signal import savgol_filter
+from scipy.special import erf
+from scipy.stats import trim_mean, t, sigmaclip
+from astropy.stats import mad_std, sigma_clip
+
+from ..config import get_array_config
+
+__all__ = ["average_overscan", "smooth_overscan", "correct_overscan", ]
+
+
+def average_overscan(image: np.ndarray,
+                     axis: int | None = None,
+                     cen_method: str | Callable = 'trim_mean',
+                     std_method: str | Callable = 'mad_std',
+                     clip_flag: bool = False,
+                     clip_sigma: float = 3.5,
+                     clip_cenfunc: str | Callable = 'median',
+                     clip_stdfunc: str | Callable = 'std',
+                     ) -> tuple[float | np.ndarray, float | np.ndarray]:
+    """
+    Average overscan data along one direction or overall.
+
+    Average overscan data along one direction or overall.
+
+    Parameters
+    ----------
+    image : numpy.ndarray
+        Input image: overscan region cutout.
+    axis : int, optional
+        Axis along which the operation is performed. If None (default), use the full array.
+    cen_method : str or Callable, optional
+        Algorithm to calculate average value: "trim_mean" (default), "median", "mean" or Callable function.
+    std_method : str or Callable, optional
+        Algorithm to calculate standard deviation: "mad_std" (default), "std" or Callable function.
+    clip_flag : bool, optional
+        If True, use sigma-clipping.
+    clip_sigma : float, optional
+        Clipping threshold for sigma-clipping. Default is 3.5.
+    clip_cenfunc : str or Callable, optional
+        Center value function for sigma-clipping: "median" (default), "mean" or Callable function.
+    clip_stdfunc : str or Callable, optional
+        Standard deviation function for sigma-clipping: "std" (default), "mad_std" or Callable function.
+
+    Returns
+    -------
+    avg : numpy.ndarray
+        Averaged overscan data (one value or 1D array).
+    err : numpy.ndarray
+        The uncertainty.
+    """
+
+    # sigma-clipping
+    if clip_flag and cen_method != 'trim_mean':
+        data = sigma_clip(image, sigma=clip_sigma, axis=axis, masked=True,
+                          cenfunc=clip_cenfunc, stdfunc=clip_stdfunc)
+
+        if cen_method == 'median':
+            avg = np.ma.median(data)
+        elif cen_method == 'mean':
+            avg = np.ma.mean(data)
+        else:
+            avg = cen_method(data)
+
+        if std_method == 'mad_std':
+            err = mad_std(data, axis=axis)
+        elif std_method == 'std':
+            err = np.ma.std(data, axis=axis)
+        else:
+            err = std_method(data, axis=axis)
+
+    # without clipping
+    else:
+        if cen_method == 'trim_mean':
+            avg = trim_mean(image, 0.1, axis=axis)
+        elif cen_method == 'median':
+            avg = np.median(image, axis=axis).astype(float)
+        elif cen_method == 'mean':
+            avg = np.mean(image, axis=axis)
+        else:
+            avg = cen_method(image, axis=axis)
+
+        if std_method == 'mad_std':
+            err = mad_std(image, axis=axis)
+        elif std_method == 'std':
+            err = np.std(image, axis=axis)
+        else:
+            err = std_method(image, axis=axis)
+
+    # uncertainty of the average
+    ny, nx = image.shape
+    if axis is None:
+        n = ny * nx
+    elif axis == 0:
+        n = ny
+    else:
+        n = nx
+    if cen_method == 'trim_mean':
+        n = n * 0.8
+    err = err / np.sqrt(n) * t.ppf(erf(1), df=n - 1)  # uncertainty of the mean
+    if cen_method == 'median' or cen_method is np.median or cen_method is np.ma.median:
+        err *= 1.2533
+
+    return avg, err
+
+
+def smooth_overscan(data: np.ndarray,
+                    filter_width: int,
+                    filter_method: str = 'savgol',
+                    filter_polyorder: int = 2,
+                    ) -> np.ndarray:
+    """
+    Smooth 1D averaged overscan data.
+
+    Smooth 1D averaged overscan data.
+
+    Parameters
+    ----------
+    data : numpy.ndarray
+        Input 1-D overscan data.
+    filter_width : int
+        The width of filter. Must be positive.
+    filter_method : str, optional
+        Process algorithm: "savgol" (default) or "median".
+    filter_polyorder : int, optional
+        The order of polynomial used for Savitzky-Golay filter.
+
+    Returns
+    -------
+    numpy.ndarray
+        Smoothed 1-D overscan data.
+    """
+
+    if filter_method == 'savgol':
+        dd = savgol_filter(data, filter_width * 2, filter_polyorder)
+        diff = data - dd
+        index = np.abs(diff) < sigmaclip(diff, 5, 5)[0].std() * 5  # clip outliers
+        x = np.arange(len(data))
+        dd = np.interp(x, x[index], data[index])
+        output = savgol_filter(dd, filter_width, filter_polyorder)
+    else:
+        output = median_filter(data, size=filter_width)
+
+    return output
+
+
+def correct_overscan(image: np.ndarray,
+                     array_config_name: str = 'CCDArrayConfig',
+                     region_type: str = 'x',
+                     correct_type: str = 'line',
+                     gap: int = 10,
+                     cen_method: str | Callable = 'trim_mean',
+                     std_method: str | Callable = 'mad_std',
+                     clip_flag: bool = False,
+                     clip_sigma: float = 3.5,
+                     clip_cenfunc: str | Callable = 'median',
+                     clip_stdfunc: str | Callable = 'std',
+                     filter_width: int = 0,
+                     filter_method: str = 'savgol',
+                     filter_polyorder: int = 2,
+                     clean_inactive: bool = True,
+                     ) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Overscan correction.
+
+    Overscan correction. All parameters are optional inputs.
+
+    Parameters
+    ----------
+    image : numpy.ndarray
+        Input 2-D image.
+    array_config_name : str, optional
+        The name of the pixel array configuration class. Default is "CCDArrayConfig".
+    region_type : str, optional
+        The type of overscan region. "x" (default) is serial overscan, "y" is parallel overscan.
+    correct_type : str, optional
+        Correction type. "line" (default) is correct by row or column. "all" is correct with all pixels.
+    gap : int, optional
+        Number of pixels after active pixels not to use for overscan correction.
+    cen_method : str or Callable, optional
+        Algorithm to calculate average value: "trim_mean" (default), "median", "mean" or Callable function.
+    std_method : str or Callable, optional
+        Algorithm to calculate standard deviation: "mad_std" (default), "std" or Callable function.
+    clip_flag : bool, optional
+        If True, use sigma-clipping.
+    clip_sigma : float, optional
+        Clipping threshold for sigma-clipping. Default is 3.5.
+    clip_cenfunc : str or Callable, optional
+        Center value function for sigma-clipping: "median" (default), "mean" or Callable function.
+    clip_stdfunc : str or Callable, optional
+        Standard deviation function for sigma-clipping: "std" (default), "mad_std" or Callable function.
+    filter_width : int, optional
+        The width of filter. If smaller than 2, then smoothing will be skipped.
+    filter_method : str, optional
+        Filtering algorithm: "savgol" (default) or "median".
+    filter_polyorder : int, optional
+        The order of polynomial used for Savitzky-Golay filter.
+    clean_inactive : bool, optional
+        If True, all pixels out of active region are set to zero. Default is True.
+
+    Returns
+    -------
+    output : numpy.ndarray
+        Overscan corrected image.
+    output_err : numpy.ndarray
+        Uncertainty image.
+
+    Examples
+    --------
+    image_corr, err_corr = correct_overscan(image, region_type='x', correct_type='line')
+    image_corr, err_corr = correct_overscan(image, region_type='y', correct_type='all', clean_inactive=True)
+    """
+
+    aconf = get_array_config(array_config_name)('raw')
+    if correct_type == 'line':
+        if region_type == 'x':
+            axis = 1
+        else:
+            axis = 0
+    else:
+        axis = None
+
+    if clean_inactive:
+        output = np.zeros(image.shape)
+    else:
+        output = image.astype(float)
+    output_err = np.zeros(image.shape)
+
+    for i in range(aconf.nchan):
+        overscan = aconf.get_cutout(i+1, image, rtype='{}_overscan'.format(region_type))
+
+        # remove the first pixels (to avoid potential CTI contamination)
+        if gap > 0:
+            if region_type == 'x':
+                overscan = overscan[:, gap:]
+            else:
+                overscan = overscan[gap:, :]
+
+        # calculate average and uncertainty
+        avg_value, err_value = average_overscan(overscan, axis=axis,
+                                                cen_method=cen_method, std_method=std_method,
+                                                clip_flag=clip_flag, clip_sigma=clip_sigma,
+                                                clip_cenfunc=clip_cenfunc, clip_stdfunc=clip_stdfunc)
+
+        # smooth if required
+        if filter_width > 1 and axis is not None:
+            avg_value = smooth_overscan(avg_value, filter_width,
+                                        filter_method=filter_method, filter_polyorder=filter_polyorder)
+            err_value = err_value / np.sqrt(filter_width) * 1.2533
+
+        # correct the active region
+        x1, x2, y1, y2 = aconf.get_boundary(i+1, rtype='active')
+        nx, ny = x2 - x1, y2 - y1
+        if axis is None:
+            output[y1:y2, x1:x2] = image[y1:y2, x1:x2] - avg_value
+            output_err[y1:y2, x1:x2] = err_value
+        elif axis == 0:
+            output[y1:y2, x1:x2] = image[y1:y2, x1:x2] - np.tile(avg_value, ny).reshape((ny, nx))
+            output_err[y1:y2, x1:x2] = np.tile(err_value, ny).reshape((ny, nx))
+        else:
+            output[y1:y2, x1:x2] = image[y1:y2, x1:x2] - np.repeat(avg_value, nx).reshape((ny, nx))
+            output_err[y1:y2, x1:x2] = np.repeat(err_value, nx).reshape((ny, nx))
+
+    return output, output_err

measurement_pipeline/L1_pipeline/csst_msc_instrument/instrument_cor_demo.head

+V_INST  = '0.0.2   '           / version of instrument correction               T_INST  = '2023-12-29T04:50:46.710' / timestamp of instrument correction        S_INST  =                    0 / status of instrument correction                S_OVSCAN=                    0 / status of overscan correction                  S_GAIN  =                    0 / status of gain correction                      R_GAIN  = 'file    '           / reference gain                                 S_BIAS  =                    0 / status of bias frame correction                R_BIAS  = 'file    '           / reference bias                                 S_CROSST=                    0 / status of crosstalk correction                 R_CROSST= 'file    '           / reference crosstalk                            S_CTI   =                    1 / status of CTI correction                       R_CTI   = 'file    '           / reference CTI                                  S_DARK  =                    0 / status of dark frame correction                R_DARK  = 'file    '           / reference dark                                 S_NLIN  =                    1 / status of non-linear correction                R_NLIN  = 'file    '           / reference non-linear                           S_SHUT  =                    1 / status of shutter effect correction            R_SHUT  = 'file    '           / reference shutter effect                       S_FLAT  =                    0 / status of flat frame correction                R_FLAT  = 'file    '           / reference flat                                 S_CRS   =                    1 / status of cosmic rays mask                     R_CRS   = 'deepCR_model'       / method and config of cosmic rays mask          CRCOUNT =               240466 / cosmic rays pixel counts                       S_SAT   =                    1 / status of satellite correction                 R_SAT   = 'file    '           / reference satellite correction                 S_FRINGE=                    1 / status of fringe correction                    R_FRINGE= 'file    '           / reference fringe                               SKY_BKG =               0.0323 / estimated sky background (e-/s per pixel)      SKY_RMS =               0.0374 / standard dev of frame background (e-/s)        SATURATE=             365.9759 / flux limit of saturated pixel (e-/s)           END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             

measurement_pipeline/L1_pipeline/ref_combine.py

+import numpy as np
+from astropy.io import fits
+
+
+def array_combine(ndarray, mode="mean") -> np.ndarray:
+    """ Function to combine 3-D data array
+
+    Parameters
+    ----------
+    ndarray: array, input data cube (3D)
+    mode: mean, median, sum, mean_clip, median_clip, default is mean
+    """
+    if mode == "median":
+        array = np.median(ndarray, axis=0)
+    elif mode == "median_clip":
+        ndarray = np.sort(ndarray, axis=0)[1:-1]
+        array = np.median(ndarray, axis=0)
+    elif mode == "sum":
+        array = np.sum(ndarray, axis=0)
+    elif mode == "mean":
+        array = np.mean(ndarray, axis=0)
+    elif mode == "mean_clip":
+        ndarray = np.sort(ndarray, axis=0)[1:-1]
+        array = np.mean(ndarray, axis=0)
+    return array
+
+
+def load_bias(path: str) -> np.ndarray:
+    with fits.open(path) as hdul:
+        du = hdul[1].data
+    du = du.astype(int)
+    return du
+
+
+def load_dark(path: str, bias) -> np.ndarray:
+    with fits.open(path) as hdul:
+        du = hdul[1].data
+        hu = hdul[0].header
+    du = du.astype(int)
+    du = du - bias
+    du = du / hu["EXPTIME"]
+    return du
+
+
+def load_flat(path: str, bias, dark) -> np.ndarray:
+    with fits.open(path) as hdul:
+        du = hdul[1].data
+        hu = hdul[0].header
+    du = du.astype(int)
+    du = du - bias - dark * hu["EXPTIME"]
+    du = du / hu["EXPTIME"]
+    du = du / np.median(du)
+    return du
+
+
+def combine(func, mode: str, path_list, *args) -> np.ndarray:
+    du_list = [func(path, *args) for path in path_list]
+    du = array_combine(du_list, mode)
+    return du
+
+
+def combine_images(b_p_lst, d_p_lst, f_p_lst, mode_list=["median", "median", "median", ]):
+    """
+
+    Parameters
+    ----------
+    b_p_lst:
+        List of currently ccd number bias file path
+    d_p_lst:
+        List of currently ccd number dark file path
+    f_p_lst:
+        List of currently ccd number flat file path
+    mode_list:
+        [0] bias combine mode
+        [1] dark combine mode
+        [2] flat combine mode
+        mean, median, sum, mean_clip, median_clip
+    """
+    bias = combine(load_bias, mode_list[0], b_p_lst)
+    dark = combine(load_dark, mode_list[1], d_p_lst, bias)
+    flat = combine(load_flat, mode_list[2], f_p_lst, bias, dark)
+    return bias.copy(), dark.copy(), flat.copy()

measurement_pipeline/run_ref_combine.py

+import os
+from glob import glob
+from astropy.io import fits
+from L1_pipeline.ref_combine import combine_images
+
+ref_path = "/public/share/yangxuliu/CSSOSDataProductsSims/outputs_cali/"
+output_path = "/public/home/fangyuedong/project/calib_data"
+
+def combine_ref_func(ref_path, output_path, num="01"):
+    bias_path_list = glob(ref_path + '*/CSST_MSC_MS_BIAS_*_' + num + '_*')
+    dark_path_list = glob(ref_path + '*/CSST_MSC_MS_DARK_*_' + num + '_*')
+    flat_path_list = glob(ref_path + '*/CSST_MSC_MS_FLAT_*_' + num + '_*')
+    bias, dark, flat = combine_images(b_p_lst=bias_path_list,
+                                    d_p_lst=dark_path_list,
+                                    f_p_lst=flat_path_list, )
+    bias_out_path = os.path.join(output_path, "bias_" + str(num) + ".fits")
+    dark_out_path = os.path.join(output_path, "dark_" + str(num) + ".fits")
+    flat_out_path = os.path.join(output_path, "flat_" + str(num) + ".fits")
+    hdu = fits.PrimaryHDU(bias)
+    hdu.writeto(bias_out_path)
+    hdu = fits.PrimaryHDU(dark)
+    hdu.writeto(dark_out_path)
+    hdu = fits.PrimaryHDU(flat)
+    hdu.writeto(flat_out_path)
+
+if __name__ == "__main__":
+    ref_path = "/public/share/yangxuliu/CSSOSDataProductsSims/outputs_cali/"
+    output_path = "/public/home/fangyuedong/project/calib_data"
+    num = '08'
+    combine_ref_func(
+        ref_path=ref_path,
+        output_path=output_path,
+        num=num
+    )
\ No newline at end of file