Merge pull request #1860 from pypeit/spatflex_redoshift

Spatial flexure reorganisation

pypeit/images/buildimage.py

@@ -160,7 +160,7 @@ def construct_file_name(cls, calib_key, calib_dir=None, basename=None):
 
 def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=None, dark=None,
                         scattlight=None, flatimages=None, maxiters=5, ignore_saturation=True,
-                        slits=None, mosaic=None, calib_dir=None, setup=None, calib_id=None):
+                        slits=None, mosaic=None, manual_spat_flexure=None, calib_dir=None, setup=None, calib_id=None):
     """
     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
@@ -222,6 +222,9 @@ def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=
             Flag processed image will be a mosaic of multiple detectors.  By
             default, this is determined by the format of ``det`` and whether or
             not this is a bias or dark frame.  *Only used for testing purposes.*
+        manual_spat_flexure (:obj:`list`, `numpy.ndarray`_, optional):
+            A list of the spatial flexures for each image in file_list.  This is only
+            used to manually correct the slit traces for spatial flexure of each image.
         calib_dir (:obj:`str`, `Path`_, optional):
             The directory for processed calibration files.  Required for
             elements of :attr:`frame_image_classes`, ignored otherwise.
@@ -247,20 +250,22 @@ def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=
         # NOTE: This should not be necessary because FrameGroupPar explicitly
         # requires frametype to be valid
         msgs.error(f'{frame_par["frametype"]} is not a valid PypeIt frame type.')
+    manual_spatflex = manual_spat_flexure if manual_spat_flexure is not None else [np.ma.masked]*len(file_list)
 
     # 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:
+    for ii, ifile in enumerate(file_list):
         # Load raw image
         rawImage = rawimage.RawImage(ifile, spectrograph, det)
         # Process
         rawImage_list.append(rawImage.process(
             frame_par['process'], scattlight=scattlight, bias=bias,
-            bpm=bpm, dark=dark, flatimages=flatimages, slits=slits, mosaic=mosaic))
+            bpm=bpm, dark=dark, flatimages=flatimages, slits=slits, mosaic=mosaic,
+            manual_spat_flexure=manual_spatflex[ii]))
 
     # Do it
     combineImage = combineimage.CombineImage(rawImage_list, frame_par['process'])

pypeit/images/rawimage.py

@@ -401,7 +401,7 @@ def build_rn2img(self, units='e-', digitization=False):
         return np.array(rn2)
 
     def process(self, par, bpm=None, scattlight=None, flatimages=None, bias=None, slits=None, dark=None,
-                mosaic=False, debug=False):
+                mosaic=False, manual_spat_flexure=None, debug=False):
         """
         Process the data.
 
