Source code for process.imag.io.outputs

'''fits outputs for imag

Created on 21 february 2019

@author: Catalano Camille, APC
'''


from astropy.io import fits
import astropy.table as astt

from common.io import fits_tools
from common.instru.array_convention import yzegp_to_ijdet_array
        
        
def update_images_headers(hdu_count, hdu_var, hdu_snr, hdu_exp, hdu_mod, time_exp_mn):
    """ add main general keywords to the images headers
    
    :param hdu_count: hdu of the cleaned sky count
    :type hdu_count: fits.ImageHDU
    :param hdu_var: hdu of the cleaned sky var
    :type hdu_var: fits.ImageHDU
    :param hdu_snr: hdu of the cleaned sky snr
    :type hdu_snr: fits.ImageHDU
    :param hdu_exp: hdu of the exposure
    :type hdu_exp: fits.ImageHDU
    :param hdu_mod: hdu of the source model image
    :type hdu_mod: fits.ImageHDU
    :param time_exp_mn: time of the observation
    :type time_exp_mn: float
    
    :return: updated hdu_count, hdu_var, hdu_snr, hdu_exp, hdu_mod
    :rtype: fits.ImageHDU * 5
    """
    hdu_count.name = 'Intensity'
    hdu_count.header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_count.header['comment'] = "Intensity Image"
    hdu_count.header['BUNIT'] = 'ct.s^-1.cm^-2'
    hdu_count.header['IMATYPE'] = ('Intensity', "Type of image")
    hdu_count.header['HDUCLAS1'] = ('IMAGE', 'Dataset contains a sky image')
    hdu_var.name = 'Variance'
    hdu_var.header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_var.header['comment'] = "Variance Image"
    hdu_var.header['BUNIT'] = 'ct^2.s^-2.cm^-4'
    hdu_var.header['IMATYPE'] = ('Variance', "Type of image")
    hdu_var.header['HDUCLAS1'] = ('IMAGE', 'Dataset contains a sky image')
    hdu_snr.name = 'SNR'
    hdu_snr.header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_snr.header['comment'] = "Signal to Noise Image"
    hdu_snr.header['IMATYPE'] = ('SNR', "Type of image")
    hdu_snr.header['HDUCLAS1'] = ('IMAGE', 'Dataset contains a sky image')
    hdu_exp.name = 'Exposure'
    hdu_exp.header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_exp.header['comment'] = "Exposure Image"
    hdu_exp.header['BUNIT'] = 'seconds'
    hdu_exp.header['IMATYPE'] = ('Exposure', "Type of image")
    hdu_mod.name = 'Model'
    hdu_mod.header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_mod.header['comment'] = "Source models Image"
    hdu_mod.header['BUNIT'] = 'ct.s^-1.cm^-2'
    hdu_mod.header['IMATYPE'] = ('Model', "Type of image")
    hdu_mod.header['HDUCLAS1'] = ('IMAGE', 'Dataset contains a sky image')
    return hdu_count, hdu_var, hdu_snr, hdu_exp, hdu_mod
        
        
