Minor improvements to quicklook, coadd2d as well as a bug fix

pypeit/coadd2d.py

@@ -27,13 +27,15 @@
 from pypeit.core import findobj_skymask
 from pypeit.core.wavecal import wvutils
 from pypeit.core import coadd
+from pypeit.core import parse 
 #from pypeit.core import parse
 from pypeit import calibrations
 from pypeit import spec2dobj
 from pypeit.core.moment import moment1d
 from pypeit.manual_extract import ManualExtractionObj
 
-
+#TODO We should decide which parameters go in through the parset 
+# and which parameters are passed in to the method as arguments
 class CoAdd2D:
 
     """
@@ -835,7 +837,8 @@ def reduce(self, pseudo_dict, show=False, clear_ginga=True, show_peaks=False, sh
         # maskdef stuff
         if parcopy['reduce']['slitmask']['assign_obj'] and slits.maskdef_designtab is not None:
             # Get plate scale
-            platescale = sciImage.detector.platescale * self.spat_samp_fact
+            platescale = parse.parse_binning(sciImage.detector.binning)[1]*sciImage.detector.platescale*self.spat_samp_fact
+            #platescale = sciImage.detector.platescale * self.spat_samp_fact
 
             # Assign slitmask design information to detected objects
             slits.assign_maskinfo(sobjs_obj, platescale, None, TOLER=parcopy['reduce']['slitmask']['obj_toler'])
@@ -870,51 +873,30 @@ def reduce(self, pseudo_dict, show=False, clear_ginga=True, show_peaks=False, sh
                objmodel_pseudo, ivarmodel_pseudo, outmask_pseudo, sobjs, sciImage.detector, slits, \
                pseudo_dict['tilts'], pseudo_dict['waveimg']
 
-    def snr_report(self, snr_bar, slitid=None):
-        """
-        ..todo.. I need a doc string
-
-        Args:
-            snr_bar:
-            slitid:
-
-        Returns:
-
-        """
-
-        # Print out a report on the SNR
-        msg_string = msgs.newline() + '-------------------------------------'
-        msg_string += msgs.newline() + '  Summary for highest S/N object'
-        if slitid is not None:
-            msg_string += msgs.newline() + '      found on slitid = {:d}            '.format(slitid)
-        msg_string += msgs.newline() + '-------------------------------------'
-        msg_string += msgs.newline() + '           exp#        S/N'
-        for iexp, snr in enumerate(snr_bar):
-            msg_string += msgs.newline() + '            {:d}         {:5.2f}'.format(iexp, snr)
-
-        msg_string += msgs.newline() + '-------------------------------------'
-        msgs.info(msg_string)
 
-    def offsets_report(self, offsets, offsets_method):
+    def offsets_report(self, offsets, platescale, offsets_method):
         """
         Print out a report on the offsets
 
         Args:
-            offsets:
-            offsets_method:
-
+            offsets (`numpy.ndarray`_)
+                Array of offsets
+            platescale (float):
+                the plate scale for this detector
+            offsets_method (str):
+                A string describing the method used to determine the offsets
         Returns:
 
         """
 
         if offsets_method is not None and offsets is not None:
-            msg_string = msgs.newline() + '---------------------------------------------'
+            msg_string = msgs.newline() + '---------------------------------------------------------------------------------'
             msg_string += msgs.newline() + ' Summary of offsets from {}     '.format(offsets_method)
-            msg_string += msgs.newline() + '---------------------------------------------'
-            msg_string += msgs.newline() + '           exp#      offset                  '
+            msg_string += msgs.newline() + '---------------------------------------------------------------------------------'
+            msg_string += msgs.newline() + '           exp#      offset (pixels)    offset (arcsec)'
             for iexp, off in enumerate(offsets):
-                msg_string += msgs.newline() + '            {:d}        {:5.2f}'.format(iexp, off)
-            msg_string += msgs.newline() + '-----------------------------------------------'
+                msg_string += msgs.newline() + '            {:2d}            {:6.2f}              {:6.3f}'.format(iexp, off, off*platescale)
+            msg_string += msgs.newline() + '---------------------------------------------------------------------------------'
             msgs.info(msg_string)
 
     def offset_slit_cen(self, slitid, offsets):
@@ -1126,17 +1108,18 @@ def compute_offsets(self, offsets):
 
         """
         msgs.info('Get Offsets')
