Merge branch 'develop' into tilts_fwhm

README.rst

@@ -46,7 +46,7 @@ PypeIt |forks| |stars|
 
 |github| |pypi| |pypi_downloads| |License|
 
-|docs| |CITests| |Coverage| 
+|docs| |CITests|
 
 |DOI_latest| |JOSS| |arxiv|
 

doc/help/run_pypeit.rst

@@ -4,7 +4,7 @@
     usage: run_pypeit [-h] [-v VERBOSITY] [-r REDUX_PATH] [-m] [-s] [-o] [-c]
                       pypeit_file
     
-    ##  [1;37;42mPypeIt : The Python Spectroscopic Data Reduction Pipeline v1.16.1.dev468+g832ee84e4[0m
+    ##  [1;37;42mPypeIt : The Python Spectroscopic Data Reduction Pipeline v1.16.1.dev95+g35e103c10[0m
     ##  
     ##  Available spectrographs include:
     ##   aat_uhrf, bok_bc, gemini_flamingos1, gemini_flamingos2,

doc/qa.rst

@@ -9,9 +9,11 @@ PypeIt QA
 =========
 
 As part of the standard reduction, PypeIt generates a series
-of Quality Assurance (QA) files. This document describes
+of fixed-format Quality Assurance (QA) figures. This document describes
 the typical outputs, in the typical order that they appear.
 
+*This page is still a work in progress.*
+
 The basic arrangement is that individual PNG files are created
 and then a set of HTML files are generated to organize
 viewing of the PNGs.
@@ -20,16 +22,10 @@ viewing of the PNGs.
 HTML
 ====
 
-When the code completes (or crashes out), a set of
-HTML files are generated in the ``QA/`` folder.  There
-is one HTML file per calibration frame set and one
-HTML file per science exposure.  Example names are
-``MF_A.html``.
-
-Open in your browser and have at 'em.
-Quick links are provided to allow one to jump between
-the various files.
-
+When the code completes (or crashes out), an HTML file is generated in the
+``QA/`` folder, one per setup that has been reduced (typically one).  An example
+filename is ``MF_A.html``.  These HTML files are out of date, so you're better
+off opening the PNG files in the ``PNGs`` directory directly.
 
 Calibration QA
 ==============
@@ -39,9 +35,38 @@ to calibration processing.  There is a unique
 one generated for each setup and detector and
 (possibly) calibration set.
 
-Generally, the title describes the type of QA and the
-sub-title indicates the user who ran PypeIt and the
-date of the processing.
+Generally, the title describes the type of QA plotted.
+
+.. _qa-order-predict:
+
+Echelle Order Prediction
+------------------------
+
+When reducing echelle observations and inserting missing orders, a QA plot is
+produced to assess the success of the predicted locations.  The example below is
+for Keck/HIRES.
+
+.. figure:: figures/Edges_A_0_MSC01_orders_qa.png
+   :width: 60%
+
+   Example QA plot showing the measured order spatial widths (blue) and gaps
+   (green) in pixels.  The widths should be nearly constant as a function of
+   position, whereas the gaps should change monotonically with spatial pixel.
+
+In the figure above, measured values that are included in the polynomial fit are
+shown as filled points.  The colored lines show the best fit polynomial model
+used for the predicted order locations.  The fit allows for an iterative
+rejection of points; measured widths and gaps that are rejected during the fit
+are shown as orange and purple crosses, respectively.  The measurements that are
+rejected during the fit are not necessarily *removed* as invalid traces, but the
+code allows you to identify outlier traces that *will be* removed.  None of the
+traces in the example image above are identified as outliers; if they exist,
+they will be plotted as orange and purple triangles for widths and gaps,
+respectively.  Missing orders that will be added are included as open squares;
+gaps are green, widths are blue.  To deal with overlap, "bracketing" orders are
+added for the overlap calculation but are removed in the final set of traces;
+the title of the plot indicates if bracketing orders are included and the
+vertical dashed lines shows the edges of the detector/mosaic.
 
 .. _qa-wave-fit:
 
@@ -52,7 +77,7 @@ PypeIt produces plots like the one below showing the result of the wavelength
 calibration.
 
 .. figure:: figures/deimos_arc1d.png