@@ -519,6 +519,9 @@ def process(self, par, bpm=None, scattlight=None, flatimages=None, bias=None, sl
                 mosaic.  If flats or slits are provided (and used), this *must*
                 be true because these objects are always defined in the mosaic
                 frame.
+            manual_spat_flexure (:obj:`float`, optional):
+                The spatial flexure of the image. This is only set if the user wishes to
+                manually correct the slit traces of this image for spatial flexure.
             debug (:obj:`bool`, optional):
                 Run in debug mode.
 
@@ -671,9 +674,11 @@ def process(self, par, bpm=None, scattlight=None, flatimages=None, bias=None, sl
         # bias and dark subtraction) and before field flattening.  Also the
         # function checks that the slits exist if running the spatial flexure
         # correction, so no need to do it again here.
-        self.spat_flexure_shift = self.spatial_flexure_shift(slits, method=self.par['spat_flexure_method'],
-                                                                 maxlag=self.par['spat_flexure_maxlag']) \
-            if self.par['spat_flexure_method'] != "skip" else None
+        self.spat_flexure_shift = None
+        if self.par['spat_flexure_method'] != "skip" or not np.ma.is_masked(manual_spat_flexure):
+            self.spat_flexure_shift = self.spatial_flexure_shift(slits, method=self.par['spat_flexure_method'],
+                                                                 manual_spat_flexure=manual_spat_flexure,
+                                                                 maxlag=self.par['spat_flexure_maxlag'])
 
         #   - Subtract scattered light... this needs to be done before flatfielding.
         if self.par['subtract_scattlight']:
@@ -766,19 +771,26 @@ def _squeeze(self):
         return _det, self.image, self.ivar, self.datasec_img, self.det_img, self.rn2img, \
                 self.base_var, self.img_scale, self.bpm
 
-    def spatial_flexure_shift(self, slits, force=False, method="detector", maxlag=20):
+    def spatial_flexure_shift(self, slits, force=False, manual_spat_flexure=np.ma.masked, method="detector", maxlag=20):
         """
         Calculate a spatial shift in the edge traces due to flexure.
 
         This is a simple wrapper for
         :func:`~pypeit.core.flexure.spat_flexure_shift`.
 
         Args:
-            slits (:class:`~pypeit.slittrace.SlitTraceSet`, optional):
+            slits (:class:`~pypeit.slittrace.SlitTraceSet`):
                 Slit edge traces
             force (:obj:`bool`, optional):
                 Force the image to be field flattened, even if the step log
                 (:attr:`steps`) indicates that it already has been.
+            manual_spat_flexure (:obj:`float`, optional):
+                Manually set the spatial flexure shift. If provided, this
+                value is used instead of calculating the shift. The default
+                value is `np.ma.masked`, which means the shift is calculated
+                from the image data. The only way this value is used is if
+                the user sets the `shift` parameter in their pypeit file to
+                be a float.
             method (:obj:`str`, optional):
                 Method to use to calculate the spatial flexure shift. Options
                 are 'detector' (default), 'slit', and 'edge'. The 'detector'
@@ -801,10 +813,34 @@ def spatial_flexure_shift(self, slits, force=False, method="detector", maxlag=20
         if self.nimg > 1:
             msgs.error('CODING ERROR: Must use a single image (single detector or detector '
                        'mosaic) to determine spatial flexure.')
-        self.spat_flexure_shift = flexure.spat_flexure_shift(self.image[0], slits, method=method, maxlag=maxlag)
+
+        # Check if the slits are provided
+        if slits is None:
+            if not np.ma.is_masked(manual_spat_flexure):
+                msgs.warn('Manual spatial flexure provided without slits - assuming no spatial flexure.')
+            else:
+                msgs.warn('Cannot calculate spatial flexure without slits - assuming no spatial flexure.')
+            return
+
+        # First check for manual flexure
+        if not np.ma.is_masked(manual_spat_flexure):
+            msgs.info(f'Adopting a manual spatial flexure of {manual_spat_flexure} pixels')
+            spat_flexure = np.full((slits.nslits, 2), np.float64(manual_spat_flexure))
+        else:
+            spat_flexure = flexure.spat_flexure_shift(self.image[0], slits, method=method, maxlag=maxlag)
+
+        # Print the flexure values
+        if np.all(spat_flexure == spat_flexure[0, 0]):
+            msgs.info(f'Spatial flexure is: {spat_flexure[0, 0]} pixels')
+        else:
+            # Print the flexure values for each slit separately
+            for slit in range(spat_flexure.shape[0]):
+                msgs.info(
+                    f'Spatial flexure for slit {slits.spat_id[slit]} is: left={spat_flexure[slit, 0]} pixels; right={spat_flexure[slit, 1]} pixels')
+
         self.steps[step] = True
         # Return
-        return self.spat_flexure_shift
+        return spat_flexure
 
     def flatfield(self, flatimages, slits=None, force=False, debug=False):
         """

pypeit/inputfiles.py

@@ -481,20 +481,21 @@ def path_and_files(self, key:str, skip_blank=False, include_commented_out=False,
         ## Build full paths to file and set frame types
         data_files = []
         for row in self.data:
+            rowkey = '' if np.ma.is_masked(row[key]) else row[key]
 
             # Skip Empty entries?
-            if skip_blank and row[key].strip() in ['', 'none', 'None']:
+            if skip_blank and rowkey.strip() in ['', 'none', 'None']:
                 continue
 
             # Skip commented out entries
-            if row[key].strip().startswith("#"):
+            if rowkey.strip().startswith("#"):
                 if not include_commented_out:
                     continue
                 # Strip the comment character and any whitespace following it
                 # from the filename
-                name = row[key].strip("# ")
+                name = rowkey.strip("# ")
             else:
-                name = row[key]
+                name = rowkey
 
             # Searching..
             if len(self.file_paths) > 0:
@@ -503,7 +504,7 @@ def path_and_files(self, key:str, skip_blank=False, include_commented_out=False,
                     if os.path.isfile(filename):
                         break
             else:
-                filename = row[key]
+                filename = rowkey
 
             # Check we got a good hit
             if check_exists and not os.path.isfile(filename):

pypeit/metadata.py

@@ -497,8 +497,9 @@ def construct_basename(self, row, obstime=None):
         _obstime = self.construct_obstime(row) if obstime is None else obstime
         tiso = time.Time(_obstime, format='isot')
         dtime = datetime.datetime.strptime(tiso.value, '%Y-%m-%dT%H:%M:%S.%f')
+        this_target = "TargetName" if np.ma.is_masked(self['target'][row]) else self['target'][row].replace(" ", "")
         return '{0}-{1}_{2}_{3}{4}'.format(self['filename'][row].split('.fits')[0],
-                                           self['target'][row].replace(" ", ""),
+                                           this_target,
                                            self.spectrograph.camera,
                                            datetime.datetime.strftime(dtime, '%Y%m%dT'),
                                            tiso.value.split("T")[1].replace(':',''))
@@ -1309,10 +1310,31 @@ def frame_paths(self, indx):
         Returns:
             list: List of the full paths of one or more frames.
         """
-        if isinstance(indx, (int,np.integer)):
+        if isinstance(indx, (int, np.integer)):
             return os.path.join(self['directory'][indx], self['filename'][indx])
         return [os.path.join(d,f) for d,f in zip(self['directory'][indx], self['filename'][indx])]
 
+    def get_shifts(self, indx):
+        """
+        Return the shifts for the provided rows.
+
+        Args:
+            indx (:obj:`int`, array-like):
+                One or more 0-indexed rows in the table with the frames
+                to return.  Can be an array of indices or a boolean
+                array of the correct length.
+
+        Returns:
+            `numpy.ndarray`_: Array with the shifts.
+        """
+        # Make indx an array
+        _indx = np.atleast_1d(indx)
+        # Check if shifts are defined, if not, return a masked array
+        if 'shift' not in self.keys():
+            return np.ma.array(np.zeros(_indx.shape), mask=np.ones(_indx.shape, dtype=bool))
+        # Otherwise, return the shifts
+        return self['shift'][indx]
+
     def set_frame_types(self, type_bits, merge=True):
         """
         Set and return a Table with the frame types and bits.
@@ -1574,7 +1596,7 @@ def set_user_added_columns(self):
 
         """
         if 'manual' not in self.keys():
-            self['manual'] = ''
+            self['manual'] = np.ma.array(np.zeros(len(self)), mask=np.ones(len(self), dtype=bool))
         if 'shift' not in self.keys():
             # Instantiate a masked array
             self['shift'] = np.ma.array(np.zeros(len(self)), mask=np.ones(len(self), dtype=bool))

pypeit/pypeit.py

@@ -776,14 +776,16 @@ def objfind_one(self, frames, det, bg_frames=None, std_outfile=None):
 
         # Build Science image
         sci_files = self.fitstbl.frame_paths(frames)
+        manual_spat_flexure = self.fitstbl.get_shifts(frames)
         sciImg = buildimage.buildimage_fromlist(
             self.spectrograph, det, frame_par,
             sci_files, bias=self.caliBrate.msbias, bpm=self.caliBrate.msbpm,
             dark=self.caliBrate.msdark,
             scattlight=self.caliBrate.msscattlight,
             flatimages=self.caliBrate.flatimages,
             slits=self.caliBrate.slits,  # For flexure correction
-            ignore_saturation=False)
+            ignore_saturation=False,
+            manual_spat_flexure=manual_spat_flexure)
 
         # get no bkg subtracted sciImg to generate a global sky model without bkg subtraction.
         # it's a dictionary with only `image` and `ivar` keys if bkg_redux=False, otherwise it's None
@@ -795,6 +797,7 @@ def objfind_one(self, frames, det, bg_frames=None, std_outfile=None):
             bkg_redux_sciimg = sciImg
             # Build the background image
             bg_file_list = self.fitstbl.frame_paths(bg_frames)
+            bg_manual_spat_flexure = self.fitstbl.get_shifts(bg_frames)
             # TODO I think we should create a separate self.par['bkgframe'] parameter set to hold the image
             # processing parameters for the background frames.  This would allow the user to specify different
             # parameters for the background frames than for the science frames.  
@@ -805,50 +808,22 @@ def objfind_one(self, frames, det, bg_frames=None, std_outfile=None):
                                                    scattlight=self.caliBrate.msscattlight,
                                                    flatimages=self.caliBrate.flatimages,
                                                    slits=self.caliBrate.slits,
-                                                   ignore_saturation=False)
+                                                   ignore_saturation=False,
+                                                   manual_spat_flexure=bg_manual_spat_flexure)
 
             # NOTE: If the spatial flexure exists for sciImg, the subtraction
             # function propagates that to the subtracted image, ignoring any
             # spatial flexure determined for the background image.
             sciImg = bkg_redux_sciimg.sub(bgimg)
 
-        # Flexure
-        spat_flexure = np.zeros((self.caliBrate.slits.nslits, 2))  # No spatial flexure, unless we find it below
-        # use the flexure correction in the "shift" column
-        manual_flexure = self.fitstbl[frames[0]]['shift']
-        if (self.objtype == 'science' and self.par['scienceframe']['process']['spat_flexure_method'] != "skip") or \
-                (self.objtype == 'standard' and self.par['calibrations']['standardframe']['process']['spat_flexure_method'] != "skip") or \
-                not np.ma.is_masked(manual_flexure):
-            # First check for manual flexure
-            if not np.ma.is_masked(manual_flexure):
-                msgs.info(f'Implementing manual spatial flexure of {manual_flexure} pixels')
-                spat_flexure = np.full((self.caliBrate.slits.nslits, 2), np.float64(manual_flexure))
-                sciImg.spat_flexure = spat_flexure
-            else:
-                if sciImg.spat_flexure is not None:
-                    msgs.info(f'Using auto-computed spatial flexure')
-                    spat_flexure = sciImg.spat_flexure
-                else:
-                    msgs.info('Assuming no spatial flexure correction')
-        else:
-            msgs.info('Assuming no spatial flexure correction')
-
-        # Print the flexure values
-        if np.all(spat_flexure == spat_flexure[0, 0]):
-            msgs.info(f'Spatial flexure is: {spat_flexure[0, 0]}')
-        else:
-            # Print the flexure values for each slit separately
-            for slit in range(spat_flexure.shape[0]):
-                msgs.info(
-                    f'Spatial flexure for slit {self.caliBrate.slits.spat_id[slit]} is: left={spat_flexure[slit, 0]} right={spat_flexure[slit, 1]}')
         # Build the initial sky mask
         initial_skymask = self.load_skyregions(initial_slits=self.spectrograph.pypeline != 'SlicerIFU',
-                                               scifile=sciImg.files[0], frame=frames[0], spat_flexure=spat_flexure)
+                                               scifile=sciImg.files[0], frame=frames[0], spat_flexure=sciImg.spat_flexure)
 
         # Deal with manual extraction
         row = self.fitstbl[frames[0]]
         manual_obj = ManualExtractionObj.by_fitstbl_input(
-            row['filename'], row['manual'], self.spectrograph) if len(row['manual'].strip()) > 0 else None
+            row['filename'], row['manual'], self.spectrograph) if not np.ma.is_masked(row['manual']) else None
 
         # Instantiate Reduce object
         # Required for pypeline specific object

pypeit/spec2dobj.py

@@ -580,7 +580,15 @@ def build_primary_hdr(self, raw_header, spectrograph, calib_dir=None,
         # Add the spectrograph-specific sub-header
         if subheader is not None:
             for key in subheader.keys():
-                hdr[key.upper()] = subheader[key]
+                # Find the value and check if it is masked
+                if isinstance(subheader[key], (tuple, list)):
+                    # value + comment
+                    _value = ('', subheader[key][1]) if np.ma.is_masked(subheader[key][0]) else subheader[key]
+                else:
+                    # value only
+                    _value = '' if np.ma.is_masked(subheader[key]) else subheader[key]
+                # Update the header card with the corresponding value
+                hdr[key.upper()] = _value
 
         # PYPEIT
         # TODO Should the spectrograph be written to the header?

pypeit/specobjs.py

@@ -811,7 +811,15 @@ def write_to_fits(self, subheader, outfile, overwrite=True, update_det=None,
                     for line in str(subheader[key.upper()]).split('\n'):
                         header[key.upper()] = line
             else:
-                header[key.upper()] = subheader[key]
+                # Find the value and check if it is masked
+                if isinstance(subheader[key], (tuple, list)):
+                    # value + comment
+                    _value = ('', subheader[key][1]) if np.ma.is_masked(subheader[key][0]) else subheader[key]
+                else:
+                    # value only
+                    _value = '' if np.ma.is_masked(subheader[key]) else subheader[key]
+                # Update the header card with the corresponding value
+                header[key.upper()] = _value
                 # Also store the datetime in ISOT format
                 if key.upper() == 'MJD':
                     if isinstance(subheader[key], (list, tuple)):