Merge branch 'develop' into spatflex_options

# Conflicts:
#	doc/calibrations/flexure.rst
#	doc/releases/1.17.2dev.rst
#	pypeit/core/flexure.py
#	pypeit/images/rawimage.py
#	pypeit/par/pypeitpar.py

doc/calibrations/flexure.rst

@@ -42,6 +42,10 @@ that the spatial flexure is different for different slits or slit edges. Also no
 that some spectrographs post-process the spatial flexure values when using the ``slit``
 or ``edge`` methods (e.g. perform a linear fit to the spatial flexure values) to improve
 the result. Consult the documentation for your instrument to see if this is the case.
+(search for the value of ``spat_flexure_correct``).
+A QA plot is generated for each frame for which the spatial flexure correction
+is applied. We recommend the user to inspect these plots to ensure the flexure
+correction is reasonable. See :ref:`qa-spat-flex` for more information.
 
 Depending on what frame types you choose to correct, the
 code will behave somewhat differently.  Here we describe
@@ -325,4 +329,29 @@ HDU Name               HDU Type                        Data Type  Description
 ``FLEXURE``            `astropy.io.fits.BinTableHDU`_  ...        All data from the :class:`~pypeit.core.flexure.MultiSlitFlexure` datamodel
 =====================  ==============================  =========  ==========================================================================
 
+Inspecting
+==========
+
+PypeIt provides the ``pypeit_chk_flexure`` script to inspect both the
+spatial and spectral flexure corrections.
+
+The script usage can be displayed by calling the script with the
+``-h`` option:
+
+.. include:: ../help/pypeit_chk_flexure.rst
+
+The script takes as input one or multiple `spec2d*.fits` or `spec1d*.fits` files
+and print to screen the flexure correction applied to each file.
+
+Here is a typical call to print the spatial flexure correction:
+
+.. code-block:: console
+
+    pypeit_chk_flexure Science/spec2d_r230417_01033-frb22022_LRISr_20230417T082242.672.fits --spat
+
+and to print the spectral flexure correction:
+
+.. code-block:: console
+
+    pypeit_chk_flexure Science/spec1d_r230417_01033-frb22022_LRISr_20230417T082242.672.fits --spec
 

doc/qa.rst

@@ -111,11 +111,26 @@ Exposure QA
 For each processed, science exposure there are a series of
 PNGs generated, per detector and (sometimes) per slit.
 
+.. _qa-spat-flex:
 
-Flexure QA
-----------
+Spatial Flexure QA
+------------------
 
-If a flexure correction was performed (default), the fit to the
+If a spatial flexure correction was performed, the result of the correction
+is shown in a plot like the one below.  The plot shows a few snippets of the
+science/standard spectral image with overlaid the slit edges as traced in the
+``trace`` image (dashed lines) and after applying the spatial flexure correction
+(solid lines). The value of the shift is also reported on the top of the plot.
+
+.. figure:: figures/r230417_01033_DET01_spat_flex_corr.png
+   :align: center
+
+.. _qa-spec-flex:
+
+Spectral Flexure QA
+-------------------
+
+If a spectral flexure correction was performed (default), the fit to the
 correlation lags per object
 is shown and the adopted shift is listed.  Here is
 an example:
@@ -133,3 +148,6 @@ Here is an example:
 .. figure:: figures/qa/flex_sky_armlsd.jpg
    :align: center
 
+
+
+

doc/releases/1.17.2dev.rst

@@ -17,6 +17,11 @@ Functionality/Performance Improvements and Additions
 - The WCS for datacubes now adopts the convention of North
   is up and East is left. In previous version of PypeIt,
   East was right.
+- Improved the code to compute the spatial flexure correction and added a QA
+  plot to inspect the results.  Two new parset parameters, `spat_flexure_sigdetect`
+  and `spat_flexure_vrange`, were added to control the detection of the slit edges
+  (used to compute the flexure) and the stretching of the image in the QA plot,
+  respectively.
 - The spatial flexure can now take a constant value for every slit, independent values for each slit,
   or independent values for each slit edge.
 
