Updates to the instruments and reduction code to account for nirc2 electronics upgrade in 2024.

kai/instruments.py

@@ -5,6 +5,7 @@
 from astropy.io import fits
 import pdb
 from astropy.io.fits.hdu.image import _ImageBaseHDU
+from astropy.time import Time
 
 module_dir = os.path.dirname(__file__)
 
@@ -52,15 +53,26 @@ def get_distortion_maps(self, date):
     def get_align_type(self, errors=False):
         pass
 
-    def get_saturation_level(self):
+    def get_saturation_level(self, hdr):
         pass
 
-    def get_lin_corr_coeffs(self):
+    def get_linearity_correction_coeffs(self, hdr):
+        pass
+
+    def get_radec(self, hdr):
+        pass
+
+    def get_mjd(self, hdr):
         pass
 
 
 class NIRC2(Instrument):
     def __init__(self):
+        """
+        NIRC2 Instrument object. During initialization, the observation year
+        can be specified in order to help KAI account for NIRC2 instrument
+        upgrade(s).
+        """    
         self.name = 'NIRC2'
 
         # Define header keywords
@@ -74,7 +86,6 @@ def __init__(self):
         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'
 
         self.bad_pixel_mask = 'nirc2mask.fits'
@@ -89,13 +100,16 @@ def __init__(self):
         return
 
     def get_filter_name(self, hdr):
-        filter1 = hdr['fwiname']
-        filter2 = hdr['fwoname']
-        filt = filter1
-        if (filter1.startswith('PK')):
-            filt = filter2
-
-        return filt
+        if 'fwiname' in hdr:
+            filter1 = hdr['fwiname']
+            filter2 = hdr['fwoname']
+            filt = filter1
+            if (filter1.startswith('PK')):
+                filt = filter2
+
+            return filt
+        else:
+            return 'None'
 
     def get_plate_scale(self, hdr):
         """
@@ -179,13 +193,18 @@ def get_align_type(self, hdr, errors=False):
 
         return atype
 
-    def get_saturation_level(self):
+    def get_saturation_level(self, hdr):
         """
         Set to the 95% saturation threshold in DN.
         """
-        return 12000.0
+        date = hdr['DATE-OBS']
+
+        if (float(date[0:4]) >= 2024):
+            return 6000.0
+        else:
+            return 12000.0
 
-    def get_linearity_correction_coeffs(self):
+    def get_linearity_correction_coeffs(self, hdr):
         """
         Returns coefficients (`coeffs`, as defined below)
         in order to perform linearity correction
