Merge pull request #1693 from pypeit/bad_order_fix

Echelle improvements

doc/releases/1.14.1dev.rst

@@ -12,7 +12,10 @@ Functionality/Performance Improvements and Additions
 
 - Started the development of instrument-specific scattered light removal. In this
   release, we only model KCWI/KCRM scattered light.
+- Turned on usage of the use_2dmodel_mask parameter for echelle spectrographs
 - Added support for Keck/KCRM RL data reduction.
+- Allow one to turn off reinit of reduce_bpm in global_skysub and made this 
+  the default for the final pass
 
 Instrument-specific Updates
 ---------------------------

pypeit/core/skysub.py

@@ -1034,6 +1034,7 @@ def ech_local_skysub_extract(sciimg, sciivar, fullmask, tilts, waveimg,
                              slitmask, sobjs, spat_pix=None,
                              fit_fwhm=False,
                              min_snr=2.0, bsp=0.6, trim_edg=(3,3), std=False, prof_nsigma=None,
+                             use_2dmodel_mask=True, 
                              niter=4, sigrej=3.5, bkpts_optimal=True,
                              force_gauss=False, sn_gauss=4.0, model_full_slit=False,
                              model_noise=True, debug_bkpts=False, show_profile=False,
@@ -1119,6 +1120,8 @@ def ech_local_skysub_extract(sciimg, sciivar, fullmask, tilts, waveimg,
         light profile like elliptical galaxies. If prof_nsigma is
         set then the profiles will no longer be apodized by an
         exponential at large distances from the trace.
+    use_2dmodel_mask : bool, optional
+        Use the mask made from profile fitting when extracting?
     niter : int, optional
         Number of iterations for successive profile fitting and local sky-subtraction
     sigrej : :obj:`float`, optional
@@ -1377,6 +1380,7 @@ def ech_local_skysub_extract(sciimg, sciivar, fullmask, tilts, waveimg,
                                        spat_pix=spat_pix, ingpm=inmask, std=std, bsp=bsp,
                                        trim_edg=trim_edg, prof_nsigma=prof_nsigma, niter=niter,
                                        sigrej=sigrej,
+                                       use_2dmodel_mask=use_2dmodel_mask, 
                                        bkpts_optimal=bkpts_optimal, force_gauss=force_gauss,
                                        sn_gauss=sn_gauss, model_full_slit=model_full_slit,
                                        model_noise=model_noise, debug_bkpts=debug_bkpts,

pypeit/extraction.py

@@ -849,6 +849,7 @@ def local_skysub_extract(self, global_sky, sobjs,
         sn_gauss = self.par['reduce']['extraction']['sn_gauss']
         model_full_slit = self.par['reduce']['extraction']['model_full_slit']
         force_gauss = self.par['reduce']['extraction']['use_user_fwhm']
+        use_2dmodel_mask = self.par['reduce']['extraction']['use_2dmodel_mask']
         self.skymodel, self.objmodel, self.ivarmodel, self.outmask, self.sobjs \
             = skysub.ech_local_skysub_extract(self.sciImg.image, self.sciImg.ivar,
                                               self.sciImg.fullmask, self.tilts, self.waveimg,
@@ -858,6 +859,7 @@ def local_skysub_extract(self, global_sky, sobjs,
                                               std=self.std_redux, fit_fwhm=fit_fwhm,
                                               min_snr=min_snr, bsp=bsp, sigrej=sigrej,
                                               force_gauss=force_gauss, sn_gauss=sn_gauss,
+                                              use_2dmodel_mask=use_2dmodel_mask,
                                               model_full_slit=model_full_slit,
                                               model_noise=model_noise,
                                               show_profile=show_profile,

pypeit/find_objects.py

@@ -361,6 +361,7 @@ def run(self, std_trace=None, show_peaks=False, show_skysub_fit=False):
         skymask_init = self.create_skymask(sobjs_obj)
         # Global sky subtract now using the skymask defined by object positions
         initial_sky = self.global_skysub(skymask=skymask_init, show_fit=show_skysub_fit)
+
         # Second pass object finding on sky-subtracted image with updated sky
         # created after masking objects
         if not self.std_redux and not self.par['reduce']['findobj']['skip_second_find']:
@@ -482,39 +483,47 @@ def get_platescale(self, slitord_id=None):
 
 
     def global_skysub(self, skymask=None, update_crmask=True, trim_edg = (0, 0),
-                      previous_sky=None, show_fit=False, show=False, show_objs=False, objs_not_masked=False):
+                      previous_sky=None, show_fit=False, show=False, 
+                      show_objs=False, objs_not_masked=False,
+                      reinit_bpm:bool=True):
         """
         Perform global sky subtraction, slit by slit
 
         Wrapper to skysub.global_skysub
 
         Args:
-            skymask (`numpy.ndarray`_, None):
+            skymask (`numpy.ndarray`_, optional):
                 A 2D image indicating sky regions (1=sky)
             update_crmask (bool, optional):
-                ??
+                Update the crmask in the science image
             trim_edg (tuple, optional):
                  A two tuple of ints that specify the number of pixels to trim from the slit edges
+                 Only used by the IFU child
             show_fit (bool, optional):
-                ??
+                Show the sky fits?
             show (bool, optional):
-                ??
+                Show the sky image generated?
             show_objs (bool, optional):
-                ??
+                If show=True, show the objects on the sky image?
             previous_sky (`numpy.ndarray`_, optional):
                 Sky model estimate from a previous run of global_sky
                 Used to generate an improved estimated of the variance
             objs_not_masked (bool, optional):
                 Set this to be True if there are objects on the slit/order that are not being masked
                 by the skymask. This is typically the case for the first pass sky-subtraction
                 before object finding, since a skymask has not yet been created.
+            reinit_bpm (:obj:`bool`, optional):
+                If True (default), the bpm is reinitialized to the initial bpm 
+                Should be False on the final run in case there was a failure
+                upstream and no sources were found in the slit/order
 
         Returns:
             `numpy.ndarray`_: image of the the global sky model
 
         """
         # reset bpm since global sky is run several times and reduce_bpm is here updated.
-        self.reduce_bpm = self.reduce_bpm_init.copy()
+        if reinit_bpm:
+            self.reduce_bpm = self.reduce_bpm_init.copy()
         # Prep
         global_sky = np.zeros_like(self.sciImg.image)
         # Parameters for a standard star
@@ -1238,12 +1247,21 @@ def joint_skysub(self, skymask=None, update_crmask=True, trim_edg=(0,0),
         return _global_sky
 
     def global_skysub(self, skymask=None, update_crmask=True, trim_edg=(0,0),
-                      previous_sky=None, show_fit=False, show=False, show_objs=False, objs_not_masked=False):
+                      previous_sky=None, show_fit=False, show=False, show_objs=False, objs_not_masked=False,
+                      reinit_bpm:bool=True):
         """
         Perform global sky subtraction. This SlicerIFU-specific routine ensures that the
         edges of the slits are not trimmed, and performs a spatial and spectral
         correction using the sky spectrum, if requested. See Reduce.global_skysub()
         for parameter definitions.
+
+        See base class method for description of parameters.
+
+        Args:
+            reinit_bpm (:obj:`bool`, optional):
+                If True (default), the bpm is reinitialized to the initial bpm 
+                Should be False on the final run in case there was a failure
+                upstream and no sources were found in the slit/order
         """
         if self.par['reduce']['findobj']['skip_skysub']:
             msgs.info("Skipping global sky sub as per user request")
@@ -1252,7 +1270,7 @@ def global_skysub(self, skymask=None, update_crmask=True, trim_edg=(0,0),
         # Generate a global sky sub for all slits separately
         global_sky_sep = super().global_skysub(skymask=skymask, update_crmask=update_crmask,
                                                trim_edg=trim_edg, show_fit=show_fit, show=show,
-                                               show_objs=show_objs)
+                                               show_objs=show_objs, reinit_bpm=reinit_bpm)
         # Check if any slits failed
         if np.any(global_sky_sep[self.slitmask >= 0] == 0) and not self.bkg_redux:
             # Cannot continue without a sky model for all slits

pypeit/pypeit.py

@@ -578,6 +578,7 @@ def reduce_exposure(self, frames, bg_frames=None, std_outfile=None):
             sciImg_list.append(sciImg)
             objFind_list.append(objFind)
 
+
         # slitmask stuff
         if len(calibrated_det) > 0 and self.par['reduce']['slitmask']['assign_obj']:
             # get object positions from slitmask design and slitmask offsets for all the detectors
@@ -953,7 +954,9 @@ def extract_one(self, frames, det, sciImg, objFind, initial_sky, sobjs_obj):
         else:
             # Update the skymask
             skymask = objFind.create_skymask(sobjs_obj)
-            final_global_sky = objFind.global_skysub(previous_sky=initial_sky, skymask=skymask, show=self.show)
+            final_global_sky = objFind.global_skysub(previous_sky=initial_sky, 
+                                                     skymask=skymask, show=self.show,
+                                                     reinit_bpm=False)
         scaleImg = objFind.scaleimg
 
         # Each spec2d file includes the slits object with unique flagging