Merge pull request #1527 from pypeit/fire_fixes

Enable Echelle Extraction to run through even when a bad order occurs

CHANGES.rst

@@ -31,6 +31,7 @@
 - Parse Keck/NIRES dither patterns, similar to MOSFIRE
 - Introduce BitMaskArray class to ease use of bitmasks
 - Fixed memory hogging by matplotlib when using version >= 3.6.1
+- Allow for bad orders during extraction  (without crashing)
 
 1.11.0 (21 Oct 2022)
 --------------------

pypeit/core/skysub.py

@@ -1210,40 +1210,53 @@ def ech_local_skysub_extract(sciimg, sciivar, fullmask, tilts, waveimg,
 
     # Find the spat IDs
     gdslit_spat = np.unique(slitmask[slitmask >= 0]).astype(int)  # Unique sorts
-    if gdslit_spat.size != norders:
-        msgs.error("You have not dealt with masked orders properly")
+    #if gdslit_spat.size != norders:
+    #    msgs.error("You have not dealt with masked orders properly")
 
-    if (np.sum(sobjs.sign > 0) % norders) == 0:
-        nobjs = int((np.sum(sobjs.sign > 0)/norders))
-    else:
-        msgs.error('Number of specobjs in sobjs is not an integer multiple of the number or ordres!')
+    #if (np.sum(sobjs.sign > 0) % norders) == 0:
+    #    nobjs = int((np.sum(sobjs.sign > 0)/norders))
+    #else:
+    #    msgs.error('Number of specobjs in sobjs is not an integer multiple of the number or ordres!')
 
-    order_snr = np.zeros((norders, nobjs))
+    # Set bad obj to -nan
     uni_objid = np.unique(sobjs[sobjs.sign > 0].ECH_OBJID)
+    nobjs = len(uni_objid)
+    order_snr = np.zeros((norders, nobjs))
+    order_snr_gpm = np.ones_like(order_snr) 
     for iord in range(norders):
         for iobj in range(nobjs):
             ind = (sobjs.ECH_ORDERINDX == iord) & (sobjs.ECH_OBJID == uni_objid[iobj])
-            order_snr[iord,iobj] = sobjs[ind].ech_snr
+            # Allow for missed/bad order
+            if np.sum(ind) == 0:
+                order_snr_gpm[iord,iobj] = False
+            else:
+                order_snr[iord,iobj] = sobjs[ind].ech_snr
 
     # Compute the average SNR and find the brightest object
-    snr_bar = np.mean(order_snr,axis=0)
+    snr_bar = np.sum(order_snr,axis=0) / np.sum(order_snr_gpm,axis=0)
     srt_obj = snr_bar.argsort()[::-1]
     ibright = srt_obj[0] # index of the brightest object
+
     # Now extract the orders in descending order of S/N for the brightest object
     srt_order_snr = order_snr[:,ibright].argsort()[::-1]
     fwhm_here = np.zeros(norders)
     fwhm_was_fit = np.zeros(norders,dtype=bool)
+
     # Print out a status message
     str_out = ''
     for iord in srt_order_snr:
-        str_out += '{:<8d}{:<8d}{:>10.2f}'.format(slit_vec[iord], order_vec[iord], order_snr[iord,ibright]) + msgs.newline()
+        if order_snr_gpm[iord,ibright]:
+            str_out += '{:<8d}{:<8d}{:>10.2f}'.format(slit_vec[iord], order_vec[iord], order_snr[iord,ibright]) + msgs.newline()
     dash = '-'*27
     dash_big = '-'*40
     msgs.info(msgs.newline() + 'Reducing orders in order of S/N of brightest object:' + msgs.newline() + dash +
               msgs.newline() + '{:<8s}{:<8s}{:>10s}'.format('slit','order','S/N') + msgs.newline() + dash +
               msgs.newline() + str_out)
     # Loop over orders in order of S/N ratio (from highest to lowest) for the brightest object
     for iord in srt_order_snr:
+        # Is this a bad slit?
+        if not np.any(order_snr_gpm[iord,:]):
+            continue
         order = order_vec[iord]
         msgs.info("Local sky subtraction and extraction for slit/order: {:d}/{:d}".format(iord,order))
         other_orders = (fwhm_here > 0) & np.invert(fwhm_was_fit)
@@ -1279,7 +1292,10 @@ def ech_local_skysub_extract(sciimg, sciivar, fullmask, tilts, waveimg,
                         spec.FWHM = fwhm_this_ord
 
                     str_out = ''
-                    for slit_now, order_now, snr_now, fwhm_now in zip(slit_vec[other_orders], order_vec[other_orders],order_snr[other_orders,ibright], fwhm_here[other_orders]):
+                    for slit_now, order_now, snr_now, fwhm_now in zip(
+                        slit_vec[other_orders], order_vec[other_orders],
+                        order_snr[other_orders,ibright], 
+                        fwhm_here[other_orders]):
                         str_out += '{:<8d}{:<8d}{:>10.2f}{:>10.2f}'.format(slit_now, order_now, snr_now, fwhm_now) + msgs.newline()
                     msgs.info(msgs.newline() + 'Using' +  fwhm_str + ' for FWHM of object={:d}'.format(uni_objid[iobj]) +
                               ' on slit/order: {:d}/{:d}'.format(iord,order) + msgs.newline() + dash_big +

pypeit/extraction.py

@@ -55,6 +55,8 @@ class Extract:
             Global spectral flexure correction for each slit (in pixels)
         vel_corr (float):
             Relativistic reference frame velocity correction (e.g. heliocentyric/barycentric/topocentric)
+        extract_bpm (`numpy.ndarray`_):
+            Bad pixel mask for extraction
 
     """
 
@@ -231,13 +233,17 @@ def __init__(self, sciImg, slits, sobjs_obj, spectrograph, par, objtype, global_
 
         # Now apply a global flexure correction to each slit provided it's not a standard star
         if self.par['flexure']['spec_method'] != 'skip' and not self.std_redux:
+            # Update slitshift values
             self.spec_flexure_correct(mode='global')
-            if self.nobj_to_extract > 0:
-                for iobj in range(self.sobjs_obj.nobj):
-                    islit = self.slits.spatid_to_zero(self.sobjs_obj[iobj].SLITID)
-                    self.sobjs_obj[iobj].update_flex_shift(self.slitshift[islit], flex_type='global')
-
-        # remove objects found in `BOXSLIT` (we don't want to extract those)
+            # Apply?
+            for iobj in range(self.sobjs_obj.nobj):
+                # Ignore negative sources
+                if self.sobjs_obj[iobj].sign < 0 and not self.return_negative:
+                    continue
+                islit = self.slits.spatid_to_zero(self.sobjs_obj[iobj].SLITID)
+                self.sobjs_obj[iobj].update_flex_shift(self.slitshift[islit], flex_type='global')
+
+        # Remove objects rejected for one or another reasons (we don't want to extract those)
         remove_idx = []
         for i, sobj in enumerate(self.sobjs_obj):
             if sobj.SLITID in list(self.slits.spat_id[self.extract_bpm]):
@@ -260,15 +266,6 @@ def nsobj_to_extract(self):
         else:
             return 0
 
-    @property
-    def nobj_to_extract(self):
-        """
-        Number of objects to extract. Defined in children
-        Returns:
-
-        """
-        return None
-
     def initialise_slits(self, slits, initial=False):
         """
         Gather all the :class:`SlitTraceSet` attributes
@@ -400,9 +397,8 @@ def run(self, model_noise=None, spat_pix=None):
                See main doc string for description
 
         """
-
         # Do we have any detected objects to extract?
-        if self.nobj_to_extract > 0:
+        if self.nsobj_to_extract > 0:
             # Extract + Return
             self.skymodel, self.objmodel, self.ivarmodel, self.outmask, self.sobjs \
                 = self.extract(self.global_sky, model_noise=model_noise, spat_pix=spat_pix)
@@ -436,9 +432,8 @@ def run(self, model_noise=None, spat_pix=None):
         # TODO: change slits.mask > 2 to use named flags.
         reduce_masked = np.where(np.invert(self.extract_bpm_init) & self.extract_bpm & (self.slits.mask > 2))[0]
         if len(reduce_masked) > 0:
-            # TODO Change BADREDUCE to BADEXTRACT
             self.slits.mask[reduce_masked] = self.slits.bitmask.turn_on(
-                self.slits.mask[reduce_masked], 'BADREDUCE')
+                self.slits.mask[reduce_masked], 'BADEXTRACT')
 
         # Return
         return self.skymodel, self.objmodel, self.ivarmodel, self.outmask, self.sobjs, self.waveimg, self.tilts
@@ -684,16 +679,6 @@ def __init__(self, sciImg, slits, sobjs_obj, spectrograph, par, objtype, **kwarg
         super().__init__(sciImg, slits, sobjs_obj, spectrograph, par, objtype, **kwargs)
 
 
-    @property
-    def nobj_to_extract(self):
-        """
-        See parent method docs.
-
-        Returns:
-
-        """
-        return self.nsobj_to_extract
-
     # TODO: JFH Should we reduce the number of iterations for standards or
     # near-IR redux where the noise model is not being updated?
     def local_skysub_extract(self, global_sky, sobjs, spat_pix=None, model_noise=True,
@@ -832,22 +817,6 @@ def __init__(self, sciImg, slits, sobjs_obj, spectrograph, par, objtype, **kwarg
             msgs.error('Unable to set Echelle orders, likely because they were incorrectly '
                        'assigned in the relevant SlitTraceSet.')
 
-
-    @property
-    def nobj_to_extract(self):
-        """
-        See parent method docs.
-
-        Returns:
-
-        """
-
-        norders = self.order_vec.size
-        if (self.nsobj_to_extract % norders) == 0:
-            return int(self.nsobj_to_extract/norders)
-        else:
-            msgs.error('Number of specobjs in sobjs is not an integer multiple of the number or ordres!')
-
     # JFH TODO Should we reduce the number of iterations for standards or near-IR redux where the noise model is not
     # being updated?
     def local_skysub_extract(self, global_sky, sobjs,

pypeit/find_objects.py

@@ -337,6 +337,7 @@ def run(self, std_trace=None, show_peaks=False, show_skysub_fit=False):
 
         # If the skip_skysub is set (i.e. image is already sky-subtracted), simply find objects
         if self.par['reduce']['findobj']['skip_skysub']:
+            msgs.info("Skipping global sky sub as per user request")
             sobjs_obj, self.nobj = self.find_objects(self.sciImg.image, self.sciImg.ivar,
                                                      std_trace=std_trace, show=self.findobj_show,
                                                      show_peaks=show_peaks)

pypeit/pypeit.py

@@ -577,9 +577,8 @@ def reduce_exposure(self, frames, bg_frames=None, std_outfile=None):
         detectors = self.spectrograph.select_detectors(subset=subset)
         msgs.info(f'Detectors to work on: {detectors}')
 
-        # Loop on Detectors
+        # Loop on Detectors -- Calibrate, process image, find objects
         # TODO: Attempt to put in a multiprocessing call here?
-        # objfind
         for self.det in detectors:
             msgs.info(f'Reducing detector {self.det}')
             # run calibration
@@ -597,8 +596,8 @@ def reduce_exposure(self, frames, bg_frames=None, std_outfile=None):
             # in the slitmask stuff in between the two loops
             calib_slits.append(self.caliBrate.slits)
             # global_sky, skymask and sciImg are needed in the extract loop
-            initial_sky, sobjs_obj, sciImg, objFind = self.objfind_one(frames, self.det, bg_frames,
-                                                                       std_outfile=std_outfile)
+            initial_sky, sobjs_obj, sciImg, objFind = self.objfind_one(
+                frames, self.det, bg_frames, std_outfile=std_outfile)
             if len(sobjs_obj)>0:
                 all_specobjs_objfind.add_sobj(sobjs_obj)
             initial_sky_list.append(initial_sky)
@@ -711,7 +710,7 @@ def calib_one(self, frames, det):
 
         """
 
-        msgs.info(f'Building calibrations for detector {det}')
+        msgs.info(f'Building/loading calibrations for detector {det}')
         # Instantiate Calibrations class
         caliBrate = calibrations.Calibrations.get_instance(
             self.fitstbl, self.par['calibrations'], self.spectrograph,
@@ -889,11 +888,19 @@ def extract_one(self, frames, det, sciImg, objFind, initial_sky, sobjs_obj):
             final_global_sky = objFind.global_skysub(previous_sky=initial_sky, skymask=skymask, show=self.show)
         scaleImg = objFind.scaleimg
 
+        # Each spec2d file includes the slits object with unique flagging
+        #  for extraction failures.  So we make a copy here before those flags
+        #  are modified.
+        maskdef_designtab = self.caliBrate.slits.maskdef_designtab
+        slits = copy.deepcopy(self.caliBrate.slits)
+        slits.maskdef_designtab = None
+
+
         # update here slits.mask since global_skysub modify reduce_bpm and we need to propagate it into extraction
         flagged_slits = np.where(objFind.reduce_bpm)[0]
         if len(flagged_slits) > 0:
-            self.caliBrate.slits.mask[flagged_slits] = \
-                self.caliBrate.slits.bitmask.turn_on(self.caliBrate.slits.mask[flagged_slits], 'BADREDUCE')
+            slits.mask[flagged_slits] = \
+                slits.bitmask.turn_on(slits.mask[flagged_slits], 'BADSKYSUB')
 
         msgs.info("Extraction begins for {} on det={}".format(self.basename, det))
 
@@ -902,7 +909,7 @@ def extract_one(self, frames, det, sciImg, objFind, initial_sky, sobjs_obj):
         # At instantiaton, the fullmask in self.sciImg is modified
         # TODO Are we repeating steps in the init for FindObjects and Extract??
         self.exTract = extraction.Extract.get_instance(
-            sciImg, self.caliBrate.slits, sobjs_obj, self.spectrograph,
+            sciImg, slits, sobjs_obj, self.spectrograph,
             self.par, self.objtype, global_sky=final_global_sky, waveTilts=self.caliBrate.wavetilts, wv_calib=self.caliBrate.wv_calib,
             bkg_redux=self.bkg_redux, return_negative=self.par['reduce']['extraction']['return_negative'],
             std_redux=self.std_redux, basename=self.basename, show=self.show)
@@ -933,14 +940,9 @@ def extract_one(self, frames, det, sciImg, objFind, initial_sky, sobjs_obj):
 
         # Construct table of spectral flexure
         spec_flex_table = Table()
-        spec_flex_table['spat_id'] = self.caliBrate.slits.spat_id
+        spec_flex_table['spat_id'] = slits.spat_id
         spec_flex_table['sci_spec_flexure'] = self.exTract.slitshift
 
-        # pull out maskdef_designtab from caliBrate.slits
-        maskdef_designtab = self.caliBrate.slits.maskdef_designtab
-        slits = copy.deepcopy(self.caliBrate.slits)
-        slits.maskdef_designtab = None
-
         # Construct the Spec2DObj
         spec2DObj = spec2dobj.Spec2DObj(sciimg=sciImg.image,
                                         ivarraw=sciImg.ivar,

pypeit/scripts/collate_1d.py

@@ -521,7 +521,7 @@ def build_parameters(args):
         params['collate1d']['match_using'] = args.match_using
 
     if args.exclude_slit_bm is not None and len(args.exclude_slit_bm) > 0:
-        params['collate1d']['exclude_slit_trace_bm'] = args.exclude_slit_bm
+        params['collate1d']['exclude_slit_trace_bm'] = args.exclude_slit_bm.split(',')
 
     if args.exclude_serendip:
         params['collate1d']['exclude_serendip'] = True
@@ -648,8 +648,8 @@ def get_parser(cls, width=None):
         parser.add_argument('--dry_run', action='store_true', help=blank_par.descr['dry_run'])
         parser.add_argument('--ignore_flux', default=False, action='store_true', help=blank_par.descr['ignore_flux'])
         parser.add_argument('--flux', default=False, action = 'store_true', help=blank_par.descr['flux'])
-        parser.add_argument('--exclude_slit_bm', type=str, nargs='*',
-                            help=blank_par.descr['exclude_slit_trace_bm'])
+        parser.add_argument('--exclude_slit_bm', type=str, # nargs='*',
+                            help=blank_par.descr['exclude_slit_trace_bm']+' Comma separated.')
         parser.add_argument('--exclude_serendip', action='store_true',
                             help=blank_par.descr['exclude_serendip'])
         parser.add_argument("--wv_rms_thresh", type=float, default = None, help=blank_par.descr['wv_rms_thresh'])

pypeit/slittrace.py

@@ -32,9 +32,10 @@ class SlitTraceBitMask(BitMask):
     """
     Mask bits used during slit tracing.
     """
-    version = '1.0.0'
+    version = '1.0.1'
 
     def __init__(self):
+        # Only ever append new bits (and don't remove old ones)
         mask = dict([
             ('SHORTSLIT', 'Slit formed by left and right edge is too short. Not ignored for flexure'),
             ('BOXSLIT', 'Slit formed by left and right edge is valid (large enough to be a valid '
@@ -44,7 +45,9 @@ def __init__(self):
             ('BADTILTCALIB', 'Tilts analysis failed for this slit'),
             ('SKIPFLATCALIB', 'Flat field generation failed for this slit. Skip flat fielding'),
             ('BADFLATCALIB', 'Flat field generation failed for this slit. Ignore it fully.'),
-            ('BADREDUCE', 'Skysub/extraction failed for this slit'),
+            ('BADREDUCE', 'Reduction failed for this slit'), # THIS IS DEPRECATED (we may remove in v1.13) BUT STAYS HERE TO ALLOW FOR BACKWARDS COMPATIBILITY
+            ('BADSKYSUB', 'Skysub failed for this slit'),
+            ('BADEXTRACT', 'Extraction failed for this slit'),
         ])
         super(SlitTraceBitMask, self).__init__(list(mask.keys()), descr=list(mask.values()))
 
@@ -58,7 +61,7 @@ def exclude_for_flexure(self):
         # Ignore these flags when performing a flexure calculation
         #  Currently they are *all* of the flags..
         return ['SHORTSLIT', 'USERIGNORE', 'BADWVCALIB', 'BADTILTCALIB',
-                'SKIPFLATCALIB', 'BADFLATCALIB', 'BADREDUCE']
+                'SKIPFLATCALIB', 'BADFLATCALIB', 'BADSKYSUB', 'BADEXTRACT']
 
 
 
@@ -258,9 +261,14 @@ def _validate(self):
         if self.slitbitm is None:
             self.slitbitm = ','.join(list(self.bitmask.keys()))
         else:
-            # Validate
-            if self.slitbitm != ','.join(list(self.bitmask.keys())):
-                msgs.error("Input BITMASK keys differ from current data model!")
+            # Validate -- All of the keys must be present and in current order, but new ones can exist
+            bitms = self.slitbitm.split(',')
+            curbitm = list(self.bitmask.keys())
+            for kk, bit in enumerate(bitms):
+                if curbitm[kk] != bit:
+                    msgs.error("Input BITMASK keys differ from current data model!")
+            # Update to current, no matter what
+            self.slitbitm = ','.join(list(self.bitmask.keys()))
         # Mask
         if self.mask is None:
             self.mask = self.mask_init.copy()