def save_sky_images(sky_images, file_path, tstart, tstop, energy_ranges, creator='IMAG?', proc_id="01"):
    """
    save images in a fits file ECL-SKY-IMA
    
    use common keywords
    
    filename is 'ECL-SKY-IMA-' + proc_id + '.fits'
    
    file contains: 
        - Intensity: cleaned sky count image
        - Variance: cleaned sky variance image
        - SNR: cleaned sky signal to noise ratio image
        - Exposure: initial sky Exposure image
        - Model: Added sources model images (for debug)
        - Initial Sky: Initial sky intensity image (for debug)
        
    save also the models images in a fits file: see save_models
    
    :param sky_images: list of sky images and models
    :type sky_images: list(SkyImages)
    :param file_path: PATH to the directory where to write the file
    :type file_path: string
    :param tstart: start time of the observation in s from mjdref
    :type tstart: float
    :param tstop: end time of the observation in s from mjdref
    :type tstop: float
    :param energy_ranges: energy ranges (in PI) of the images [[low_energy, high_energy],[]]
    :type energy_ranges: [[int, int]]
    :param creator: program that has generated the file. Default='IMAG?'
    :type creator: string
    :param proc_id: id of the process. Default="01".
    :type proc_id: string
    """
    # construction of the fits file
    prihdr = fits.Header()
    prihdu = fits.PrimaryHDU(header=prihdr)
    hdu_list = [prihdu]
    hdu_images = []
    grp_table = astt.Table([[],[],[],[],[],[],[],[],[],[],[]],names=('member_xtension', 'member_name', 'member_version', 'member_position', 'member_location', 'member_uri_type', 'imatype', 'chanmin', 'chanmax', 'e_min', 'e_max'), dtype=('S20', 'S20', 'u1', 'u1', 'S20', 'S20', 'S20', 'f4', 'f4', 'f4', 'f4'))

    for image_idx in range(len(sky_images)):
        chanmin = energy_ranges[image_idx][0]
        chanmax = energy_ranges[image_idx][1]
        e_min = chanmin
        e_max = chanmax
        grp_table.add_row(['IMAGE   ', 'Intensity', '1', image_idx*5+2, '', '', 'intensity', chanmin, chanmax, e_min, e_max])
        grp_table.add_row(['IMAGE   ', 'Variance', '1', image_idx*5+3, '', '', 'variance', chanmin, chanmax, e_min, e_max])
        grp_table.add_row(['IMAGE   ', 'SNR', '1', image_idx*5+4, '', '', 'SNR', chanmin, chanmax, e_min, e_max])
        grp_table.add_row(['IMAGE   ', 'Exposure', '1', image_idx*5+5, '', '', 'Exposure', chanmin, chanmax, e_min, e_max])
        grp_table.add_row(['IMAGE   ', 'Model', '1', image_idx*5+6, '', '', 'Model', chanmin, chanmax, e_min, e_max])

        hdu_count = fits.ImageHDU(sky_images[image_idx].cleaned_sky_count)
        hdu_var = fits.ImageHDU(sky_images[image_idx].cleaned_sky_var)
        hdu_snr = fits.ImageHDU(sky_images[image_idx].cleaned_sky_snr)
        hdu_exp = fits.ImageHDU(sky_images[image_idx].exposure)
        hdu_mod = fits.ImageHDU(sky_images[image_idx].model_image)
        fits_tools.common_keyword(hdu_count, tstart, tstop, creator, 1.0, proc_id)
        fits_tools.energy_keywords(hdu_count, chanmin, chanmax, e_min, e_max)
        hdu_count.header.extend(sky_images[image_idx].wcs.to_header().cards)
        fits_tools.common_keyword(hdu_var, tstart, tstop, creator, 1.0, proc_id)
        fits_tools.energy_keywords(hdu_var, chanmin, chanmax, e_min, e_max)
        hdu_var.header.extend(sky_images[image_idx].wcs.to_header().cards)
        fits_tools.common_keyword(hdu_snr, tstart, tstop, creator, 1.0, proc_id)
        fits_tools.energy_keywords(hdu_snr, chanmin, chanmax, e_min, e_max)
        hdu_snr.header.extend(sky_images[image_idx].wcs.to_header().cards)
        fits_tools.common_keyword(hdu_exp, tstart, tstop, creator, 1.0, proc_id)
        fits_tools.energy_keywords(hdu_exp, chanmin, chanmax, e_min, e_max)
        hdu_exp.header.extend(sky_images[image_idx].wcs.to_header().cards)
        fits_tools.common_keyword(hdu_mod, tstart, tstop, creator, 1.0, proc_id)
        fits_tools.energy_keywords(hdu_mod, chanmin, chanmax, e_min, e_max)
        hdu_mod.header.extend(sky_images[image_idx].wcs.to_header().cards)
        hdu_count, hdu_var, hdu_snr, hdu_exp, hdu_mod = update_images_headers(hdu_count, hdu_var, hdu_snr, hdu_exp, hdu_mod, sky_images[image_idx].time_exp_mn)
        hdu_images.extend([hdu_count, hdu_var, hdu_snr, hdu_exp, hdu_mod])
        
    hdu_grp = fits.table_to_hdu(grp_table)
    hdu_grp.header['comment'] = 'ECLAIRs reconstructed and cleaned individual Sky Images'
    fits_tools.common_keyword(hdu_grp, tstart, tstop, creator, 1.0, proc_id)
    fits_tools.grouphdu_keyword(hdu_grp, 'ECL-SKY-IMA-IDX')
    hdu_list.append(hdu_grp)
    hdu_list.extend(hdu_images)
    thdulist = fits.HDUList(hdu_list)
    thdulist.writeto(file_path + '/ECL-SKY-IMA-' + proc_id + '.fits', overwrite=True, checksum=True)
        
    save_src_models(sky_images, file_path, tstart, tstop, energy_ranges, creator, proc_id)
    
    
