An error occurred while loading the file. Please try again.
ImageHeader_1.py 13.55 KiB
"""
generate image header
"""
import numpy as np
from astropy.io import fits
import astropy.wcs as pywcs
from scipy import math
import random
import os
import sys
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])
        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 Header_extention(xlen = 9232, ylen = 9216, 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+']' name = [] value = [] description = [] for k in range(1,2): # f = open("extension"+str(k)+"_image.param","w") j = row_num i = col_num ccdnum = str((j-1)*5+i) name = ['EXTNAME', 'BSCALE', 'BZERO', 'OBSID', 'CCDNAME', 'AMPNAME', 'GAIN', 'RDNOISE', 'DARK', 'SATURATE', 'RSPEED', 'CHIPTEMP', 'CCDCHIP', 'DATASEC', 'CCDSUM', 'NSUM', 'LTM1_1', 'LTM2_2', 'LTV1', 'LTV2', 'ATM1_1', 'ATM2_2', 'ATV1',
'ATV2', 'DTV1', 'DTV2', 'DTM1_1', 'DTM2_2'] value = ['IM'+str(k), 1.0, 0.0, 'CSST.20200101T000000', 'ccd' + ccdnum.rjust(2,'0'), 'ccd' + ccdnum.rjust(2,'0') + ':'+str(k), gain, readout, dark, saturation, 10.0, -100.0, 'ccd' + ccdnum.rjust(2,'0'), data_sec, '1 1', '1 1', flag_ltm_x[k], flag_ltm_y[k], flag_ltv_x[k]*(detector_size_x-20*2+1), flag_ltv_y[k]*(detector_size_y+1), flag_ltm_x[k], flag_ltm_y[k], flag_ltv_x[k]*(detector_size_x-20*2+1), flag_ltv_y[k]*(detector_size_y+1), 0, 0, 1, 1] description = ['Extension name', ' ', ' ', 'Observation ID', 'CCD name', 'Amplifier name', 'Gain (e-/ADU)', 'Readout noise (e-/pixel)', 'Dark noise (e-/pixel/s)', 'Saturation (e-)', 'Read speed', 'Chip temperature', 'CCD chip ID', 'Data section', 'CCD pixel summing', 'CCD pixel summing', 'CCD to image transformation', 'CCD to image transformation', 'CCD to image transformation', 'CCD to image transformation', 'CCD to amplifier transformation', 'CCD to amplifier transformation', 'CCD to amplifier transformation', 'CCD to amplifier transformation', 'CCD to detector transformatio', 'CCD to detector transformatio', 'CCD to detector transformatio', 'CCD to detector transformatio'] return name, value, description ##9232 9216 898 534 1309 60 -40 -23.4333 def WCS_def(xlen = 9232, ylen = 9216, gapx = 898.0, gapy1 = 534, gapy2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1): """ 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 = 5 y_num = 6 detector_num = x_num*y_num detector_size_x = xlen detector_size_y = ylen gap_x = gapx gap_y = [gapy1,gapy2] ra_ref = ra dec_ref = dec pa_aper = pa pixel_size = psize gap_y1_num = 3 gap_y2_num = 2 x_center = (detector_size_x*x_num+gap_x*(x_num-1))/2 y_center = (detector_size_y*y_num+gap_y[0]*gap_y1_num+gap_y[1]*gap_y2_num)/2 gap_y_map = np.array([[0,0,0,0,0],[gap_y[0],gap_y[1],gap_y[1],gap_y[1],gap_y[1]],[gap_y[1],gap_y[0],gap_y[0],gap_y[0],gap_y[0]],[gap_y[0],gap_y[0],gap_y[0],gap_y[0],gap_y[0]],[gap_y[0],gap_y[0],gap_y[0],gap_y[0],gap_y[1]],[gap_y[1],gap_y[1],gap_y[1],gap_y[1],gap_y[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+gap_x)*i-gap_x-detector_size_x/2, detector_size_y*j + sum(gap_y_map[0:j,i-1]) - detector_size_y/2 # print(i,j,x_ref,y_ref,ra_ref,dec_ref) 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") 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 = 9232, ylen = 9216,pointNum = '1', ra = 60, dec = -40, psize = 0.074, row_num = 1, col_num = 1): # array_size1, array_size2, flux, sigma = int(argv[1]), int(argv[2]), 1000.0, 5.0 filerParm_fn = os.path.split(os.path.realpath(__file__))[0] + '/filter.lst' f = open(filerParm_fn) s = f.readline() s = s.strip("\n") filter = s.split(' ') k = (row_num-1)*5+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) h_prim = fits.Header() date = '200930' time_obs = '120000' for i in range(len(name)): if(name[i]=='FILTER'): value[i] = filter[k-1] if(name[i]=='FILENAME'): value[i] = 'CSST_' + date + '_' +time_obs + '_' + pointNum.rjust(6,'0') + '_' +ccdnum.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('==================================================================') h_prim.add_comment('==================================================================',after='FILTER') h_prim.add_comment('Other information') h_prim.add_comment('==================================================================') return h_prim def generateExtensionHeader(xlen = 9232, ylen = 9216,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): h_ext = fits.Header() for i in range(1,2): # NAXIS1:Number of pixels per row; NAXIS2:Number of rows 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)): h_ext[name[j]] = (value[j],description[j]) name, value, description = WCS_def(xlen = xlen, ylen = ylen, gapx = 898.0, gapy1 = 534, gapy2 = 1309, ra = ra, dec = dec, pa = pa ,psize = psize, row_num = row_num, col_num = col_num) 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='DTM2_2') h_ext.add_comment('WCS information') h_ext.add_comment('==================================================================') 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)