"""
Creates output .fits files for BUBE
Created on 21 february 2019
@author: Catalano Camille, APC
"""
import logging
import numpy as np
import astropy.units as u
import astropy.table as astt
from astropy.io import fits
from astropy.time import Time
from ecpi.common.io import fits_tools
from ecpi.process.bube.core.main_bube import DetectorImages
from ecpi.common.mission.observation import EclairsObservation
from ecpi.common.instru.array_convention import yzegp_to_ijdet_array
from ecpi.common.instru.model_effect import ECLAIRsDetectorEffect
logger = logging.getLogger(__name__)
[docs]def save_ecl_det_ubc(det_image: DetectorImages, tstart, tstop,
file_path, attitude, eded: ECLAIRsDetectorEffect,
creator='BUBE?', proc_id="01"):
"""save detector corrected images in a fits file : ECL-DET-UBC
common keywords are used
..warning: it is assumed the global efficiency matrix is the same for each channel !
filename is ('ECL-DET-UBC-' + proc_id + '.fits')
:param det_image: detector image to save
:type det_image: DetectorImages
: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 file_path: PATH to the directory where to write the file
:type file_path: string
:param eded: ECLAIRsDetectorEffect object containing path to ARF file
:type eded: ECLAIRsDetectorEffect object
:param attitude: [ra, dec, ori] in degrees
:type attitude: [float, float, float]
:param creator: program that has generated the file. Default='BUBE?'.
:type creator: string
:param proc_id: id of the process. Default="01".
:type proc_id: string
"""
prihdr = fits.Header()
prihdu = fits.PrimaryHDU(header=prihdr)
obstime = tstop - tstart
fits_tools.pointing_keywords(prihdu, attitude, obstime)
fits_tools.common_keyword(prihdu, tstart, tstop, creator, 1.0, proc_id)
prihdu.header['CARD'] = ("ECL-DET-UBC", "Product type")
hdu_list = [prihdu]
hdu_images = []
grp_table = astt.Table([[]]*11,
names=('member_xtension', 'member_name', 'member_version',
'member_position', 'member_location', 'member_uri_type',
'imatype', 'chanmin', 'chanmax', 'e_min', 'e_max'),
dtype=('S8', 'S32', 'i4', 'i4', 'S256', 'S3', 'S20', 'f4', 'f4', 'f4', 'f4'))
grp_table['e_min'].unit = u.keV
grp_table['e_max'].unit = u.keV
for shadow_idx in range(len(det_image.shadowgrams)):
chanmin = det_image.energy_ranges[shadow_idx][0]
chanmax = det_image.energy_ranges[shadow_idx][1]
e_min = eded.get_energymin_val(chanmin)#
e_max = eded.get_energymax_val(chanmax)
# add row in astropy table for each field..
grp_table.add_row(['IMAGE ', 'Intensity', '1', shadow_idx*5+2, '', '',
'intensity', chanmin, chanmax, e_min, e_max])
grp_table.add_row(['IMAGE ', 'Variance', '1', shadow_idx*5+3, '', '',
'variance', chanmin, chanmax, e_min, e_max])
grp_table.add_row(['IMAGE ', 'GlobalEfficiency', '1', shadow_idx*5+4, '', '',
'GlabalEfficiency', chanmin, chanmax, e_min, e_max])
# declare elementary HDU for each filed..
shd_bg_corr = yzegp_to_ijdet_array(det_image.shadowgrams[shadow_idx])
shd_var = yzegp_to_ijdet_array(det_image.shadowgrams_var[shadow_idx])
globeff = det_image.global_efficiency
hdu_intensity = fits.ImageHDU(np.array(shd_bg_corr, dtype='float32'))
hdu_var = fits.ImageHDU(np.array(shd_var, dtype='float32'))
hdu_globeff = fits.ImageHDU(np.array(globeff, dtype='float32'))
hdu_intensity.name = 'ShadowgramCorr_{}_{}'.format(int(e_min), int(e_max))
hdu_var.name = 'ShadowgramVar_{}_{}'.format(int(e_min), int(e_max))
hdu_globeff.name = 'GlobalEfficiencyMatrix_{}_{}'.format(int(e_min), int(e_max))
# fill each elementary HDU with appropriate KW.
fits_tools.common_keyword(hdu_intensity, tstart, tstop, creator, 1.0, proc_id)
fits_tools.energy_keywords(hdu_intensity, chanmin, chanmax, e_min, e_max)
#hdu_intensity.header.extend(sky_images[image_idx].wcs.to_header().cards)
hdu_intensity.header['CARD'] = ("ECL-DET-UBC", "Product type")
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)
hdu_var.header['CARD'] = ("ECL-DET-UBC", "Product type")
fits_tools.common_keyword(hdu_globeff, tstart, tstop, creator, 1.0, proc_id)
fits_tools.energy_keywords(hdu_globeff, chanmin, chanmax, e_min, e_max)
#hdu_globeff.header.extend(sky_images[image_idx].wcs.to_header().cards)
hdu_globeff.header['CARD'] = ("ECL-DET-UBC", "Product type")
hdu_images.extend([hdu_intensity, hdu_var, hdu_globeff])
hdu_grp = fits.table_to_hdu(grp_table)
hdu_grp.header['comment'] = ("ECLAIRs detector images corrected from uniformity,",
"background and earth occultation/albedo")
hdu_grp.header['EXTNAME'] = ("ECL-DET-UBC", "Product type")
fits_tools.grouphdu_keyword(hdu_grp, 'ECL-DET-UBC-GRP', attitude)
fits_tools.common_keyword(hdu_grp, tstart, tstop, creator, 1.0, proc_id)
hdu_list.append(hdu_grp)
hdu_list.extend(hdu_images)
thdulist = fits.HDUList(hdu_list)
time_isot_tag = create_time_tag_from_tstart(tstart)
thdulist.writeto(file_path + '/ECL-DET-UBC-UTC' + time_isot_tag + '.fits',
overwrite=True, checksum=True)
[docs]def create_time_tag_from_tstart(tstart):
"""Creates time tag of scientific products.
:param tstart: starting time of the observation
:type tstart: float
"""
t_ref = Time("2017-01-01T00:00:00.000").mjd
tstart_day = tstart / 24 / 3600.
time_isot_tag = Time((tstart_day+t_ref), scale='tt', format='mjd').isot
time_isot_tag = time_isot_tag.replace("-", "")
time_isot_tag = time_isot_tag.replace(":", "")
return time_isot_tag.split(".")[0]
[docs]def read_attitude_files(l_files):
"""Load informations contained in SVO-ATT-CNV files.
..warning: works with a single SVO-ATT-CNV file for the moment.
..warning: works with a SVO-ATT-CNV file with a single line for the moment.
:param l_files: list of SVO-ATT-CNV files as specified in parameter file.
:type l_files: list
:return: dictionary of read values in SVO-ATT-CNV file
:rtype: dict
"""
att_file = l_files[0]
hdul = fits.open(att_file)
dict_att = {
'TIME_AAV': hdul[1].data[0][0],
'QPARAM': hdul[1].data[0][1],
'ANGLE_VEL': hdul[1].data[0][2]
#other columns exist that could be read.
}
logger.info(f"loading in memory attitude file {att_file}")
hdul.close()
return dict_att
[docs]def read_table_svo_att_cnv(p_files):
'''
:param p_files:
'''
return read_data_ext(p_files)
[docs]def read_table_svo_orb_cnv(p_files):
'''
:param p_files:
'''
return read_data_ext(p_files)
[docs]def read_data_ext(p_files, n_ext=1):
'''
:param p_files:
:param n_ext:
'''
hdul = fits.open(p_files)
ext_data = hdul[n_ext].data
hdul.close()
return ext_data
[docs]def read_orbit_files(l_files):
"""Load informations contained in SVO-ORB-CNV files.
..warning: works with a single SVO-ORB-CNV file for the moment.
..warning: works with a SVO-ORB-CNV file with a single line for the moment.
:param l_files: list of SVO-ORB-CNV files as specified in parameter file.
:type l_files: list
:return: dictionary of read values in SVO-ORB-CNV file
:rtype: dict
"""
orb_file = l_files[0]
hdul = fits.open(orb_file)
dict_orb = {
'TIME_PVT': hdul[1].data[0][0],
'POSITION': hdul[1].data[0][2],
'VELOCITY': hdul[1].data[0][3]
#other columns exist that could be read.
}
logger.info(f"loading in memory orbit file {orb_file}")
hdul.close()
return dict_orb
[docs]def check_time_coherence_evt_att(t_start, d_attitude):
"""Check the time consistence between event and
attitude files.
:param t_start: starting time of the observation in second.
:type t_start: float
:param d_attitude: dictionnary collecting infos from attitude files.
:type d_attitude: dict
:return: whether origin of time is coherent
:rtype: boolean
"""
date_att = d_attitude['TIME_AAV']
logger.debug(f'check time coherence: t_start={t_start} | att file={date_att}')
return np.isclose(t_start, date_att, atol=1e-1)