def save_src_models(sky_images, file_path, tstart, tstop, energy_ranges, creator, proc_id="01"):
    """save models independant shadowgrams in fits file ECL-DET-MOD
    
    use common keywords
    
    filename is 'ECL-DET-MOD-' + proc_id + '.fits'
    
    file contains: for each source model, the shadowgram of the source
    
    :param sky_images: list of sky images and models
    :type sky_images: list(SkyImages)
    :param file_path: PATH to the directory where to write the file
    :type file_path: string
    :param tstart: start time of the observation in s from mjdref
    :type tstart: float
    :param tstop: end time of the observation in s from mjdref
    :type tstop: float
    :param energy_ranges: energy ranges (in PI) of the images [[low_energy, high_energy],[]]
    :type energy_ranges: [[int, int]]
    :param creator: program that has generated the file.
    :type creator: string
    :param proc_id: id of the process. Default="01".
    :type proc_id: string
    """
    prihdr = fits.Header()
    prihdu = fits.PrimaryHDU(header=prihdr)
    hdu_list = [prihdu]
    hdu_images = []
    grp_table = astt.Table([[],[],[],[],[],[],[],[],[],[],[]],names=('member_xtension', 'member_name', 'member_version', 'member_position', 'member_location', 'member_uri_type', 'imatype', 'chanmin', 'chanmax', 'e_min', 'e_max'), dtype=('S20', 'S20', 'u1', 'u1', 'S20', 'S20', 'S20', 'f4', 'f4', 'f4', 'f4'))

    model_fits_position = 2
    for image_idx in range(len(sky_images)):
        for [model_shadowgram, model_name] in sky_images[image_idx].models_shadowgrams:
            chanmin = energy_ranges[image_idx][0]
            chanmax = energy_ranges[image_idx][1]
            e_min = chanmin
            e_max = chanmax
            grp_table.add_row(['IMAGE   ', model_name, '1', model_fits_position, '', '', 'Intensity', chanmin, chanmax, e_min, e_max])
    
            hdu_mod = fits.ImageHDU(yzegp_to_ijdet_array(model_shadowgram))
            hdu_mod.name = model_name
            fits_tools.common_keyword(hdu_mod, tstart, tstop, creator, 1.0, proc_id)
            fits_tools.energy_keywords(hdu_mod, chanmin, chanmax, e_min, e_max)
            hdu_mod.header.extend(sky_images[image_idx].wcs.to_header().cards)
            hdu_mod.header['comment'] = model_name + ' source model shadowgram'
            hdu_mod.header['BUNIT'] = 'ct.s^-1.cm^-2'
            hdu_mod.header['IMATYPE'] = ('Intensity', "Type of image")
            hdu_images.append(hdu_mod)
            model_fits_position += 1
    hdu_grp = fits.table_to_hdu(grp_table)
    hdu_grp.header['comment'] = 'ECLAIRs detector model images of point sources'
    fits_tools.common_keyword(hdu_grp, tstart, tstop, creator, 1.0, proc_id)
    fits_tools.grouphdu_keyword(hdu_grp, 'ECL-DET-MOD-IDX')
    hdu_list.append(hdu_grp)
    hdu_list.extend(hdu_images)
    thdulist = fits.HDUList(hdu_list)
    thdulist.writeto(file_path + '/ECL-DET-MOD-' + proc_id + '.fits', overwrite=True, checksum=True)
    
    