+        platescale = parse.parse_binning(self.stack_dict['detectors'][0].binning)[1]*self.stack_dict['detectors'][0].platescale
+        #platescale = self.stack_dict['detectors'][0].platescale
         # 1) offsets are provided in the header of the spec2d files
         if offsets == 'header':
             msgs.info('Using offsets from header')
-            pscale = self.stack_dict['detectors'][0].platescale
             dithoffs = [self.spectrograph.get_meta_value(f, 'dithoff') for f in self.spec2d]
             if None in dithoffs:
                 msgs.error('Dither offsets keyword not found for one or more spec2d files. '
                            'Choose another option for `offsets`')
-            dithoffs_pix = - np.array(dithoffs) / pscale
+            dithoffs_pix = - np.array(dithoffs) / platescale
             self.offsets = dithoffs_pix[0] - dithoffs_pix
-            self.offsets_report(self.offsets, 'header keyword')
+            self.offsets_report(self.offsets, platescale, 'header keyword')
 
         elif self.objid_bri is None and offsets == 'auto':
             msgs.error('Offsets cannot be computed because no unique reference object '
@@ -1147,15 +1130,15 @@ def compute_offsets(self, offsets):
             msgs.info('Using user input offsets')
             # use them
             self.offsets = self.check_input(offsets, 'offsets')
-            self.offsets_report(self.offsets, 'user input')
+            self.offsets_report(self.offsets, platescale, 'user input')
 
         # 3) parset `offsets` is = 'maskdef_offsets' (no matter if we have a bright object or not)
         elif offsets == 'maskdef_offsets':
             if self.maskdef_offset is not None:
                 # the offsets computed during the main reduction (`run_pypeit`) are used
                 msgs.info('Determining offsets using maskdef_offset recoded in SlitTraceSet')
                 self.offsets = self.maskdef_offset[0] - self.maskdef_offset
-                self.offsets_report(self.offsets, 'maskdef_offset')
+                self.offsets_report(self.offsets, platescale, 'maskdef_offset')
             else:
                 # if maskdef_offsets were not computed during the main reduction, we cannot continue
                 msgs.error('No maskdef_offset recoded in SlitTraceSet')
@@ -1211,7 +1194,7 @@ def compute_weights(self, weights):
             msgs.error('Invalid value for `weights`')
 
 
-    def get_brightest_object(self, specobjs_list, spat_ids):
+    def get_brightest_obj(self, specobjs_list, spat_ids):
         """
         Dummy method to identify the brightest object. Overloaded by child methods.
 
@@ -1312,7 +1295,7 @@ def __init__(self, spec2d_files, spectrograph, par, det=1, offsets=None, weights
             self.objid_bri, self.slitidx_bri, self.spatid_bri, self.snr_bar_bri = \
                 np.repeat(user_objid, self.nexp), _slitidx_bri, user_slit, None
         elif len(self.stack_dict['specobjs_list']) > 0 and (offsets == 'auto' or weights == 'auto'):
-            self.objid_bri, self.slitidx_bri, self.spatid_bri, self.snr_bar_bri = \
+            self.objid_bri, self.spat_pixpos_bri, self.slitidx_bri, self.spatid_bri, self.snr_bar_bri = \
                 self.get_brightest_obj(self.stack_dict['specobjs_list'], self.spat_ids)
         else:
             self.objid_bri, self.slitidx_bri, self.spatid_bri, self.snr_bar_bri = (None,)*4
@@ -1402,7 +1385,8 @@ def compute_offsets(self, offsets):
                 plt.show()
 
             self.offsets = offsets
-            self.offsets_report(self.offsets, offsets_method)
+            platescale = parse.parse_binning(self.stack_dict['detectors'][0].binning)[1]*self.stack_dict['detectors'][0].platescale
+            self.offsets_report(self.offsets, platescale, offsets_method)
 
     def compute_weights(self, weights):
         """
@@ -1427,47 +1411,55 @@ def compute_weights(self, weights):
             _, self.use_weights = self.optimal_weights(self.spatid_bri, self.objid_bri, weight_method='constant')
             if self.par['coadd2d']['user_obj'] is not None:
                 msgs.info(f'Weights computed using a unique reference object in slit={self.spatid_bri} provided by the user')
+                # TODO Should we not still printo out a report for this usage case?
             else:
                 msgs.info(f'Weights computed using a unique reference object in slit={self.spatid_bri} with the highest S/N')
-                self.snr_report(self.snr_bar_bri, slitid=self.spatid_bri)
+                self.snr_report(self.spatid_bri, self.spat_pixpos_bri, self.snr_bar_bri)
 
-    def get_brightest_obj(self, specobjs_list, spat_ids):
+    def get_brightest_obj(self, specobjs_list, slit_spat_ids):
 
         """
-        Utility routine to find the brightest object in each exposure given a specobjs_list for MultiSlit reductions.
+        Utility routine to find the brightest reference object in each exposure given a specobjs_list 
+        for MultiSlit reductions.
 
         Args:
             specobjs_list: list
                List of SpecObjs objects.
-            spat_ids (`numpy.ndarray`_):
+            slit_spat_ids (`numpy.ndarray`_):
 
         Returns:
             tuple: Returns the following:
-                - objid: ndarray, int, shape (len(specobjs_list),):
-                  Array of object ids representing the brightest object
+                - objid: ndarray, int, shape=(len(specobjs_list),):
+                  Array of object ids representing the brightest reference object
+                  in each exposure
+                - spatid_pixpos: ndarray, float, shape=(len(specobjs_list),):
+                  Array of spatial pixel positions of the brightest reference object
                   in each exposure
-                - slit_idx (int): 0-based index
-                - spat_id (int): SPAT_ID for slit that highest S/N ratio object is on
-                  (only for pypeline=MultiSlit)
+                - slit_idx (int): 
+                  A zero-based index for the slit that the brightest object is on
+                - spat_id (int): 
+                  The SPAT_ID for the slit that the highest S/N ratio object is on
                 - snr_bar: ndarray, float, shape (len(list),): Average
-                  S/N over all the orders for this object
+                  S/N computed over all the exposures for this brightest reference object
         """
         msgs.info('Finding brightest object')
         nexp = len(specobjs_list)
-        nslits = spat_ids.size
+        nslits = slit_spat_ids.size
 
         slit_snr_max = np.zeros((nslits, nexp), dtype=float)
         bpm = np.ones(slit_snr_max.shape, dtype=bool)
         objid_max = np.zeros((nslits, nexp), dtype=int)
+        spat_pixpos_max = np.zeros((nslits, nexp), dtype=float)
         # Loop over each exposure, slit, find the brightest object on that slit for every exposure
         for iexp, sobjs in enumerate(specobjs_list):
             msgs.info("Working on exposure {}".format(iexp))
-            for islit, spat_id in enumerate(spat_ids):
+            for islit, spat_id in enumerate(slit_spat_ids):
                 if len(sobjs) == 0:
                     continue
                 ithis = np.abs(sobjs.SLITID - spat_id) <= self.par['coadd2d']['spat_toler']
                 if np.any(ithis):
                     objid_this = sobjs[ithis].OBJID
+                    spat_pixpos_this = sobjs[ithis].SPAT_PIXPOS
                     fluxes, ivars, gpms = [], [], []
                     for iobj, spec in enumerate(sobjs[ithis]):
                         # check if OPT_COUNTS is available
@@ -1495,6 +1487,7 @@ def get_brightest_obj(self, specobjs_list, spat_ids):
                         imax = np.argmax(rms_sn)
                         slit_snr_max[islit, iexp] = rms_sn[imax]
                         objid_max[islit, iexp] = objid_this[imax]
+                        spat_pixpos_max[islit, iexp] = spat_pixpos_this[imax]
                         bpm[islit, iexp] = False
 
         # If a slit has bpm = True for some exposures and not for others, set bpm = True for all exposures
@@ -1517,8 +1510,41 @@ def get_brightest_obj(self, specobjs_list, spat_ids):
             snr_bar_mean = slit_snr[slitid]
             snr_bar = slit_snr_max[slitid, :]
             objid = objid_max[slitid, :]
+            spat_pixpos = spat_pixpos_max[slitid, :]
+
+        return objid, spat_pixpos, slitid, slit_spat_ids[slitid], snr_bar
+
+
+    def snr_report(self, slitid, spat_pixpos, snr_bar):
+        """
+
+        Print out a SNR report for the reference object used to compute the weights for multislit 2D coadds. 
+
+        Args:
+            slitid (:obj:`int`):
+                The SPAT_ID of the slit that the reference object is on
+            spat_pixpos (:obj:`numpy.ndarray`):
+                Array of spatial pixel position of the reference object in the slit for each exposure shape = (nexp,)
+            snr_bar (:obj:`numpy.ndarray`):
+                Array of average S/N ratios for the reference object in each exposure, shape = (nexp,)
+
+
+        Returns:
+
+        """
+
+        # Print out a report on the SNR
+        msg_string = msgs.newline() + '-------------------------------------'
+        msg_string += msgs.newline() + '  Summary for highest S/N object'
+        msg_string += msgs.newline() + '      found on slitid = {:d}            '.format(slitid)
+        msg_string += msgs.newline() + '-------------------------------------'
+        msg_string += msgs.newline() + '       exp#   spat_pixpos     S/N'
+        msg_string += msgs.newline() + '-------------------------------------'
+        for iexp, (spat,snr) in enumerate(zip(spat_pixpos, snr_bar)):
+            msg_string += msgs.newline() + '       {:2d}      {:7.1f}      {:5.2f}'.format(iexp, spat, snr)
+        msg_string += msgs.newline() + '-------------------------------------'
+        msgs.info(msg_string)
 
-        return objid, slitid, spat_ids[slitid], snr_bar
 
     # TODO add an option here to actually use the reference trace for cases where they are on the same slit and it is
     # single slit???
@@ -1799,6 +1825,30 @@ def get_brightest_obj(self, specobjs_list, nslits):
             return None, None, None
         return objid, None, snr_bar
 
+    def snr_report(self, snr_bar):
+        """
+        Printo out a SNR report for echelle 2D coadds.
+
+        Args:
+            snr_bar (:obj:`numpy.ndarray`):
+                Array of average S/N ratios for the brightest object in each exposure. Shape = (nexp,)
+
+        Returns:
+
+        """
+
+        # Print out a report on the SNR
+        msg_string = msgs.newline() + '-------------------------------------'
+        msg_string += msgs.newline() + '  Summary for highest S/N object'
+        msg_string += msgs.newline() + '-------------------------------------'
+        msg_string += msgs.newline() + '           exp#        S/N'
+        for iexp, snr in enumerate(snr_bar):
+            msg_string += msgs.newline() + '            {:d}         {:5.2f}'.format(iexp, snr)
+
+        msg_string += msgs.newline() + '-------------------------------------'
+        msgs.info(msg_string)
+
+
     def reference_trace_stack(self, slitid, offsets=None, objid=None):
         """
         Utility function for determining the reference trace about

pypeit/core/coadd.py

@@ -882,7 +882,7 @@ def sn_weights(fluxes, ivars, gpms, sn_smooth_npix=None, weight_method='auto', v
 
     # Check if relative weights input
     if verbose:
-        msgs.info('Computing weights with weight_method=%s'.format(weight_method))
+        msgs.info('Computing weights with weight_method={:s}'.format(weight_method))
 
     weights = []
 
@@ -1859,18 +1859,18 @@ def spec_reject_comb(wave_grid, wave_grid_mid, waves_list, fluxes_list, ivars_li
         # Compute the stack
         wave_stack, flux_stack, ivar_stack, gpm_stack, nused = compute_stack(
             wave_grid, waves_list, fluxes_list, ivars_list, utils.array_to_explist(this_gpms, nspec_list=nspec_list), weights_list)
-        # Interpolate the individual spectra onto the wavelength grid of the stack. Use wave_grid_mid for this
-        # since it has no masked values
+        # Interpolate the stack onto the wavelength grids of the individual spectra. This will be used to perform 
+        # the rejection of pixels in the individual spectra.
         flux_stack_nat, ivar_stack_nat, gpm_stack_nat, _ = interp_spec(
             waves, wave_grid_mid, flux_stack, ivar_stack, gpm_stack)
-        ## TESTING
-        #nused_stack_nat, _, _ = interp_spec(
-        #    waves, wave_grid_mid, nused, ivar_stack, mask_stack)
-        #embed()
-        rejivars, sigma_corrs, outchi, chigpm = update_errors(fluxes, ivars, this_gpms,
-                                                               flux_stack_nat,  ivar_stack_nat, gpm_stack_nat, sn_clip=sn_clip)
-        this_gpms, qdone = pydl.djs_reject(fluxes, flux_stack_nat, outmask=this_gpms,inmask=gpms, invvar=rejivars,
+        rejivars, sigma_corrs, outchi, chigpm = update_errors(
+            fluxes, ivars, this_gpms, flux_stack_nat,  ivar_stack_nat, gpm_stack_nat, sn_clip=sn_clip)
+        this_gpms, qdone = pydl.djs_reject(fluxes, flux_stack_nat, outmask=this_gpms, inmask=gpms, invvar=rejivars,
                                           lower=lower,upper=upper, maxrej=maxrej, sticky=False)
+        #rejivars, sigma_corrs, outchi, chigpm = update_errors(
+        #    fluxes, ivars, this_gpms, flux_stack_nat,  ivar_stack_nat, gpm_stack_nat, sn_clip=sn_clip)
+        #this_gpms, qdone = pydl.djs_reject(fluxes, flux_stack_nat, outmask=this_gpms, inmask=gpms, invvar=rejivars,
+        #                                  lower=lower,upper=upper, maxrej=maxrej, sticky=False) 
         iter += 1
 
 

pypeit/images/buildimage.py

@@ -162,7 +162,7 @@ def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=
                         scattlight=None, flatimages=None, maxiters=5, ignore_saturation=True,
                         slits=None, mosaic=None, calib_dir=None, setup=None, calib_id=None):
     """
-    Perform basic image processing on a list of images and combine the results. All 
+    Perform basic image processing on a list of images and combine the results. All
     core processing steps for each image are handled by :class:`~pypeit.images.rawimage.RawImage` and
     image combination is handled by :class:`~pypeit.images.combineimage.CombineImage`.
     This function can be used to process both single images, lists of images, and detector mosaics.
@@ -251,15 +251,15 @@ def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=
     # Should the detectors be reformatted into a single image mosaic?
     if mosaic is None:
         mosaic = isinstance(det, tuple) and frame_par['frametype'] not in ['bias', 'dark']
-        
+
     rawImage_list = []
     # Loop on the files
     for ifile in file_list:
         # Load raw image
         rawImage = rawimage.RawImage(ifile, spectrograph, det)
         # Process
         rawImage_list.append(rawImage.process(
-            frame_par['process'], scattlight=scattlight, bias=bias, 
+            frame_par['process'], scattlight=scattlight, bias=bias,
             bpm=bpm, dark=dark, flatimages=flatimages, slits=slits, mosaic=mosaic))
 
     # Do it