Commit 3f932427 authored by Xin Zhang's avatar Xin Zhang
Browse files

modify header

parent 0001f663

ObservationSim/Config/Header/ImageHeader.py

+137 −181
Original line number Diff line number Diff line
@@ -12,6 +12,8 @@ import random 


import os

import sys

import astropy.coordinates as coord

from astropy.time import Time



def chara2digit(char):

    """ Function to judge and convert characters to digitals
@@ -86,68 +88,68 @@ def rotate_CD_matrix(cd, pa_aper): 
    return cd_rot





def Header_extention(xlen = 9216, ylen = 9232, gain = 1.0, readout = 5.0, dark = 0.02,saturation=90000, row_num = 1, col_num = 1):



    """ Creat an image frame for CCST with multiple extensions



    Parameters

    ----------



    """



    flag_ltm_x = [0,1,-1,1,-1]

    flag_ltm_y = [0,1,1,-1,-1]

    flag_ltv_x = [0,0,1,0,1]

    flag_ltv_y = [0,0,0,1,1]



    detector_size_x = int(xlen)

    detector_size_y = int(ylen)



    data_x = str(int(detector_size_x))

    data_y = str(int(detector_size_y))



    data_sec = '[1:'+data_x+',1:'+data_y+']'

    e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/extension_header.param'

    name, value, description = read_header_parameter(e_header_fn)

    f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst')

    s = f.readline()

    s = s.strip("\n")

    filters = s.split(' ')

    s = f.readline()

    s = s.strip("\n")

    filterID = s.split()



    s = f.readline()

    s = s.strip("\n")

    CCDID = s.split()



    k = (row_num-1)*6+col_num



    h_iter = 0

    for n1,v1,d1 in zip(name, value, description):

        if n1=='EXTNAME':

            value[h_iter] = 'RAW,'+CCDID[k-1].rjust(2,'0')

        if n1=='CCDNAME':

            value[h_iter] = 'ccd' + CCDID[k-1].rjust(2,'0')

        if n1=='AMPNAME':

            value[h_iter] = 'ccd' + CCDID[k-1].rjust(2,'0') + ':A'

        if n1=='GAIN':

            value[h_iter] = gain

        if n1=='RDNOISE':

            value[h_iter] = readout

        if n1=='SATURATE':

            value[h_iter] = saturation

        if n1=='CCDCHIP':

            value[h_iter] = 'ccd' + CCDID[k-1].rjust(2,'0')

        if n1=='CCDLABEL':

            value[h_iter] = filters[k-1] + '-' + filterID[k-1]

        if n1=='DATASEC':

            value[h_iter] = data_sec



        h_iter = h_iter + 1



   

    return name, value, description

# def Header_extention(xlen = 9216, ylen = 9232, gain = 1.0, readout = 5.0, dark = 0.02,saturation=90000, row_num = 1, col_num = 1):

#

#     """ Creat an image frame for CCST with multiple extensions

#

#     Parameters

#     ----------

#

#     """

#

#     flag_ltm_x = [0,1,-1,1,-1]

#     flag_ltm_y = [0,1,1,-1,-1]

#     flag_ltv_x = [0,0,1,0,1]

#     flag_ltv_y = [0,0,0,1,1]

#

#     detector_size_x = int(xlen)

#     detector_size_y = int(ylen)

#

#     data_x = str(int(detector_size_x))

#     data_y = str(int(detector_size_y))

#

#     data_sec = '[1:'+data_x+',1:'+data_y+']'

#     e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/extension_header.param'

#     name, value, description = read_header_parameter(e_header_fn)

#     f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst')

#     s = f.readline()

#     s = s.strip("\n")

#     filters = s.split(' ')

#     s = f.readline()

#     s = s.strip("\n")

#     filterID = s.split()

#

#     s = f.readline()

#     s = s.strip("\n")

#     CCDID = s.split()

#

#     k = (row_num-1)*6+col_num

#

#     h_iter = 0

#     for n1,v1,d1 in zip(name, value, description):

#         if n1=='EXTNAME':

#             value[h_iter] = 'RAW,'+CCDID[k-1].rjust(2,'0')

#         if n1=='CCDNAME':

#             value[h_iter] = 'ccd' + CCDID[k-1].rjust(2,'0')

#         if n1=='AMPNAME':

#             value[h_iter] = 'ccd' + CCDID[k-1].rjust(2,'0') + ':A'

#         if n1=='GAIN':

#             value[h_iter] = gain

#         if n1=='RDNOISE':

#             value[h_iter] = readout

#         if n1=='SATURATE':

#             value[h_iter] = saturation

#         if n1=='CCDCHIP':

#             value[h_iter] = 'ccd' + CCDID[k-1].rjust(2,'0')

#         if n1=='CCDLABEL':

#             value[h_iter] = filters[k-1] + '-' + filterID[k-1]

