Flat calibration

corgidrp/data.py

@@ -158,6 +158,31 @@ def add_error_term(self, input_error, err_name):
         self.all_err = np.array([frame.err for frame in self.frames])   
         for i, frame in enumerate(self.frames):
             frame.err = self.all_err[i]
+
+    def rescale_error(self, input_error, err_name):
+        """
+        Calls Image.rescale_errors() for each frame.
+        Updates Dataset.all_err
+
+        Args:
+          input_error (np.array): 2-d error layer or 3-d layer
+          err_name (str): name of the uncertainty layer
+        """
+        if input_error.ndim == 3:
+            for i,frame in enumerate(self.frames):
+                frame.rescale_error(input_error[i], err_name)
+
+        elif input_error.ndim ==2:
+            for frame in self.frames:
+                frame.rescale_error(input_error, err_name)
+
+        else:
+            raise ValueError("input_error is not either a 2D or 3D array.")
+
+        # Preserve pointer links between Dataset.all_err and Image.err
+        self.all_err = np.array([frame.err for frame in self.frames])   
+        for i, frame in enumerate(self.frames):
+            frame.err = self.all_err[i]               
 
 class Image():
     """
@@ -440,6 +465,39 @@ def add_error_term(self, input_error, err_name):
 
         # record history since 2-D error map doesn't track individual terms
         self.err_hdr['HISTORY'] = "Added error term: {0}".format(err_name)
+
+    def rescale_error(self, input_error, err_name):
+        """
+        Add a layer of a specific additive uncertainty on the 3-dim error array extension
+        and update the combined uncertainty in the first layer.
+        Update the error header and assign the error name. 
+
+        Only tracks individual errors if the "track_individual_errors" setting is set to True
+        in the configuration file
+
+        Args:
+          input_error (np.array): 2-d error layer
+          err_name (str): name of the uncertainty layer
+        """
+        if input_error.ndim != 2 or input_error.shape != self.data.shape:
+            raise ValueError("we expect a 2-dimensional error layer with dimensions {0}".format(self.data.shape))
+
+        #first layer is always the updated combined error
+        self.err = self.err*input_error
+        self.err_hdr["Layer_1"] = "combined_error"
+
+        if corgidrp.track_individual_errors:
+            #append new error as layer on 3D cube
+            self.err=np.append(self.err, [input_error], axis=0)
+
+            layer = str(self.err.shape[0])
+            #self.err_hdr["Layer_1"] = "combined_error"
+            self.err_hdr["Layer_" + layer] = err_name    
+
+        # record history since 2-D error map doesn't track individual terms
+        self.err_hdr['HISTORY'] = "Errors rescaled by: {0}".format(err_name)    
+
+
 
     def get_hash(self):
         """
@@ -497,6 +555,45 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, input_dataset=N
         if 'DATATYPE' not in self.ext_hdr or self.ext_hdr['DATATYPE'] != 'Dark':
             raise ValueError("File that was loaded was not a Dark file.")
 