@@ -201,8 +220,56 @@ def get_linearity_correction_coeffs(self):
         """
 
         lin_corr_coeffs = np.array([1.001, -6.9e-6, -0.70e-10])
+
+        # Post 2024 upgrade to NIRC2 electronics, gain is ~2x compared to the
+        # values calculated by Metchev+ (i.e. gain is now ~8 e-/DN rather ~4
+        # e-/DN, and so nonlinearity now starts at ~4000 DN/coadd rather than
+        # ~8000 DN/coadd)
+        # Therefore, x can be multiplied by 2 for the same coeffs,
+        # or equivalently: coeffs can be multiplied by 2^n: [2^0, 2^1, 2^2]
+
+        date = hdr['DATE-OBS']
+        if (float(date[0:4]) >= 2024):
+            lin_corr_coeffs = np.array([
+                1.001 * 2**0,
+                -6.9e-6 * 2**1,
+                -0.70e-10 * 2**2,
+            ])
+
         return lin_corr_coeffs
+
+    def get_radec(self, hdr):
+        """Return list of RA and Dec as decimal degrees"""
 
+        date = hdr['DATE-OBS']
+
+        if (float(date[0:4]) < 2024):
+            # Previous to 2023/2024 electronics upgrade, RA and DEC stored
+            # as decimal degrees
+            return [float(hdr['RA']), float(hdr['DEC'])]
+        else:
+            # New coordinates are in HH:MM:SS.SSS or DD:MM:SS.SSS, so have
+            # to convert to decimals degrees
+            RA_split = hdr['RA'].split(':')
+            RA_decimal_hrs = float(RA_split[0]) +\
+                float(RA_split[1])/60. + float(RA_split[2])/3600.
+            RA_decimal_degs = (RA_decimal_hrs / 24.) * 360.
+
+            DEC_split = hdr['RA'].split(':')
+            DEC_decimal_degs = float(DEC_split[0]) +\
+                float(DEC_split[1])/60. + float(DEC_split[2])/3600.
+
+            return [RA_decimal_degs, DEC_decimal_degs]
+
+    def get_mjd(self, hdr):
+        """Return observation time in MJD"""
+
+        date = hdr['DATE-OBS']
+
+        if (float(date[0:4]) < 2024):
+            return float(hdr['MJD-OBS'])
+        else:
+            return float(hdr['MJD'])
 
 
 class OSIRIS(Instrument):
@@ -414,11 +481,21 @@ def get_align_type(self, hdr, errors=False):
 
         return atype
 
-    def get_saturation_level(self):
+    def get_saturation_level(self, hdr):
         """
         Set to the 95% saturation threshold in DN.
         """
         return 17000.0
+
+    def get_radec(self, hdr):
+        """Return list of RA and Dec as decimal degrees"""
+
+        return [float(hdr['RA']), float(hdr['DEC'])]
+
+    def get_mjd(self, hdr):
+        """Return observation time in MJD"""
+
+        return float(hdr['MJD-OBS'])
 
 ##################################################
 #

kai/reduce/data.py

@@ -184,7 +184,7 @@ def clean(
         # This is the image for which refSrc is relevant.
         firstFile = instrument.make_filenames([files[0]], rootDir=rawDir)[0]
         hdr1 = fits.getheader(firstFile, ignore_missing_end=True)
-        radecRef = [float(hdr1['RA']), float(hdr1['DEC'])]
+        radecRef = instrument.get_radec(hdr1)
         aotsxyRef = kai_util.getAotsxy(hdr1)
 
         if ref_offset_method == 'pcu':
@@ -320,7 +320,9 @@ def clean(
             clean_drizzle(distXgeoim, distYgeoim, _bp, _ce, _wgt, _dlog,
                           fixDAR=fixDAR, instrument=instrument,
                           use_koa_weather=use_koa_weather)
-
+
+            hdr = fits.getheader(_raw, ignore_missing_end=True)
+
             ### Make .max file ###
             # Determine the non-linearity level. Raw data level of
             # non-linearity is 12,000 but we subtracted
@@ -330,7 +332,7 @@ def clean(
             nonlinSky = skyObj.getNonlinearCorrection(sky)
 
             coadds = fits.getval(_ss, instrument.hdr_keys['coadds'])
-            satLevel = (coadds*instrument.get_saturation_level()) - nonlinSky - bkg
+            satLevel = (coadds*instrument.get_saturation_level(hdr)) - nonlinSky - bkg
             file(_max, 'w').write(str(satLevel))
 
             ### Rename and clean up files ###
@@ -339,7 +341,6 @@ def clean(
 
             ### Make the *.coo file and update headers ###
             # First check if PA is not zero
-            hdr = fits.getheader(_raw, ignore_missing_end=True)
             phi = instrument.get_position_angle(hdr)
 
             clean_makecoo(_ce, _cc, refSrc, strSrc, aotsxyRef, radecRef,
@@ -1148,7 +1149,7 @@ def combine_drizzle(imgsize, cleanDir, roots, outroot, weights, shifts,
 
 
         # Read in MJD of current file from FITS header
-        mjd = float(hdr['MJD-OBS'])
+        mjd = instrument.get_mjd(hdr)
         mjd_weightedSum += weights[i] * mjd
 
         # Drizzle this file ontop of all previous ones.
@@ -1231,6 +1232,10 @@ def combine_drizzle(imgsize, cleanDir, roots, outroot, weights, shifts,
         'MJD-OBS', mjd_weightedMean,
         'Weighted modified julian date of combined observations'
     )
+    fits_f[0].header.set(
+        'MJD', mjd_weightedMean,
+        'Weighted modified julian date of combined observations'
+    )
 
     ## Also update date field in header
     fits_f[0].header.set(
@@ -1364,7 +1369,7 @@ def combine_submaps(
             ir.drizzle.ygeoim = distYgeoim
 
         # Read in MJD of current file from FITS header
-        mjd = float(hdr['MJD-OBS'])
+        mjd = instrument.get_mjd(hdr)
         mjd_weightedSums[sub] += weights[i] * mjd
 
         # Drizzle this file ontop of all previous ones.
@@ -1439,10 +1444,23 @@ def combine_submaps(
                               'Y Distortion Image')
 
         # Store weighted MJDs in header
-        fits_f[0].header.set('MJD-OBS', mjd_weightedMeans[s], 'Weighted modified julian date of combined observations')
-
+        fits_f[0].header.set(
+            'MJD-OBS',
+            mjd_weightedMeans[s],
+            'Weighted modified julian date of combined observations',
+        )
+        fits_f[0].header.set(
+            'MJD',
+            mjd_weightedMeans[s],
+            'Weighted modified julian date of combined observations',
+        )
+
         ## Also update date field in header
-        fits_f[0].header.set('DATE', '{0}'.format(submaps_time_obs[s].fits), 'Weighted observation date')
+        fits_f[0].header.set(
+            'DATE',
+            '{0}'.format(submaps_time_obs[s].fits),
+            'Weighted observation date',
+        )
 
         # Write out final submap fits file
         fits_f[0].writeto(_fits[s], output_verify=outputVerify)
@@ -1981,7 +1999,7 @@ def clean_makecoo(_ce, _cc, refSrc, strSrc, aotsxyRef, radecRef,
 
     hdr = fits.getheader(_ce, ignore_missing_end=True)
 
-    radec = [float(hdr['RA']), float(hdr['DEC'])]
+    radec = instrument.get_radec(hdr)
     aotsxy = kai_util.getAotsxy(hdr)
     if offset_method == 'pcu':
         pcuxy = [float(hdr['PCSFX']), float(hdr['PCSFY'])]  #New version may be PCUX and PCUY

kai/reduce/kai_util.py

@@ -58,7 +58,10 @@ def makelog(directory, outfile='image_log.txt',
             f.write('%16s  ' % frame)
 
             # Second column is object name
-            f.write('%-16s  ' % hdr[ihk['object_name']].replace(' ', ''))
+            if ihk['object_name'] in hdr:
+                f.write('%-16s  ' % hdr[ihk['object_name']].replace(' ', ''))
+            else:
+                f.write('%-16s  ' % '[No object name]')
 
             # Next is integration time, coadds, sampmode, multisam
             f.write('%8.3f  %3d  ' % (hdr[ihk['itime']], hdr[ihk['coadds']]))
@@ -70,10 +73,16 @@ def makelog(directory, outfile='image_log.txt',
             f.write('%-10s ' % filt)
 
             # Camera name
-            f.write('%-6s ' % hdr[ihk['camera']])
+            if ihk['camera'] in hdr:
+                f.write('%-6s ' % hdr[ihk['camera']])
+            else:
+                f.write('%-6s ' % 'None')
 
             # Shutter state
-            f.write('%-6s ' % hdr[ihk['shutter']])
+            if ihk['shutter'] in hdr:
+                f.write('%-6s ' % hdr[ihk['shutter']])
+            else:
+                f.write('%-6s ' % 'None')
 
             # TRICK dichroic (only for OSIRIS)
             if isinstance(instrument, instruments.OSIRIS):

kai/reduce/lin_correction.py

@@ -31,7 +31,7 @@ def lin_correction(file, instrument=instruments.default_inst):
     num_coadds = im_header['COADDS']
 
     x = im_data / num_coadds
-    coeffs = instrument.get_linearity_correction_coeffs()
+    coeffs = instrument.get_linearity_correction_coeffs(im_header)
 
     # Determine order of polynomial correction
     norm_poly_order = len(coeffs)