ChipUtils.py

import os
import galsim
import ctypes
import numpy as np
from astropy.io import fits
from datetime import datetime

from ObservationSim.Instrument.Chip import Effects as effects
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader

try:
    import importlib.resources as pkg_resources
except ImportError:
    # Try backported to PY<37 'importlib_resources'
    import importlib_resources as pkg_resources


def log_info(msg, logger=None):
    if logger:
        logger.info(msg)
    else:
        print(msg, flush=True)

def getChipSLSConf(chipID):
    confFile = ''
    if chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf']
    if chipID == 2: confFile = ['CSST_GV4.conf', 'CSST_GV3.conf']
    if chipID == 3: confFile = ['CSST_GU2.conf', 'CSST_GU1.conf']
    if chipID == 4: confFile = ['CSST_GU4.conf', 'CSST_GU3.conf']
    if chipID == 5: confFile = ['CSST_GV2.conf', 'CSST_GV1.conf']
    if chipID == 10: confFile = ['CSST_GI4.conf', 'CSST_GI3.conf']
    if chipID == 21: confFile = ['CSST_GI6.conf', 'CSST_GI5.conf']
    if chipID == 26: confFile = ['CSST_GV8.conf', 'CSST_GV7.conf']
    if chipID == 27: confFile = ['CSST_GU6.conf', 'CSST_GU5.conf']
    if chipID == 28: confFile = ['CSST_GU8.conf', 'CSST_GU7.conf']
    if chipID == 29: confFile = ['CSST_GV6.conf', 'CSST_GV5.conf']
    if chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf']
    return confFile

def generateHeader(chip, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime=150., timestamp = 1621915200):
    datetime_obs = datetime.utcfromtimestamp(timestamp)
    date_obs = datetime_obs.strftime("%y%m%d")
    time_obs = datetime_obs.strftime("%H%M%S")
    h_prim = generatePrimaryHeader(
        xlen=chip.npix_x, 
        ylen=chip.npix_y, 
        pointNum = str(pointing_ID),
        ra=ra_cen, 
        dec=dec_cen, 
        pixel_scale=chip.pix_scale,
        date=date_obs,
        time_obs=time_obs,
        im_type = im_type,
        exptime=exptime,
        chip_name=str(chip.chipID).rjust(2, '0')
        )
    h_ext = generateExtensionHeader(
        chip=chip,
        xlen=chip.npix_x, 
        ylen=chip.npix_y, 
        ra=ra_cen, 
        dec=dec_cen,
        pa=img_rot.deg, 
        gain=chip.gain, 
        readout=chip.read_noise, 
        dark=chip.dark_noise, 
        saturation=90000, 
        pixel_scale=chip.pix_scale, 
        pixel_size=chip.pix_size,
        xcen=chip.x_cen,
        ycen=chip.y_cen,
        extName='SCI',
        timestamp = timestamp,
        exptime = exptime,
        readoutTime = chip.readout_time)
    return h_prim, h_ext

def outputCal(chip, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, output_dir, exptime=150., timestamp = 1621915200):
    h_prim, h_ext = generateHeader(
        chip=chip,
        ra_cen=ra_cen,
        dec_cen=dec_cen,
        img_rot=img_rot,
        im_type=im_type,
        pointing_ID=pointing_ID,
        exptime=exptime,
        timestamp=timestamp)
    hdu1 = fits.PrimaryHDU(header=h_prim)
    hdu1.add_checksum()
    hdu1.header.comments['CHECKSUM'] = 'HDU checksum'
    hdu1.header.comments['DATASUM'] = 'data unit checksum'
    hdu2 = fits.ImageHDU(img.array, header=h_ext)
    hdu2.add_checksum()
    hdu2.header.comments['XTENSION'] = 'extension type'
    hdu2.header.comments['CHECKSUM'] = 'HDU checksum'
    hdu2.header.comments['DATASUM'] = 'data unit checksum'
    hdu1 = fits.HDUList([hdu1, hdu2])
    fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
    hdu1.writeto(fname, output_verify='ignore', overwrite=True)