-   :width: 60 %
+   :width: 60%
 
    An example QA plot for Keck/DEIMOS wavelength calibration.  The extracted arc
    spectrum is shown to the left with arc lines used for the wavelength solution

doc/releases/1.16.1dev.rst

@@ -21,9 +21,18 @@ Dependency Changes
 Functionality/Performance Improvements and Additions
 ----------------------------------------------------
 
+- Added the ``max_overlap`` parameter, which limits the set of new order traces
+  added, to compensate for orders missed during automated edge-tracing, to those
+  that have less than a given fractional overlap with adjacent orders.
+- Added the ``order_fitrej`` and ``order_outlier`` parameters used to set the
+  sigma-clipping threshold used when fitting Legendre functions to the order
+  widths and gaps.
 - Added the possibility to decide if the extracted standard star spectrum should be
   used as a crutch for tracing the object in the science frame (before it was done as default).
   This is done by setting the parameter ``use_std_trace`` in FindObjPar.
+- Now PypeIt can handle the case where "Standard star trace does not match the
+  number of orders in the echelle data" both in `run_pypeit` and in
+  `pypeit_coadd_1dspec`.
 - Now PypeIt can handle the case where "Standard star trace does not match the number of orders in the echelle data"
   both in `run_pypeit` and in `pypeit_coadd_1dspec`.
 - Added the functionality to use slitless flats to create pixelflats. Note: new frametype
@@ -109,6 +118,14 @@ Under-the-hood Improvements
 - Introduced :class:`~pypeit.pypeitdata.PypeItDataPaths` to handle all
   interactions with the ``pypeit/data`` directory, which provides a unified
   interface for accessing on-disk and cached files.
+- When adding missing orders, the full syncing procedure is no longer performed.
+  The code now only checks that the edges are still synced after the missed
+  orders are added.
+- When detecting overlapping orders/slits, the code now forces each edge used to
+  have been directly detected; i.e., if an edge is inserted, the fact that the
+  resulting slit is abnormally short should not trigger the overlap detection.
+- Improved the QA plot resulting from fitting order widths and gaps as a
+  function of spatial position.
 - Updated general raw image reader so that it correctly accounts for
   spectrographs that read the data and overscan sections directly from the file
   headers.

pypeit/core/datacube.py

