Commit 72c1f3f5 authored by GZhao's avatar GZhao
Browse files

toml data model; add dataset keyword; update fits header

parent 66d8b23f
Loading
Loading
Loading
Loading
Loading
+20 −10
Original line number Diff line number Diff line
@@ -474,7 +474,6 @@ class CpicVisEmccd(object):
        np.ndarray
            bias frame
        """

        shape = self.bias_shape
        TPI = np.pi * 2

@@ -681,15 +680,17 @@ class CpicVisEmccd(object):

        self.volt = volt_func(em_set)

        volt_3ff = volt_func(int('3ff', 16))
        volt_190 = volt_func(int('190', 16))
        # volt_3ff = volt_func(int('3ff', 16))
        # volt_190 = volt_func(int('190', 16))

        em_coe_c = 0.24
        em_coe_c = 0.2465
        # using the expression of em_b = ln(g199) / constant to make fitting easier
        constant = (np.exp(em_coe_c * volt_190) - np.exp(em_coe_c * volt_3ff))
        # fitting from the ccd test result
        ln_g190 = (-ccd_temp - 7) * 0.0325
        em_coe_b = ln_g190 / constant
        # constant = (np.exp(em_coe_c * volt_190) - np.exp(em_coe_c * volt_3ff))
        # # fitting from the ccd test result
        # ln_g190 = (-ccd_temp - 7) * 0.0325
        # em_coe_b = ln_g190 / constant
        em_coe_b = -0.31041118 * ccd_temp + 8.61774308
        em_coe_b = em_coe_b * 1e-5

        emgain = np.exp(em_coe_b * np.exp(em_coe_c * self.volt))
        emgain = np.maximum(1, emgain)
@@ -759,7 +760,7 @@ class CpicVisEmccd(object):
    #     img_line = np.convolve(image.flatten(), self.cti_trail, mode='full')
    #     return img_line[:shape[0]*shape[1]].reshape(shape)

    def readout(self, image_focal, em_set, expt_set, image_cosmic_ray=False, emgain=None):
    def readout(self, image_focal, em_set, expt_set, image_cosmic_ray=False, emgain=None, seed=-1):
        """From focal planet image to ccd output.
        Interface function for emccd. Simulate the readout process.

@@ -775,12 +776,16 @@ class CpicVisEmccd(object):
            cosmic ray image. Unit: photon-electron/pixel
        emgain: float, optional
            if not None, use the given emgain. Else, calculate the emgain using em_set
        seed: int, optional
            random seed. Default is -1, which means no random seed.

        Returns
        -------
        np.ndarray
            ccd output image. Unit: ADU
        """
        if seed != -1:
            np.random.seed(seed)
            
        expt = expt_set
        if expt_set == 0:
@@ -844,6 +849,7 @@ class CpicVisEmccd(object):
            image_shutter = np.random.poisson(image_shutter)
            image[:, self.pscan1+self.ldark:-self.oscan1-self.rdark] += image_shutter


        if self.switch['cic']:
            cic_frame = np.zeros((self.dark_shape[0], self.bias_shape[1])) + self.cic
            image[self.pscan2:-self.oscan2, :] += np.random.poisson(cic_frame)
@@ -855,9 +861,12 @@ class CpicVisEmccd(object):
        # >>>     image += np.random.binomial(image, pEM)
        # This code is too slow, so we used a modified gamma



        em_fix = self.em_fix_fuc_fit(emgain) * emgain
        image = np.random.gamma(image, em_fix) + image * (emgain - em_fix)
        

        if self.switch['em_blooming']:
            image = self.emregester_blooming(image)

@@ -870,6 +879,7 @@ class CpicVisEmccd(object):
                residual = np.random.binomial(image, 1-big_cte)
                image[:, 1:] += residual[:, :-1] - residual[:, 1:]


        bias = self.bias_frame()

        image = image / self.ph_per_adu + bias
+2 −3
Original line number Diff line number Diff line
@@ -81,7 +81,7 @@ config['bands'] = {
}
config['diameter'] = 2  # in meters
config['platescale'] = 0.016153
config['datamodel'] = f'{cpism_refdata}/io/csst-cpic-l0.yaml'
config['datamodel'] = f'{cpism_refdata}/io/csst-cpic-l0.toml'

config['log_dir'] = f'{cpism_refdata}/log'
config['log_level'] = 'info'
@@ -91,7 +91,7 @@ config['dm_pickle'] = f'{cpism_refdata}/optics/dm_model.pkl'
config['pysyn_refdata'] = f'{cpism_refdata}/starmodel/grp/redcat/trds'
config['catalog_folder'] = f'{cpism_refdata}/demo_catalog'
config['csst_format'] = True
config['nsample'] = 5
config['nsample'] = 1

update_able_keys = [
    'apm_file', 'actuator_file', 'aberration',
@@ -99,7 +99,6 @@ update_able_keys = [
    'nsample', 'csst_format', 'output', 'check_fits_header'
]


def replace_cpism_refdata(
        config: dict,
        output: str = '$') -> None:
+36 −10
Original line number Diff line number Diff line
import yaml
import os
from datetime import datetime
import numpy as np
@@ -6,6 +5,7 @@ import pandas as pd
from astropy.io import fits
from astropy.coordinates import SkyCoord
import astropy.units as u
import toml

from .config import __version__, which_focalplane
from .utils import Logger
@@ -44,21 +44,25 @@ def check_and_update_fits_header(header):
    # if hdu == 'primary':
    #     model_file = f"{cpism_refdata}/io/primary_header.json"
    model_file = config['datamodel']
    with open(model_file, 'r', encoding='utf-8') as fid:
        data_model = yaml.load(fid, Loader=yaml.FullLoader)
    data_model = toml.load(model_file)

    if 'FILETYPE' in header.keys():
    if 'SIMPLE' in header.keys():
        header_model = data_model['HDU0']
        hdu = 'hdu0'
    else:
        header_model = data_model['HDU1']
        header_model = data_model['HDUN']
        hdu = 'hdu1'
    dm_comment = {}


    def print_warning(info):
        if header_check:
            log.warning(info)

    has_all_keywords = True
    has_unexpected_keywords = True
    has_correct_order = True

    # check existance and format of keyword in fits header
    for _, content in header_model.items():
        comment = content['comment']
@@ -76,6 +80,7 @@ def check_and_update_fits_header(header):

        if keyword not in header.keys():
            print_warning(f"Keyword {keyword} not found in header.")
            has_all_keywords = False

        elif not pd.isnull(dtype):
            value = header[keyword]
@@ -108,9 +113,28 @@ def check_and_update_fits_header(header):
        if keyword not in dm_comment.keys():
            print_warning(
                f"Keyword {keyword} not found in the [{hdu}] data model.")
            has_unexpected_keywords = False
        elif keyword != 'COMMENT':  # comment keyword is not allowed to be updated
            header[keyword] = (header[keyword], dm_comment[keyword])

    if has_unexpected_keywords and has_all_keywords:  # if there is no other keywords, then use the default header model
        for header_key, model_key in zip(header.keys(),dm_comment.keys()):
            if header_key != model_key:
                print_warning(
                    f"Keyword {header_key} not correspond to {model_key}")
                has_correct_order = False
    else:
        has_correct_order = False

    overview = "data model check: [OVERVIEW]\n"
    if not has_correct_order:
        overview += "The header is not in the correct order\n"
    if not has_all_keywords:
        overview += "The header is missing some keywords\n"
    if not has_unexpected_keywords:
        overview += "The header contains some unexpected keywords\n"
    print_warning(overview)

    return header


@@ -123,7 +147,7 @@ def obsid_parser(
    ----------
    obsid: str
        The obsid of the observation.
        Obsid must be 11 digits and start with 5 for CPIC.
        Obsid must be 11 digits and start with 4 for CPIC.

    Returns
    -------
@@ -133,7 +157,7 @@ def obsid_parser(
    Raises
    ------
    ValueError
        If the obsid is not 11 digits or does not start with 5.
        If the obsid is not 11 digits or does not start with 4.
    """
    obsid = str(obsid)
    if len(obsid) != 11:
@@ -261,9 +285,11 @@ def primary_hdu(
    if len(heaer_filename) > 68:
        heaer_filename = heaer_filename[:68]
    header['FILENAME'] = heaer_filename
    header['FILETYPE'] = obsid_parser(obsid)
    header['TELESCOP'] = 'CSST'
    header['INSTRUME'] = 'CPIC'
    header['OBSID'] = str(obsid)
    header['OBSTYPE'] = obsid_parser(obsid)
    header['DATASET'] = obs_info.get("dataset", "default")
    header['RADECSYS'] = 'ICRS'
    header['EQUINOX'] = 2000.0
    header['FITSSWV'] = f'CPISM V{__version__}'
@@ -280,7 +306,7 @@ def primary_hdu(

    header['OBJECT'] = cstar.get('name', target_name)
    header['TARGET'] = target_name
    header['OBSID'] = str(obsid)

    header['RA_OBJ'] = radec.ra.degree
    header['DEC_OBJ'] = radec.dec.degree

@@ -404,9 +430,9 @@ def frame_header(obs_info, index, primary_header, camera_dict):
    header['DETSN'] = '00000000000'
    header['DETNAME'] = camera_config['detector_name']
    header['CHIPLAB'] = camera_config['ccd_label']
    header['DEWTEMP'] = float(camera_config['cooler_temp']) + 273.15
    frame_info = obs_info['frame_info'][index]
    header['CHIPTEMP'] = float(frame_info['chiptemp']) + 273.15
    header['DEWTEMP'] = float(camera_config['cooler_temp']) + 273.15
    header['DETSIZE'] = f"{imgszx} * {imgszy}"
    header['IMGINDEX'] = index + 1

+6 −2
Original line number Diff line number Diff line
@@ -38,7 +38,8 @@ def vis_observation(
        nsample: int = 1,
        emgain=None,
        prograss_bar=None,
        output='./'
        output='./',
        dataset='default'
) -> np.ndarray:
    """Observation simulation main process on visable band using EMCCD.

@@ -110,6 +111,7 @@ def vis_observation(
        'obsid': obsid,
        'rotation': rotation,
        'shift': shift,
        'dataset': dataset
    }

    all_frame_info = []
@@ -262,6 +264,7 @@ def quick_run_v2(
        emgain_value = None

    camera = CpicVisEmccd()
    camera.switch['bias_ci'] = False
    if not camera_effect:
        camera.switch['bias_vp'] = False
        camera.switch['bias_hp'] = False
@@ -288,7 +291,8 @@ def quick_run_v2(
        rotation=rotation,
        camera=camera,
        prograss_bar=prograss_bar,
        output=output
        output=output,
        nsample=10
    )


+1 −1
Original line number Diff line number Diff line
band: f662
band: f661
emgain: 200
emset: -1
expt: 100
Loading