Implementing Calibrated Contrast Measurements

.gitignore

@@ -3,6 +3,9 @@ __pycache__/
 *.py[cod]
 *$py.class
 
+# Transmission mask files
+spaceKLIP/resources/transmissions/NRC*
+
 # C extensions
 *.so
 

requirements.txt

@@ -4,5 +4,6 @@ jdaviz
 jwst
 matplotlib
 numpy
+pysynphot
 scipy
 webbpsf

spaceKLIP/database.py

@@ -48,7 +48,7 @@
 # Load NIRCam, NIRISS, and MIRI filters
 wave_nircam, weff_nircam, do_svo = get_filter_info('NIRCAM', return_more=True)
 wave_niriss, weff_niriss = get_filter_info('NIRISS', do_svo=do_svo)
-wave_miri,   weff_miri   = get_filter_info('MIRI',   do_svo=do_svo)
+wave_miri,   weff_miri   = get_filter_info('MIRI',   do_svo=False) 
 
 class Database():
     """
@@ -1286,7 +1286,7 @@ def create_database(output_dir,
 
     Automatically searches for uncal.fits in the input directory and creates 
     a database of the JWST data. Only works for stage0, stage1, or stage2 data
-    by default; set readlevel3=True to read in stage 3 outputs.
+    by default; set readlevel=3 to read in stage 3 outputs.
 
     Parameters
     ----------
@@ -1309,6 +1309,11 @@ def create_database(output_dir,
         default is None.
     assoc_using_targname : bool, optional
         Associate observations using the TARGNAME keyword. The default is True.
+    readlevel : str or int
+        Level of data products to read in, either '012' (or an integer 0,1,2) to read in
+        individual exposures, or '3' to read in level-3 PSF-subtracted products,
+        or '4' to read in extracted PSF fitting products.
+
     verbose : bool, optional
         Print information to the screen. The default is True.
 
@@ -1362,13 +1367,10 @@ def create_database(output_dir,
     elif str(readlevel) == '3':
         # the above get_files usage won't match KLIP outputsa, so find them here
         datapaths_klip = sorted(glob.glob(os.path.join(input_dir, "*KLmodes-all.fits")))
-
         db.read_jwst_s3_data(datapaths=datapaths+datapaths_klip,
-                             verbose=verbose
                             )
     elif str(readlevel) == '4':
         db.read_jwst_s4_data(datapaths=datapaths,
-                             verbose=verbose
                             )
     else:
         raise ValueError("Invalid/unknown value for readlevel parameter")

spaceKLIP/imagetools.py

@@ -36,6 +36,8 @@
 from webbpsf_ext.coords import dist_image
 from webbpsf_ext.webbpsf_ext_core import _transmission_map
 
+from webbpsf.constants import JWST_CIRCUMSCRIBED_DIAMETER
+
 import logging
 log = logging.getLogger(__name__)
 log.setLevel(logging.INFO)
@@ -605,6 +607,7 @@ def subtract_background(self,
             set to 5, a median of every 5 background images will be subtracted from
             the corresponding 5 target images. The default is None (i.e. a median 
             across all images).
+
         subdir : str, optional
             Name of the directory where the data products shall be saved. The
             default is 'bgsub'.
@@ -659,6 +662,7 @@ def subtract_background(self,
                 sci_bg_data = np.concatenate(sci_bg_data)
                 sci_bg_erro = np.concatenate(sci_bg_erro)
                 sci_bg_pxdq = np.concatenate(sci_bg_pxdq)
+
                 sci_bg_data_split = np.array_split(sci_bg_data, split_inds, axis=0)
                 sci_bg_erro_split = np.array_split(sci_bg_erro, split_inds, axis=0)
                 sci_bg_pxdq_split = np.array_split(sci_bg_pxdq, split_inds, axis=0)
@@ -695,6 +699,7 @@ def subtract_background(self,
                 ref_bg_data = np.concatenate(ref_bg_data)
                 ref_bg_erro = np.concatenate(ref_bg_erro)
                 ref_bg_pxdq = np.concatenate(ref_bg_pxdq)
+
                 ref_bg_data_split = np.array_split(ref_bg_data, split_inds, axis=0)
                 ref_bg_erro_split = np.array_split(ref_bg_erro, split_inds, axis=0)
                 ref_bg_pxdq_split = np.array_split(ref_bg_pxdq, split_inds, axis=0)
@@ -1391,7 +1396,7 @@ def blur_frames(self,
                         if self.database.obs[key]['EXP_TYPE'][j] in ['NRC_CORON']:
                             diam = 5.2
                         else:
-                            diam = 6.5
+                            diam = JWST_CIRCUMSCRIBED_DIAMETER
                     else:
                         raise UserWarning('Data originates from unknown telescope')
                     if fact_temp is not None:
@@ -2094,8 +2099,8 @@ def align_frames(self,
                         log.warning('  --> The following frames might not be properly aligned: '+str(ww))
 
                 # Write FITS file and PSF mask.
-                head_pri['XOFFSET'] = xoffset
-                head_pri['YOFFSET'] = yoffset
+                head_pri['XOFFSET'] = xoffset #arcseconds
+                head_pri['YOFFSET'] = yoffset #arcseconds 
                 head_sci['CRPIX1'] = crpix1
                 head_sci['CRPIX2'] = crpix2
                 fitsfile = ut.write_obs(fitsfile, output_dir, data, erro, pxdq, head_pri, head_sci, is2d, imshifts, maskoffs)
@@ -2115,6 +2120,7 @@ def align_frames(self,
                            label='PA = %.0f deg' % self.database.obs[key]['ROLL_REF'][j])
             ax.axhline(0., color='gray', lw=1, zorder=-1)  # set zorder to ensure lines are drawn behind all the scatter points
             ax.axvline(0., color='gray', lw=1, zorder=-1)
+
             ax.set_aspect('equal')
             xlim = ax.get_xlim()
             ylim = ax.get_ylim()
@@ -2186,3 +2192,4 @@ def align_frames(self,
                 plt.savefig(output_file)
                 log.info(f" Plot saved in {output_file}")
                 plt.close()
+

spaceKLIP/make_psfmasks.py

@@ -109,7 +109,9 @@
 
         # Write PSF mask.
         hdul['SCI'].data = psfmask
-        hdul.writeto('resources/transmissions/' + apername + '_' + filt + '.fits', output_verify='fix', overwrite=True)
+        outfile = 'resources/transmissions/' + apername + '_' + filt + '.fits'
+        print(f"Wrote to {outfile}")
+        hdul.writeto(outfile, output_verify='fix', overwrite=True)
 
 fitsfiles = ['JWST_MIRI_F1065C_transmission_webbpsf-ext_v2.fits',
              'JWST_MIRI_F1140C_transmission_webbpsf-ext_v2.fits',
@@ -122,6 +124,7 @@
     hdu.header['EXTNAME'] = 'SCI'
     hdul.append(hdu)
     hdul.writeto('resources/transmissions/' + fitsfile, output_verify='fix', overwrite=True)
+    print(f"Wrote to resources/transmissions/{fitsfile}")
 
 
 def get_bar_offset_from_siaf(siaf, filt, channel='LW'):
@@ -198,7 +201,9 @@ def get_bar_offset_from_siaf(siaf, filt, channel='LW'):
 
     # Write PSF mask.
     hdul['SCI'].data = psfmask
-    hdul.writeto('resources/transmissions/' + apername + '_' + filt.upper() + '.fits', output_verify='fix', overwrite=True)
+    outfile = 'resources/transmissions/' + apername + '_' + filt.upper() + '.fits'
+    hdul.writeto(outfile, output_verify='fix', overwrite=True)
+    print(f"Wrote to {outfile}")
 
 # Loop through apertures and filters.
 fitsfile = 'jwst_nircam_psfmask_0003.fits'