def update_images_debug_headers(hdu_images_tmp, time_exp_mn):
    """add main general keywords to the debug images headers
    
    :param hdu_images_tmp: list of hdu of debug images
    :type hdu_images_tmp: list(hdu)
    :param time_exp_mn: time of the observation
    :type time_exp_mn: float
    
    :return: updated hdus
    :rtype: list(hdu)
    """
    hdu_images_tmp[0].name = 'Initial Sky'
    hdu_images_tmp[0].header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_images_tmp[0].header['comment'] = "Initial Sky Intensity Image"
    hdu_images_tmp[0].header['BUNIT'] = 'ct.s^-1.cm^-2'
    hdu_images_tmp[0].header['IMATYPE'] = ("Intensity", "Type of image")
    hdu_images_tmp[0].header['HDUCLAS1'] = ('IMAGE', 'Dataset contains a sky image')
    hdu_images_tmp[1].name = 'Initial Variance'
    hdu_images_tmp[1].header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_images_tmp[1].header['comment'] = "Initial Variance Image"
    hdu_images_tmp[1].header['BUNIT'] = 'ct^2.s^-2.cm^-4'
    hdu_images_tmp[1].header['IMATYPE'] = ("Variance", "Type of image")
    hdu_images_tmp[1].header['HDUCLAS1'] = ('IMAGE', 'Dataset contains a sky image')
    hdu_images_tmp[2].name = 'Initial SNR'
    hdu_images_tmp[2].header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_images_tmp[2].header['comment'] = "Initial Signal to Noise Image"
    hdu_images_tmp[2].header['IMATYPE'] = ("SNR", "Type of image")
    hdu_images_tmp[2].header['HDUCLAS1'] = ('IMAGE', 'Dataset contains a sky image')
    hdu_images_tmp[3].name = 'Remaining Sky'
    hdu_images_tmp[3].header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_images_tmp[3].header['comment'] = "Remaining Sky Intensity Image"
    hdu_images_tmp[3].header['BUNIT'] = 'ct.s^-1.cm^-2'
    hdu_images_tmp[3].header['IMATYPE'] = ("Intensity", "Type of image")
    hdu_images_tmp[3].header['HDUCLAS1'] = ('IMAGE', 'Dataset contains a sky image')
    hdu_images_tmp[4].name = 'Remaining Variance'
    hdu_images_tmp[4].header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_images_tmp[4].header['comment'] = "Remaining Variance Image"
    hdu_images_tmp[4].header['BUNIT'] = 'ct^2.s^-2.cm^-4'
    hdu_images_tmp[4].header['IMATYPE'] = ("Variance", "Type of image")
    hdu_images_tmp[4].header['HDUCLAS1'] = ('IMAGE', 'Dataset contains a sky image')
    hdu_images_tmp[5].name = 'Remaining SNR'
    hdu_images_tmp[5].header['timeexp'] = str(time_exp_mn) + ' minutes'
    hdu_images_tmp[5].header['comment'] = "Remaining Signal to Noise Image"
    hdu_images_tmp[5].header['IMATYPE'] = ("SNR", "Type of image")
    hdu_images_tmp[5].header['HDUCLAS1'] = ('IMAGE', 'Dataset contains a sky image')
    return hdu_images_tmp
        
        