#         if n1=='DATASEC':

#             value[h_iter] = data_sec

#

#         h_iter = h_iter + 1

#

#

#     return name, value, description





##9232 9216  898 534 1309 60 -40  -23.4333
@@ -229,48 +231,9 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r 
        return r_dat





    #     name = ['EQUINOX',

    #             'WCSDIM',

    #             'CTYPE1',

    #             'CTYPE2',

    #             'CRVAL1',

    #             'CRVAL2',

    #             'CRPIX1',

    #             'CRPIX2',

    #             'CD1_1',

    #             'CD1_2',

    #             'CD2_1',

    #             'CD2_2']

    #     value = [2000.0,

    #              2.0,

    #              'RA---TAN',

    #              'DEC--TAN',

    #              ra_ref,

    #              dec_ref,

    #              flag_ext_x[k]*((x_ref+flag_ext_x[k]*detector_size_x/2)-x_center),

    #              flag_ext_y[k]*((y_ref+flag_ext_y[k]*detector_size_y/2)-y_center),

    #              cd_rot[0,0],

    #              cd_rot[0,1],

    #              cd_rot[1,0],

    #              cd_rot[1,1]]

    #     description = ['Equinox of WCS',

    #                    'WCS Dimensionality',

    #                    'Coordinate type',

    #                    'Coordinate typ',

    #                    'Coordinate reference value',

    #                    'Coordinate reference value',

    #                    'Coordinate reference pixel',

    #                    'Coordinate reference pixel',

    #                    'Coordinate matrix',

    #                    'Coordinate matrix',

    #                    'Coordinate matrix',

    #                    'Coordinate matrix']



    # return name, value, description







def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, psize = 0.074, row_num = 1, col_num = 1, date='200930', time_obs='120000'):

def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, psize = 0.074, row_num = 1, col_num = 1, date='200930', time_obs='120000', im_type = 'MS'):



    # array_size1, array_size2, flux, sigma = int(argv[1]), int(argv[2]), 1000.0, 5.0
@@ -278,8 +241,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec 
    k = (row_num-1)*6+col_num

    # ccdnum = str(k)



    g_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/global_header.param'

    name, value, description = read_header_parameter(g_header_fn)

    g_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/global_header.header'

    f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst')

    s = f.readline()

    s = s.strip("\n")
@@ -293,87 +255,81 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec 
    CCDID = s.split()



    h_prim = fits.Header()

    h_prim = fits.Header.fromfile(g_header_fn)



    h_prim['PIXSIZE1'] = xlen

    h_prim['PIXSIZE2'] = ylen

    h_prim['FILENAME'] = 'MSC_' + im_type + '_' + date + time_obs + '_1' + pointNum.rjust(8,'0') + '_' +CCDID[k-1].rjust(2,'0')+'_raw'

    h_prim['DATE'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6]

    h_prim['TIME'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]

    h_prim['DETECTOR'] = 'CCD'+CCDID[k-1].rjust(2,'0')

    h_prim['RA_OBJ'] = ra

    h_prim['DEC_OBJ'] = dec

    h_prim['OBJECT'] = '1'+ pointNum.rjust(8,'0')

    h_prim['OBSID'] = '1'+ pointNum.rjust(8,'0')

    h_prim['TELFOCUS'] = 'f/14'



    co = coord.SkyCoord(ra, dec, unit='deg')



    ra_hms = format(co.ra.hms.h, '02.0f') + ':' + format(co.ra.hms.m, '02.0f') + ':' + format(co.ra.hms.s, '05.2f')

    dec_hms = format(co.dec.dms.d, '02.0f') + ':' + format(abs(co.dec.dms.m), '02.0f') + ':' + format(abs(co.dec.dms.s),

                                                                                                      '05.2f')



    h_prim['RA_NOM'] = ra_hms

    h_prim['DEC_NOM'] = dec_hms

    h_prim['PIXSCAL1'] = str(psize)

    h_prim['PIXSCAL2'] = str(psize)



    ttt = h_prim['DATE'] + 'T' + h_prim['TIME']

    tstart = Time(ttt)

    h_prim['EXPSTART'] = round(tstart.mjd, 5)

    tend = Time(tstart.cxcsec + h_prim['EXPTIME'], format="cxcsec")

    h_prim['EXPEND'] = round(tend.mjd, 5)



    for i in range(len(name)):

        

        if(name[i]=='FILENAME'):

            # value[i] = 'CSST_' + date + '_' +time_obs + '_' + pointNum.rjust(6,'0') + '_' +ccdnum.rjust(2,'0')+'_raw'

            # value[i] = 'CSST_' + date + '_' +time_obs + '_' + pointNum.rjust(6,'0') + '_' +CCDID[k-1].rjust(2,'0')+'_raw'

            value[i] = 'MSC_' + date + time_obs + '_' + pointNum.rjust(7,'0') + '_' +CCDID[k-1].rjust(2,'0')+'_raw'

        

        if(name[i]=='DETSIZE'):

            value[i] = '[1:' + str(int(xlen)) + ',1:'+  str(int(ylen)) + ']'



        if(name[i]=='PIXSCAL1'):

            value[i] = str(psize)



        if(name[i]=='PIXSCAL2'):

            value[i] = str(psize)

        



        h_prim[name[i]] = (value[i],description[i])



    h_prim.add_comment('==================================================================',after='FILETYPE')

    h_prim.add_comment('Target information')

    h_prim.add_comment('==================================================================')



    h_prim.add_comment('==================================================================',after='EQUINOX')

    h_prim.add_comment('Exposure information')

    h_prim.add_comment('==================================================================')



    h_prim.add_comment('==================================================================',after='MJDEND')

    h_prim.add_comment('Telescope information')

    h_prim.add_comment('==================================================================')



    h_prim.add_comment('==================================================================',after='REFFRAME')

    h_prim.add_comment('Detector information')

    h_prim.add_comment('==================================================================')

    return h_prim



    h_prim.add_comment('==================================================================',after='DETETEMP')

    h_prim.add_comment('Other information')

    h_prim.add_comment('==================================================================')

