Merge pull request #1646 from pypeit/coadd1d_badspec1d

Coadd1d bad spec1d

pypeit/coadd1d.py

@@ -12,6 +12,7 @@
 import numpy as np
 
 from astropy.io import fits
+from astropy import stats
 
 from pypeit.spectrographs.util import load_spectrograph
 from pypeit.onespec import OneSpec
@@ -83,6 +84,7 @@ def __init__(self, spec1dfiles, objids, spectrograph=None, par=None, sensfuncfil
         self.show = show
         self.nexp = len(self.spec1dfiles) # Number of exposures
         self.coaddfile = None
+        self.gpm_exp = np.ones(self.nexp, dtype=bool).tolist()  # list of bool indicating the exposures that have been coadded
 
     def run(self):
         """
@@ -133,7 +135,7 @@ def save(self, coaddfile, telluric=None, obj_model=None, overwrite=True):
 
         # Add history entries for coadding.
         history = History()
-        history.add_coadd1d(self.spec1dfiles, self.objids)
+        history.add_coadd1d(self.spec1dfiles, self.objids, gpm_exp=self.gpm_exp)
 
         # Add on others
         if telluric is not None:
@@ -207,9 +209,97 @@ def load(self):
                 header_out['DEC_OBJ'] = sobjs[indx][0]['DEC']
             headers.append(header_out)
 
-
         return waves, fluxes, ivars, gpms, headers
 
+    def check_exposures(self):
+        """
+        Check if there are bad exposures.
+        Exposures with flux masked everywhere are always removed.
+        Exposures that are considered bad based on their S/N compared to the
+        average S/N among all the exposures, are removed only if self.par['sigrej_exp'] is set.
+        The attributes self.waves, self.fluxes, self.ivars, self.gpms need to be defined.
+
+        Returns
+        -------
+        gpm_exp: list of bool
+            List of boolean that indicates which exposures
+            have been coadded. The length of the list is nexp.
+        _waves : list of float `numpy.ndarray`_
+            Updated list of wavelength arrays.
+        _fluxes : list of float `numpy.ndarray`_
+            Updated list of flux arrays.
+        _ivars : list of float `numpy.ndarray`_
+            Updated list of inverse variance arrays.
+        _gpms : list of bool `numpy.ndarray`_
+            Updated list of good pixel mask variance arrays.
+        """
+
+        # initialize the exposures lists
+        _waves = [wave for wave in self.waves]
+        _fluxes = [flux for flux in self.fluxes]
+        _ivars = [ivar for ivar in self.ivars]
+        _gpms = [gpm for gpm in self.gpms]
+        _spec1dfiles = [spec1dfile for spec1dfile in self.spec1dfiles]
+        _objids = [objid for objid in self.objids]
+
+        # good exposures index
+        goodindx_exp = np.arange(self.nexp)
+
+        # check if there are exposures that are completely masked out, i.e., gpms = False for all spectral pixels
+        masked_exps = [np.all(np.logical_not(gpm)) for gpm in _gpms]
+        if np.any(masked_exps):
+            msgs.warn(f'The following exposure(s) is/are completely masked out. It/They will not be coadded.')
+            [msgs.warn(f"Exposure {i}: {fname.split('/')[-1]}  {obj}")
+             for i, (fname, obj, masked_exp) in enumerate(zip(_spec1dfiles, _objids, masked_exps)) if masked_exp]
+            # remove masked out exposure
+            _waves = [wave for (wave, masked_exp) in zip(_waves, masked_exps) if not masked_exp]
+            _fluxes = [flux for (flux, masked_exp) in zip(_fluxes, masked_exps) if not masked_exp]
+            _ivars = [ivar for (ivar, masked_exp) in zip(_ivars, masked_exps) if not masked_exp]
+            _gpms = [gpm for (gpm, masked_exp) in zip(_gpms, masked_exps) if not masked_exp]
+            _spec1dfiles = [spec1dfile for (spec1dfile, masked_exp) in zip(_spec1dfiles, masked_exps) if not masked_exp]
+            _objids = [objid for (objid, masked_exp) in zip(_objids, masked_exps) if not masked_exp]
+            # update good exposures index
+            goodindx_exp = goodindx_exp[np.logical_not(masked_exps)]
+
+        # check if there is still more than 1 exposure left
+        if len(_fluxes) < 2:
+            msgs.error('At least 2 unmasked exposures are required for coadding.')
+
+        # check if there is any bad exposure by comparing the rms_sn with the median rms_sn among all exposures
+        if len(_fluxes) > 2:
+            # Evaluate the sn_weights.
+            rms_sn, weights = coadd.sn_weights(_fluxes, _ivars, _gpms, const_weights=True)
+            # some stats
+            mean, med, sigma = stats.sigma_clipped_stats(rms_sn, sigma_lower=2., sigma_upper=2.)
+            _sigrej = self.par['sigrej_exp'] if self.par['sigrej_exp'] is not None else 10.0
+            # we set thresh_value to never be less than 0.2
+            thresh_value = round(0.2 + med + _sigrej * sigma, 2)
+            bad_exps = rms_sn > thresh_value
+            if np.any(bad_exps):
+                warn_msg = f'The following exposure(s) has/have S/N > {thresh_value:.2f} ' \
+                           f'({_sigrej} sigma above the median S/N in the stack).'
+                if self.par['sigrej_exp'] is not None:
+                        warn_msg += ' It/They WILL NOT BE COADDED.'
+                msgs.warn(warn_msg)
+                [msgs.warn(f"Exposure {i}: {fname.split('/')[-1]}  {obj}")
+                 for i, (fname, obj, bad_exp) in enumerate(zip(_spec1dfiles, _objids, bad_exps)) if bad_exp]
+                if self.par['sigrej_exp'] is not None:
+                    # remove bad exposure
+                    _waves = [wave for (wave, bad_exp) in zip(_waves, bad_exps) if not bad_exp]
+                    _fluxes = [flux for (flux, bad_exp) in zip(_fluxes, bad_exps) if not bad_exp]
+                    _ivars = [ivar for (ivar, bad_exp) in zip(_ivars, bad_exps) if not bad_exp]
+                    _gpms = [gpm for (gpm, bad_exp) in zip(_gpms, bad_exps) if not bad_exp]
+                    _spec1dfiles = [spec1dfile for (spec1dfile, bad_exp) in zip(_spec1dfiles, bad_exps) if not bad_exp]
+                    _objids = [objid for (objid, bad_exp) in zip(_objids, bad_exps) if not bad_exp]
+                    # update good exposures index
+                    goodindx_exp = goodindx_exp[np.logical_not(bad_exps)]
+
+        # gpm for the exposures, i.e., which exposures have been coadded
+        gpm_exp = np.zeros(self.nexp, dtype=bool)
+        gpm_exp[goodindx_exp] = True
+
+        return gpm_exp.tolist(), _waves, _fluxes, _ivars, _gpms
+
     def coadd(self):
         """
         Perform coadd for for Multi/Longslit data using multi_combspec
@@ -220,9 +310,10 @@ def coadd(self):
         """
         # Load the data
         self.waves, self.fluxes, self.ivars, self.gpms, self.headers = self.load()
+        # check if there are bad exposures and remove them
+        self.gpm_exp, _waves, _fluxes, _ivars, _gpms = self.check_exposures()
         # Perform and return the coadd
-        return coadd.multi_combspec(
-            self.waves, self.fluxes, self.ivars, self.gpms,
+        return coadd.multi_combspec(_waves, _fluxes, _ivars, _gpms,
             sn_smooth_npix=self.par['sn_smooth_npix'], wave_method=self.par['wave_method'],
             dv=self.par['dv'], dwave=self.par['dwave'], dloglam=self.par['dloglam'],
             wave_grid_min=self.par['wave_grid_min'], wave_grid_max=self.par['wave_grid_max'],

pypeit/coadd2d.py

@@ -503,14 +503,14 @@ def optimal_weights(self, slitorderid, objid, const_weights=False):
             if not np.any(ithis):
                 msgs.error('Slit/order or OBJID provided not valid. Optimal weights cannot be determined.')
             # check if OPT_COUNTS is available
-            if sobjs[ithis][0].has_opt_ext():
+            if sobjs[ithis][0].has_opt_ext() and np.any(sobjs[ithis][0].OPT_MASK):
                 wave_iexp, flux_iexp, ivar_iexp, gpm_iexp = sobjs[ithis][0].get_opt_ext()
                 waves.append(wave_iexp)
                 fluxes.append(flux_iexp)
                 ivars.append(ivar_iexp)
                 gpms.append(gpm_iexp)
             # check if BOX_COUNTS is available
-            elif sobjs[ithis][0].has_box_ext():
+            elif sobjs[ithis][0].has_box_ext() and np.any(sobjs[ithis][0].BOX_MASK):
                 wave_iexp, flux_iexp, ivar_iexp, gpm_iexp = sobjs[ithis][0].get_box_ext()
                 waves.append(wave_iexp)
                 fluxes.append(flux_iexp)
@@ -959,19 +959,22 @@ def get_wave_grid(self, wave_method):
                     wave_box = moment1d(waveimg * mask, trace_spat, 2 * box_radius,
                                     row=row)[0] / (box_denom + (box_denom == 0.0))
                     gpm_box = box_denom > 0.
-                    waves += [wave for wave in wave_box.T]
-                    gpms  += [(wave > 0.) & gpm for (wave, gpm) in zip(wave_box.T, gpm_box.T)]
+                    waves += [wave for (wave, gpm) in zip(wave_box.T, gpm_box.T) if np.any(gpm)]
+                    gpms += [(wave > 0.) & gpm for (wave, gpm) in zip(wave_box.T, gpm_box.T) if np.any(gpm)]
 
         else:
             waves, gpms = [], []
             for iexp, spec_this in enumerate(self.stack_dict['specobjs_list']):
                 for spec in spec_this:
                     # NOTE: BOX extraction usage needed for quicklook
-                    waves.append(spec.OPT_WAVE if spec.OPT_WAVE is not None else spec.BOX_WAVE)
-                    gpms.append(spec.OPT_MASK if spec.OPT_MASK is not None else spec.BOX_MASK)
-                    # TODO -- OPT_MASK is likely to become a bpm with int values
-                    #gpm[:self.nspec_array[iexp], indx] = spec.OPT_MASK
-                    #indx += 1
+                    good_opt_ext = spec.has_opt_ext() and np.any(spec.OPT_MASK)
+                    good_box_ext = spec.has_box_ext() and np.any(spec.BOX_MASK)
+                    if good_opt_ext or good_box_ext:
+                        waves.append(spec.OPT_WAVE if good_opt_ext else spec.BOX_WAVE)
+                        gpms.append(spec.OPT_MASK if good_opt_ext else spec.BOX_MASK)
+                        # TODO -- OPT_MASK is likely to become a bpm with int values
+                        #gpm[:self.nspec_array[iexp], indx] = spec.OPT_MASK
+                        #indx += 1
 
         return wvutils.get_wave_grid(waves=waves, gpms=gpms, wave_method=wave_method,
                                                                 spec_samp_fact=self.spec_samp_fact)
@@ -1433,19 +1436,21 @@ def get_brightest_obj(self, specobjs_list, spat_ids):
         for iexp, sobjs in enumerate(specobjs_list):
             msgs.info("Working on exposure {}".format(iexp))
             for islit, spat_id in enumerate(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
                     fluxes, ivars, gpms = [], [], []
                     for iobj, spec in enumerate(sobjs[ithis]):
                         # check if OPT_COUNTS is available
-                        if spec.has_opt_ext():
+                        if spec.has_opt_ext() and np.any(spec.OPT_MASK):
                             _, flux_iobj, ivar_iobj, gpm_iobj = spec.get_opt_ext()
                             fluxes.append(flux_iobj)
                             ivars.append(ivar_iobj)
                             gpms.append(gpm_iobj)
                         # check if BOX_COUNTS is available
-                        elif spec.has_box_ext():
+                        elif spec.has_box_ext() and np.any(spec.BOX_MASK):
                             _, flux_iobj, ivar_iobj, gpm_iobj = spec.get_box_ext()
                             fluxes.append(flux_iobj)
                             ivars.append(ivar_iobj)
@@ -1549,9 +1554,10 @@ def get_maskdef_dict(self, slit_idx, ref_trace_stack):
             objpos_dspat_vec = np.zeros(self.nexp)
             for iexp in range(self.nexp):
                 # get maskdef_slitcen
-                maskdef_slitcen_pixpos = \
-                    self.stack_dict['slits_list'][iexp].maskdef_slitcen[self.nspec_array[0]//2, slit_idx] \
-                    + self.maskdef_offset[iexp]
+                mslitcen_pixpos = self.stack_dict['slits_list'][iexp].maskdef_slitcen
+                if mslitcen_pixpos.ndim < 2:
+                    mslitcen_pixpos = mslitcen_pixpos[:, None]
+                maskdef_slitcen_pixpos = mslitcen_pixpos[self.nspec_array[0]//2, slit_idx] + self.maskdef_offset[iexp]
 
                 # get maskdef_objpos
                 # find left edge
@@ -1730,10 +1736,10 @@ def get_brightest_obj(self, specobjs_list, nslits):
                     flux = None
                     ind = (sobjs.ECH_ORDERINDX == iord) & (sobjs.ECH_OBJID == uni_objid[iobj])
                     # check if OPT_COUNTS is available
-                    if sobjs[ind][0].has_opt_ext():
+                    if sobjs[ind][0].has_opt_ext() and np.any(sobjs[ind][0].OPT_MASK):
                         _, flux, ivar, mask = sobjs[ind][0].get_opt_ext()
                     # check if BOX_COUNTS is available
-                    elif sobjs[ind][0].has_box_ext():
+                    elif sobjs[ind][0].has_box_ext() and np.any(sobjs[ind][0].BOX_MASK):
                         _, flux, ivar, mask = sobjs[ind][0].get_box_ext()
                         msgs.warn(f'Optimal extraction not available for object {sobjs[ind][0].ECH_OBJID} '
                                   f'in order {sobjs[ind][0].ECH_ORDER}. Using box extraction.')

pypeit/core/wavecal/wvutils.py

@@ -229,7 +229,7 @@ def get_wave_grid(waves=None, gpms=None, wave_method='linear', iref=0, wave_grid
             gpms = [wave > 1.0 for wave in waves]
 
         if wave_grid_min is None:
-            wave_grid_min = np.min([wave[gpm].min()for wave, gpm in zip(waves, gpms)])
+            wave_grid_min = np.min([wave[gpm].min() for wave, gpm in zip(waves, gpms)])
         if wave_grid_max is None:
             wave_grid_max = np.max([wave[gpm].max() for wave, gpm in zip(waves, gpms)])
 

pypeit/history.py

@@ -6,6 +6,8 @@
 """
 
 import os.path
+import numpy as np
+from IPython import embed
 
 from astropy.time import Time
 from astropy.io import fits
@@ -85,7 +87,7 @@ def add_reduce(self, calib_id, metadata, frames, bg_frames):
             for frame in calib_frames:
                 self.append(f'{frame["frametype"]} "{frame["filename"]}"', add_date=False)
 
-    def add_coadd1d(self, spec1d_files, objids):
+    def add_coadd1d(self, spec1d_files, objids, gpm_exp=None):
         """
         Add history entries for 1D coadding.
         
@@ -104,27 +106,62 @@ def add_coadd1d(self, spec1d_files, objids):
         Args:
             spec1d_files (:obj:`list`): List of the spec1d files used for coadding.
             objids (:obj:`list`): List of the PypeIt object ids used in coadding.
+            gpm_exp (:obj:`list`, optional): List of boolean indicating which exposures were coadded.
         """
 
-        combined_files_objids = list(zip(spec1d_files, objids))
-        self.append(f'PypeIt Coadded {len(combined_files_objids)} objects from {len(set(spec1d_files))} spec1d files')
-
-        current_spec1d = ""
-        for (spec1d, objid) in combined_files_objids:
-            if spec1d != current_spec1d:
-                current_spec1d = spec1d
-
-                self.append(f'From "{os.path.basename(spec1d)}"', add_date=False)
-                header = fits.getheader(spec1d)
-                additional_info = None
-                if 'SEMESTER' in header:
-                    additional_info = f"Semester: {header['SEMESTER']}"
-                if 'PROGID' in header:
-                    additional_info += f" Program ID: {header['PROGID']}"
-                if additional_info is not None:
-                    self.append(additional_info, add_date=False)
-            self.append(objid, add_date=False)
+        if gpm_exp is not None:
+            # Not coadded files and objids
+            notcoadded_spec1d_files = [spec1d_file for (spec1d_file, gpm_exp) in zip(spec1d_files, gpm_exp) if not gpm_exp]
+            notcoadded_objids = [objid for (objid, gpm_exp) in zip(objids, gpm_exp) if not gpm_exp]
+            combined_notcoadd_files_objids = list(zip(notcoadded_spec1d_files, notcoadded_objids))
 
+            # Coadded files and objids
+            coadded_spec1d_files = [spec1d_file for (spec1d_file, gpm_exp) in zip(spec1d_files, gpm_exp) if gpm_exp]
+            coadded_objids = [objid for (objid, gpm_exp) in zip(objids, gpm_exp) if gpm_exp]
+            combined_files_objids = list(zip(coadded_spec1d_files, coadded_objids))
+        else:
+            combined_files_objids = list(zip(spec1d_files, objids))
+            combined_notcoadd_files_objids = None
+
+        files_objids = [combined_files_objids, combined_notcoadd_files_objids]
+        # add history
+        for file_objid in files_objids:
+            if file_objid is None:
+                continue
+            elif file_objid == combined_files_objids:
+                self.append(f'PypeIt Coadded {len(file_objid)} objects '
+                            f'from {np.unique([f[0] for f in file_objid]).size} spec1d files')
+            elif file_objid == combined_notcoadd_files_objids and len(file_objid) > 0:
+                self.append(f'PypeIt DID NOT COADD {len(file_objid)} objects '
+                            f'from {np.unique([f[0] for f in file_objid]).size} spec1d files', add_date=False)
+
+            current_spec1d = ""
+            for (spec1d, objid) in file_objid:
+                if spec1d != current_spec1d:
+                    current_spec1d = spec1d
+
+                    self.append(f'From "{os.path.basename(spec1d)}"', add_date=False)
+                    header = fits.getheader(spec1d)
+                    additional_info = None
+                    if 'SEMESTER' in header:
+                        additional_info = f"Semester: {header['SEMESTER']}"
+                    if 'PROGID' in header:
+                        additional_info += f" Program ID: {header['PROGID']}"
+                    if additional_info is not None:
+                        self.append(additional_info, add_date=False)
+                obj_info = objid
+                # get extension names
+                hnames = [h.name for h in fits.open(spec1d)]
+                # find the extension name that include objid
+                ind_ext = np.where([objid in h for h in hnames])[0]
+                if ind_ext.size > 0:
+                    # get the header for this extension
+                    this_ext_header = fits.getheader(spec1d, ext=ind_ext[0])
+                    if 'MASKDEF_ID' in this_ext_header:
+                        obj_info += f" {this_ext_header['MASKDEF_ID']}"
+                    if 'MASKDEF_OBJNAME' in this_ext_header:
+                        obj_info += f" {this_ext_header['MASKDEF_OBJNAME']}"
+                self.append(obj_info, add_date=False)
 
     def append(self, history, add_date=True):
         """Append a new history entry.