@@ -43,4 +48,5 @@ Bug Fixes
   in both spatial dimensions equal to half the field of view.
 - Fix the code for the extraction of the 1D flat spectrum, so that
   the spectrum is extracted even when `pixelflat_model` does not exist.
-
+- The variance cube previously had a cubic term, and this has been changed
+  to a quadratic term. This has been cross-checked with simple subpixel calculations.

doc/scripts.rst

@@ -586,6 +586,30 @@ The script usage can be displayed by calling the script with the
 
 .. include:: help/pypeit_show_pixflat.rst
 
+.. _pypeit_chk_flexure:
+
+pypeit_chk_flexure
+--------------------
+
+This script takes as input one or multiple `spec2d*.fits` or `spec1d*.fits` files
+and print to screen the spatial or spectral flexure correction applied to each file.
+Standard command-line calls are:
+
+.. code-block:: console
+
+    pypeit_chk_flexure Science/spec2d_r230417_01033-frb22022_LRISr_20230417T082242.672.fits --spat
+
+or:
+
+.. code-block:: console
+
+    pypeit_chk_flexure Science/spec1d_r230417_01033-frb22022_LRISr_20230417T082242.672.fits --spec
+
+The script usage can be displayed by calling the script with the
+``-h`` option:
+
+.. include:: help/pypeit_chk_flexure.rst
+
 pypeit_show_2dspec
 ------------------
 

pypeit/core/arc.py

@@ -1001,10 +1001,10 @@ def detect_lines(censpec, sigdetect=5.0, fwhm=4.0, fit_frac_fwhm=1.25, input_thr
         (mean, med, stddev) = stats.sigma_clipped_stats(arc[cont_mask & np.logical_not(bpm_out)], 
                                                         sigma_lower=3.0, sigma_upper=3.0, cenfunc= np.nanmedian,
                                                         stdfunc = np.nanstd)
-        thresh = med + sigdetect*stddev
         if stddev == 0.0:
-            msgs.warn('stddev = 0.0, so resetting to 1.0')
-            stddev = 1.0
+            msgs.warn('stddev = 0.0, so resetting to 0.1')
+            stddev = 0.1
+        thresh = med + sigdetect * stddev
     else:
         med = 0.0
         if isinstance(input_thresh,(float, int)):
@@ -1055,6 +1055,7 @@ def detect_lines(censpec, sigdetect=5.0, fwhm=4.0, fit_frac_fwhm=1.25, input_thr
             tampl_true = tampl_true[ikeep]
             tampl = tampl[ikeep]
             tcent = tcent[ikeep]
+            pixt = pixt[ikeep]
             twid = twid[ikeep]
             centerr = centerr[ikeep]
             ww = np.where(good[ikeep])[0]

pypeit/core/datacube.py

@@ -1822,7 +1822,7 @@ def subpixellate(output_wcs, bins, sciImg, ivarImg, waveImg, slitid_img_gpm, wgh
             vox_coord = vox_coord.reshape(numpix * num_all_subpixels, 3)
             # Use the "fast histogram" algorithm, that assumes regular bin spacing
             flxcube += histogramdd(vox_coord, bins=outshape, range=binrng, weights=np.repeat(this_sci[this_sl] * this_wght_subpix[this_sl], num_all_subpixels) * subpix_wght)
-            varcube += histogramdd(vox_coord, bins=outshape, range=binrng, weights=np.repeat(this_var[this_sl] * this_wght_subpix[this_sl]**2, num_all_subpixels) * subpix_wght**3)
+            varcube += histogramdd(vox_coord, bins=outshape, range=binrng, weights=np.repeat(this_var[this_sl] * this_wght_subpix[this_sl]**2, num_all_subpixels) * subpix_wght**2)  # NOTE :: This was changed from subpix_wght**3 to subpix_wght**2 by RJC on 2024-12-18
             normcube += histogramdd(vox_coord, bins=outshape, range=binrng, weights=np.repeat(this_wght_subpix[this_sl], num_all_subpixels) * subpix_wght)
 
     # Normalise the datacube and variance cube