def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, psize = 0.074, row_num = 1, col_num = 1, extName='SCI'):



    return h_prim

    e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/extension_header.header'

    f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst')

    s = f.readline()

    s = s.strip("\n")

    filters = s.split(' ')

    s = f.readline()

    s = s.strip("\n")

    filterID = s.split()



def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, psize = 0.074, row_num = 1, col_num = 1):

    s = f.readline()

    s = s.strip("\n")

    CCDID = s.split()



    h_ext = fits.Header()

    k = (row_num - 1) * 6 + col_num



    for i in range(1,2):

    h_ext = fits.Header.fromfile(e_header_fn)



        # NAXIS1:Number of pixels per row;  NAXIS2:Number of rows

    h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1]

    h_ext['FILTER'] = filters[k-1]

    h_ext['NAXIS1'] = xlen

    h_ext['NAXIS2'] = ylen

        name, value, description = Header_extention(xlen = xlen, ylen = ylen, gain = gain, readout = readout, dark = dark, saturation=saturation, row_num = row_num, col_num = col_num)



        for j in range(len(name)):

            # print(name[j],value[j],description[j])

            h_ext[name[j]] = (value[j],description[j])



        header_wcs = WCS_def(xlen = xlen, ylen = ylen, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = ra, dec = dec, pa = pa ,psize = psize, row_num = row_num, col_num = col_num)



        for k,v in header_wcs.items():

            h_ext.set(k,v)





        # for j in range(len(name)):

        #     h_ext[name[j]] = (value[j],description[j])



        h_ext.add_comment('==================================================================',after='OBSID')

        h_ext.add_comment('Readout information')

        h_ext.add_comment('==================================================================')



        h_ext.add_comment('==================================================================',after='CHIPTEMP')

        h_ext.add_comment('Chip information')

        h_ext.add_comment('==================================================================')



        h_ext.add_comment('==================================================================',after='TRIMSEC')

        h_ext.add_comment('WCS information')

        h_ext.add_comment('==================================================================')



        h_ext.add_comment('==================================================================',after='CD2_2')

        h_ext.add_comment('Other information')

        h_ext.add_comment('==================================================================')

    h_ext['EXTNAME'] = extName

    h_ext['GAIN1'] = gain

    h_ext['RDNOISE1'] = readout

    h_ext['CCDCHIP'] =  'ccd' + CCDID[k-1].rjust(2,'0')

    header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra=ra, dec=dec, pa=pa, psize=psize,

                         row_num=row_num, col_num=col_num)



    h_ext['CRPIX1'] = header_wcs['CRPIX1']

    h_ext['CRPIX2'] = header_wcs['CRPIX2']

    h_ext['CRVAL1'] = header_wcs['CRVAL1']

    h_ext['CRVAL2'] = header_wcs['CRVAL2']

    h_ext['CD1_1'] = header_wcs['CD1_1']

    h_ext['CD1_2'] = header_wcs['CD1_2']

    h_ext['CD2_1'] = header_wcs['CD2_1']

    h_ext['CD2_2'] = header_wcs['CD2_2']

    h_ext['EQUINOX'] = header_wcs['EQUINOX']

    h_ext['WCSDIM'] = header_wcs['WCSDIM']

    h_ext['CTYPE1'] = header_wcs['CTYPE1']

    h_ext['CTYPE2'] = header_wcs['CTYPE2']



    return h_ext

ObservationSim/Config/Header/extension_header.header