+class FlatField(Image):
+    """
+    Master flat generated from raster scan of uranus or Neptune.
+
+     Args:
+        data_or_filepath (str or np.array): either the filepath to the FITS file to read in OR the 2D image data
+        pri_hdr (astropy.io.fits.Header): the primary header (required only if raw 2D data is passed in)
+        ext_hdr (astropy.io.fits.Header): the image extension header (required only if raw 2D data is passed in)
+        input_dataset (corgidrp.data.Dataset): the Image files combined together to make this flat file (required only if raw 2D data is passed in)
+    """
+    def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, input_dataset=None):
+        # run the image class contructor
+        super().__init__(data_or_filepath, pri_hdr=pri_hdr, ext_hdr=ext_hdr)
+
+        # if this is a new master flat, we need to bookkeep it in the header
+        # b/c of logic in the super.__init__, we just need to check this to see if it is a new masterflat
+        if ext_hdr is not None:
+            if input_dataset is None:
+                # error check. this is required in this case
+                raise ValueError("This appears to be a master flat. The dataset of input files needs to be passed in to the input_dataset keyword to record history of this flat")
+            self.ext_hdr['DATATYPE'] = 'FlatField' # corgidrp specific keyword for saving to disk
+
+            # log all the data that went into making this flat
+            self._record_parent_filenames(input_dataset)
+
+            # add to history
+            self.ext_hdr['HISTORY'] = "Flat with exptime = {0} s created from {1} frames".format(self.ext_hdr['EXPTIME'], self.ext_hdr['DRPNFILE'])
+
+            # give it a default filename using the first input file as the base
+            orig_input_filename = input_dataset[0].filename.split(".fits")[0]
+            self.filename = "{0}_flatfield.fits".format(orig_input_filename)
+
+
+        # double check that this is actually a masterflat file that got read in
+        # since if only a filepath was passed in, any file could have been read in
+        if 'DATATYPE' not in self.ext_hdr or self.ext_hdr['DATATYPE'] != 'FlatField':
+            raise ValueError("File that was loaded was not a FlatField file.")
+
+
 class NonLinearityCalibration(Image):
     """
     Class for non-linearity calibration files. Although it's not stricly an image that you might look at, it is a 2D array of data
@@ -834,6 +931,7 @@ def get_hash(self):
               "NonLinearityCalibration" : NonLinearityCalibration,
               "KGain" : KGain, 
               "BadPixelMap" : BadPixelMap,
+              "FlatField" : FlatField,
               "DetectorParams" : DetectorParams }
 
 def autoload(filepath):
@@ -868,4 +966,4 @@ def autoload(filepath):
     # use the class constructor to load in the data
     frame = data_class(filepath)
 
-    return frame
+    return frame

corgidrp/detector.py

@@ -28,6 +28,33 @@ def create_dark_calib(dark_dataset):
 
     return new_dark
 
+def create_flatfield(flat_dataset):
+
+    """
+    Turn this dataset of image frames that were taken for performing the flat calibration and
+    to make one master flat image
+    this is currently a placeholder, until the final flat fielding calibration package is completed. 
+
+    Args:
+       flat_dataset (corgidrp.data.Dataset): a dataset of Image frames (L2a-level)
+
+    Returns:
+        data.masterflat: a master flat for flat calibration
+    """  
+
+
+    combined_frame = np.nanmean(flat_dataset.all_data, axis=0)
+
+    flat_field = data.FlatField(combined_frame, pri_hdr=flat_dataset[0].pri_hdr.copy(),
+                         ext_hdr=flat_dataset[0].ext_hdr.copy(), input_dataset=flat_dataset)
+
+    #determine the standard error of the mean: stddev/sqrt(n_frames)
+    flat_field.err = np.nanstd(flat_dataset.all_data, axis=0)/np.sqrt(len(flat_dataset))
+    flat_field.err=flat_field.err.reshape((1,)+flat_field.err.shape) # Get it into the right dimension
+
+
+    return flat_field
+
 def get_relgains(frame, em_gain, non_lin_correction):
     """
     For a given bias subtracted frame of dn counts, return a same sized

corgidrp/l2a_to_l2b.py

@@ -57,8 +57,42 @@ def dark_subtraction(input_dataset, dark_frame):
     darksub_dataset.update_after_processing_step(history_msg, new_all_data=darksub_cube, header_entries = {"BUNIT":"photoelectrons"})
 
     return darksub_dataset
+
+def flat_division(input_dataset, flat_field):
+    """
+
+    Divide the dataset by the master flat field. 
+
+    Args:
+        input_dataset (corgidrp.data.Dataset): a dataset of Images (L2a-level)
+        flat_field (corgidrp.data.FlatField): a master flat field to divide by
+
+    Returns:
+        corgidrp.data.Dataset: a version of the input dataset with the flat field divided out
+    """
+
+     # copy of the dataset
+    flatdiv_dataset = input_dataset.copy()
+
+    #Divide by the master flat
+    flatdiv_cube = flatdiv_dataset.all_data /  flat_field.data
+
+    # propagate the error of the master flat frame  
+    if hasattr(flat_field, "err"):
+        flatdiv_dataset.rescale_error(1/flat_field.data, "FlatField") 
+        flatdiv_dataset.add_error_term(flatdiv_dataset.all_data*flat_field.err[0]/(flat_field.data**2), "FlatField_error")
+    else:
+        raise Warning("no error attribute in the FlatField")
+
+    history_msg = "Flat calibration done using Flat field {0}".format(flat_field.filename)
+
+    # update the output dataset with this new flat calibrated data and update the history
+    flatdiv_dataset.update_after_processing_step(history_msg,new_all_data=flatdiv_cube)
+
+    return flatdiv_dataset
 
 def frame_select(input_dataset, bpix_frac=100., overexp=False, tt_thres=None):
