Merge branch 'release/4.x' into feature/gnirsxd_tutorials

Author: KathleenLabrie

gemini_instruments/gnirs/lookup.py

@@ -0,0 +1,133 @@
+import numpy as np
+
+from gempy.gemini import gemini_tools as gt
+
+from ..gemini.primitives_gemini import Gemini
+from . import parameters_irdc
+
+from recipe_system.utils.decorators import parameter_override, capture_provenance
+
+
+@parameter_override
+@capture_provenance
+class IRDC(Gemini):
+    """
+    This is the class containing all of the preprocessing primitives
+    for instruments that use the Gemini Infrared Detector Controller (IRDC).
+    """
+    tagset = {"GEMINI"}
+
+    def _initialize(self, adinputs, **kwargs):
+        super()._initialize(adinputs, **kwargs)
+        self._param_update(parameters_irdc)
+
+    def nonlinearityCorrect(self, adinputs=None, suffix=None):
+        """
+        Run on raw or nprepared Gemini NIRI data, this script calculates and
+        applies a per-pixel linearity correction based on the counts in the
+        pixel, the exposure time, the read mode, the bias level and the ROI.
+        Pixels over the maximum correctable value are set to BADVAL unless
+        given the force flag. Note that you may use glob expansion in infile,
+        however, any pattern matching characters (*,?) must be either quoted
+        or escaped with a backslash. Do we need a badval parameter that defines
+        a value to assign to uncorrectable pixels, or do we want to just add
+        those pixels to the DQ plane with a specific value?
+
+        Parameters
+        ----------
+        suffix: str
+            suffix to be added to output files
+        """
+        # Instantiate the log
+        log = self.log
+        log.debug(gt.log_message("primitive", self.myself(), "starting"))
+        timestamp_key = self.timestamp_keys[self.myself()]
+
+        def linearize(counts, coeffs):
+            """Return a linearized version of the counts in electrons per coadd"""
+            # The coefficients are in the form of a {region: coeffs} dictionary.
+            corrected_counts = np.empty_like(counts)
+            for _slice, coeff in coeffs.items():
+                log.debug("Coefficients for {} rows = {:.6f} "
+                             "{:.9e} {:.9e}".format(
+                    _slice, coeff.time_delta, coeff.gamma, coeff.eta))
+                corrected_counts[_slice] = (
+                        counts[_slice] * (1 + counts[_slice] *
+                                          (np.float32(coeff.gamma) +
+                                           counts[_slice] * np.float32(coeff.eta))))
+            return corrected_counts
+
+        for ad in adinputs:
+            if ad.phu.get(timestamp_key):
+                log.warning("No changes will be made to {}, since it has "
+                            "already been processed by nonlinearityCorrect".
+                            format(ad.filename))
+                continue
+
+            total_exptime = ad.exposure_time()
+            coadds = ad.coadds()
+            # Check the raw exposure time (i.e., per coadd). First, convert
+            # the total exposure time returned by the descriptor back to
+            # the raw exposure time
+            exptime = np.float32(total_exptime /
+                                 (coadds if ad.is_coadds_summed() else 1))
+            if exptime > 600:
+                log.warning(f"Exposure time {exptime} for {ad.filename} is "
+                            "outside the range used to derive correction.")
+
+            for ext, gain, coeffs in zip(ad, ad.gain(), self._nonlinearity_coeffs(ad)):
+                if coeffs is None:
+                    log.warning("No nonlinearity coefficients found for "
+                                f"{ad.filename} extension {ext.id} - "
+                                "no correction applied")
+                    continue
+                elif not isinstance(coeffs, dict):
+                    # coeffs apply to entire array
+                    coeffs = {None: coeffs}
+
+                raw_mean_value = np.mean(ext.data) / coadds
+                log.fullinfo("The mean value of the raw pixel data in " \
+                             "{} is {:.8f}".format(ext.filename, raw_mean_value))
+
+                # Create a new array that contains the corrected pixel data.
+                # Remember that gain() returns 1.0 after ADUToElectrons
+                raw_pixel_data = ext.data * np.float32(gain / coadds)
+                corrected_pixel_data = (linearize(raw_pixel_data, coeffs) *
+                                        np.float32(coadds / gain))
+
+                # Try to do something useful with the VAR plane, if it exists
+                # Since the data are fairly pristine, VAR will simply be the
+                # Poisson noise (divided by gain if in ADU), possibly plus RN**2
+                # So making an additive correction will sort this out,
+                # irrespective of whether there's read noise
+                if ext.variance is not None:
+                    ext.variance += (corrected_pixel_data - ext.data) / np.float32(gain)
+                # Now update the SCI extension
+                ext.data = corrected_pixel_data
+
+                # Correct for the exposure time issue by scaling the counts
+                # to the nominal exposure time
+                time_delta = np.float32(np.mean([v.time_delta for v in coeffs.values()]))
+                ext.multiply(exptime / (exptime + time_delta))
+
+                # Determine the mean of the corrected pixel data
+                corrected_mean_value = np.mean(ext.data) / coadds
+                log.fullinfo("The mean value of the corrected pixel data in "
+                             "{} is {:.8f}".format(ext.filename, corrected_mean_value))
+
+                # Correct the exposure time by adding coeff1 * coadds
+                total_exptime = total_exptime + time_delta * coadds
+
+                # Update descriptors for saturation and nonlinear thresholds
+                log.fullinfo(f"The true total exposure time = {total_exptime}")
+                for desc in ('saturation_level', 'non_linear_level'):
+                    current_value = getattr(ext, desc)()
+                    new_value = linearize(
+                        np.array([current_value * gain / coadds]),
+                        coeffs)[0] * coadds / gain
+                    ext.hdr[ad._keyword_for(desc)] = np.round(new_value, 3)
+
+            # Timestamp and update filename
+            gt.mark_history(ad, primname=self.myself(), keyword=timestamp_key)
+            ad.update_filename(suffix=suffix, strip=True)
+        return adinputs

