ImageHeader.py

"""
generate image header
"""
import numpy as np

from astropy.io import fits
import astropy.wcs as pywcs
from collections import OrderedDict

# from scipy import math
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

    Parameters
    ----------

    """

    try:
        float(char) # for int, long and float
    except ValueError:
        pass
        return char
    else:
        data = float(char)
        return data


def read_header_parameter(filename='global_header.param'):
    """ Function to read the header parameters

    Parameters
    ----------

    """

    name = []
    value = []
    description = []
    for line in open(filename):
        line = line.strip("\n")
        arr = line.split('|')
#        csvReader = csv.reader(csvDataFile)
#        for arr in csvReader:
        name.append(arr[0])
        # print(arr[0],arr[1])
        value.append(chara2digit(arr[1]))
        description.append(arr[2])

#    print(value)
    return name, value, description

def rotate_CD_matrix(cd, pa_aper):
    """Rotate CD matrix
    
    Parameters
    ----------
    cd: (2,2) array
        CD matrix
    
    pa_aper: float
        Position angle, in degrees E from N, of y axis of the detector
    
    Returns
    -------
    cd_rot: (2,2) array
        Rotated CD matrix
    
    Comments
    --------
    `astropy.wcs.WCS.rotateCD` doesn't work for non-square pixels in that it
    doesn't preserve the pixel scale!  The bug seems to come from the fact
    that `rotateCD` assumes a transposed version of its own CD matrix.
    
    """
    rad = np.deg2rad(-pa_aper)
    mat = np.zeros((2,2))
    mat[0,:] = np.array([np.cos(rad),-np.sin(rad)])
    mat[1,:] = np.array([np.sin(rad),np.cos(rad)])
    cd_rot = np.dot(mat, cd)
    return cd_rot

def calcaluteSLSRotSkyCoor(pix_xy = None,rot_angle = 1, xlen = 9216, ylen = 9232, w = None):
    rad = np.deg2rad(rot_angle)
    mat = np.zeros((2,2))
    mat[0,:] = np.array([np.cos(rad),-np.sin(rad)])
    mat[1,:] = np.array([np.sin(rad),np.cos(rad)])
    center = np.array([xlen/2, ylen/2])
    rot_pix = np.dot(mat, pix_xy-center) + center
    skyCoor = w.wcs_pix2world(np.array([rot_pix]), 1)

    return skyCoor


# 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
def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1, filter = 'GI'):

    """ Creat a wcs frame for CCST with multiple extensions

    Parameters
    ----------

    """

    flag_x = [0, 1, -1, 1, -1]
    flag_y = [0, 1, 1, -1, -1]
    flag_ext_x = [0,-1,1,-1,1]
    flag_ext_y = [0,-1,-1,1,1]
    x_num = 6
    y_num = 5
    detector_num = x_num*y_num


    detector_size_x = xlen
    detector_size_y = ylen
    gap_y = gapy
    gap_x = [gapx1,gapx2]
    ra_ref = ra
    dec_ref = dec

    pa_aper = pa

    pixel_size = psize

    gap_x1_num = 3
    gap_x2_num = 2

    y_center = (detector_size_y*y_num+gap_y*(y_num-1))/2

    x_center = (detector_size_x*x_num+gap_x[0]*gap_x1_num+gap_x[1]*gap_x2_num)/2

    gap_x_map = np.array([[0,0,0,0,0],[gap_x[0],gap_x[1],gap_x[1],gap_x[1],gap_x[1]],[gap_x[1],gap_x[0],gap_x[0],gap_x[0],gap_x[0]],[gap_x[0],gap_x[0],gap_x[0],gap_x[0],gap_x[0]],[gap_x[0],gap_x[0],gap_x[0],gap_x[0],gap_x[1]],[gap_x[1],gap_x[1],gap_x[1],gap_x[1],gap_x[0]]])


    # frame_array = np.empty((5,6),dtype=np.float64)
    # print(x_center,y_center)
    
    j = row_num
    i = col_num
    # ccdnum = str((j-1)*5+i)

    x_ref, y_ref = detector_size_x*i + sum(gap_x_map[0:i,j-1]) - detector_size_x/2. , (detector_size_y+gap_y)*j-gap_y-detector_size_y/2

    # print(i,j,x_ref,y_ref,ra_ref,dec_ref)

    r_dat = OrderedDict()

    # name = []
    # value = []
    # description = []
    
    for k in range(1,2):
        
        cd = np.array([[ pixel_size,  0], [0, pixel_size]])/3600.*flag_x[k]
        cd_rot = rotate_CD_matrix(cd, pa_aper)

        # f = open("CCD"+ccdnum.rjust(2,'0')+"_extension"+str(k)+"_wcs.param","w")

        r_dat['EQUINOX'] = 2000.0
        r_dat['WCSDIM'] = 2.0
        r_dat['CTYPE1'] = 'RA---TAN'
        r_dat['CTYPE2'] = 'DEC--TAN'
        r_dat['CRVAL1'] = ra_ref
        r_dat['CRVAL2'] = dec_ref
        r_dat['CRPIX1'] = flag_ext_x[k]*((x_ref+flag_ext_x[k]*detector_size_x/2)-x_center)
        r_dat['CRPIX2'] = flag_ext_y[k]*((y_ref+flag_ext_y[k]*detector_size_y/2)-y_center)
        r_dat['CD1_1'] = cd_rot[0,0]
        r_dat['CD1_2'] = cd_rot[0,1]
        r_dat['CD2_1'] = cd_rot[1,0]
        r_dat['CD2_2'] = cd_rot[1,1]

        if filter in ['GU', 'GV', 'GI']:
            from astropy import wcs

            w = wcs.WCS(naxis=2)
            w.wcs.crpix = [r_dat['CRPIX1'], r_dat['CRPIX2']]
            w.wcs.cd = cd_rot
            w.wcs.crval = [ra_ref, dec_ref]
            w.wcs.ctype = [r_dat['CTYPE1'], r_dat['CTYPE2']]

            # test_center_o = w.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1)

            sls_rot = 1
            if i > 2:
                sls_rot = -sls_rot

            sn_x = 30
            sn_y = 30
            x_pixs = np.zeros(sn_y * sn_x)
            y_pixs = np.zeros(sn_y * sn_x)
            xpixs_line = np.linspace(1, xlen, sn_x)
            ypixs_line = np.linspace(1, ylen, sn_y)

            sky_coors = []

            for n1, y in enumerate(ypixs_line):
                for n2, x in enumerate(xpixs_line):
                    i_pix = n1 * sn_x + n2
                    x_pixs[i_pix] = x
                    y_pixs[i_pix] = y

                    pix_coor = np.array([x, y])
                    sc1 = calcaluteSLSRotSkyCoor(pix_xy=pix_coor, rot_angle=sls_rot, w=w)
                    # print(sc1[0,0],sc1[0,1])
                    sky_coors.append((sc1[0, 0], sc1[0, 1]))

            from astropy.coordinates import SkyCoord
            from astropy.wcs.utils import fit_wcs_from_points

            wcs_new = fit_wcs_from_points(xy=np.array([x_pixs, y_pixs]),
                                          world_coords=SkyCoord(sky_coors, frame="icrs", unit="deg"), projection='TAN')

            # print(wcs_new)
            # test_center = wcs_new.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1)
            #
            # print(test_center - test_center_o)

            r_dat['CD1_1'] = wcs_new.wcs.cd[0, 0]
            r_dat['CD1_2'] = wcs_new.wcs.cd[0, 1]
            r_dat['CD2_1'] = wcs_new.wcs.cd[1, 0]
            r_dat['CD2_2'] = wcs_new.wcs.cd[1, 1]
            r_dat['CRPIX1'] = wcs_new.wcs.crpix[0]
            r_dat['CRPIX2'] = wcs_new.wcs.crpix[1]

            r_dat['CRVAL1'] = wcs_new.wcs.crval[0]
            r_dat['CRVAL2'] = wcs_new.wcs.crval[1]

        return r_dat



#TODO project_cycle is temporary， is not in header defined, delete in future
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', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.], project_cycle=6):

    # array_size1, array_size2, flux, sigma = int(argv[1]), int(argv[2]), 1000.0, 5.0

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

    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")
    filters = s.split(' ')
    s = f.readline()
    s = s.strip("\n")
    filterID = s.split()

    s = f.readline()
    s = s.strip("\n")
    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['DATE'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6] + 'T' + time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
    # h_prim['TIME'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
    h_prim['DATE-OBS'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6] + 'T' + time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
    # h_prim['TIME-OBS'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
    # h_prim['DETECTOR'] = 'CHIP'+CCDID[k-1].rjust(2,'0')
    h_prim['OBJ_RA'] = ra
    h_prim['OBJ_DEC'] = dec
    h_prim['OBJECT'] = '1'+ str(int(project_cycle)) + pointNum.rjust(7,'0')
    h_prim['OBSID'] = '1'+ str(int(project_cycle)) + pointNum.rjust(7,'0')
    # h_prim['TELFOCUS'] = 'f/14'
    h_prim['EXPTIME'] = exptime
    
    # Define file types
    file_type = {'SCI':'sci', 'BIAS':'bias', 'DARK':'dark', 'FLAT':'flat', 'CRS':'cosmic_ray', 'CRD':'cosmic_ray'}
    h_prim['OBSTYPE'] = file_type[im_type]

    # 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['RA_PNT0'] = ra
    h_prim['DEC_PNT0'] = dec
    h_prim['RA_PNT1'] = ra
    h_prim['DEC_PNT1'] = dec



    # h_prim['PIXSCAL1'] = psize
    # h_prim['PIXSCAL2'] = psize

    ttt = h_prim['DATE']
    tstart = Time(ttt)
    h_prim['EXPSTART'] = round(tstart.mjd, 5)
    h_prim['CABSTART'] = h_prim['EXPSTART']
    # tend = Time(tstart.cxcsec + h_prim['EXPTIME'], format="cxcsec")
    tend = Time(tstart.mjd + h_prim['EXPTIME']/86400., format="mjd")
    h_prim['EXPEND'] = round(tend.mjd, 5)
    h_prim['CABEND'] = h_prim['EXPEND']

    file_start_time = '20' + date[0:6] + time_obs[0:6]
    end_time_str = str(tend.datetime)
    file_end_time = end_time_str[0:4] + end_time_str[5:7]+end_time_str[8:10] + end_time_str[11:13] + end_time_str[14:16] + end_time_str[17:19]
    h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_1' + pointNum.rjust(8, '0') + '_' + CCDID[
        k - 1].rjust(2, '0') + '_L0_1'


    h_prim['POSI0_X'] = sat_pos[0]
    h_prim['POSI0_Y'] = sat_pos[1]
    h_prim['POSI0_Z'] = sat_pos[2]

    h_prim['VELO0_X'] = sat_vel[0]
    h_prim['VELO0_Y'] = sat_vel[1]
    h_prim['VELO0_Z'] = sat_vel[2]
    # h_prim['RA_PNT0'] = ra_hms
    # h_prim['DEC_PNT0'] = dec_hms
    
    # Get version of CSSTSim Package
    from pkg_resources import get_distribution
    # h_prim['SIM_VER'] = (get_distribution("CSSTSim").version, "Version of CSST MSC simulation software")
    h_prim['FITSCREA'] = get_distribution("CSSTSim").version

    return h_prim

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'):

    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()

    s = f.readline()
    s = s.strip("\n")
    CCDID = s.split()

    k = (row_num - 1) * 6 + col_num

    h_ext = fits.Header.fromfile(e_header_fn)

    h_ext['CCDCHIP'] = CCDID[k - 1].rjust(2, '0')
    h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1]
    h_ext['FILTER'] = filters[k-1]
    h_ext['NAXIS1'] = xlen
    h_ext['NAXIS2'] = ylen
    h_ext['EXTNAME'] = extName
    h_ext['GAIN1'] = gain
    h_ext['GAIN2'] = gain
    h_ext['GAIN3'] = gain
    h_ext['GAIN4'] = gain
    h_ext['GAIN5'] = gain
    h_ext['GAIN6'] = gain
    h_ext['GAIN7'] = gain
    h_ext['GAIN8'] = gain
    h_ext['GAIN9'] = gain
    h_ext['GAIN10'] = gain
    h_ext['GAIN11'] = gain
    h_ext['GAIN12'] = gain
    h_ext['GAIN13'] = gain
    h_ext['GAIN14'] = gain
    h_ext['GAIN15'] = gain
    h_ext['GAIN16'] = gain
    h_ext['RDNOIS1'] = readout
    h_ext['RDNOIS2'] = readout
    h_ext['RDNOIS3'] = readout
    h_ext['RDNOIS4'] = readout
    h_ext['RDNOIS5'] = readout
    h_ext['RDNOIS6'] = readout
    h_ext['RDNOIS7'] = readout
    h_ext['RDNOIS8'] = readout
    h_ext['RDNOIS9'] = readout
    h_ext['RDNOIS10'] = readout
    h_ext['RDNOIS11'] = readout
    h_ext['RDNOIS12'] = readout
    h_ext['RDNOIS13'] = readout
    h_ext['RDNOIS14'] = readout
    h_ext['RDNOIS15'] = readout
    h_ext['RDNOIS16'] = readout

    h_ext['PIXSCAL1'] = psize
    h_ext['PIXSCAL2'] = psize

    # h_ext['POS_ANG'] = pa
    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, filter = h_ext['FILTER'])

    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



def main(argv):

    xlen = int(argv[1])
    ylen = int(argv[2])
    pointingNum = argv[3]
    ra = float(argv[4])
    dec = float(argv[5])
    pSize = float(argv[6])
    ccd_row_num = int(argv[7])
    ccd_col_num = int(argv[8])
    pa_aper = float(argv[9])
    gain = float(argv[10])
    readout = float(argv[11])
    dark = float(argv[12])
    fw = float(argv[13])



    h_prim = generatePrimaryHeader(xlen = xlen, ylen = ylen,ra = ra, dec = dec, psize = pSize, row_num = ccd_row_num, col_num = ccd_col_num, pointNum = pointingNum)

    h_ext = generateExtensionHeader(xlen = xlen, ylen = ylen,ra = ra, dec = dec, pa = pa_aper, gain = gain, readout = readout, dark = dark, saturation=fw, psize = pSize, row_num = ccd_row_num, col_num = ccd_col_num)
    hdu1 = fits.PrimaryHDU(header=h_prim)
    hdu2 = fits.ImageHDU(np.zeros([ylen,xlen]),header = h_ext)

    hdul = fits.HDUList([hdu1,hdu2])

    hdul.writeto(h_prim['FILENAME']+'.fits',output_verify='ignore')

# if __name__ == "__main__":
#     main(sys.argv)