better e2e test agreement with same file order, cosmic detection fix from latest II&T

corgidrp/darks.py

@@ -139,7 +139,6 @@ def build_trad_dark(dataset, detector_params, detector_regions=None, full_frame=
     dark current is so high that it stands far above other noise that is
     not smeared upon readout, such as clock-induced charge, 
     fixed-pattern noise, and read noise.
-    - have been divided by EM gain.
 
     Also, add_photon_noise() should NOT have been applied to the frames in
     dataset.  And note that creation of the
@@ -275,7 +274,6 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None):
     dark current is so high that it stands far above other noise that is
     not smeared upon readout, such as clock-induced charge, 
     fixed-pattern noise, and read noise.
-    - have been divided by EM gain.
 
     Also, add_photon_noise() should NOT have been applied to the frames in
     dataset.  And note that creation of the
@@ -737,7 +735,7 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None):
                            input_err, input_dq, err_hdr=err_hdr)
 
     l2a_data_filename = dataset.copy()[0].filename
-    noise_maps.filename =  l2a_data_filename[:24] + '_DetectorNoiseMaps.fits'
+    noise_maps.filename =  l2a_data_filename[:-5] + '_DetectorNoiseMaps.fits'
 
     return noise_maps
 

corgidrp/detector.py

@@ -377,7 +377,7 @@ def imaging_slice(obstype, frame, detector_regions=None):
     return sl
 
 def flag_cosmics(cube, fwc, sat_thresh, plat_thresh, cosm_filter, cosm_box,
-                   cosm_tail, mode='image'):
+                   cosm_tail, mode='image', detector_regions=None, obstype='SCI'):
     """Identify and remove saturated cosmic ray hits and tails.
 
     Use sat_thresh (interval 0 to 1) to set the threshold above which cosmics
@@ -422,6 +422,11 @@ def flag_cosmics(cube, fwc, sat_thresh, plat_thresh, cosm_filter, cosm_box,
             If 'full', a full-frame input is assumed, and if the input tail length
             is longer than the length to the end of the full-frame row, the masking
             continues onto the next row.  Defaults to 'image'.
+        detector_regions: (dict):  
+            A dictionary of detector geometry properties.  Keys should be as 
+            found in detector_areas in detector.py. Defaults to detector_areas in detector.py.
+        obstype (string):
+            'ARRTYPE' from the header associated with the input data cube.  Defaults to 'SCI'.
 
     Returns:
         array_like, int:
@@ -440,17 +445,35 @@ def flag_cosmics(cube, fwc, sat_thresh, plat_thresh, cosm_filter, cosm_box,
      ......|<-plateau->|<------------------tail---------->|.........
 
     B Nemati and S Miller - UAH - 02-Oct-2018
+    Kevin Ludwick - UAH - 2024
 
     """
     mask = np.zeros(cube.shape, dtype=int)
 
+    if detector_regions is None:
+        detector_regions = detector_areas
+
+    if mode=='full':
+        im_num_rows = detector_regions[obstype]['image']['rows']
+        im_num_cols = detector_regions[obstype]['image']['cols']
+        im_starting_row = detector_regions[obstype]['image']['r0c0'][0]
+        im_ending_row = im_starting_row + im_num_rows
+        im_starting_col = detector_regions[obstype]['image']['r0c0'][1]
+        im_ending_col = im_starting_col + im_num_cols
+    else:
+        im_starting_row = 0
+        im_ending_row = mask.shape[1] - 1 # - 1 to get the index, not size
+        im_starting_col = 0
+        im_ending_col = mask.shape[2] - 1 # - 1 to get the index, not size
+
     # Do a cheap prefilter for rows that don't have anything bright
     max_rows = np.max(cube, axis=-1,keepdims=True)
     ji_streak_rows = np.transpose(np.array((max_rows >= sat_thresh*fwc).nonzero()[:-1]))
 
     for j,i in ji_streak_rows:
         row = cube[j,i]
-
+        if i < im_starting_row or i > im_ending_row:
+            continue
         # Find if and where saturated plateaus start in streak row
         i_begs = find_plateaus(row, fwc, sat_thresh, plat_thresh, cosm_filter)
 
