diff --git a/doc/calibrations/flexure.rst b/doc/calibrations/flexure.rst index a53f8b462..467444e60 100644 --- a/doc/calibrations/flexure.rst +++ b/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 diff --git a/doc/figures/r230417_01033_DET01_spat_flex_corr.png b/doc/figures/r230417_01033_DET01_spat_flex_corr.png new file mode 100644 index 000000000..337435e95 Binary files /dev/null and b/doc/figures/r230417_01033_DET01_spat_flex_corr.png differ diff --git a/doc/qa.rst b/doc/qa.rst index f895fc600..b96843fa2 100644 --- a/doc/qa.rst +++ b/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 + + + diff --git a/doc/releases/1.17.2dev.rst b/doc/releases/1.17.2dev.rst index 1f647d987..fdb1426a1 100644 --- a/doc/releases/1.17.2dev.rst +++ b/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. diff --git a/doc/scripts.rst b/doc/scripts.rst index e5c4aee3d..392529f84 100644 --- a/doc/scripts.rst +++ b/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 ------------------ diff --git a/pypeit/core/arc.py b/pypeit/core/arc.py index e41fc7fc8..810f20764 100644 --- a/pypeit/core/arc.py +++ b/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] diff --git a/pypeit/core/datacube.py b/pypeit/core/datacube.py index bcca07d36..917a6b871 100644 --- a/pypeit/core/datacube.py +++ b/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 diff --git a/pypeit/core/flexure.py b/pypeit/core/flexure.py index 1cae726fc..0ad726b8d 100644 --- a/pypeit/core/flexure.py +++ b/pypeit/core/flexure.py @@ -12,11 +12,13 @@ from matplotlib import pyplot as plt from matplotlib import gridspec +from matplotlib.lines import Line2D import matplotlib from astropy import stats from astropy import units from astropy.io import ascii +from astropy.stats import sigma_clipped_stats import scipy.signal import scipy.optimize as opt from scipy import interpolate @@ -26,13 +28,12 @@ from pypeit import msgs from pypeit import dataPaths from pypeit import io -from pypeit import utils -from pypeit.display import display from pypeit.core.wavecal import autoid from pypeit.core import arc from pypeit.core import extract from pypeit.core import fitting from pypeit.core import qa +from pypeit.core import trace from pypeit.datamodel import DataContainer from pypeit.images.detector_container import DetectorContainer from pypeit.images.mosaic import Mosaic @@ -42,7 +43,7 @@ from IPython import embed -def spat_flexure_shift(sciimg, slits, method="detector", maxlag=20, debug=False): +def spat_flexure_shift(sciimg, slits, method="detector", bpm=None, maxlag=20, sigdetect=10., debug=False, qa_outfile=None, qa_vrange=None): """ Calculate a rigid flexure shift in the spatial dimension between the slitmask and the science image. @@ -64,10 +65,21 @@ def spat_flexure_shift(sciimg, slits, method="detector", maxlag=20, debug=False) 'slit' method calculates the shift for each slit independently, and the 'edge' method calculates the shift for each slit edge independently. + Slits object + bpm (`numpy.ndarray`_, optional): + Bad pixel mask (True = Bad) maxlag (:obj:`int`, optional): Maximum flexure searched for + sigdetect (:obj:`float`, optional): + Sigma threshold above fluctuations for the slit detection + in the collapsed sobel image debug (:obj:`bool`, optional): Run in debug mode + qa_outfile (:obj:`str`, optional): + Path to the output file where the QA is saved. If None, the QA is not generated. + qa_vrange (:obj:`tuple`, optional): + Tuple with the vmin and vmax values for the imshow plot in the QA. If None, the + vmin and vmax values are calculated from the data. Returns: float: The spatial flexure shift relative to the initial slits @@ -75,59 +87,212 @@ def spat_flexure_shift(sciimg, slits, method="detector", maxlag=20, debug=False) """ # TODO :: Need to implement different methods + msgs.info("Measuring spatial flexure") # Mask -- Includes short slits and those excluded by the user (e.g. ['rdx']['slitspatnum']) slitmask = slits.slit_img(initial=True, exclude_flag=slits.bitmask.exclude_for_flexure) _sciimg = sciimg if slitmask.shape == sciimg.shape \ - else arc.resize_mask2arc(slitmask.shape, sciimg) - onslits = slitmask > -1 - corr_slits = onslits.astype(float).flatten() - - # Compute - mean_sci, med_sci, stddev_sci = stats.sigma_clipped_stats(_sciimg[onslits]) - thresh = med_sci + 5.0*stddev_sci - corr_sci = np.fmin(_sciimg.flatten(), thresh) - lags, xcorr = utils.cross_correlate(corr_sci, corr_slits, maxlag) - xcorr_denom = np.sqrt(np.sum(corr_sci*corr_sci)*np.sum(corr_slits*corr_slits)) - xcorr_norm = xcorr / xcorr_denom - # TODO -- Generate a QA plot - - tampl_true, tampl, pix_max, twid, centerr, ww, arc_cont, nsig \ - = arc.detect_lines(xcorr_norm, sigdetect=3.0, fit_frac_fwhm=1.5, fwhm=5.0, - cont_frac_fwhm=1.0, cont_samp=30, nfind=1, debug=debug) + else arc.resize_mask2arc(slitmask.shape, sciimg) + + # mask (as much as possible) the objects on the slits to help the cross-correlation + # need to copy the bpm to avoid changing the input bpm + _bpm = np.zeros_like(_sciimg, dtype=int) if bpm is None else copy.deepcopy(bpm) + for i in range(slits.nslits): + left_edge = np.round(slits.left_init[:, i]).astype(int) + right_edge = np.round(slits.right_init[:, i]).astype(int) + for j in range(_sciimg.shape[0]): + # mask the region between the left and right edges leaving a margin of maxlag pixels + if left_edge[j]+maxlag < right_edge[j]-maxlag: + _bpm[j, left_edge[j]+maxlag:right_edge[j]-maxlag] = 1 + + # # create sobel images of both slitmask and the science image + sci_sobel, sci_edges = trace.detect_slit_edges(_sciimg, bpm=_bpm, sigdetect=sigdetect) + slits_sobel, slits_edges = trace.detect_slit_edges(slitmask, bpm=bpm, sigdetect=1.) + corr = scipy.signal.fftconvolve(sci_edges, np.fliplr(slits_edges), mode='same', axes=1) + xcorr = np.sum(corr, axis=0) + lags = scipy.signal.correlation_lags(sci_edges.shape[1], slits_edges.shape[1], mode='same') + lag0 = np.where(lags == 0)[0][0] + xcorr_max = xcorr[lag0 - maxlag:lag0 + maxlag] + lags_max = lags[lag0 - maxlag:lag0 + maxlag] + + # detect the highest peak in the cross-correlation + _, _, pix_max, _, _, _, _, _ = arc.detect_lines(xcorr_max, cont_subtract=False, input_thresh=0., nfind=1, debug=debug) # No peak? -- e.g. data fills the entire detector - if len(tampl) == 0: - msgs.warn('No peak found in spatial flexure. Assuming there is none...') + if (len(pix_max) == 0) or pix_max[0] == -999.0: + msgs.warn('No peak found in the x-correlation between the traced slits and the science/calib image.' + ' Assuming there is NO SPATIAL FLEXURE.'+msgs.newline() + 'If a flexure is expected, ' + 'consider either changing the maximum lag for the cross-correlation, ' + 'or the "spat_flexure_sigdetect" parameter, or use the manual flexure correction.') + return np.zeros((slits.nslits, 2), dtype=float) - # Find the peak - xcorr_max = np.interp(pix_max, np.arange(lags.shape[0]), xcorr_norm) - lag_max = np.interp(pix_max, np.arange(lags.shape[0]), lags) - msgs.info('Spatial flexure measured: {}'.format(lag_max[0])) + lag0_max = np.where(lags_max == 0)[0][0] + shift = round(pix_max[0] - lag0_max, 3) + msgs.info('Spatial flexure measured: {} pixels'.format(shift)) if debug: - plt.figure(figsize=(14, 6)) - plt.plot(lags, xcorr_norm, color='black', drawstyle='steps-mid', lw=3, label='x-corr') - plt.plot(lag_max[0], xcorr_max[0], 'g+', markersize=6.0, label='peak') - plt.title('Best shift = {:5.3f}'.format(lag_max[0]) + ', corr_max = {:5.3f}'.format(xcorr_max[0])) + # 1D plot of the cross-correlation + plt.figure(figsize=(10, 6)) + plt.minorticks_on() + plt.tick_params(axis='both', direction='in', top=True, right=True, which='both') + # plot xcorr_max but add a buffer of 20 pixels on each side + pad = 20 + _xcorr_max = xcorr[lag0 - (maxlag+pad):lag0 + (maxlag+pad)] + _lags_max = lags[lag0 - (maxlag+pad):lag0 + (maxlag+pad)] + plt.plot(_lags_max, _xcorr_max, 'k-', lw=1) + plt.axvline(shift, color='r', linestyle='--', label=f'Measured shift = {shift:.1f} pixels') + plt.axvline(maxlag, color='g', linestyle='--', label='Max lag') + plt.axvline(-maxlag, color='g', linestyle='--') + plt.xlabel('Lag (pixels)') + plt.ylabel('Cross-correlation') + plt.title('Spatial Flexure Cross-correlation') plt.legend() + plt.tight_layout() plt.show() - #tslits_shift = trace_slits.shift_slits(tslits_dict, lag_max) - # Now translate the tilts + # 2D plot + spat_flexure_qa(sciimg, slits, shift, gpm=np.logical_not(bpm), vrange=qa_vrange) + + if qa_outfile is not None: + # Generate the QA plot + msgs.info("Generating QA plot for spatial flexure") + spat_flexure_qa(sciimg, slits, shift, gpm=np.logical_not(bpm), vrange=qa_vrange, outfile=qa_outfile) + + return np.full((slits.nslits, 2), shift) + + +def spat_flexure_qa(img, slits, shift, gpm=None, vrange=None, outfile=None): + """ + Generate QA for the spatial flexure + Args: + img (`numpy.ndarray`_): + Image of the detector + slits (:class:`pypeit.slittrace.SlitTraceSet`): + Slits object + shift (:obj:`float`): + Shift in pixels + gpm (`numpy.ndarray`_, optional): + Good pixel mask (True = Bad) + vrange (:obj:`tuple`, optional): + Tuple with the min and max values for the imshow plot + outfile (:obj:`str`, optional): + Path to the output file where the QA is saved. If None, the QA is + shown on screen and not saved. + + + """ + debug = True if outfile is None else False + + # check that vrange is a tuple + if vrange is not None and not isinstance(vrange, tuple): + msgs.warn('vrange must be a tuple with the min and max values for the imshow plot. Ignoring vrange.') + vrange = None - #slitmask_shift = pixels.tslits2mask(tslits_shift) - #slitmask_shift = slits.slit_img(flexure=lag_max[0]) + # TODO: should we use initial or tweaked slits in this plot? + left_slits, right_slits, mask_slits = slits.select_edges(initial=True, flexure=None) + left_flex, right_flex, mask = slits.select_edges(initial=True, flexure=shift) if debug: - # Now translate the slits in the tslits_dict - all_left_flexure, all_right_flexure, mask = slits.select_edges(spat_flexure=lag_max[0]) - gpm = mask == 0 - viewer, ch = display.show_image(_sciimg) - #display.show_slits(viewer, ch, left_flexure[:,gpm], right_flexure)[:,gpm]#, slits.id) #, args.det) - #embed(header='83 of flexure.py') - - return np.full((slits.nslits, 2), lag_max[0]) + # where to start and end the plot in the spatial&spectral direction + nxsnip = 1 + spat_starts = [0] + spat_ends = [img.shape[1]] + upper_ystart = 0 + upper_yend = img.shape[0] + + else: + # where to start and end the plot in the spatial direction + xstart = int(np.floor(np.min([left_slits, left_flex]) - 20)) + xend = int(np.ceil(np.max([right_slits, right_flex]) + 20)) + + # how many snippets to plot in the spatial direction + if slits.nslits == 1: + # if longslit plot 2 snippets, one for the left edge and one for the right edge + nxsnip = 2 + snippet = int((xend - xstart) // nxsnip) + spat_starts = [xstart, xstart + snippet] + spat_ends = [xend - snippet, xend] + elif slits.nslits <= 12: + # if 12 or less slits plot 3-4 snippets equally spaced + nxsnip = 3 if slits.nslits <= 6 else 4 + snippet = int((xend - xstart) // nxsnip) + spat_starts = [xstart, xstart + snippet, xstart + 2*snippet] + spat_ends = [xend - 2*snippet, xend - snippet, xend] + if slits.nslits > 6: + # add the 4th snippet + spat_starts.append(xstart + 3*snippet) + spat_ends.insert(0, xend - 3*snippet) + else: + # if more than 12 slits plot 4 snippets + nxsnip = 4 + # approximately, we want 3 slits in each snippet + snippet = int(3 * (xend - xstart)/slits.nslits) + # this would give nx many snippets + nx = int((xend - xstart) // snippet) + # but we want to plot only nxsnip of those snippets + spat_starts = [xstart + i * snippet for i in np.linspace(0, nx - 1, nxsnip, dtype=int)] + spat_ends = [xstart + i * snippet for i in np.linspace(1, nx, nxsnip, dtype=int)] + + # where to start and end the plot in the spectral direction for both the upper and lower sections + lower_ystart = 0 + lower_yend = int(snippet) + upper_ystart = int(img.shape[0] - snippet) + upper_yend = img.shape[0] + + # plot the spatial flexure + rows = 1 if debug else 2 + fig = plt.figure(figsize=(9, 8) if debug else (nxsnip*4, 8)) + gs = gridspec.GridSpec(rows, nxsnip, figure=fig) + # spectral vector for plotting the slits + spec = np.tile(np.arange(slits.nspec), (slits.nslits, 1)).T + thin = 10 + # legend elements + legend_elements = [Line2D([0], [0], color='C3', lw=1, ls='--', label='initial left edges'), + Line2D([0], [0], color='C1', lw=1, ls='--', label='initial right edges'), + Line2D([0], [0], color='C3', lw=1, label='shifted left edges'), + Line2D([0], [0], color='C1', lw=1, label='shifted right edges')] + # loop over the 2 rows if we save the plot in the output directory, otherwise plot the whole detector + for r in range(rows): + _ystar, _yend = (upper_ystart, upper_yend) if r == 0 else (lower_ystart, lower_yend) + # loop over the snippets + for s in range(nxsnip): + ax = fig.add_subplot(gs[r, s]) + if vrange is None: + # get vmin and vmax for imshow + _xstart = spat_starts[s] if spat_starts[s] >= 0 else 0 + _xend = spat_ends[s] if spat_ends[s] <= img.shape[1] else img.shape[1] + _img = img[_ystar:_yend, _xstart:_xend] + _gpm = gpm[_ystar:_yend, _xstart:_xend] if gpm is not None else np.ones_like(_img, dtype=bool) + m, med, sig = sigma_clipped_stats(_img[_gpm], sigma_lower=5.0, sigma_upper=5.0) + vmin = m - 1.0 * sig + vmax = m + 4.0 * sig + else: + vmin, vmax = vrange + # imshow img instead of _img to show the actual pixel values in each snippet + ax.imshow(img, origin='lower', vmin=vmin, vmax=vmax) + ax.set_ylim(_ystar, _yend) + ax.set_xlim(spat_starts[s], spat_ends[s]) + + # plot the slits + for i in range(slits.nslits): + plt.plot(left_slits[::thin, i], spec[::thin, i], color='C3', lw=1, ls='--', zorder=5) + plt.plot(right_slits[::thin, i], spec[::thin, i], color='C1', lw=1, ls='--', zorder=5) + plt.plot(left_flex[::thin, i], spec[::thin, i], color='C3', lw=1, zorder=6) + plt.plot(right_flex[::thin, i], spec[::thin, i], color='C1', lw=1, zorder=6) + ax.tick_params(axis='both', labelsize=6) + if r == 0 and s == 0: + plt.suptitle(f'Shift={shift:.1f} pixels', fontsize=18) + ax.legend(handles=legend_elements, fontsize=7) + if not debug: + ax.set_ylabel('Upper snippets', fontsize=18) + elif r == 1 and s == 0: + ax.set_ylabel('Lower snippets', fontsize=18) + plt.tight_layout() + if debug: + plt.show() + else: + fig.savefig(outfile, dpi=200) + plt.close(fig) def spec_flex_shift(obj_skyspec, sky_file=None, arx_skyspec=None, arx_fwhm_pix=None, diff --git a/pypeit/core/qa.py b/pypeit/core/qa.py index 8c8e498e1..a7929f2cb 100644 --- a/pypeit/core/qa.py +++ b/pypeit/core/qa.py @@ -79,6 +79,8 @@ def set_qa_filename(root, method, det=None, slit=None, prefix=None, out_dir=None outfile = 'QA/PNGs/{:s}_{:s}_obj_trace.png'.format(root, det) elif method == 'obj_profile_qa': outfile = 'PNGs/{:s}_{:s}_S{:04d}_obj_prof.png'.format(root, det, slit) + elif method == 'spat_flexure_qa_corr': + outfile = 'QA/PNGs/{:s}_spat_flex_corr.png'.format(root) elif method == 'spec_flexure_qa_corr': # outfile = 'QA/PNGs/{:s}_D{:02d}_S{:04d}_spec_flex_corr.png'.format(root, det, slit) outfile = 'PNGs/{:s}_S{:04d}_spec_flex_corr.png'.format(root, slit) diff --git a/pypeit/core/wavecal/wvutils.py b/pypeit/core/wavecal/wvutils.py index d5b44af57..d100d5c38 100644 --- a/pypeit/core/wavecal/wvutils.py +++ b/pypeit/core/wavecal/wvutils.py @@ -458,8 +458,8 @@ def zerolag_shift_stretch(theta, y1, y2, stretch_func = 'quadratic'): return -corr_norm - -def get_xcorr_arc(inspec1, sigdetect=5.0, sig_ceil=10.0, percent_ceil=50.0, use_raw_arc=False, fwhm = 4.0, debug=False): +def get_xcorr_arc(inspec1, sigdetect=5.0, input_thresh=None, sig_ceil=10.0, percent_ceil=50.0, use_raw_arc=False, + fwhm=4.0, cont_sub=True, debug=False): """ Utility routine to create a synthetic arc spectrum for cross-correlation using the location of the peaks in the input spectrum. @@ -468,7 +468,10 @@ def get_xcorr_arc(inspec1, sigdetect=5.0, sig_ceil=10.0, percent_ceil=50.0, use_ inspec1 (`numpy.ndarray`_): Input spectrum, shape = (nspec,) sigdetect (float, optional, default=3.0): - Peak finding threshold for lines that will be used to create the synthetic xcorr_arc + Sigma threshold above fluctuations for finding peaks that will be used to create the synthetic xcorr_arc + input_thresh (float, optional): + Input threshold for finding peaks that will be used to create the synthetic xcorr_arc. If set, sigdetect + will be ignored. sig_ceil (float, optional, default = 10.0): Significance threshold for peaks that will be used to determine the line amplitude clipping threshold. For peaks with significance > sig_ceil, the code will find the amplitude corresponding to @@ -480,6 +483,8 @@ def get_xcorr_arc(inspec1, sigdetect=5.0, sig_ceil=10.0, percent_ceil=50.0, use_ If True, use amplitudes from the raw arc, i.e. do not continuum subtract. Default = False fwhm (float, optional): Fwhm of arc lines. Used for peak finding and to assign a fwhm in the xcorr_arc. + cont_sub (bool, optional): + Perform a simple continuum subtraction when detecting the peaks. Default is True. debug (bool, optional): Show plots for line detection debugging. Default = False @@ -492,7 +497,9 @@ def get_xcorr_arc(inspec1, sigdetect=5.0, sig_ceil=10.0, percent_ceil=50.0, use_ # Run line detection to get the locations and amplitudes of the lines tampl1, tampl1_cont, tcent1, twid1, centerr1, w1, arc1, nsig1 = arc.detect_lines(inspec1, sigdetect=sigdetect, - fwhm=fwhm, debug=debug) + input_thresh=input_thresh, + fwhm=fwhm, cont_subtract=cont_sub, + debug=debug) ampl = tampl1 if use_raw_arc else tampl1_cont @@ -506,7 +513,7 @@ def get_xcorr_arc(inspec1, sigdetect=5.0, sig_ceil=10.0, percent_ceil=50.0, use_ ampl_clip = np.clip(ampl, None, ceil_upper) if ampl_clip.size == 0: msgs.warn('No lines were detected in the arc spectrum. Cannot create a synthetic arc spectrum for cross-correlation.') - return None + return np.zeros_like(inspec1) # Make a fake arc by plopping down Gaussians at the location of every centroided line we found xcorr_arc = np.zeros_like(inspec1) @@ -519,7 +526,6 @@ def get_xcorr_arc(inspec1, sigdetect=5.0, sig_ceil=10.0, percent_ceil=50.0, use_ continue xcorr_arc += ampl_clip[ind]*np.exp(-0.5*((spec_vec - tcent1[ind])/sigma)**2) - return xcorr_arc @@ -750,7 +756,7 @@ def xcorr_shift_stretch(inspec1, inspec2, cc_thresh=-1.0, percent_ceil=50.0, use y2 = get_xcorr_arc(inspec2, percent_ceil=percent_ceil, use_raw_arc=use_raw_arc, sigdetect=sigdetect, sig_ceil=sig_ceil, fwhm=fwhm) - if y1 is None or y2 is None: + if np.all(y1 == 0) or np.all(y2 == 0): msgs.warn('No lines detected punting on shift/stretch') return 0, None, None, None, None, None, None diff --git a/pypeit/images/combineimage.py b/pypeit/images/combineimage.py index 31980b00b..26c0a61a1 100644 --- a/pypeit/images/combineimage.py +++ b/pypeit/images/combineimage.py @@ -160,9 +160,12 @@ def run(self, ignore_saturation=False, maxiters=5): basev_stack = np.zeros(shape, dtype=float) gpm_stack = np.zeros(shape, dtype=bool) exptime = np.zeros(self.nimgs, dtype=float) + spat_flex = np.zeros(self.nimgs, dtype=float) # Save the exposure time to check if it's consistent for all images. exptime[kk] = rawImage.exptime + # Save the spatial flexure to check if it's consistent for all images and propagate it to the combined image + spat_flex[kk] = rawImage.spat_flexure # Processed image img_stack[kk] = rawImage.image # Get the count scaling @@ -194,6 +197,19 @@ def run(self, ignore_saturation=False, maxiters=5): else: comb_texp = exptime[0] + # Check that all spatial flexure values are consistent + comb_spat_flex = None + # remove nan (None) values. Since spat_flex is a float array, + # if rawImage.spat_flexure is None, it will be converted to nan + no_nan = np.logical_not(np.isnan(spat_flex)) + if np.sum(no_nan) > 0: + if np.any(np.absolute(np.diff(spat_flex[no_nan])) > 0.1): + msgs.warn(f'Spatial flexure is not consistent for all images being combined: {spat_flex}.') + comb_spat_flex = np.round(np.mean(spat_flex[no_nan]),3) + msgs.warn(f'Using the average: {comb_spat_flex}.') + else: + comb_spat_flex = spat_flex[no_nan][0] + # scale the images to their mean, if requested, before combining if self.par['scale_to_mean']: msgs.info("Scaling images to have the same mean before combining") @@ -274,6 +290,7 @@ def run(self, ignore_saturation=False, maxiters=5): # NOTE: The detector is needed here so # that we can get the dark current later. detector=rawImage.detector, + spat_flexure=comb_spat_flex, PYP_SPEC=rawImage.PYP_SPEC, units='e-' if self.par['apply_gain'] else 'ADU', exptime=comb_texp, noise_floor=self.par['noise_floor'], diff --git a/pypeit/images/rawimage.py b/pypeit/images/rawimage.py index 4e5bcd225..1dcf24712 100644 --- a/pypeit/images/rawimage.py +++ b/pypeit/images/rawimage.py @@ -20,11 +20,13 @@ from pypeit.core import flat from pypeit.core import flexure from pypeit.core import scattlight +from pypeit.core import qa from pypeit.core.mosaic import build_image_mosaic from pypeit.images import pypeitimage from pypeit import utils from pypeit.display import display - +from pypeit import io +from pathlib import Path # TODO: I don't understand why we have some of these attributes. E.g., why do # we need both hdu and headarr? @@ -678,7 +680,8 @@ def process(self, par, bpm=None, scattlight=None, flatimages=None, bias=None, sl 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']) + maxlag=self.par['spat_flexure_maxlag'], + debug=debug) # - Subtract scattered light... this needs to be done before flatfielding. if self.par['subtract_scattlight']: @@ -771,7 +774,7 @@ 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, manual_spat_flexure=np.ma.masked, method="detector", maxlag=20): + def spatial_flexure_shift(self, slits, force=False, manual_spat_flexure=np.ma.masked, method="detector", debug=False): """ Calculate a spatial shift in the edge traces due to flexure. @@ -798,8 +801,8 @@ def spatial_flexure_shift(self, slits, force=False, manual_spat_flexure=np.ma.ma 'slit' method calculates the shift for each slit independently, and the 'edge' method calculates the shift for each slit edge independently. - maxlag (:obj:'float', optional): - Maximum range of lag values over which to compute the CCF. + debug (:obj:`bool`, optional): + Run in debug mode. Return: `numpy.ndarray`_: The calculated flexure correction for the edge of each slit shape is (nslits, 2) @@ -827,7 +830,16 @@ def spatial_flexure_shift(self, slits, force=False, manual_spat_flexure=np.ma.ma 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) + # get filename for QA + basename = f'{io.remove_suffix(self.filename)}_{self.spectrograph.get_det_name(self.det)}' + outdir = str(Path(slits.calib_dir).parent) if slits.calib_dir is not None else None + qa_outfile = qa.set_qa_filename(basename, 'spat_flexure_qa_corr', out_dir=outdir) + + spat_flexure = flexure.spat_flexure_shift(self.image[0], slits, method=method, bpm=self._bpm[0], + maxlag=self.par['spat_flexure_maxlag'], + sigdetect=self.par['spat_flexure_sigdetect'], + debug=debug, qa_outfile=qa_outfile, + qa_vrange=self.par['spat_flexure_vrange']) # Print the flexure values if np.all(spat_flexure == spat_flexure[0, 0]): diff --git a/pypeit/par/pypeitpar.py b/pypeit/par/pypeitpar.py index 447162005..4df12bdd9 100644 --- a/pypeit/par/pypeitpar.py +++ b/pypeit/par/pypeitpar.py @@ -222,7 +222,8 @@ def __init__(self, trim=None, apply_gain=None, orient=None, empirical_rn=None, shot_noise=None, noise_floor=None, use_pixelflat=None, use_illumflat=None, use_specillum=None, use_pattern=None, subtract_scattlight=None, scattlight=None, subtract_continuum=None, - spat_flexure_method=None, spat_flexure_maxlag=None): + spat_flexure_method=None, spat_flexure_maxlag=None, + spat_flexure_sigdetect=None, spat_flexure_vrange=None): # Grab the parameter names and values from the function # arguments @@ -370,6 +371,19 @@ def __init__(self, trim=None, apply_gain=None, orient=None, dtypes['spat_flexure_maxlag'] = int descr['spat_flexure_maxlag'] = 'Maximum of possible spatial flexure correction, in pixels' + defaults['spat_flexure_sigdetect'] = 5. + dtypes['spat_flexure_sigdetect'] = [int, float] + descr['spat_flexure_sigdetect'] = 'Sigma threshold above fluctuations in the ' \ + 'Sobel-filtered significance image, used for ' \ + 'finding slit edges in the spectral image, ' \ + 'for which the spatial flexure is computed.' + + defaults['spat_flexure_vrange'] = None + dtypes['spat_flexure_vrange'] = tuple + descr['spat_flexure_vrange'] = 'This parameter is used when generating the QA plot for the spatial flexure. ' \ + 'It sets the data range (vmin,vmax) used by the colormap when showing the ' \ + 'spectral image. If None, the range is set automatically.' + defaults['combine'] = 'mean' options['combine'] = ProcessImagesPar.valid_combine_methods() dtypes['combine'] = str @@ -467,7 +481,8 @@ def from_dict(cls, cfg): parkeys = ['trim', 'apply_gain', 'orient', 'use_biasimage', 'subtract_continuum', 'subtract_scattlight', 'scattlight', 'use_pattern', 'use_overscan', 'overscan_method', 'overscan_par', 'use_darkimage', 'dark_expscale', - 'spat_flexure_method', 'spat_flexure_maxlag', 'use_illumflat', 'use_specillum', + 'spat_flexure_method', 'spat_flexure_maxlag', 'spat_flexure_sigdetect', + 'spat_flexure_vrange', 'use_illumflat', 'use_specillum', 'empirical_rn', 'shot_noise', 'noise_floor', 'use_pixelflat', 'combine', 'scale_to_mean', 'correct_nonlinear', 'satpix', #'calib_setup_and_bit', 'n_lohi', 'mask_cr', 'lamaxiter', 'grow', 'clip', 'comb_sigrej', 'rmcompact', diff --git a/pypeit/spectrographs/keck_lris.py b/pypeit/spectrographs/keck_lris.py index 2b63aa1c7..2ce95ea22 100644 --- a/pypeit/spectrographs/keck_lris.py +++ b/pypeit/spectrographs/keck_lris.py @@ -150,6 +150,10 @@ def config_specific_par(self, scifile, inp_par=None): # Arc lamps list from header par['calibrations']['wavelengths']['lamps'] = ['use_header'] + # spatial flexure maxshift for non-longslit + if 'long' not in self.get_meta_value(scifile, 'decker'): + par['scienceframe']['process']['spat_flexure_maxlag'] = 10 + return par def init_meta(self):