@@ -223,7 +223,7 @@ def extract_point_source(wave, flxcube, ivarcube, bpmcube, wcscube, exptime,
 
     Returns
     -------
-    sobjs : :class:`pypeit.specobjs.SpecObjs`
+    sobjs : :class:`~pypeit.specobjs.SpecObjs`
         SpecObjs object containing the extracted spectrum
     """
     if whitelight_range is None:

pypeit/core/flexure.py

@@ -37,7 +37,7 @@
 from pypeit.datamodel import DataContainer
 from pypeit.images.detector_container import DetectorContainer
 from pypeit.images.mosaic import Mosaic
-from pypeit import specobj, specobjs, spec2dobj
+from pypeit import specobj, specobjs
 from pypeit import wavemodel
 
 from IPython import embed
@@ -1395,83 +1395,6 @@ def sky_em_residuals(wave:np.ndarray, flux:np.ndarray,
     return dwave[m], diff[m], diff_err[m], los[m], los_err[m]
 
 
-def flexure_diagnostic(file, file_type='spec2d', flexure_type='spec', chk_version=False):
-    """
-    Print the spectral or spatial flexure of a spec2d or spec1d file
-
-    Args:
-        file (:obj:`str`, `Path`_):
-            Filename of the spec2d or spec1d file to check
-        file_type (:obj:`str`, optional):
-            Type of the file to check. Options are 'spec2d' or 'spec1d'. Default
-            is 'spec2d'.
-        flexure_type (:obj:`str`, optional):
-            Type of flexure to check. Options are 'spec' or 'spat'. Default is
-            'spec'.
-        chk_version (:obj:`bool`, optional):
-            If True, raise an error if the datamodel version or type check
-            failed.  If False, throw a warning only. Default is False.
-
-    Returns:
-        :obj:`astropy.table.Table`, :obj:`float`: If file_type is 'spec2d' and
-        flexure_type is 'spec', return a table with the spectral flexure.  If
-        file_type is 'spec2d' and flexure_type is 'spat', return the spatial
-        flexure.  If file_type is 'spec1d', return a table with the spectral
-        flexure.
-    """
-
-    # value to return
-    return_flex = None
-
-    if file_type == 'spec2d':
-        # load the spec2d file
-        allspec2D = spec2dobj.AllSpec2DObj.from_fits(file, chk_version=chk_version)
-        # Loop on Detectors
-        for det in allspec2D.detectors:
-            print('')
-            print('=' * 50 + f'{det:^7}' + '=' * 51)
-            # get and print the spectral flexure
-            if flexure_type == 'spec':
-                spec_flex = allspec2D[det].sci_spec_flexure
-                spec_flex.rename_column('sci_spec_flexure', 'global_spec_shift')
-                if np.all(spec_flex['global_spec_shift'] != None):
-                    spec_flex['global_spec_shift'].format = '0.3f'
-                # print the table
-                spec_flex.pprint_all()
-                # return the table
-                return_flex = spec_flex
-            # get and print the spatial flexure
-            if flexure_type == 'spat':
-                spat_flex = allspec2D[det].sci_spat_flexure
-                # print the value
-                print(f'Spat shift: {spat_flex}')
-                # return the value
-                return_flex = spat_flex
-    elif file_type == 'spec1d':
-        # no spat flexure in spec1d file
-        if flexure_type == 'spat':
-            msgs.error("Spat flexure not available in the spec1d file, try with a spec2d file")
-        # load the spec1d file
-        sobjs = specobjs.SpecObjs.from_fitsfile(file, chk_version=chk_version)
-        spec_flex = Table()
-        spec_flex['NAME'] = sobjs.NAME
-        spec_flex['global_spec_shift'] = sobjs.FLEX_SHIFT_GLOBAL
-        if np.all(spec_flex['global_spec_shift'] != None):
-            spec_flex['global_spec_shift'].format = '0.3f'
-        spec_flex['local_spec_shift'] = sobjs.FLEX_SHIFT_LOCAL
-        if np.all(spec_flex['local_spec_shift'] != None):
-            spec_flex['local_spec_shift'].format = '0.3f'
-        spec_flex['total_spec_shift'] = sobjs.FLEX_SHIFT_TOTAL
-        if np.all(spec_flex['total_spec_shift'] != None):
-            spec_flex['total_spec_shift'].format = '0.3f'
-        # print the table
-        spec_flex.pprint_all()
-        # return the table
-        return_flex = spec_flex
-
-    return return_flex
-
-
 # TODO -- Consider separating the methods from the DataContainer as per calibrations
 class MultiSlitFlexure(DataContainer):
     """

pypeit/core/skysub.py

@@ -1090,7 +1090,7 @@ def ech_local_skysub_extract(sciimg, sciivar, fullmask, tilts, waveimg,
     Parameters
     ----------
     sciimg : `numpy.ndarray`_
-        science image, usually with a global sky subtracted.
+        Science image, usually with a global sky subtracted.
         shape = (nspec, nspat)
     sciivar : `numpy.ndarray`_
         inverse variance of science image.
@@ -1181,9 +1181,9 @@ def ech_local_skysub_extract(sciimg, sciivar, fullmask, tilts, waveimg,
             fullbkpt = bset.breakpoints
 
     force_gauss : bool, default = False
-        If True, a Gaussian profile will always be assumed for the
-        optimal extraction using the FWHM determined from object finding (or provided by the user) for the spatial
-        profile.
+        If True, a Gaussian profile will always be assumed for the optimal
+        extraction using the FWHM determined from object finding (or provided by
+        the user) for the spatial profile.
     sn_gauss : int or float, default = 4.0
         The signal to noise threshold above which optimal extraction
         with non-parametric b-spline fits to the objects spatial

pypeit/core/trace.py

@@ -1647,7 +1647,41 @@ def predicted_center_difference(lower_spat, spat, width_fit, gap_fit):
     return np.absolute(spat - test_spat)
 
 
-def extrapolate_orders(cen, width_fit, gap_fit, min_spat, max_spat, tol=0.01):
+def extrapolate_prev_order(cen, width_fit, gap_fit, tol=0.01):
+    """
+    Predict the location of the order to the spatial left (lower spatial pixel
+    numbers) of the order center provided.
+
+    This is largely a support function for :func:`extrapolate_orders`.
+
+    Args:
+        cen (:obj:`float`):
+            Center of the known order.
+        width_fit (:class:`~pypeit.core.fitting.PypeItFit`):
+            Model of the order width as a function of the order center.
+        gap_fit (:class:`~pypeit.core.fitting.PypeItFit`):
+            Model of the order gap *after* each order as a function of the order
+            center.
+        tol (:obj:`float`, optional):
+            Tolerance used when optimizing the order locations predicted toward
+            lower spatial pixels.
+
+    Returns:
+        :obj:`float`: Predicted center of the previous order.
+    """
+    # Guess the position of the previous order
+    w = width_fit.eval(cen)
+    guess = np.array([cen - w - gap_fit.eval(cen)])
+    # Set the bounds based on this guess and the expected order width
+    bounds = optimize.Bounds(lb=guess - w/2, ub=guess + w/2)
+    # Optimize the spatial position
+    res = optimize.minimize(predicted_center_difference, guess,
+                            args=(cen, width_fit, gap_fit),
+                            method='trust-constr', jac='2-point', bounds=bounds, tol=tol)    
+    return res.x[0]
+
+
+def extrapolate_orders(cen, width_fit, gap_fit, min_spat, max_spat, tol=0.01, bracket=False):
     """
     Predict the locations of additional orders by extrapolation.
 
@@ -1676,6 +1710,9 @@ def extrapolate_orders(cen, width_fit, gap_fit, min_spat, max_spat, tol=0.01):
         tol (:obj:`float`, optional):
             Tolerance used when optimizing the order locations predicted toward
             lower spatial pixels.
+        bracket (:obj:`bool`, optional):
+            Bracket the added orders with one additional order on either side.
+            This can be useful for dealing with predicted overlap.
 
     Returns:
         :obj:`tuple`: Two arrays with orders centers (1) below the first and (2)
@@ -1686,25 +1723,20 @@ def extrapolate_orders(cen, width_fit, gap_fit, min_spat, max_spat, tol=0.01):
     # Extrapolate toward lower spatial positions
     lower_spat = [cen[0]]
     while lower_spat[-1] > min_spat:
-        # Guess the position of the previous order
-        l = width_fit.eval(lower_spat[-1])
-        guess = np.array([lower_spat[-1] - l - gap_fit.eval(lower_spat[-1])])
-        # Set the bounds based on this guess and the expected order width
-        bounds = optimize.Bounds(lb=guess - l/2, ub=guess + l/2)
-        # Optimize the spatial position
-        res = optimize.minimize(predicted_center_difference, guess,
-                                args=(lower_spat[-1], width_fit, gap_fit),
-                                method='trust-constr', jac='2-point', bounds=bounds, tol=tol)
-        lower_spat += [res.x[0]]
+        lower_spat += [extrapolate_prev_order(lower_spat[-1], width_fit, gap_fit, tol=tol)]
 
     # Extrapolate toward larger spatial positions
     upper_spat = [cen[-1]]
     while upper_spat[-1] < max_spat:
         upper_spat += [upper_spat[-1] + width_fit.eval(upper_spat[-1]) 
                         + gap_fit.eval(upper_spat[-1])]
 
-    # Return arrays after removing the first and last spatial position (which
-    # are either repeats of values in `cen` or outside the spatial range)
+    # Return arrays after removing the first spatial position because it is a
+    # repeat of the input.  If not bracketting, also remove the last point
+    # because it will be outside the minimum or maximum spatial position (set by
+    # min_spat, max_spat).
+    if bracket:
+        return np.array(lower_spat[-1:0:-1]), np.array(upper_spat[1:])
     return np.array(lower_spat[-2:0:-1]), np.array(upper_spat[1:-1])