+
     """
 
     Selects the frames that we want to use for further processing.
@@ -155,21 +189,6 @@ def cti_correction(input_dataset):
 
     return input_dataset.copy()
 
-def flat_division(input_dataset, master_flat):
-    """
-
-    Divide the dataset by the master flat field. 
-
-    Args:
-        input_dataset (corgidrp.data.Dataset): a dataset of Images (L2a-level)
-        master_flat (corgidrp.data.Flat): a master flat field to divide by
-
-    Returns:
-        corgidrp.data.Dataset: a version of the input dataset with the flat field divided out
-    """
-
-    return input_dataset.copy()
-
 def desmear(input_dataset, detector_params):
     """
 

corgidrp/mocks.py

@@ -37,6 +37,65 @@ def create_dark_calib_files(filedir=None, numfiles=10):
     dataset = data.Dataset(frames)
     return dataset
 
+def create_simflat_dataset(filedir=None, numfiles=10):
+    """
+    Create simulated data to check the flat division
+
+    Args:
+        filedir (str): (Optional) Full path to directory to save to.
+        numfiles (int): Number of files in dataset.  Defaults to 10.
+
+    Returns:
+        corgidrp.data.Dataset:
+        The simulated dataset
+    """
+    # Make filedir if it does not exist
+    if (filedir is not None) and (not os.path.exists(filedir)):
+        os.mkdir(filedir)
+
+    filepattern = "sim_flat_{0:04d}.fits"
+    frames = []
+    for i in range(numfiles):
+        prihdr, exthdr = create_default_headers()
+        # generate images in normal distribution with mean 1 and std 0.01
+        np.random.seed(456+i); sim_data = np.random.poisson(lam=150., size=(1024, 1024)).astype(np.float64)
+        frame = data.Image(sim_data, pri_hdr=prihdr, ext_hdr=exthdr)
+        if filedir is not None:
+            frame.save(filedir=filedir, filename=filepattern.format(i))
+        frames.append(frame)
+    dataset = data.Dataset(frames)
+    return dataset
+
+def create_flatfield_dummy(filedir=None, numfiles=2):
+
+    """
+    Turn this flat field dataset of image frames that were taken for performing the flat calibration and
+    to make one master flat image
+
+    Args:
+        filedir (str): (Optional) Full path to directory to save to.
+        numfiles (int): Number of files in dataset.  Defaults to 1 to create the dummy flat can be changed to any number
+
+    Returns:
+        corgidrp.data.Dataset: 
+        a set of flat field images 
+    """
+    ## Make filedir if it does not exist
+    if (filedir is not None) and (not os.path.exists(filedir)):
+        os.mkdir(filedir)
+
+    filepattern= "flat_field_{0:01d}.fits"
+    frames=[]
+    for i in range(numfiles):
+        prihdr, exthdr = create_default_headers()
+        np.random.seed(456+i); sim_data = np.random.normal(loc=1.0, scale=0.01, size=(1024, 1024))
+        frame = data.Image(sim_data, pri_hdr=prihdr, ext_hdr=exthdr)
+        if filedir is not None:
+            frame.save(filedir=filedir, filename=filepattern.format(i))
+        frames.append(frame)
+    flatfield = data.Dataset(frames)
+    return flatfield
+
 def create_nonlinear_dataset(filedir=None, numfiles=2,em_gain=2000):
     """
     Create simulated data to non-linear data to test non-linearity correction.