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.coordinates import SkyCoord
from astropy.wcs.utils import fit_wcs_from_points
from astropy.time import Time
from astropy import wcs

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_ref = 60, dec_ref = -40, pa = -23.433, pixel_scale = 0.074, pixel_size=1e-2, 
    rotate_chip=0., filter = 'GI', row_num = None, col_num = None, xcen = None, ycen = None):

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

    Parameters
    ----------

    """
    r_dat = OrderedDict()
    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
    
    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]
    pa_aper = pa

    if (row_num is not None) and (col_num is not None):

        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]

        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]]])
        
        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
        
        for k in range(1,2):
            
            cd = np.array([[ pixel_scale,  0], [0, pixel_scale]])/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['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']:

                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, xlen=xlen, ylen=ylen, w=w)
                        # print(sc1[0,0],sc1[0,1])
                        sky_coors.append((sc1[0, 0], sc1[0, 1]))

                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]

    elif (xcen is not None) and (ycen is not None):
        xcen, ycen = xcen/pixel_size, ycen/pixel_size
        x1, y1 = xcen - xlen/2., ycen - ylen/2.
        r_dat['CRPIX1'] = -x1
        r_dat['CRPIX2'] = -y1

        cd = np.array([[ pixel_scale,  0], [0, pixel_scale]])/3600.*flag_x[1]
        cd_rot = rotate_CD_matrix(cd, pa_aper)
        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]

        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']]

        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=rotate_chip, xlen=xlen, ylen=ylen, w=w)
                sky_coors.append((sc1[0, 0], sc1[0, 1]))
        wcs_new = fit_wcs_from_points(xy=np.array([x_pixs, y_pixs]), 
                                      world_coords=SkyCoord(sky_coors, frame="icrs", unit="deg"), 
                                      projection='TAN')
        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]

    else:
        raise ValueError('In function WCS_def(): Either (row_num, col_num) or (xcen, ycen, pixel_size) should be given')
    
    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, pixel_scale = 0.074, date='200930', time_obs='120000', im_type = 'MS', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.], project_cycle=6, chip_name="01"):

    # 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

    obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'}

    OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + pointNum.rjust(7,'0')

    h_prim['OBJECT'] = str(int(project_cycle)) + pointNum.rjust(7,'0')
    h_prim['OBSID'] = OBS_id
    # h_prim['TELFOCUS'] = 'f/14'
    h_prim['EXPTIME'] = exptime
    
    # Define file types
    file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'cosmic_ray', 'CRD':'cosmic_ray','CALS':'CALS','CALF':'CALF'}
    h_prim['FILETYPE'] = 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, '02.0f')
    dec_hms = format(co.dec.dms.d, '02.0f') + format(abs(co.dec.dms.m), '02.0f') + format(abs(co.dec.dms.s), '02.0f')
    h_prim['TARGET'] = ra_hms + '+' + dec_hms
    #
    # 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'] = pixel_scale
    # h_prim['PIXSCAL2'] = pixel_scale

    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 + '_' + OBS_id + '_' + CCDID[
    #     k - 1].rjust(2, '0') + '_L0_V01'
    h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'


    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(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, pixel_scale = 0.074, pixel_size=1e-2, 
                            extName='SCI', row_num = None, col_num = None, xcen=None, ycen=None):

    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['CCDCHIP'] = chip.chip_name
    h_ext['CCDLABEL'] = chip.filter_type + '-' + str(chip.chipID).rjust(2, '0')
    h_ext['FILTER'] = chip.filter_type
    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'] = pixel_scale
    h_ext['PIXSCAL2'] = pixel_scale

    # h_ext['POS_ANG'] = pa
    header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra_ref=ra, dec_ref=dec, pa=pa, pixel_scale=pixel_scale, pixel_size=pixel_size, 
                         rotate_chip=chip.rotate_angle, filter = h_ext['FILTER'], row_num=row_num, col_num=col_num, xcen = xcen, ycen = ycen)

    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)