def save_sky_images_debug(sky_images, file_path, tstart, tstop, energy_ranges, creator, proc_id="01"):
    """save auxiliary images in fits file ECL-SKY-IMA_DEBUG
    
    use common keywords
    
    filename is 'ECL-SKY-IMA-' + proc_id + '_DEBUG.fits'
    
    file contains: for each source model, the shadowgram of the source
    
    :param sky_images: list of sky images and models
    :type sky_images: list(SkyImages)
    :param file_path: PATH to the directory where to write the file
    :type file_path: string
    :param tstart: start time of the observation in s from mjdref
    :type tstart: float
    :param tstop: end time of the observation in s from mjdref
    :type tstop: float
    :param energy_ranges: energy ranges (in PI) of the images [[low_energy, high_energy],[]]
    :type energy_ranges: [[int, int]]
    :param creator: program that has generated the file.
    :type creator: string
    :param proc_id: id of the process. Default="01".
    :type proc_id: string
    """
    prihdr = fits.Header()
    prihdu = fits.PrimaryHDU(header=prihdr)
    hdu_list = [prihdu]
    hdu_images = []
    grp_table = astt.Table([[],[],[],[],[],[],[],[],[],[],[]],names=('member_xtension', 'member_name', 'member_version', 'member_position', 'member_location', 'member_uri_type', 'imatype', 'chanmin', 'chanmax', 'e_min', 'e_max'), dtype=('S20', 'S20', 'u1', 'u1', 'S20', 'S20', 'S20', 'f4', 'f4', 'f4', 'f4'))

    for image_idx in range(len(sky_images)):
        chanmin = energy_ranges[image_idx][0]
        chanmax = energy_ranges[image_idx][1]
        e_min = chanmin
        e_max = chanmax
        grp_table.add_row(['IMAGE   ', 'Initial Sky', '1', image_idx*6+2, '', '', 'intensity', chanmin, chanmax, e_min, e_max])
        grp_table.add_row(['IMAGE   ', 'Initial Variance', '1', image_idx*6+3, '', '', 'variance', chanmin, chanmax, e_min, e_max])
        grp_table.add_row(['IMAGE   ', 'Initial SNR', '1', image_idx*6+4, '', '', 'SNR', chanmin, chanmax, e_min, e_max])
        grp_table.add_row(['IMAGE   ', 'Remaining Sky', '1', image_idx*6+5, '', '', 'Intensity', chanmin, chanmax, e_min, e_max])
        grp_table.add_row(['IMAGE   ', 'Remaining Variance', '1', image_idx*6+6, '', '', 'Variance', chanmin, chanmax, e_min, e_max])
        grp_table.add_row(['IMAGE   ', 'Remaining SNR', '1', image_idx*6+7, '', '', 'SNR', chanmin, chanmax, e_min, e_max])

        hdu_init_sky = fits.ImageHDU(sky_images[image_idx].initial_sky_count)
        hdu_init_var = fits.ImageHDU(sky_images[image_idx].initial_sky_var)
        hdu_init_snr = fits.ImageHDU(sky_images[image_idx].initial_sky_snr)
        hdu_remain_sky = fits.ImageHDU(sky_images[image_idx].remains_sky_count)
        hdu_remain_var = fits.ImageHDU(sky_images[image_idx].remains_sky_var)
        hdu_remain_snr = fits.ImageHDU(sky_images[image_idx].remains_sky_snr)
        hdu_images_tmp = [hdu_init_sky, hdu_init_var, hdu_init_snr, hdu_remain_sky, hdu_remain_var, hdu_remain_snr]
        for hdu in hdu_images_tmp:
            fits_tools.common_keyword(hdu, tstart, tstop, creator, 1.0, proc_id)
            fits_tools.energy_keywords(hdu, chanmin, chanmax, e_min, e_max)
            hdu.header.extend(sky_images[image_idx].wcs.to_header().cards)
        hdu_images_tmp = update_images_debug_headers(hdu_images_tmp, sky_images[image_idx].time_exp_mn)
        hdu_images.extend(hdu_images_tmp)
    hdu_grp = fits.table_to_hdu(grp_table)
    fits_tools.common_keyword(hdu_grp, tstart, tstop, creator, 1.0, proc_id)
    fits_tools.grouphdu_keyword(hdu_grp, 'ECL-SKY-IMA-IDX')
    hdu_list.append(hdu_grp)
    hdu_list.extend(hdu_images)
    thdulist = fits.HDUList(hdu_list)
    thdulist.writeto(file_path + '/ECL-SKY-IMA-' + proc_id + '_DEBUG.fits', overwrite=True, checksum=True)

    
