import numpy as np
import pylab as plt
import os
import collections
from astropy.io import fits
import pdb
from astropy.io.fits.hdu.image import _ImageBaseHDU
[docs]
module_dir = os.path.dirname(__file__)
[docs]
class Instrument(object):
def __init__(self):
# Must define the following parameters as a
# Define
self.hdr_keys['filename'] = 'filename'
return
[docs]
def get_bad_pixel_mask_name(self):
return self.bad_pixel_mask
[docs]
def get_filter_name(self, hdr):
pass
[docs]
def get_plate_scale(self, hdr):
pass
[docs]
def get_position_angle(self, hdr):
pass
[docs]
def get_parallactic_angle(self,hdr):
pass
[docs]
def get_instrument_angle(self, hdr):
pass
[docs]
def get_central_wavelength(self, hdr):
pass
[docs]
def get_gain(self, hdr):
pass
[docs]
def make_filenames(self, files, rootDir='', prefix='n'):
pass
[docs]
def get_distortion_maps(self, date):
pass
[docs]
def get_align_type(self, errors=False):
pass
[docs]
def get_saturation_level(self):
pass
[docs]
def get_lin_corr_coeffs(self):
pass
[docs]
class NIRC2(Instrument):
def __init__(self):
# Define header keywords
self.hdr_keys['filename'] = 'FILENAME'
self.hdr_keys['object_name'] = 'OBJECT'
self.hdr_keys['itime'] = 'ITIME'
self.hdr_keys['coadds'] = 'COADDS'
self.hdr_keys['sampmode'] = 'SAMPMODE'
self.hdr_keys['nfowler'] = 'MULTISAM'
self.hdr_keys['camera'] = 'CAMNAME'
self.hdr_keys['shutter'] = 'SHRNAME'
self.hdr_keys['mjd'] = 'MJD-OBS'
self.hdr_keys['elevation'] = 'EL'
[docs]
self.bad_pixel_mask = 'nirc2mask.fits'
[docs]
self.distXgeoim = module_dir + '/reduce//distortion/nirc2_narrow_xgeoim.fits'
[docs]
self.distYgeoim = module_dir + '/reduce//distortion/nirc2_narrow_ygeoim.fits'
[docs]
self.telescope = 'Keck'
[docs]
self.telescope_diam = 10.5 # telescope diameter in meters
return
[docs]
def get_filter_name(self, hdr):
filter1 = hdr['fwiname']
filter2 = hdr['fwoname']
filt = filter1
if (filter1.startswith('PK')):
filt = filter2
return filt
[docs]
def get_plate_scale(self, hdr):
"""
Return the plate scale in arcsec/pixel.
"""
# Setup NIRC2 plate scales
# Units are arcsec/pixel
scales = {"narrow": 0.009952,
"medium": 0.019829,
"wide": 0.039686}
scale = scales[hdr['CAMNAME']]
return scale
[docs]
def get_position_angle(self, hdr):
"""
Get the sky PA in degrees East of North.
"""
return float(hdr['ROTPOSN']) - float(hdr['INSTANGL'])
[docs]
def get_parallactic_angle(self,hdr):
"""
Get the parallactic angle in degrees East of North
"""
q = hdr['PARANG']
return q
[docs]
def get_instrument_angle(self, hdr):
return float(hdr['INSTANGL'])
[docs]
def get_central_wavelength(self, hdr):
"""
Return the central wavelength of the filter for
this observation in microns.
"""
return float(hdr['CENWAVE'])
[docs]
def get_gain(self, hdr):
return hdr['GAIN']
[docs]
def make_filenames(self, files, rootDir='', prefix='n'):
file_names = [rootDir + prefix + str(i).zfill(4) + '.fits' for i in files]
return file_names
[docs]
def get_distortion_maps(self, hdr):
"""
Inputs
----------
date : str
Date in string format such as '2015-10-02'.
"""
date = hdr['DATE-OBS']
if (float(date[0:4]) < 2015):
distXgeoim = module_dir + '/reduce/distortion/nirc2_narrow_xgeoim.fits'
distYgeoim = module_dir + '/reduce/distortion/nirc2_narrow_ygeoim.fits'
if (float(date[0:4]) == 2015) & (float(date[5:7]) < 0o5):
distXgeoim = module_dir + '/reduce/distortion/nirc2_narrow_xgeoim.fits'
distYgeoim = module_dir + '/reduce/distortion/nirc2_narrow_ygeoim.fits'
if (float(date[0:4]) == 2015) & (float(date[5:7]) >= 0o5):
distXgeoim = module_dir + '/reduce/distortion/nirc2_narrow_xgeoim_post20150413.fits'
distYgeoim = module_dir + '/reduce/distortion/nirc2_narrow_ygeoim_post20150413.fits'
if (float(date[0:4]) > 2015):
distXgeoim = module_dir + '/reduce/distortion/nirc2_narrow_xgeoim_post20150413.fits'
distYgeoim = module_dir + '/reduce/distortion/nirc2_narrow_ygeoim_post20150413.fits'
return distXgeoim, distYgeoim
[docs]
def get_align_type(self, hdr, errors=False):
# Setup NIRC2 plate scales
# Units are arcsec/pixel
atypes = {"narrow": 8,
"medium": 14,
"wide": 12}
atype = atypes[hdr['CAMNAME']]
if errors == True:
atype += 1
return atype
[docs]
def get_saturation_level(self):
"""
Set to the 95% saturation threshold in DN.
"""
return 12000.0
[docs]
def get_linearity_correction_coeffs(self):
"""
Returns coefficients (`coeffs`, as defined below)
in order to perform linearity correction
x = (FITS_orig) / (No. of coadds)
norm = coeffs[0] + coeffs[1] * x + coeffs[2] * x^2
FITS_corrected = FITS_orig / norm
From Stanimir Metchev's linearity correction code
(http://www.astro.sunysb.edu/metchev/ao.html)
"""
lin_corr_coeffs = np.array([1.001, -6.9e-6, -0.70e-10])
return lin_corr_coeffs
[docs]
class OSIRIS(Instrument):
"""
OSIRIS Imager - after 2019
"""
def __init__(self):
# Define
self.hdr_keys['filename'] = 'datafile'
self.hdr_keys['object_name'] = 'object'
self.hdr_keys['itime'] = 'truitime'
self.hdr_keys['coadds'] = 'coadds'
self.hdr_keys['sampmode'] = 'sampmode'
self.hdr_keys['nfowler'] = 'numreads'
self.hdr_keys['camera'] = 'instr'
self.hdr_keys['shutter'] = 'ifilter'
self.hdr_keys['mjd'] = 'MJD-OBS'
self.hdr_keys['elevation'] = 'EL'
[docs]
self.bad_pixel_mask = 'osiris_img_mask.fits'
[docs]
self.distXgeoim = module_dir + '/reduce/distortion/OSIRIS_im_x_2021.fits'
[docs]
self.distYgeoim = module_dir + '/reduce/distortion/OSIRIS_im_y_2021.fits'
[docs]
self.telescope = 'Keck'
[docs]
self.telescope_diam = 10.5 # telescope diameter in meters
return
[docs]
def get_filter_name(self, hdr):
f = hdr['ifilter']
return f.split('-')[0]
[docs]
def get_plate_scale(self, hdr):
"""
Return the plate scale in arcsec/pix.
"""
# scale = 0.00995
date = hdr['DATE-OBS']
if (float(date[0:4]+date[5:7]+date[8:10]) < 20201116):
scale = 0.0099418
elif (float(date[0:4]+date[5:7]+date[8:10]) >= 20201116):
scale = 0.0099576
return scale
[docs]
def get_pcu_scale(self,hdr):
"""
The scale on the PCU pinhole mask, in mm (at the front focus) per pixel
"""
scale = 1/138.5
return scale
[docs]
def get_pcuxyRef(self,hdr):
"""
The reference position, when the PCU is on-axis
"""
pcuxyRef = [90,185]
return pcuxyRef
[docs]
def get_position_angle(self, hdr):
"""
Get the sky PA in degrees East of North.
"""
# pa = float(hdr['ROTPOSN']) - self.get_instrument_angle(hdr)
pa = hdr['PA_IMAG']
#if in PCU mode, read the PCU rotation angle instead
if 'PCUZ' in hdr.keys():
if float(hdr['PCUZ']) > 20.0:
pcu_angle = float(hdr['PCUR'])
pinhole_angle = 65.703 #the angle at which the pihole mask is horizontal.
# rotator_angle = hdr['ROTPPOSN']
# default_rotator_angle = 90
pa = pcu_angle - pinhole_angle
return pa
[docs]
def get_instrument_angle(self, hdr):
"""
Get the angle of the instrument w.r.t. to the telescope or
AO bench in degrees.
"""
inst_angle = (hdr['INSTANGL'] - 42.5)
return inst_angle
[docs]
def get_parallactic_angle(self,hdr):
"""
Get the parallactic angle in degrees East of North
"""
q = hdr['PARANG']
return q
[docs]
def get_central_wavelength(self, hdr):
"""
Return the central wavelength of the filter for
this observation in microns.
"""
filt_name = hdr['IFILTER']
# These are very approximate for now.
wave_dict = {'Kp-LHex': 2.12,
'Kn3-LHex': 2.12,
'Kcont-LHex': 2.270,
'Kcont': 2.270,
'Hbb-LHex': 1.65,
'Drk': 0.00,
'Kp': 2.12,
'Kp-sHex': 2.12,
'Kn3': 2.12,
'Kn3-sHex': 2.12,
'Hbb': 1.65,
'Hbb-LAnn':1.65,
'Hn3': 1.635,
'Hcont':1.5832,
'BrGamma':2.169,
'BrGamma-sAnn':2.169,
}
if filt_name not in wave_dict.keys():
print('NO information available on this filter: ' + filt_name)
return 2.12
else:
return wave_dict[filt_name]
[docs]
def get_gain(self, hdr):
return hdr['DETGAIN']
[docs]
def make_filenames(self, files, rootDir='', prefix=''):
file_names = [rootDir + prefix + i + '.fits' for i in files]
return file_names
[docs]
def flip_images(self, files, rootDir=''):
"""
Flip images (as they come from the detector flipped) and
subtract reference pixels.
"""
for ff in range(len(files)):
old_file = files[ff]
new_file = files[ff].replace('.fits', '_flip.fits')
hdu_list = fits.open(old_file)
# Fetch the date and figure out how to
# best flip the images.
year = int(hdu_list[0].header['DATE-OBS'].split('-')[0])
month = int(hdu_list[0].header['DATE-OBS'].split('-')[1])
date = int(hdu_list[0].header['DATE-OBS'].split('-')[2])
for hh in range(len(hdu_list)):
if isinstance(hdu_list[hh], _ImageBaseHDU):
# Subtract the reference pixels
new_data = self.subtract_reference_pixels(hdu_list[hh].data)
# if year == 2019:
# hdu_list[hh].data = new_data[:, ::-1]
# else:
# hdu_list[hh].data = new_data[::-1, :]
hdu_list[hh].data = new_data[::-1, :]
# Need to modify PA_IMAG to account for the flip. Added 2023-10-30 by J. Lu.
pa_orig = hdu_list[0].header['PA_IMAG']
hdu_list[0].header['PA_IMAG'] = 360.0 - pa_orig
hdu_list.writeto(new_file, overwrite=True)
# Add header values.
wave = self.get_central_wavelength(hdu_list[0].header)
fits.setval(new_file, 'EFFWAVE', value= wave)
fits.setval(new_file, 'CENWAVE', value= wave)
fits.setval(new_file, 'CAMNAME', value = 'narrow') # from NIRC2
return
[docs]
def subtract_reference_pixels(self, img):
horiz_ref_pixels = np.concatenate([img[:, 0:4], img[:, -4:]], axis=1)
ref_pix_median = np.median(horiz_ref_pixels, axis=1)
new_img = img - np.array([ref_pix_median]).T
return new_img
[docs]
def get_distortion_maps(self, hdr):
"""
Inputs
----------
date : str
Date in string format such as '2015-10-02'.
"""
date = hdr['DATE-OBS']
if (float(date[0:4]+date[5:7]+date[8:10]) < 20201116):
self.distXgeoim = module_dir + '/reduce/distortion/OSIRIS_im_x_2020.fits'
self.distYgeoim = module_dir + '/reduce/distortion/OSIRIS_im_y_2020.fits'
elif (float(date[0:4]+date[5:7]+date[8:10]) >= 20201116):
self.distXgeoim = module_dir + '/reduce/distortion/OSIRIS_im_x_2021.fits'
self.distYgeoim = module_dir + '/reduce/distortion/OSIRIS_im_y_2021.fits'
return self.distXgeoim, self.distYgeoim
[docs]
def get_align_type(self, hdr, errors=False):
atype = 14
if errors == True:
atype += 1
return atype
[docs]
def get_saturation_level(self):
"""
Set to the 95% saturation threshold in DN.
"""
return 17000.0
##################################################
#
# SET DEFAULT INSTRUMENT FOR MODULE.
#
##################################################