@@ -459,6 +482,8 @@ def flag_cosmics(cube, fwc, sat_thresh, plat_thresh, cosm_filter, cosm_box,
         ex_l = np.array([])
         if i_begs is not None:
             for i_beg in i_begs:
+                if i_beg < im_starting_col or i_beg > im_ending_col:
+                    continue
                 # implement cosm_tail
                 if i_beg+cosm_filter+cosm_tail+1 > mask.shape[2]:
                     ex_l = np.append(ex_l,
@@ -468,10 +493,11 @@ def flag_cosmics(cube, fwc, sat_thresh, plat_thresh, cosm_filter, cosm_box,
                                 mask.shape[2]))
                 mask[j, i, i_beg:streak_end] = 1
                 # implement cosm_box
-                st_row = max(i-cosm_box, 0)
-                end_row = min(i+cosm_box+1, mask.shape[1])
-                st_col = max(i_beg-cosm_box, 0)
-                end_col = min(i_beg+cosm_box+1, mask.shape[2])
+                # can't have cosm_box appear in non-image pixels
+                st_row = max(i-cosm_box, im_starting_row)
+                end_row = min(i+cosm_box+1, im_ending_row+1)
+                st_col = max(i_beg-cosm_box, im_starting_col)
+                end_col = min(i_beg+cosm_box+1, im_ending_col+1)
                 mask[j, st_row:end_row, st_col:end_col] = 1
                 pass
 

corgidrp/l1_to_l2a.py

@@ -156,7 +156,7 @@ def prescan_biassub(input_dataset, noise_maps=None, return_full_frame=False,
 
 def detect_cosmic_rays(input_dataset, detector_params, sat_thresh=0.7,
                        plat_thresh=0.7, cosm_filter=1, cosm_box=3, cosm_tail=10,
-                       mode='image'):
+                       mode='image', detector_regions=None):
     """
     Detects cosmic rays in a given dataset. Updates the DQ to reflect the pixels that are affected.
     TODO: (Eventually) Decide if we want to invest time in improving CR rejection (modeling and subtracting the hit
@@ -192,6 +192,9 @@ def detect_cosmic_rays(input_dataset, detector_params, sat_thresh=0.7,
             If 'full', a full-frame input is assumed, and if the input tail length
             is longer than the length to the end of the full-frame row, the masking
             continues onto the next row.  Defaults to 'image'.
+        detector_regions: (dict):  
+            A dictionary of detector geometry properties.  Keys should be as 
+            found in detector_areas in detector.py. Defaults to detector_areas in detector.py.
 
     Returns:
         corgidrp.data.Dataset:
@@ -200,12 +203,14 @@ def detect_cosmic_rays(input_dataset, detector_params, sat_thresh=0.7,
     sat_dqval = 32 # DQ value corresponding to full well saturation
     cr_dqval = 128 # DQ value corresponding to CR hit
 
+    if detector_regions is None:
+        detector_regions = detector_areas
+
     # you should make a copy the dataset to start
     crmasked_dataset = input_dataset.copy()
 
     crmasked_cube = crmasked_dataset.all_data
 
-
     # Calculate the full well capacity for every frame in the dataset
     kgain = np.array([detector_params.params['kgain'] for frame in crmasked_dataset])
     emgain_list = []
@@ -246,14 +251,17 @@ def detect_cosmic_rays(input_dataset, detector_params, sat_thresh=0.7,
     m2 = np.zeros_like(crmasked_cube)
 
     for i in range(len(crmasked_cube)):
+        obstype = crmasked_dataset.frames[i].ext_hdr['ARRTYPE']
         m2[i,:,:] = flag_cosmics(cube=crmasked_cube[i:i+1,:,:],
                         fwc=fwcem_dn_arr[i],
                         sat_thresh=sat_thresh,
                         plat_thresh=plat_thresh,
                         cosm_filter=cosm_filter,
                         cosm_box=cosm_box,
                         cosm_tail=cosm_tail,
-                        mode=mode
+                        mode=mode,
+                        detector_regions=detector_regions,
+                        obstype=obstype
                         ) * cr_dqval
 
     # add the two masks to the all_dq mask

corgidrp/recipe_templates/build_trad_dark_full_frame.json

@@ -40,9 +40,6 @@
                 "KGain" : "AUTOMATIC"
             }
         },
-        {
-            "name" : "em_gain_division"
-        },
         {
             "name" : "build_trad_dark",
             "calibs" : {

corgidrp/recipe_templates/build_trad_dark_image.json

@@ -40,9 +40,6 @@
                 "KGain" : "AUTOMATIC"
             }
         },
-        {
-            "name" : "em_gain_division"
-        },
         {
             "name" : "build_trad_dark",
             "calibs" : {

corgidrp/recipe_templates/l1_to_l2a_noisemap.json

@@ -17,6 +17,9 @@
             "name" : "detect_cosmic_rays",
             "calibs" : {
                 "DetectorParams" : "AUTOMATIC"
+            },
+            "keywords" : {
+                "mode" : "full"
             }
         },
         {

tests/e2e_tests/bp_map_e2e.py

@@ -255,6 +255,7 @@ def test_bp_map_simulated_dark_e2e(tvacdata_path, e2eoutput_path):
     # Create a dark object and save it
     master_dark = data.Dark(simple_dark_data, pri_hdr=pri_hdr, ext_hdr=ext_hdr,
                             input_dataset=mock_input_dataset)
+
     master_dark.save(filedir=bp_map_outputdir, filename="dark_mock.fits")
     this_caldb.create_entry(master_dark)
     master_dark_ref = master_dark.filepath
@@ -314,7 +315,8 @@ def test_bp_map_simulated_dark_e2e(tvacdata_path, e2eoutput_path):
 
 if __name__ == "__main__":
     # Set default paths and parse command-line arguments
-    tvacdata_dir = "/Users/jmilton/Documents/CGI/CGI_TVAC_Data"
+    # tvacdata_dir = "/Users/jmilton/Documents/CGI/CGI_TVAC_Data"
+    tvacdata_dir = "/Users/kevinludwick/Library/CloudStorage/Box-Box/CGI_TVAC_Data/"
     outputdir = thisfile_dir
 
     # Argument parser setup

tests/e2e_tests/kgain_e2e.py

@@ -10,13 +10,86 @@
 import corgidrp.walker as walker
 import corgidrp.caldb as caldb
 
+try:
+    from cal.kgain.calibrate_kgain import calibrate_kgain
+    import cal
+except:
+    pass
+
 thisfile_dir = os.path.dirname(__file__) # this file's folder
 
-tvac_kgain = 8.8145 #e/DN, result from new iit drp code
-tvac_readnoise = 130.12 #e, result from new iit drp code
+# tvac_kgain = 8.49404981510777 #8.8145 #e/DN, result from new iit code with specified file input order
+# tvac_readnoise = 121.76070832489948 # 130.12 #e, result from new iit code with specified file input order
 
 @pytest.mark.e2e
 def test_l1_to_kgain(tvacdata_path, e2eoutput_path):
+
+    # run II&T code 
+    default_config_file = os.path.join(cal.lib_dir, 'kgain', 'config_files', 'kgain_parms.yaml')
+    stack_arr2_f = []
+    stack_arr_f = []
+    box_data = os.path.join(tvacdata_path, 'TV-20_EXCAM_noise_characterization', 'kgain') 
+    for f in os.listdir(box_data):
+        file = os.path.join(box_data, f)
+        if not file.endswith('.fits'):
+            continue
+        for i in range(51841, 51871):
+            if str(i) in file:
+                stack_arr2_f.append(file)
+                with fits.open(file, mode='update') as hdus:
+                    try:
+                        hdus[0].header['OBSTYPE'] = 'MNFRAME'
+                    except:
+                        pass
+                    try:
+                        hdus[1].header['OBSTYPE'] = 'MNFRAME'
+                    except:
+                        pass
+                exit
+        for i in range(51731, 51841):
+            if str(i) in file:
+                stack_arr_f.append(file)
+                with fits.open(file, mode='update') as hdus:
+                    try:
+                        hdus[0].header['OBSTYPE'] = 'KGAIN'
+                    except:
+                        pass
+                    try:
+                        hdus[1].header['OBSTYPE'] = 'KGAIN'
+                    except:
+                        pass
+                exit
+    #stack_arr2 = np.stack(stack_arr2)
+    # fileorder_filepath = os.path.join(os.path.split(box_data)[0], 'results', 'TVAC_kgain_file_order.npy')
+    #np.save(fileorder_filepath, stack_arr_f+stack_arr2_f)
+    stack_arr_f = sorted(stack_arr_f)
+    stack_dat = data.Dataset(stack_arr_f)
+    stack2_dat = data.Dataset(stack_arr2_f)
+    stack_arr2 = stack2_dat.all_data
+
+    split, _ = stack_dat.split_dataset(exthdr_keywords=['EXPTIME'])
+    stack_arr = []
+    for dset in split:
+        # Breaking up the one set that has 10 frames at the same exptime instead of 5 like all the rest
+        # Making the set 2 separate sets of 5 each perhaps unfairly weights this exptime 
+        # which is doubly represented now, but doing this results in the same 8.8145 number from before
+        if dset.all_data.shape[0] == 10:
+            stack_arr.append(dset.all_data[:5])
+            # for ind in [5,6,7,8,9]:
+            #     fp = dset.frames[ind].filepath
+            #     stack_arr_f.remove(fp)
+
+            stack_arr.append(dset.all_data[5:])
+            continue
+        stack_arr.append(dset.all_data)
+    stack_arr = np.stack(stack_arr)
+    pass
+    ordered_filelist = stack_arr_f+stack_arr2_f
+
+    (tvac_kgain, tvac_readnoise, mean_rn_std_e, ptc) = calibrate_kgain(stack_arr, stack_arr2, emgain=1, min_val=800, max_val=3000, 
+                    binwidth=68, config_file=default_config_file, 
+                    mkplot=None, verbose=None)
+
     # figure out paths, assuming everything is located in the same relative location
     l1_datadir = os.path.join(tvacdata_path, "TV-20_EXCAM_noise_characterization", "kgain")
 
@@ -26,29 +99,11 @@ def test_l1_to_kgain(tvacdata_path, e2eoutput_path):
         shutil.rmtree(kgain_outputdir)
     os.mkdir(kgain_outputdir)
 
-    # define the raw science data to process
-
-    l1_data_filelist_same_exp = [os.path.join(l1_datadir, "CGI_EXCAM_L1_00000{0}.fits".format(i)) for i in np.arange(51841,51871)]#[51841, 51851]] # just grab some files of it
-    l1_data_filelist_range_exp = [os.path.join(l1_datadir, "CGI_EXCAM_L1_00000{0}.fits".format(i)) for i in np.arange(51731, 51841)]#[51731, 51761]]
-
     ####### Run the walker on some test_data
-    for file in l1_data_filelist_range_exp:
-        image = data.Image(file)
-        # This should not be necessary anymore after the updates of the OBSTYPE keyword, up to now it is only "SCI"
-        if image.pri_hdr['OBSTYPE'] != 'KGAIN':
-            image.pri_hdr['OBSTYPE'] = 'KGAIN'
-            image.save(filename = file)
-
-    for file in l1_data_filelist_same_exp:
-        image = data.Image(file)
-        # This should not be necessary anymore after the updates of the OBSTYPE keyword, up to now it is only "SCI"
-        if image.pri_hdr['OBSTYPE'] != 'MNFRAME':
-            image.pri_hdr['OBSTYPE'] = 'MNFRAME'
-            image.save(filename = file)
-        l1_data_filelist_range_exp.append(file)
-
-    walker.walk_corgidrp(l1_data_filelist_range_exp, "", kgain_outputdir)
-    kgain_file = os.path.join(kgain_outputdir, "CGI_EXCAM_L1_0000051731_kgain.fits")
+
+    walker.walk_corgidrp(ordered_filelist, "", kgain_outputdir, template="l1_to_kgain.json")
+    kgain_file = os.path.join(kgain_outputdir, os.path.split(ordered_filelist[0])[1][:-5]+'_kgain.fits') #"CGI_EXCAM_L1_0000051731_kgain.fits")
+
     kgain = data.KGain(kgain_file)
 
     ##### Check against TVAC kgain, readnoise
@@ -64,8 +119,8 @@ def test_l1_to_kgain(tvacdata_path, e2eoutput_path):
     print ("error of kgain:", kgain.error)
     print ("error of readnoise:", kgain.ext_hdr["RN_ERR"])
 
-    assert np.abs(diff_kgain) < 0.01
-    assert np.abs(diff_readnoise) < 3
+    assert np.abs(diff_kgain) == 0 #< 0.01
+    assert np.abs(diff_readnoise) == 0 #< 3
 
     this_caldb = caldb.CalDB()
     this_caldb.remove_entry(kgain)
@@ -77,7 +132,7 @@ def test_l1_to_kgain(tvacdata_path, e2eoutput_path):
     # to edit the file. The arguments use the variables in this file as their
     # defaults allowing the use to edit the file if that is their preferred
     # workflow.
-    tvacdata_dir = "/home/schreiber/DataCopy/corgi/CGI_TVAC_Data/"
+    tvacdata_dir = '/Users/kevinludwick/Library/CloudStorage/Box-Box/CGI_TVAC_Data/Working_Folder'  #"/home/schreiber/DataCopy/corgi/CGI_TVAC_Data/"
     outputdir = thisfile_dir
 
     ap = argparse.ArgumentParser(description="run the l1->kgain end-to-end test")

tests/e2e_tests/l1_to_l2b_e2e.py

@@ -157,7 +157,7 @@ def test_l1_to_l2b(tvacdata_path, e2eoutput_path):
     # to edit the file. The arguments use the variables in this file as their
     # defaults allowing the use to edit the file if that is their preferred
     # workflow.
-    tvacdata_dir = "/Users/jmilton/Documents/CGI/CGI_TVAC_Data"
+    tvacdata_dir =  "/Users/kevinludwick/Library/CloudStorage/Box-Box/CGI_TVAC_Data/Working_Folder/"
     outputdir = thisfile_dir
 
     ap = argparse.ArgumentParser(description="run the l1->l2a end-to-end test")