def save_catalogs(catalogs_output, file_path, tstart, tstop, energy_ranges, creator, proc_id="01"):
    """save output sources catalogs in fits file ECL-SOP-IMA
    
    use common keywords
    
    filename is 'ECL-SOP-IMA-' + proc_id + '.fits'
    
    file contains: for each energy bands, the catalog of identified sources
    
    :param catalogs_output: list of sources catalogs
    :type catalogs_output: list(CatalogIdentifiedSources)
    :param file_path: PATH to the directory where to write the file
    :type file_path: string
    :param tstart: start time of the observation in s from mjdref
    :type tstart: float
    :param tstop: end time of the observation in s from mjdref
    :type tstop: float
    :param energy_ranges: energy ranges (in PI) of the catalogs [[low_energy, high_energy],[]]
    :type energy_ranges: [[int, int]]
    :param creator: program that has generated the file.
    :type creator: string
    :param proc_id: id of the process. Default="01".
    :type proc_id: string
    """
    prihdr = fits.Header()
    prihdu = fits.PrimaryHDU(header=prihdr)
    hdu_list = [prihdu]
    hdu_catalogs = []
    grp_table = astt.Table([[],[],[],[],[],[],[],[],[],[]],names=('member_xtension', 'member_name', 'member_version', 'member_position', 'member_location', 'member_uri_type', 'chanmin', 'chanmax', 'e_min', 'e_max'), dtype=('S20', 'S20', 'u1', 'u1', 'S20', 'S20', 'f4', 'f4', 'f4', 'f4'))

    for catalog_idx in range(len(catalogs_output)):
        chanmin = energy_ranges[catalog_idx][0]
        chanmax = energy_ranges[catalog_idx][1]
        e_min = chanmin
        e_max = chanmax
        grp_table.add_row(['BINTABLE   ', "ECL-SKY-RES", '1', catalog_idx+2, '', '', chanmin, chanmax, e_min, e_max])

        hdu = fits.table_to_hdu(catalogs_output[catalog_idx]._catalog)
        hdu.name = "ECL-SKY-RES"
        fits_tools.common_keyword(hdu, tstart, tstop, creator, 1.0, proc_id)
        fits_tools.energy_keywords(hdu, chanmin, chanmax, e_min, e_max)
        hdu.header['comment'] = 'source parameters from sky images'
        hdu.header['BASETYPE'] = ('DAL_TABLE', 'Data Access Layer base type')
        hdu_catalogs.append(hdu)
    hdu_grp = fits.table_to_hdu(grp_table)
    hdu_grp.header['comment'] = 'Source Parameters derived from single reconstructed sky image'
    fits_tools.common_keyword(hdu_grp, tstart, tstop, creator, 1.0, proc_id)
    fits_tools.grouphdu_keyword(hdu_grp, 'ECL-SOP-IMA-IDX')
    hdu_list.append(hdu_grp)
    hdu_list.extend(hdu_catalogs)
    thdulist = fits.HDUList(hdu_list)
    thdulist.writeto(file_path + '/ECL-SOP-IMA-' + proc_id + '.fits', overwrite=True, checksum=True)