0 → 100644
+1 −0
Original line number Diff line number Diff line 
XTENSION= 'IMAGE   '           / Image extension                                BITPIX  =                   16 / array data type                                NAXIS   =                    2 / number of array dimensions                     NAXIS1  =                 9216 / length of first array axis                     NAXIS2  =                 9232 / length of second array axis                    PCOUNT  =                    0                                                  GCOUNT  =                    1                                                  EXTNAME = 'SCI     '                                                            EXTVER  =                    1                                                  BSCALE  =                    1                                                  BZERO   =                32768                                                  BUNIT   = 'ADU     '                                                            COMMENT ==================================================================      COMMENT CCD chip information                                                    COMMENT ==================================================================      NCHAN   =                   16 / Number of readout channels                     NCHAN1  =                    8 / Number of horizintal channels                  NCHAN2  =                    2 / Number of virtical channels                    DETSIZE = '[1:9216:9232]'      / detector size                                  BIASSEC = '[9216:9217,9232:9234]' / bias section                                CCDCHIP =                    6 / CCD chip ID                                    CCDLABEL= 'y-1     '           / CCD chip label                                 PSCAN1  =                   27 / horizontal prescan width                       PSCAN2  =                    8 / virtical prescan height                        OSCAN1  =                   16 / horizontal overscan width                      OSCAN2  =                    8 / virtical overscan height                       FILTER  = 'y       '           / filter name                                    COMMENT ==================================================================      COMMENT WORLD COORDINATE SYSTEM AND RELATED PARAMETERS                          COMMENT ==================================================================      WCSDIM  =                    2 / Number of World Coordinate System axes         EQUINOX =               2000.0 / Epoch (year)                                   CRPIX1  =             -10017.0 / Coordinate reference pixel of x                CRPIX2  =              24876.0 / Coordinate reference pixel of y                CRVAL1  =            62.228226 / Coordinate reference value of x                CRVAL2  =           -42.316932 / Coordinate reference value of y                CTYPE1  = 'RA---TAN'           / the coordinate type                            CTYPE2  = 'DEC--TAN'           / the coordinate type                            CD1_1   = 1.88602083707394E-05 / partial of first axis coordinate of x          CD2_1   = -8.1745583617600E-06 / partial of first axis coordinate of y          CD1_2   = 8.17455836176000E-06 / partial of second axis coordinate of x         CD2_2   = 1.88602083707394E-05 / partial of second axis coordinate of y         POS_ANG =                 30.0 / Positon angle of detector array                COMMENT ==================================================================      COMMENT Readout information                                                     COMMENT ==================================================================      GAIN1   =                  1.1 / CCD gain                                       RDNOISE1=                  5.0 / read noise                                     READTIME=                 40.0 / read time (sec)                                RDSPEED =                 10.0 / read speed (in MHz)                            CHIPTEMP=               -100.0 / chip temperature (in K)                        COMMENT ==================================================================      COMMENT Other information                                                       COMMENT ==================================================================      CHECKSUM= '''abcde'''          / SHA256 checksum of global headers              END

ObservationSim/Config/Header/global_header.header

0 → 100644
+1 −0

File added.

Preview size limit exceeded, changes collapsed.

ObservationSim/ObservationSim.py

+4 −2
Original line number Diff line number Diff line
@@ -212,7 +212,8 @@ class Observation(object): 
            row_num=chip.rowID, 

            col_num=chip.colID,

            date=self.config["date_obs"],

            time_obs=self.config["time_obs"])

            time_obs=self.config["time_obs"],

            im_type = 'MS')

        h_ext = generateExtensionHeader(

            xlen=chip.npix_x, 

            ylen=chip.npix_y,
@@ -225,7 +226,8 @@ class Observation(object): 
            saturation=90000, 

            psize=chip.pix_scale, 

            row_num=chip.rowID, 

            col_num=chip.colID)

            col_num=chip.colID,

            extName='SCI')

        chip.img = galsim.Image(chip.img.array, dtype=np.uint16)

        hdu1 = fits.PrimaryHDU(header=h_prim)

        hdu2 = fits.ImageHDU(chip.img.array, header=h_ext)

test/ObservationSim.cfg

+1 −1
Original line number Diff line number Diff line
@@ -10,7 +10,7 @@ time_obs 093000 # Observation time [hhmmss] 
ra_center     60.0                    # Telesscope pointing center [degree] (Default)

dec_center   -40.0                     

psf_rcont     0.15,0.8                # Radius of the 80% flux concentration (for a Gaussian PSF)

image_rot     23.4333                 # Image orientation [degree]

image_rot    -113.4333                # Image orientation [degree]

sigma_spin    0.0                     # Image spin [degree]

rotateEll     0.0                     # Rotate the intrinsic ellipticity by a given angle  

                                      # (currently only support 0.0, 45.0, 90.0, 135.0 [degree])