geminidr/gemini/parameters_irdc.py

@@ -71,5 +71,3 @@ class tracePinholeAperturesConfig(parameters_spect.tracePinholeAperturesConfig):
     """
     max_shift = config.RangeField("Maximum shift per pixel in line position",
                                   float, 0.4, min=0.001, max=0.5, inclusiveMax=True)
-    def setDefaults(self):
-        self.start_pos = 255

geminidr/gemini/primitives_irdc.py

@@ -7,7 +7,7 @@
 
 from gemini_instruments.gnirs import lookup as adlookup
 from ..core import NearIR
-from ..gemini.primitives_gemini import Gemini
+from ..gemini.primitives_irdc import IRDC
 from . import parameters_gnirs
 
 from recipe_system.utils.decorators import parameter_override, capture_provenance
@@ -16,11 +16,10 @@
 # ------------------------------------------------------------------------------
 @parameter_override
 @capture_provenance
-class GNIRS(Gemini, NearIR):
+class GNIRS(IRDC, NearIR):
     """
-    This is the class containing all of the preprocessing primitives
-    for the F2 level of the type hierarchy tree. It inherits all
-    the primitives from the level above
+    This is the class containing all of the primitives used by all GNIRS
+    modes. It inherits all the primitives from the level above.
     """
     tagset = {"GEMINI", "GNIRS"}
 
@@ -98,3 +97,27 @@ def standardizeInstrumentHeaders(self, adinputs=None, suffix=None):
             gt.mark_history(ad, primname=self.myself(), keyword=timestamp_key)
             ad.update_filename(suffix=suffix, strip=True)
         return adinputs
+
+    def _nonlinearity_coeffs(self, ad):
+        """
+        Returns a list of namedtuples containing the necessary information to
+        perform a nonlinearity correction. The list contains one namedtuple
+        per extension (although normal GNIRS data only has a single extension).
+
+        Returns
+        -------
+        list of namedtuples
+            nonlinearity info (max counts, exptime correction, gamma, eta)
+        """
+        array_name = set(ad.array_name())
+        assert len(array_name) == 1, ("Multiple array names found in {}".
+                                      format(ad.filename))
+        try:
+            nonlin_coeffs = self.inst_adlookup.nonlin_coeffs[array_name.pop()]
+        except KeyError:
+            self.log.warning("No nonlinearity coefficients are available for "
+                             "this array/detector controller.")
+            return [None] * len(ad)
+        read_mode = ad.read_mode()
+        well_depth = ad.well_depth_setting()
+        return [nonlin_coeffs.get((read_mode,well_depth))] * len(ad)

geminidr/gnirs/parameters_gnirs_crossdispersed.py

@@ -13,6 +13,7 @@
 from .primitives_gnirs import GNIRS
 from . import parameters_gnirs_spect
 
+from gempy.gemini import gemini_tools as gt
 from gempy.library import wavecal
 
 from recipe_system.utils.decorators import parameter_override, capture_provenance
@@ -32,6 +33,66 @@ def _initialize(self, adinputs, **kwargs):
         super()._initialize(adinputs, **kwargs)
         self._param_update(parameters_gnirs_spect)
 
+    def normalizeFlat(self, adinputs=None, **params):
+        """
+        A GNIRS-specific primitive to normalize the spectroscopic flatfield.
+        Because of the odd/even row effect, the trace of the brightness has
+        two loci and can be difficult to fit, sometimes jumping from one locus
+        to the other. This primitive will attempt to remove the odd/even
+        offset temporarily while the fit takes place, by scaling the odd rows
+        by a constance value to match the even rows.
+
+        Parameters
+        ----------
+        suffix : str/None
+            suffix to be added to output files
+        center : int/None
+            central row/column for 1D extraction (None => use middle)
+        nsum : int
+            number of rows/columns to average (about "center")
+        function : str
+            type of function to fit (splineN or polynomial types)
+        order : int
+            Order of the spline fit to be performed
+        lsigma : float/None
+            lower rejection limit in standard deviations
+        hsigma : float/None
+            upper rejection limit in standard deviations
+        niter : int
+            maximum number of rejection iterations
+        grow : float/False
+            growth radius for rejected pixels
+        interactive : bool
+            set to activate an interactive preview to fine tune the input parameters
+        """
+        log = self.log
+        log.debug(gt.log_message("primitive", self.myself(), "starting"))
+
+        for ad in adinputs:
+            # XD will have multiple extensions at this point
+            all_scaling_data = []
+            for ext in ad:
+                row_values = np.ma.median(np.ma.masked_array(
+                    ext.data, mask=ext.mask), axis=1)
+                bright_enough = row_values[1::2] > 50 * ext.read_noise()
+                if bright_enough.any():
+                    scaling = np.ma.median((row_values[::2] /
+                                            row_values[1::2])[bright_enough])
+                else:
+                    scaling = 1.0
+                log.debug(f"Scaling for {ad.filename}:{ext.id}: {scaling:8.6f}")
+                scaling_data = np.ones_like(ext.data)
+                scaling_data[1::2] *= scaling
+                ext.multiply(scaling_data)
+                all_scaling_data.append(scaling_data)
+
+            # normalizeFlat() operates in-place; we don't need a return value
+            super().normalizeFlat([ad], **params)
+            for ext, scaling_data in zip(ad, all_scaling_data):
+                ext.divide(scaling_data)
+
+        return adinputs
+
     def standardizeWCS(self, adinputs=None, **params):
         """
         This primitive updates the WCS attribute of each NDAstroData extension