def add_sky_background(img, filt, exptime, sky_map=None, tel=None):
    # Add sky background
    if sky_map is None:
        sky_map = filt.getSkyNoise(exptime=exptime)
        sky_map = sky_map * np.ones_like(img.array)
        sky_map = galsim.Image(array=sky_map)
        # Apply Poisson noise to the sky map
        # # (NOTE): only for photometric chips if it utilizes the photon shooting to draw stamps
        # if self.survey_type == "photometric":
        #     sky_map.addNoise(poisson_noise)
    elif img.array.shape != sky_map.shape:
        raise ValueError("The shape img and sky_map must be equal.")
    elif tel is not None: # If sky_map is given in flux
        sky_map = sky_map * tel.pupil_area * exptime
    img += sky_map
    return img, sky_map

def get_flat(img, seed):
    flat_img = effects.MakeFlatSmooth(
                GSBounds=img.bounds, 
                seed=seed)
    flat_normal = flat_img / np.mean(flat_img.array)
    return flat_img, flat_normal

def add_cosmic_rays(img, chip, exptime=150, seed=0):
    cr_map, cr_event_num = effects.produceCR_Map(
        xLen=chip.npix_x, yLen=chip.npix_y, 
        exTime=exptime+0.5*chip.readout_time, 
        cr_pixelRatio=0.003*(exptime+0.5*chip.readout_time)/600.,
        gain=chip.gain, 
        attachedSizes=chip.attachedSizes,
        seed=seed)   # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
    img += cr_map
    cr_map[cr_map > 65535] = 65535
    cr_map[cr_map < 0] = 0
    crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
    del cr_map
    return img, crmap_gsimg, cr_event_num

def add_PRNU(img, chip, seed=0):
    prnu_img = effects.PRNU_Img(
        xsize=chip.npix_x, 
        ysize=chip.npix_y, 
        sigma=0.01, 
        seed=seed)
    img *= prnu_img
    return img, prnu_img

def get_poisson(seed=0, sky_level=0.):
    rng_poisson = galsim.BaseDeviate(seed)
    poisson_noise = galsim.PoissonNoise(rng_poisson, sky_level=sky_level)
    return rng_poisson, poisson_noise

def get_base_img(img, read_noise, readout_time, dark_noise, exptime=150.):
    base_level = read_noise**2 + dark_noise*(exptime+0.5*readout_time)
    base_img = base_level * np.ones_like(img.array)
    return base_img

def add_poisson(img, chip, exptime=150., seed=0, sky_level=0., poisson_noise=None, dark_noise=None):
    if poisson_noise is None:
        _, poisson_noise = get_poisson(seed=seed, sky_level=sky_level)
    read_noise = chip.read_noise
    if dark_noise is None:
        dark_noise = chip.dark_noise
    base_img = get_base_img(img=img, read_noise=read_noise, readout_time=chip.readout_time, dark_noise=dark_noise, exptime=exptime)
    img += base_img
    img.addNoise(poisson_noise)
    img -= read_noise**2
    return img, base_img

def add_brighter_fatter(img):
    #Inital dynamic lib
    try:
        with pkg_resources.files('ObservationSim.Instrument.Chip.lib_bf').joinpath("libmoduleBF.so") as lib_path:
            lib_bf = ctypes.CDLL(lib_path)
    except AttributeError:
        with pkg_resources.path('ObservationSim.Instrument.Chip.lib_bf', "libmoduleBF.so") as lib_path:
            lib_bf = ctypes.CDLL(lib_path)
    lib_bf.addEffects.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_int]
    
    # Set bit flag
    bit_flag = 1
    bit_flag = bit_flag | (1 << 2)

    nx, ny = img.array.shape
    nn = nx * ny
    arr_ima= (ctypes.c_float*nn)()
    arr_imc= (ctypes.c_float*nn)()

    arr_ima[:]= img.array.reshape(nn)
    arr_imc[:]= np.zeros(nn)

    lib_bf.addEffects(nx, ny, arr_ima, arr_imc, bit_flag)
    img.array[:, :] = np.array(arr_imc[:]).reshape([nx, ny])
    del arr_ima, arr_imc
    return img