geminidr/gnirs/primitives_gnirs.py

@@ -20,9 +20,9 @@ def makeProcessedFlat(p):
 
     p.prepare()
     p.addDQ(add_illum_mask=False)
-    #p.nonlinearityCorrect()
     p.ADUToElectrons()
     p.addVAR(poisson_noise=True, read_noise=True)
+    p.nonlinearityCorrect()
     p.makeLampFlat()
     p.addIllumMaskToDQ()
     p.normalizeFlat()

geminidr/gnirs/primitives_gnirs_spect.py

@@ -18,6 +18,7 @@ def makeProcessedFlat(p):
     p.addDQ()
     p.ADUToElectrons()
     p.addVAR(poisson_noise=True, read_noise=True)
+    p.nonlinearityCorrect()
     p.stackFlats()
     p.determineSlitEdges()
     p.maskBeyondSlit()

geminidr/gnirs/recipes/sq/recipes_FLAT_IMAGE.py

@@ -28,6 +28,7 @@ def makeProcessedFlat(p):
     p.addDQ()
     p.ADUToElectrons()
     p.addVAR(poisson_noise=True, read_noise=True)
+    p.nonlinearityCorrect()
     p.selectFromInputs(tags='GCAL_IR_ON,LAMPON', outstream='IRHigh')
     p.removeFromInputs(tags='GCAL_IR_ON,LAMPON')
     p.stackFlats(stream='main')

geminidr/gnirs/recipes/sq/recipes_FLAT_LS_SPECT.py

@@ -22,6 +22,7 @@ def reduce(p):
     p.addDQ()
     p.ADUToElectrons()
     p.addVAR(read_noise=True, poisson_noise=True)
+    p.nonlinearityCorrect()
     p.darkCorrect()
     p.flatCorrect()
     p.applyDQPlane()