diff --git a/doc/calibrations/flexure.rst b/doc/calibrations/flexure.rst index 55bac8e218..467444e60e 100644 --- a/doc/calibrations/flexure.rst +++ b/doc/calibrations/flexure.rst @@ -26,11 +26,22 @@ Spatial Flexure The code has a simple yet relatively robust method to cross-correlate the slits against any input image to determine a rigid, spatial offset. This algorithm is performed for any frame type by setting -``spat_flexure_correct = True`` in the ``process`` block -of :ref:`processimagespar`. +``spat_flexure_method = `` in the ``process`` block +of :ref:`processimagespar`, where you should replace ```` with one of the +following valid methods: +``skip`` - this will skip the spatial flexure correction +``detector`` - this will calculate one spatial flexure value for all slits on the detector +``slit`` - this will calculate the spatial flexure for each slit separately +``edge`` - this will calculate the spatial flexure for each slit edge separately We have made our own determinations for which instruments to enable this as the default. Inspect the :ref:`instr_par` list to see if your instrument is included +(search for the value of ``spat_flexure_method``). We recommend that you use the +``detector`` method for most instruments, unless you have good reason to believe +that the spatial flexure is different for different slits or slit edges. Also note +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 @@ -52,11 +63,11 @@ add this to your PypeIt file: [scienceframe] [[process]] - spat_flexure_correct = True + spat_flexure_method = detector [calibrations] [[standardframe]] [[[process]]] - spat_flexure_correct = True + spat_flexure_method = detector This will: @@ -78,7 +89,7 @@ following to your :doc:`../pypeit_file`: [calibrations] [[tiltframe]] [[[process]]] - spat_flexure_correct = True + spat_flexure_method = detector This will: diff --git a/doc/calibrations/image_proc.rst b/doc/calibrations/image_proc.rst index b80caecf80..b6785573e0 100644 --- a/doc/calibrations/image_proc.rst +++ b/doc/calibrations/image_proc.rst @@ -288,8 +288,8 @@ Spatial Flexure Shift A spatial shift in the slit positions due to instrument flexure is calculated using :func:`~pypeit.core.flexure.spat_flexure_shift` if the -``spat_flexure_correct`` parameter is true. See :ref:`flexure` for additional -discussion. +``spat_flexure_method`` parameter is set to an option other than ``skip``. +See :ref:`flexure` for additional discussion. Flat-fielding ------------- diff --git a/doc/include/class_datamodel_pypeitimage.rst b/doc/include/class_datamodel_pypeitimage.rst index 9008ab7bc3..ecc344bfab 100644 --- a/doc/include/class_datamodel_pypeitimage.rst +++ b/doc/include/class_datamodel_pypeitimage.rst @@ -1,24 +1,24 @@ -**Version**: 1.3.0 +**Version**: 1.3.1 -================ =================================================================================================== ================= ======================================================================================================================================================================================== -Attribute Type Array Type Description -================ =================================================================================================== ================= ======================================================================================================================================================================================== -``PYP_SPEC`` str PypeIt spectrograph name -``amp_img`` `numpy.ndarray`_ `numpy.integer`_ Provides the amplifier that contributed to each pixel. If this is a detector mosaic, this must be used in combination with ``det_img`` to select pixels for a given detector amplifier. -``base_var`` `numpy.ndarray`_ `numpy.floating`_ Base-level image variance, excluding count shot-noise -``det_img`` `numpy.ndarray`_ `numpy.integer`_ If a detector mosaic, this image provides the detector that contributed to each pixel. -``detector`` :class:`~pypeit.images.detector_container.DetectorContainer`, :class:`~pypeit.images.mosaic.Mosaic` The detector (see :class:`~pypeit.images.detector_container.DetectorContainer`) or mosaic (see :class:`~pypeit.images.mosaic.Mosaic`) parameters -``exptime`` int, float Effective exposure time (s) -``filename`` str Filename for the image -``fullmask`` :class:`~pypeit.images.imagebitmask.ImageBitMaskArray` Image mask -``image`` `numpy.ndarray`_ `numpy.floating`_ Primary image data -``img_scale`` `numpy.ndarray`_ `numpy.floating`_ Image count scaling applied (e.g., 1/flat-field) -``ivar`` `numpy.ndarray`_ `numpy.floating`_ Inverse variance image -``nimg`` `numpy.ndarray`_ `numpy.integer`_ If a combination of multiple images, this is the number of images that contributed to each pixel -``noise_floor`` float Noise floor included in variance -``rn2img`` `numpy.ndarray`_ `numpy.floating`_ Read noise squared image -``shot_noise`` bool Shot-noise included in variance -``spat_flexure`` float Shift, in spatial pixels, between this image and SlitTrace -``units`` str (Unscaled) Pixel units (e- or ADU) -================ =================================================================================================== ================= ======================================================================================================================================================================================== +================ =================================================================================================== ================= ================================================================================================================================================================================================================================ +Attribute Type Array Type Description +================ =================================================================================================== ================= ================================================================================================================================================================================================================================ +``PYP_SPEC`` str PypeIt spectrograph name +``amp_img`` `numpy.ndarray`_ `numpy.integer`_ Provides the amplifier that contributed to each pixel. If this is a detector mosaic, this must be used in combination with ``det_img`` to select pixels for a given detector amplifier. +``base_var`` `numpy.ndarray`_ `numpy.floating`_ Base-level image variance, excluding count shot-noise +``det_img`` `numpy.ndarray`_ `numpy.integer`_ If a detector mosaic, this image provides the detector that contributed to each pixel. +``detector`` :class:`~pypeit.images.detector_container.DetectorContainer`, :class:`~pypeit.images.mosaic.Mosaic` The detector (see :class:`~pypeit.images.detector_container.DetectorContainer`) or mosaic (see :class:`~pypeit.images.mosaic.Mosaic`) parameters +``exptime`` int, float Effective exposure time (s) +``filename`` str Filename for the image +``fullmask`` :class:`~pypeit.images.imagebitmask.ImageBitMaskArray` Image mask +``image`` `numpy.ndarray`_ `numpy.floating`_ Primary image data +``img_scale`` `numpy.ndarray`_ `numpy.floating`_ Image count scaling applied (e.g., 1/flat-field) +``ivar`` `numpy.ndarray`_ `numpy.floating`_ Inverse variance image +``nimg`` `numpy.ndarray`_ `numpy.integer`_ If a combination of multiple images, this is the number of images that contributed to each pixel +``noise_floor`` float Noise floor included in variance +``rn2img`` `numpy.ndarray`_ `numpy.floating`_ Read noise squared image +``shot_noise`` bool Shot-noise included in variance +``spat_flexure`` `numpy.ndarray`_ `numpy.floating`_ Shift, in spatial pixels, between this image and SlitTrace. Shape is (nslits, 2), wherespat_flexure[i,0] is the spatial shift of the left edge of slit i and spat_flexure[i,1] is the spatial shift of the right edge of slit i. +``units`` str (Unscaled) Pixel units (e- or ADU) +================ =================================================================================================== ================= ================================================================================================================================================================================================================================ diff --git a/doc/include/class_datamodel_spec2dobj.rst b/doc/include/class_datamodel_spec2dobj.rst index b68107e926..c0803b9f40 100644 --- a/doc/include/class_datamodel_spec2dobj.rst +++ b/doc/include/class_datamodel_spec2dobj.rst @@ -1,28 +1,28 @@ -**Version**: 1.1.1 +**Version**: 1.1.2 -====================== =================================================================================================== ================= ================================================================================================================================================================================ -Attribute Type Array Type Description -====================== =================================================================================================== ================= ================================================================================================================================================================================ -``bkg_redux_skymodel`` `numpy.ndarray`_ `numpy.floating`_ 2D sky model image without the background subtraction (float32) -``bpmmask`` :class:`~pypeit.images.imagebitmask.ImageBitMaskArray` 2D bad-pixel mask for the image -``det`` int Detector index -``detector`` :class:`~pypeit.images.detector_container.DetectorContainer`, :class:`~pypeit.images.mosaic.Mosaic` Detector or Mosaic metadata -``ivarmodel`` `numpy.ndarray`_ `numpy.floating`_ 2D ivar model image (float32) -``ivarraw`` `numpy.ndarray`_ `numpy.floating`_ 2D processed inverse variance image (float32) -``maskdef_designtab`` `astropy.table.table.Table`_ Table with slitmask design and object info -``med_chis`` `numpy.ndarray`_ `numpy.floating`_ Median of the chi image for each slit/order -``objmodel`` `numpy.ndarray`_ `numpy.floating`_ 2D object model image (float32) -``scaleimg`` `numpy.ndarray`_ `numpy.floating`_ 2D multiplicative scale image [or a single scalar as an array] that has been applied to the science image (float32) -``sci_spat_flexure`` float Shift, in spatial pixels, between this image and SlitTrace -``sci_spec_flexure`` `astropy.table.table.Table`_ Global shift of the spectrum to correct for spectralflexure (pixels). This is based on the sky spectrum atthe center of each slit -``sciimg`` `numpy.ndarray`_ `numpy.floating`_ 2D processed science image (float32) -``skymodel`` `numpy.ndarray`_ `numpy.floating`_ 2D sky model image (float32) -``slits`` :class:`~pypeit.slittrace.SlitTraceSet` SlitTraceSet defining the slits -``std_chis`` `numpy.ndarray`_ `numpy.floating`_ std of the chi image for each slit/order -``tilts`` `numpy.ndarray`_ `numpy.floating`_ 2D tilts image (float64) -``vel_corr`` float Relativistic velocity correction for wavelengths -``vel_type`` str Type of reference frame correction (if any). Options are listed in the routine: WavelengthSolutionPar.valid_reference_frames() Current list: observed, heliocentric, barycentric -``waveimg`` `numpy.ndarray`_ `numpy.floating`_ 2D wavelength image in vacuum (float64) -``wavesol`` `astropy.table.table.Table`_ Table with WaveCalib diagnostic info -====================== =================================================================================================== ================= ================================================================================================================================================================================ +====================== =================================================================================================== ================= ================================================================================================================================================================================================================================= +Attribute Type Array Type Description +====================== =================================================================================================== ================= ================================================================================================================================================================================================================================= +``bkg_redux_skymodel`` `numpy.ndarray`_ `numpy.floating`_ 2D sky model image without the background subtraction (float32) +``bpmmask`` :class:`~pypeit.images.imagebitmask.ImageBitMaskArray` 2D bad-pixel mask for the image +``det`` int Detector index +``detector`` :class:`~pypeit.images.detector_container.DetectorContainer`, :class:`~pypeit.images.mosaic.Mosaic` Detector or Mosaic metadata +``ivarmodel`` `numpy.ndarray`_ `numpy.floating`_ 2D ivar model image (float32) +``ivarraw`` `numpy.ndarray`_ `numpy.floating`_ 2D processed inverse variance image (float32) +``maskdef_designtab`` `astropy.table.table.Table`_ Table with slitmask design and object info +``med_chis`` `numpy.ndarray`_ `numpy.floating`_ Median of the chi image for each slit/order +``objmodel`` `numpy.ndarray`_ `numpy.floating`_ 2D object model image (float32) +``scaleimg`` `numpy.ndarray`_ `numpy.floating`_ 2D multiplicative scale image [or a single scalar as an array] that has been applied to the science image (float32) +``sci_spat_flexure`` `numpy.ndarray`_ `numpy.floating`_ Shift, in spatial pixels, between this image and SlitTrace. Shape is (nslits, 2), where spat_flexure[i,0] is the spatial shift of the left edge of slit i and spat_flexure[i,1] is the spatial shift of the right edge of slit i. +``sci_spec_flexure`` `astropy.table.table.Table`_ Global shift of the spectrum to correct for spectral flexure (pixels). This is based on the sky spectrum at the center of each slit +``sciimg`` `numpy.ndarray`_ `numpy.floating`_ 2D processed science image (float32) +``skymodel`` `numpy.ndarray`_ `numpy.floating`_ 2D sky model image (float32) +``slits`` :class:`~pypeit.slittrace.SlitTraceSet` SlitTraceSet defining the slits +``std_chis`` `numpy.ndarray`_ `numpy.floating`_ std of the chi image for each slit/order +``tilts`` `numpy.ndarray`_ `numpy.floating`_ 2D tilts image (float64) +``vel_corr`` float Relativistic velocity correction for wavelengths +``vel_type`` str Type of reference frame correction (if any). Options are listed in the routine: WavelengthSolutionPar.valid_reference_frames() Current list: observed, heliocentric, barycentric +``waveimg`` `numpy.ndarray`_ `numpy.floating`_ 2D wavelength image in vacuum (float64) +``wavesol`` `astropy.table.table.Table`_ Table with WaveCalib diagnostic info +====================== =================================================================================================== ================= ================================================================================================================================================================================================================================= diff --git a/doc/include/class_datamodel_wavetilts.rst b/doc/include/class_datamodel_wavetilts.rst index f071e3f17d..ada33ea306 100644 --- a/doc/include/class_datamodel_wavetilts.rst +++ b/doc/include/class_datamodel_wavetilts.rst @@ -1,19 +1,19 @@ -**Version**: 1.2.0 +**Version**: 1.2.1 -==================== ============================ ================= =============================================================================================================================================================== -Attribute Type Array Type Description -==================== ============================ ================= =============================================================================================================================================================== -``PYP_SPEC`` str PypeIt spectrograph name -``bpmtilts`` `numpy.ndarray`_ `numpy.integer`_ Bad pixel mask for tilt solutions. Keys are taken from SlitTraceSetBitmask -``coeffs`` `numpy.ndarray`_ `numpy.floating`_ 2D coefficents for the fit on the initial slits. One set per slit/order (3D array). -``func2d`` str Function used for the 2D fit -``nslit`` int Total number of slits. This can include masked slits -``slits_filename`` str Path to SlitTraceSet file. This helps to find the Slits calibration file when running pypeit_chk_tilts() -``spat_flexure`` float Flexure shift from the input TiltImage -``spat_id`` `numpy.ndarray`_ `numpy.integer`_ Slit spat_id -``spat_order`` `numpy.ndarray`_ `numpy.integer`_ Order for spatial fit (nslit) -``spec_order`` `numpy.ndarray`_ `numpy.integer`_ Order for spectral fit (nslit) -``tilt_traces`` `astropy.table.table.Table`_ Table with the positions of the traced and fitted tilts for all the slits. see :func:`~pypeit.wavetilts.BuildWaveTilts.make_tbl_tilt_traces` for more details. -``tiltimg_filename`` str Path to Tiltimg file. This helps to find Tiltimg file when running pypeit_chk_tilts() -==================== ============================ ================= =============================================================================================================================================================== +==================== ============================ ================= ================================================================================================================================================================================================================================================ +Attribute Type Array Type Description +==================== ============================ ================= ================================================================================================================================================================================================================================================ +``PYP_SPEC`` str PypeIt spectrograph name +``bpmtilts`` `numpy.ndarray`_ `numpy.integer`_ Bad pixel mask for tilt solutions. Keys are taken from SlitTraceSetBitmask +``coeffs`` `numpy.ndarray`_ `numpy.floating`_ 2D coefficents for the fit on the initial slits. One set per slit/order (3D array). +``func2d`` str Function used for the 2D fit +``nslit`` int Total number of slits. This can include masked slits +``slits_filename`` str Path to SlitTraceSet file. This helps to find the Slits calibration file when running pypeit_chk_tilts() +``spat_flexure`` `numpy.ndarray`_ `numpy.floating`_ Spatial flexure shift, in spatial pixels, between TiltImage and SlitTrace. Shape is (nslits, 2), where spat_flexure[i,0] is the spatial shift of the left edge of slit i and spat_flexure[i,1] is the spatial shift of the right edge of slit i. +``spat_id`` `numpy.ndarray`_ `numpy.integer`_ Slit spat_id +``spat_order`` `numpy.ndarray`_ `numpy.integer`_ Order for spatial fit (nslit) +``spec_order`` `numpy.ndarray`_ `numpy.integer`_ Order for spectral fit (nslit) +``tilt_traces`` `astropy.table.table.Table`_ Table with the positions of the traced and fitted tilts for all the slits. see :func:`~pypeit.wavetilts.BuildWaveTilts.make_tbl_tilt_traces` for more details. +``tiltimg_filename`` str Path to Tiltimg file. This helps to find Tiltimg file when running pypeit_chk_tilts() +==================== ============================ ================= ================================================================================================================================================================================================================================================ diff --git a/doc/include/datamodel_arcimage.rst b/doc/include/datamodel_arcimage.rst index 513d279776..221c8887a8 100644 --- a/doc/include/datamodel_arcimage.rst +++ b/doc/include/datamodel_arcimage.rst @@ -1,5 +1,5 @@ -Version 1.3.0 +Version 1.3.1 ================ ============================== ========= ================================================================================================================================================ HDU Name HDU Type Data Type Description diff --git a/doc/include/datamodel_biasimage.rst b/doc/include/datamodel_biasimage.rst index fcc54ad6c9..a9710f5d4b 100644 --- a/doc/include/datamodel_biasimage.rst +++ b/doc/include/datamodel_biasimage.rst @@ -1,5 +1,5 @@ -Version 1.3.0 +Version 1.3.1 ================= ============================== ========= ================================================================================================================================================ HDU Name HDU Type Data Type Description diff --git a/doc/include/datamodel_darkimage.rst b/doc/include/datamodel_darkimage.rst index 787919fc75..e56057a0fb 100644 --- a/doc/include/datamodel_darkimage.rst +++ b/doc/include/datamodel_darkimage.rst @@ -1,5 +1,5 @@ -Version 1.3.0 +Version 1.3.1 ================= ============================== ========= ================================================================================================================================================ HDU Name HDU Type Data Type Description diff --git a/doc/include/datamodel_spec2dobj.rst b/doc/include/datamodel_spec2dobj.rst index 01417f2214..02ce84d137 100644 --- a/doc/include/datamodel_spec2dobj.rst +++ b/doc/include/datamodel_spec2dobj.rst @@ -1,29 +1,29 @@ -Version: 1.1.1 +Version: 1.1.2 -====================== ========================= ========== ================================================================================================================================================================================ -Obj Key Obj Type Array Type Description -====================== ========================= ========== ================================================================================================================================================================================ -``bkg_redux_skymodel`` ndarray floating 2D sky model image without the background subtraction (float32) -``bpmmask`` ImageBitMaskArray 2D bad-pixel mask for the image -``det`` int Detector index -``detector`` DetectorContainer, Mosaic Detector or Mosaic metadata -``ivarmodel`` ndarray floating 2D ivar model image (float32) -``ivarraw`` ndarray floating 2D processed inverse variance image (float32) -``maskdef_designtab`` Table Table with slitmask design and object info -``med_chis`` ndarray floating Median of the chi image for each slit/order -``objmodel`` ndarray floating 2D object model image (float32) -``scaleimg`` ndarray floating 2D multiplicative scale image [or a single scalar as an array] that has been applied to the science image (float32) -``sci_spat_flexure`` float Shift, in spatial pixels, between this image and SlitTrace -``sci_spec_flexure`` Table Global shift of the spectrum to correct for spectralflexure (pixels). This is based on the sky spectrum atthe center of each slit -``sciimg`` ndarray floating 2D processed science image (float32) -``skymodel`` ndarray floating 2D sky model image (float32) -``slits`` SlitTraceSet SlitTraceSet defining the slits -``std_chis`` ndarray floating std of the chi image for each slit/order -``tilts`` ndarray floating 2D tilts image (float64) -``vel_corr`` float Relativistic velocity correction for wavelengths -``vel_type`` str Type of reference frame correction (if any). Options are listed in the routine: WavelengthSolutionPar.valid_reference_frames() Current list: observed, heliocentric, barycentric -``waveimg`` ndarray floating 2D wavelength image in vacuum (float64) -``wavesol`` Table Table with WaveCalib diagnostic info -====================== ========================= ========== ================================================================================================================================================================================ +====================== ========================= ========== ================================================================================================================================================================================================================================= +Obj Key Obj Type Array Type Description +====================== ========================= ========== ================================================================================================================================================================================================================================= +``bkg_redux_skymodel`` ndarray floating 2D sky model image without the background subtraction (float32) +``bpmmask`` ImageBitMaskArray 2D bad-pixel mask for the image +``det`` int Detector index +``detector`` DetectorContainer, Mosaic Detector or Mosaic metadata +``ivarmodel`` ndarray floating 2D ivar model image (float32) +``ivarraw`` ndarray floating 2D processed inverse variance image (float32) +``maskdef_designtab`` Table Table with slitmask design and object info +``med_chis`` ndarray floating Median of the chi image for each slit/order +``objmodel`` ndarray floating 2D object model image (float32) +``scaleimg`` ndarray floating 2D multiplicative scale image [or a single scalar as an array] that has been applied to the science image (float32) +``sci_spat_flexure`` ndarray floating Shift, in spatial pixels, between this image and SlitTrace. Shape is (nslits, 2), where spat_flexure[i,0] is the spatial shift of the left edge of slit i and spat_flexure[i,1] is the spatial shift of the right edge of slit i. +``sci_spec_flexure`` Table Global shift of the spectrum to correct for spectral flexure (pixels). This is based on the sky spectrum at the center of each slit +``sciimg`` ndarray floating 2D processed science image (float32) +``skymodel`` ndarray floating 2D sky model image (float32) +``slits`` SlitTraceSet SlitTraceSet defining the slits +``std_chis`` ndarray floating std of the chi image for each slit/order +``tilts`` ndarray floating 2D tilts image (float64) +``vel_corr`` float Relativistic velocity correction for wavelengths +``vel_type`` str Type of reference frame correction (if any). Options are listed in the routine: WavelengthSolutionPar.valid_reference_frames() Current list: observed, heliocentric, barycentric +``waveimg`` ndarray floating 2D wavelength image in vacuum (float64) +``wavesol`` Table Table with WaveCalib diagnostic info +====================== ========================= ========== ================================================================================================================================================================================================================================= diff --git a/doc/include/datamodel_tiltimage.rst b/doc/include/datamodel_tiltimage.rst index 55189a8428..abf90e427f 100644 --- a/doc/include/datamodel_tiltimage.rst +++ b/doc/include/datamodel_tiltimage.rst @@ -1,5 +1,5 @@ -Version 1.3.0 +Version 1.3.1 ================= ============================== ========= ================================================================================================================================================ HDU Name HDU Type Data Type Description diff --git a/doc/include/datamodel_wavetilts.rst b/doc/include/datamodel_wavetilts.rst index 6e9f122e53..4a95b30aa0 100644 --- a/doc/include/datamodel_wavetilts.rst +++ b/doc/include/datamodel_wavetilts.rst @@ -1,5 +1,5 @@ -Version 1.2.0 +Version 1.2.1 =========== ============================== ========= ===================================================== HDU Name HDU Type Data Type Description diff --git a/doc/include/imgproc_defaults_table.rst b/doc/include/imgproc_defaults_table.rst index 2fc2ffd84d..bbeb19e0f6 100644 --- a/doc/include/imgproc_defaults_table.rst +++ b/doc/include/imgproc_defaults_table.rst @@ -1,19 +1,19 @@ -======================== ========= ========= ========= ========= ========= ========= ========= ============= ============= ======== ============ =========== -Parameter Default ``bias`` ``dark`` ``trace`` ``arc`` ``tilt`` ``align`` ``pixelflat`` ``illumflat`` ``sky`` ``standard`` ``science`` -======================== ========= ========= ========= ========= ========= ========= ========= ============= ============= ======== ============ =========== -``apply_gain`` ``True`` -``use_pattern`` ``False`` -``empirical_rn`` ``False`` -``use_overscan`` ``True`` -``trim`` ``True`` -``orient`` ``True`` -``use_biasimage`` ``True`` ``False`` -``use_darkimage`` ``False`` -``spat_flexure_correct`` ``False`` -``use_pixelflat`` ``True`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` -``use_illumflat`` ``True`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` -``use_specillum`` ``False`` -``shot_noise`` ``True`` ``False`` -``noise_floor`` ``0.0`` ``0.01`` ``0.01`` ``0.01`` -``mask_cr`` ``False`` ``True`` ``True`` ``True`` ``True`` -======================== ========= ========= ========= ========= ========= ========= ========= ============= ============= ======== ============ =========== +======================= ========= ========= ========= ========= ========= ========= ========= ============= ============= ======== ============ =========== +Parameter Default ``bias`` ``dark`` ``trace`` ``arc`` ``tilt`` ``align`` ``pixelflat`` ``illumflat`` ``sky`` ``standard`` ``science`` +======================= ========= ========= ========= ========= ========= ========= ========= ============= ============= ======== ============ =========== +``apply_gain`` ``True`` +``use_pattern`` ``False`` +``empirical_rn`` ``False`` +``use_overscan`` ``True`` +``trim`` ``True`` +``orient`` ``True`` +``use_biasimage`` ``True`` ``False`` +``use_darkimage`` ``False`` +``spat_flexure_method`` ``skip`` +``use_pixelflat`` ``True`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` +``use_illumflat`` ``True`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` ``False`` +``use_specillum`` ``False`` +``shot_noise`` ``True`` ``False`` +``noise_floor`` ``0.0`` ``0.01`` ``0.01`` ``0.01`` +``mask_cr`` ``False`` ``True`` ``True`` ``True`` ``True`` +======================= ========= ========= ========= ========= ========= ========= ========= ============= ============= ======== ============ =========== diff --git a/doc/scripts/build_imgproc_defaults.py b/doc/scripts/build_imgproc_defaults.py index 5ab91bf955..e52f349e31 100644 --- a/doc/scripts/build_imgproc_defaults.py +++ b/doc/scripts/build_imgproc_defaults.py @@ -25,7 +25,7 @@ def write_imgproc_def_table(ofile, spec=None): 'orient', 'use_biasimage', 'use_darkimage', - 'spat_flexure_correct', + 'spat_flexure_method', 'use_pixelflat', 'use_illumflat', 'use_specillum', diff --git a/pypeit/calibrations.py b/pypeit/calibrations.py index 6154ac7a15..f54d913a37 100644 --- a/pypeit/calibrations.py +++ b/pypeit/calibrations.py @@ -258,6 +258,9 @@ def find_calibrations(self, frametype, frameclass): raw_files : :obj:`list` The list of raw files in :attr:`fitstbl` with the provided frametype. + rows : :obj:`numpy.ndarray` + The indices of the raw files in :attr:`fitstbl` with the provided + frametype. cal_file : `Path`_ The path with/for the processed calibration frame calib_key : :obj:`str` @@ -285,12 +288,12 @@ def find_calibrations(self, frametype, frameclass): setup = self.fitstbl['setup'][self.frame] cal_file = frameclass.glob(self.calib_dir, setup, self.calib_ID, detname=detname) if cal_file is None or len(cal_file) > 1: - return [], None, None, setup, None, detname + return [], rows, None, None, setup, None, detname cal_file = cal_file[0] calib_key = frameclass.parse_key_dir(str(cal_file), from_filename=True)[0] calib_id = frameclass.parse_calib_key(calib_key)[1] - return [], cal_file, calib_key, setup, frameclass.ingest_calib_id(calib_id), detname + return [], rows, cal_file, calib_key, setup, frameclass.ingest_calib_id(calib_id), detname # Otherwise, use the metadata for the raw frames to set the name of # the processed calibration frame. @@ -300,7 +303,7 @@ def find_calibrations(self, frametype, frameclass): # Construct the expected calibration frame file name cal_file = Path(frameclass.construct_file_name(calib_key, calib_dir=self.calib_dir)) - return self.fitstbl.frame_paths(rows), cal_file, calib_key, setup, \ + return self.fitstbl.frame_paths(rows), rows, cal_file, calib_key, setup, \ frameclass.ingest_calib_id(calib_id), detname def set_config(self, frame, det, par=None): @@ -354,7 +357,7 @@ def get_arc(self, force:str=None): # Find the calibrations frame = {'type': 'arc', 'class': buildimage.ArcImage} - raw_files, cal_file, calib_key, setup, calib_id, detname \ + raw_files, raw_index, cal_file, calib_key, setup, calib_id, detname \ = self.find_calibrations(frame['type'], frame['class']) if len(raw_files) == 0 and cal_file is None: @@ -375,13 +378,17 @@ def get_arc(self, force:str=None): # Perform a check on the files self.check_calibrations(raw_files) + # If a manual spatial flexure is requested for any of the frames, collect it now + manual_spat_flexure = self.fitstbl.get_shifts(raw_index) + # Otherwise, create the processed file. msgs.info(f'Preparing a {frame["class"].calib_type} calibration frame.') self.msarc = buildimage.buildimage_fromlist(self.spectrograph, self.det, self.par['arcframe'], raw_files, bias=self.msbias, bpm=self.msbpm, - dark=self.msdark, calib_dir=self.calib_dir, - setup=setup, calib_id=calib_id) + dark=self.msdark, slits=self.slits, calib_dir=self.calib_dir, + setup=setup, calib_id=calib_id, + manual_spat_flexure=manual_spat_flexure) # Save the result self.msarc.to_file() # Return it @@ -406,7 +413,7 @@ def get_tiltimg(self, force:str=None): # Find the calibrations frame = {'type': 'tilt', 'class':buildimage.TiltImage} - raw_files, cal_file, calib_key, setup, calib_id, detname \ + raw_files, raw_index, cal_file, calib_key, setup, calib_id, detname \ = self.find_calibrations(frame['type'], frame['class']) if len(raw_files) == 0 and cal_file is None: @@ -427,6 +434,9 @@ def get_tiltimg(self, force:str=None): # Perform a check on the files self.check_calibrations(raw_files) + # If a manual spatial flexure is requested for any of the frames, collect it now + manual_spat_flexure = self.fitstbl.get_shifts(raw_index) + # Otherwise, create the processed file. msgs.info(f'Preparing a {frame["class"].calib_type} calibration frame.') self.mstilt = buildimage.buildimage_fromlist(self.spectrograph, self.det, @@ -434,7 +444,8 @@ def get_tiltimg(self, force:str=None): bias=self.msbias, bpm=self.msbpm, dark=self.msdark, slits=self.slits, calib_dir=self.calib_dir, setup=setup, - calib_id=calib_id) + calib_id=calib_id, + manual_spat_flexure=manual_spat_flexure) # Save the result self.mstilt.to_file() # Return it @@ -463,7 +474,7 @@ def get_align(self, force:str=None): # Find the calibrations frame = {'type': 'align', 'class': alignframe.Alignments} - raw_files, cal_file, calib_key, setup, calib_id, detname \ + raw_files, raw_index, cal_file, calib_key, setup, calib_id, detname \ = self.find_calibrations(frame['type'], frame['class']) if len(raw_files) == 0 and cal_file is None: @@ -487,13 +498,17 @@ def get_align(self, force:str=None): # Perform a check on the files self.check_calibrations(raw_files) + # If a manual spatial flexure is requested for any of the frames, collect it now + manual_spat_flexure = self.fitstbl.get_shifts(raw_index) + # Otherwise, create the processed file. msgs.info(f'Preparing a {frame["class"].calib_type} calibration frame.') msalign = buildimage.buildimage_fromlist(self.spectrograph, self.det, self.par['alignframe'], raw_files, bias=self.msbias, bpm=self.msbpm, dark=self.msdark, calib_dir=self.calib_dir, - setup=setup, calib_id=calib_id) + setup=setup, calib_id=calib_id, + manual_spat_flexure=manual_spat_flexure) # Instantiate # TODO: From JFH: Do we need the bpm here? Check that this was in the previous code. @@ -525,7 +540,7 @@ def get_bias(self, force:str=None): # Find the calibrations frame = {'type': 'bias', 'class': buildimage.BiasImage} - raw_files, cal_file, calib_key, setup, calib_id, detname \ + raw_files, raw_index, cal_file, calib_key, setup, calib_id, detname \ = self.find_calibrations(frame['type'], frame['class']) # If no raw files are available and no processed calibration frame @@ -574,7 +589,7 @@ def get_dark(self, force:str=None): # Find the calibrations frame = {'type': 'dark', 'class': buildimage.DarkImage} - raw_files, cal_file, calib_key, setup, calib_id, detname \ + raw_files, raw_index, cal_file, calib_key, setup, calib_id, detname \ = self.find_calibrations(frame['type'], frame['class']) if len(raw_files) == 0 and cal_file is None: @@ -665,7 +680,7 @@ def get_scattlight(self, force:str=None): # Prep frame = {'type': 'scattlight', 'class': scattlight.ScatteredLight} - raw_scattlight_files, cal_file, calib_key, setup, calib_id, detname = \ + raw_scattlight_files, raw_scattlight_index, cal_file, calib_key, setup, calib_id, detname = \ self.find_calibrations(frame['type'], frame['class']) scatt_idx = self.fitstbl.find_frames(frame['type'], calib_ID=self.calib_ID, index=True) @@ -692,17 +707,21 @@ def get_scattlight(self, force:str=None): # Perform a check on the files self.check_calibrations(raw_scattlight_files) + # If a manual spatial flexure is requested for any of the frames, collect it now + manual_spat_flexure = self.fitstbl.get_shifts(raw_scattlight_index) + binning = self.fitstbl[scatt_idx[0]]['binning'] dispname = self.fitstbl[scatt_idx[0]]['dispname'] scattlightImage = buildimage.buildimage_fromlist(self.spectrograph, self.det, self.par['scattlightframe'], raw_scattlight_files, bias=self.msbias, bpm=self.msbpm, dark=self.msdark, calib_dir=self.calib_dir, - setup=setup, calib_id=calib_id) + setup=setup, calib_id=calib_id, + manual_spat_flexure=manual_spat_flexure) spatbin = parse.parse_binning(binning)[1] pad = self.par['scattlight_pad'] // spatbin - offslitmask = self.slits.slit_img(pad=pad, flexure=None) == -1 + offslitmask = self.slits.slit_img(pad=pad, spat_flexure=None) == -1 # Get starting parameters for the scattered light model x0, bounds = self.spectrograph.scattered_light_archive(binning, dispname) @@ -786,20 +805,21 @@ def get_flats(self, force:str=None): # get illumination flat frames illum_frame = {'type': 'illumflat', 'class': flatfield.FlatImages} - raw_illum_files, illum_cal_file, illum_calib_key, illum_setup, illum_calib_id, detname \ + raw_illum_files, raw_illum_index, illum_cal_file, illum_calib_key, illum_setup, illum_calib_id, detname \ = self.find_calibrations(illum_frame['type'], illum_frame['class']) # get pixel flat frames pixel_frame = {'type': 'pixelflat', 'class': flatfield.FlatImages} - raw_pixel_files, pixel_cal_file, pixel_calib_key, pixel_setup, pixel_calib_id, detname \ - = [], None, None, illum_setup, None, detname + raw_pixel_files, raw_pixel_index, pixel_cal_file, pixel_calib_key, pixel_setup, pixel_calib_id, detname \ + = [], [], None, None, illum_setup, None, detname # read in the raw pixelflat frames only if the user has not provided a pixelflat_file if self.par['flatfield']['pixelflat_file'] is None: - raw_pixel_files, pixel_cal_file, pixel_calib_key, pixel_setup, pixel_calib_id, detname \ + raw_pixel_files, raw_pixel_index, pixel_cal_file, pixel_calib_key, pixel_setup, pixel_calib_id, detname \ = self.find_calibrations(pixel_frame['type'], pixel_frame['class']) # get lamp off flat frames raw_lampoff_files = self.fitstbl.find_frame_files('lampoffflats', calib_ID=self.calib_ID) + raw_lampoff_index = self.fitstbl.find_frames('lampoffflats', calib_ID=self.calib_ID, index=True) # Check if we have any calibration frames to work with if len(raw_pixel_files) == 0 and pixel_cal_file is None \ @@ -859,6 +879,9 @@ def get_flats(self, force:str=None): # Perform a check on the files self.check_calibrations(raw_pixel_files) + # If a manual spatial flexure is requested for any of the frames, collect it now + manual_spat_flexure = self.fitstbl.get_shifts(raw_pixel_index) + msgs.info('Creating pixel-flat calibration frame using files: ') for f in raw_pixel_files: msgs.prindent(f'{Path(f).name}') @@ -867,7 +890,8 @@ def get_flats(self, force:str=None): raw_pixel_files, dark=self.msdark, slits=self.slits, bias=self.msbias, bpm=self.msbpm, - scattlight=self.msscattlight) + scattlight=self.msscattlight, + manual_spat_flexure=manual_spat_flexure) if len(raw_lampoff_files) > 0: # Reset the BPM self.get_bpm(frame=raw_lampoff_files[0]) @@ -875,6 +899,9 @@ def get_flats(self, force:str=None): # Perform a check on the files self.check_calibrations(raw_lampoff_files) + # If a manual spatial flexure is requested for any of the frames, collect it now + manual_spat_flexure = self.fitstbl.get_shifts(raw_lampoff_index) + msgs.info('Subtracting lamp off flats using files: ') for f in raw_lampoff_files: msgs.prindent(f'{Path(f).name}') @@ -883,7 +910,8 @@ def get_flats(self, force:str=None): raw_lampoff_files, slits=self.slits, dark=self.msdark, bias=self.msbias, - bpm=self.msbpm, scattlight=self.msscattlight) + bpm=self.msbpm, scattlight=self.msscattlight, + manual_spat_flexure=manual_spat_flexure) pixel_flat = pixel_flat.sub(lampoff_flat) # Initialise the pixel flat @@ -905,6 +933,9 @@ def get_flats(self, force:str=None): # Perform a check on the files self.check_calibrations(raw_illum_files) + # If a manual spatial flexure is requested for any of the frames, collect it now + manual_spat_flexure = self.fitstbl.get_shifts(raw_illum_index) + msgs.info('Creating slit-illumination flat calibration frame using files: ') for f in raw_illum_files: msgs.prindent(f'{Path(f).name}') @@ -912,7 +943,8 @@ def get_flats(self, force:str=None): illum_flat = buildimage.buildimage_fromlist(self.spectrograph, self.det, self.par['illumflatframe'], raw_illum_files, dark=self.msdark, bias=self.msbias, scattlight=self.msscattlight, - slits=self.slits, flatimages=self.flatimages, bpm=self.msbpm) + slits=self.slits, flatimages=self.flatimages, bpm=self.msbpm, + manual_spat_flexure=manual_spat_flexure) if len(raw_lampoff_files) > 0: msgs.info('Subtracting lamp off flats using files: ') for f in raw_lampoff_files: @@ -921,6 +953,9 @@ def get_flats(self, force:str=None): # Perform a check on the files self.check_calibrations(raw_lampoff_files) + # If a manual spatial flexure is requested for any of the frames, collect it now + manual_spat_flexure = self.fitstbl.get_shifts(raw_lampoff_index) + # Build the image lampoff_flat = buildimage.buildimage_fromlist(self.spectrograph, self.det, self.par['lampoffflatsframe'], @@ -929,7 +964,8 @@ def get_flats(self, force:str=None): bias=self.msbias, slits=self.slits, scattlight=self.msscattlight, - bpm=self.msbpm) + bpm=self.msbpm, + manual_spat_flexure=manual_spat_flexure) illum_flat = illum_flat.sub(lampoff_flat) # Initialise the illum flat @@ -980,6 +1016,12 @@ def get_slits(self, force:str=None): """ Load or generate the definition of the slit boundaries. + Args: + force (:obj:`str`, optional): + 'remake' -- Force the frame to be remade. + 'reload' -- Reload the frame if it exists. + None -- Load the existing frame if it exists and reuse_calibs=True + Returns: :class:`~pypeit.slittrace.SlitTraceSet`: Traces of the slit edges; also kept internally as :attr:`slits`. @@ -993,9 +1035,10 @@ def get_slits(self, force:str=None): # Prep frame = {'type': 'trace', 'class': slittrace.SlitTraceSet} - raw_trace_files, cal_file, calib_key, setup, calib_id, detname \ - = self.find_calibrations(frame['type'], frame['class']) + raw_trace_files, raw_trace_index, cal_file, calib_key, setup, calib_id, detname \ + = self.find_calibrations(frame['type'], frame['class']) raw_lampoff_files = self.fitstbl.find_frame_files('lampoffflats', calib_ID=self.calib_ID) + raw_lampoff_index = self.fitstbl.find_frames('lampoffflats', calib_ID=self.calib_ID, index=True) if len(raw_trace_files) == 0 and cal_file is None: msgs.warn(f'No raw {frame["type"]} frames found and unable to identify a relevant ' @@ -1040,6 +1083,9 @@ def get_slits(self, force:str=None): # Perform a check on the files self.check_calibrations(raw_trace_files) + # If a manual spatial flexure is requested for any of the frames, collect it now + manual_spat_flexure = self.fitstbl.get_shifts(raw_trace_index) + # NOTE: self.msscattlight is *always* created after identifying the # slits, meaning that it is redundant to pass the scattlight argument # here. @@ -1047,7 +1093,8 @@ def get_slits(self, force:str=None): self.par['traceframe'], raw_trace_files, bias=self.msbias, bpm=self.msbpm, dark=self.msdark, calib_dir=self.calib_dir, - setup=setup, calib_id=calib_id) + setup=setup, calib_id=calib_id, + manual_spat_flexure=manual_spat_flexure) if len(raw_lampoff_files) > 0: msgs.info('Subtracting lamp off flats using files: ') for f in raw_lampoff_files: @@ -1059,10 +1106,14 @@ def get_slits(self, force:str=None): # Perform a check on the files self.check_calibrations(raw_lampoff_files) + # If a manual spatial flexure is requested for any of the frames, collect it now + manual_spat_flexure = self.fitstbl.get_shifts(raw_lampoff_index) + lampoff_flat = buildimage.buildimage_fromlist(self.spectrograph, self.det, self.par['lampoffflatsframe'], raw_lampoff_files, dark=self.msdark, - bias=self.msbias, bpm=self.msbpm) + bias=self.msbias, bpm=self.msbpm, + manual_spat_flexure=manual_spat_flexure) traceImage = traceImage.sub(lampoff_flat) edges = edgetrace.EdgeTraceSet(traceImage, self.spectrograph, self.par['slitedges'], @@ -1126,7 +1177,7 @@ def get_wv_calib(self, force:str=None): # Find the calibrations frame = {'type': 'arc', 'class': wavecalib.WaveCalib} - raw_files, cal_file, calib_key, setup, calib_id, detname \ + raw_files, raw_index, cal_file, calib_key, setup, calib_id, detname \ = self.find_calibrations(frame['type'], frame['class']) if len(raw_files) == 0 and cal_file is None: @@ -1212,7 +1263,7 @@ def get_tilts(self, force:str=None): # Find the calibrations frame = {'type': 'tilt', 'class': wavetilts.WaveTilts} - raw_files, cal_file, calib_key, setup, calib_id, detname \ + raw_files, raw_index, cal_file, calib_key, setup, calib_id, detname \ = self.find_calibrations(frame['type'], frame['class']) if len(raw_files) == 0 and cal_file is None: @@ -1232,8 +1283,8 @@ def get_tilts(self, force:str=None): return self.wavetilts # Get flexure - _spat_flexure = self.mstilt.spat_flexure \ - if self.par['tiltframe']['process']['spat_flexure_correct'] else None + _spat_flexure = self.mstilt.spat_flexure if self.par['tiltframe']['process']['spat_flexure_method'] != "skip" \ + else None # get measured fwhm from wv_calib measured_fwhms = [wvfit.fwhm for wvfit in self.wv_calib.wv_fits] @@ -1267,7 +1318,7 @@ def process_load_selection(self, frame, cal_file, force): force : :obj:`str` Defines how to treat a pre-existing calibration file. Must be one of the following options: - + - ``'remake'``: Force the calibration be remade. - ``'reload'``: Reload the frame if it exists. @@ -1290,7 +1341,7 @@ def process_load_selection(self, frame, cal_file, force): f"{frame['class'].__name__} calibration.") self.success = False return None - if force == 'reload' or (self.reuse_calibs and _cal_file.exists()): + if force == 'reload' or (self.reuse_calibs and _cal_file.exists()): return frame['class'].from_file(_cal_file, chk_version=self.chk_version) def run_the_steps(self, stop_at_step:str=None): @@ -1301,7 +1352,7 @@ def run_the_steps(self, stop_at_step:str=None): for step in self.steps: if stop_at_step is not None and step == stop_at_step: force = 'remake' - msgs.info(f"Calibrations will stop at {stop_at_step}") + msgs.info(f"Calibrations will stop at {stop_at_step}") else: force = None getattr(self, f'get_{step}')(force=force) @@ -1309,7 +1360,7 @@ def run_the_steps(self, stop_at_step:str=None): self.failed_step = f'get_{step}' return if stop_at_step is not None and step == stop_at_step: - msgs.info(f"Calibrations stopping at {stop_at_step}") + msgs.info(f"Calibrations stopping at {stop_at_step}") return msgs.info("Calibration complete and/or fully loaded!") msgs.info("#######################################################################") @@ -1617,7 +1668,7 @@ def default_steps(): # Order matters! And the name must match a viable "get_{step}" method # in Calibrations. # TODO: Does the bpm need to be done after the dark? - return ['bias', 'dark', 'bpm', 'slits', 'arc', 'tiltimg', + return ['bias', 'dark', 'bpm', 'slits', 'arc', 'tiltimg', 'wv_calib', 'tilts', 'scattlight', 'flats'] diff --git a/pypeit/coadd2d.py b/pypeit/coadd2d.py index c82f8bf236..eb9f713fc5 100644 --- a/pypeit/coadd2d.py +++ b/pypeit/coadd2d.py @@ -26,7 +26,7 @@ from pypeit.core import findobj_skymask from pypeit.core.wavecal import wvutils from pypeit.core import coadd -from pypeit.core import parse +from pypeit.core import parse from pypeit import spec2dobj from pypeit.core.moment import moment1d from pypeit.manual_extract import ManualExtractionObj @@ -66,13 +66,13 @@ def get_instance(cls, spec2dfiles, spectrograph, par, det=1, return next(c for c in cls.__subclasses__() if c.__name__ == (spectrograph.pypeline + 'CoAdd2D'))( - spec2dfiles, spectrograph, par, det=det, + spec2dfiles, spectrograph, par, det=det, only_slits=only_slits, exclude_slits=exclude_slits, sn_smooth_npix=sn_smooth_npix, bkg_redux=bkg_redux, find_negative=find_negative, show=show, show_peaks=show_peaks, debug_offsets=debug_offsets, debug=debug) - def __init__(self, spec2d, spectrograph, par, det=1, - only_slits=None, exclude_slits=None, + def __init__(self, spec2d, spectrograph, par, det=1, + only_slits=None, exclude_slits=None, sn_smooth_npix=None, bkg_redux=False, find_negative=False, show=False, show_peaks=False, debug_offsets=False, debug=False): """ @@ -212,11 +212,11 @@ def __init__(self, spec2d, spectrograph, par, det=1, # effective exposure time self.exptime_coadd = self.stack_dict['exptime_coadd'] - + # define the wavelength grid for the 2d coadd self.wave_grid, self.wave_grid_mid, self.dsamp = self.get_wave_grid() - - + + # Handle the reference object self.handle_reference_obj() @@ -800,7 +800,7 @@ def reduce(self, pseudo_dict, show=False, clear_ginga=True, show_peaks=False, sh resampled_pixscale = parse.parse_binning(sciImage.detector.binning)[1]*sciImage.detector.platescale*self.par['coadd2d']['spat_samp_fact'] # Assign slitmask design information to detected objects - slits.assign_maskinfo(sobjs_obj, resampled_pixscale, None, TOLER=parcopy['reduce']['slitmask']['obj_toler']) + slits.assign_maskinfo(sobjs_obj, resampled_pixscale, None, tolerance=parcopy['reduce']['slitmask']['obj_toler']) if parcopy['reduce']['slitmask']['extract_missing_objs'] is True: # Set the FWHM for the extraction of missing objects @@ -998,7 +998,7 @@ def load_coadd2d_stacks(self, spec2d, chk_version=False): skymodel_stack.append(s2dobj.skymodel) sciivar_stack.append(s2dobj.ivarmodel) mask_stack.append(s2dobj.bpmmask.mask) - slitmask_stack.append(s2dobj.slits.slit_img(flexure=s2dobj.sci_spat_flexure)) + slitmask_stack.append(s2dobj.slits.slit_img(spat_flexure=s2dobj.sci_spat_flexure)) # check if exptime is consistent for all images exptime_coadd = np.percentile(exptime_stack, 50., method='higher') @@ -1121,7 +1121,7 @@ def compute_weights(self): # TODO maybe better behavior here would be to crash out to force the user to change the weight method explicitly # to 'uniform'. What I don't like here is that we are using uniform weights even though the user requested 'auto' # and they might miss the warning. Its debatable though. - + # warn if the user had put `auto` in the parset msgs.warn('Weights cannot be computed because no unique reference object ' 'with the highest S/N was found. Using uniform weights instead.') @@ -1206,12 +1206,12 @@ def get_brightest_obj(self, specobjs_list, spat_ids): """ msgs.error('The get_brightest_obj() method should be overloaded by the child class.') - + def handle_reference_obj(self): """ Dummy method to handle the reference object. Overloaded by child methods. """ - + msgs.error('The handle_reference_obj() method should be overloaded by the child class.') @@ -1263,7 +1263,7 @@ def get_maskdef_dict(self, slit_idx, ref_trace_stack): def wave_method(self): """ Get the wavelength method to be used in the coadd2d. This is a dummy method since - it is overloaded by child classes. + it is overloaded by child classes. Returns: str: The wavelength method to be used in the coadd2d. @@ -1290,7 +1290,7 @@ class MultiSlitCoAdd2D(CoAdd2D): """ - def __init__(self, spec2d_files, spectrograph, par, det=1, + def __init__(self, spec2d_files, spectrograph, par, det=1, only_slits=None, exclude_slits=None, sn_smooth_npix=None, bkg_redux=False, find_negative=False, show=False, show_peaks=False, debug_offsets=False, debug=False): @@ -1304,7 +1304,7 @@ def __init__(self, spec2d_files, spectrograph, par, det=1, # This will be an array of the object spatial pixel positions used for auto weights in each exposure self.spat_pixpos_id_weights = None - super().__init__(spec2d_files, spectrograph, det=det, + super().__init__(spec2d_files, spectrograph, det=det, only_slits=only_slits, exclude_slits=exclude_slits, sn_smooth_npix=sn_smooth_npix, bkg_redux=bkg_redux, find_negative=find_negative, par=par, show=show, show_peaks=show_peaks, debug_offsets=debug_offsets, @@ -1313,14 +1313,14 @@ def __init__(self, spec2d_files, spectrograph, par, det=1, def handle_reference_obj(self): """ - Method to handle the syntrax for the reference object to be used for offsets and weights. - - """ + Method to handle the syntrax for the reference object to be used for offsets and weights. + + """ # otherwise, find if there is a bright object we could use if len(self.stack_dict['specobjs_list']) > 0 and (self.par['coadd2d']['offsets'] == 'auto' or self.par['coadd2d']['weights'] == 'auto'): - # If the user passed in user_obj_ids, we will use these for the brighest object to - # be optionally used for offsets and weights. + # If the user passed in user_obj_ids, we will use these for the brighest object to + # be optionally used for offsets and weights. if self.par['coadd2d']['user_obj_ids'] is not None: if self.par['coadd2d']['weights'] != 'auto': msgs.error('Parameter `user_obj_ids` can only be used if weights are set to `auto`.') @@ -1350,14 +1350,14 @@ def handle_reference_obj(self): msgs.error('Not all spatial IDs are within spat_toler of each other') self.spatid_bri = int(np.rint(np.mean(spatids))) self.spat_pixpos_bri = np.array(spat_pixpos) - self.snr_bar_bri, _ = coadd.calc_snr(fluxes, ivars, gpms) - self.obj_id_bri = self.spat_pixpos_id_bri + self.snr_bar_bri, _ = coadd.calc_snr(fluxes, ivars, gpms) + self.obj_id_bri = self.spat_pixpos_id_bri else: - # Otherwise, find the brightest object in the stack and obtain the relevant information + # Otherwise, find the brightest object in the stack and obtain the relevant information self.spat_pixpos_id_bri, self.spat_pixpos_bri, self.spatid_bri, self.snr_bar_bri = self.get_brightest_obj(self.stack_dict['specobjs_list'], self.spat_ids) - self.obj_id_bri = self.spat_pixpos_id_bri + self.obj_id_bri = self.spat_pixpos_id_bri + - # TODO When we run multislit, we actually compute the rebinned images twice. Once here to compute the offsets @@ -1421,8 +1421,8 @@ def compute_offsets(self): trim_edg=self.par['reduce']['findobj']['find_trim_edge'], maxdev=self.par['reduce']['findobj']['find_maxdev'], numiterfit=self.par['reduce']['findobj']['find_numiterfit'], - ncoeff=3, snr_thresh=self.par['reduce']['findobj']['snr_thresh'], - nperslit=1 if self.par['coadd2d']['user_obj_ids'] is None else None, + ncoeff=3, snr_thresh=self.par['reduce']['findobj']['snr_thresh'], + nperslit=1 if self.par['coadd2d']['user_obj_ids'] is None else None, find_min_max=self.par['reduce']['findobj']['find_min_max'], show_trace=self.debug_offsets, show_peaks=self.debug_offsets) if len(sobjs_exp) == 0: @@ -1452,7 +1452,7 @@ def compute_offsets(self): f'to the trace for the user object {self.par["coadd2d"]["user_obj_ids"][iexp]} ' f'in exposure {iexp+1} within the specified spatial ' f'tolerance ={self.par["coadd2d"]["spat_toler"]}') - else: + else: traces_rect[:, iexp] = sobjs_exp.TRACE_SPAT if self.par['coadd2d']['user_obj_ids'] is not None: @@ -1508,7 +1508,7 @@ def compute_weights(self): def get_brightest_obj(self, specobjs_list, slit_spat_ids): """ - Utility routine to find the brightest reference object in each exposure given a specobjs_list + Utility routine to find the brightest reference object in each exposure given a specobjs_list for MultiSlit reductions. Args: @@ -1524,7 +1524,7 @@ def get_brightest_obj(self, specobjs_list, slit_spat_ids): - spat_pixpos: ndarray, float, shape=(len(specobjs_list),): Array of spatial pixel positions of the brightest reference object in each exposure - - spat_id (int): + - spat_id (int): The SPAT_ID for the slit that the highest S/N ratio object is on - snr_bar: ndarray, float, shape (len(list),): RMS S/N computed for this brightest reference object in each exposure """ @@ -1589,8 +1589,8 @@ def get_brightest_obj(self, specobjs_list, slit_spat_ids): def snr_report(self, slitid, spat_pixpos, snr_bar): """ - - Print out a SNR report for the reference object used to compute the weights for multislit 2D coadds. + + Print out a SNR report for the reference object used to compute the weights for multislit 2D coadds. Args: slitid (:obj:`int`): @@ -1599,7 +1599,7 @@ def snr_report(self, slitid, spat_pixpos, snr_bar): Array of spatial pixel position of the reference object in the slit for each exposure shape = (nexp,) snr_bar (:obj:`numpy.ndarray`): Array of average S/N ratios for the reference object in each exposure, shape = (nexp,) - + Returns: @@ -1699,7 +1699,8 @@ def get_maskdef_dict(self, slit_idx, ref_trace_stack): # get maskdef_objpos # find left edge slits_left, _, _ = \ - self.stack_dict['slits_list'][iexp].select_edges(flexure=self.stack_dict['spat_flexure_list'][iexp]) + self.stack_dict['slits_list'][iexp].select_edges( + spat_flexure=self.stack_dict['spat_flexure_list'][iexp]) # targeted object spat pix maskdef_obj_pixpos = \ self.stack_dict['slits_list'][iexp].maskdef_objpos[slit_idx] + self.maskdef_offset[iexp] \ @@ -1751,10 +1752,10 @@ class EchelleCoAdd2D(CoAdd2D): objects ``obj_id`` on each order """ - def __init__(self, spec2d_files, spectrograph, par, det=1, + def __init__(self, spec2d_files, spectrograph, par, det=1, only_slits=None, exclude_slits=None, sn_smooth_npix=None, bkg_redux=False, find_negative=False, show=False, show_peaks=False, debug_offsets=False, debug=False): - super().__init__(spec2d_files, spectrograph, det=det, + super().__init__(spec2d_files, spectrograph, det=det, only_slits=only_slits, exclude_slits=exclude_slits, sn_smooth_npix=sn_smooth_npix, bkg_redux=bkg_redux, find_negative=find_negative, par=par, show=show, show_peaks=show_peaks, debug_offsets=debug_offsets, @@ -1763,8 +1764,8 @@ def __init__(self, spec2d_files, spectrograph, par, det=1, def handle_reference_obj(self): """ - Method to handle the syntrax for the reference object to be used for offsets and weights. - + Method to handle the syntrax for the reference object to be used for offsets and weights. + """ # If a user-input object to compute offsets and weights is provided, check if it exists and get the needed info @@ -1786,7 +1787,7 @@ def handle_reference_obj(self): if flux is not None and ivar is not None and mask is not None: user_obj_exist[iexp, iord] = True - + if not np.all(user_obj_exist): msgs.error('Object provided through `user_obj_ids` does not exist in all the exposures.') @@ -1796,8 +1797,8 @@ def handle_reference_obj(self): elif len(self.stack_dict['specobjs_list']) > 0 and (self.par['coadd2d']['offsets'] == 'auto' or self.par['coadd2d']['weights'] == 'auto'): self.obj_id_bri, self.snr_bar_bri = \ self.get_brightest_obj(self.stack_dict['specobjs_list'], self.stack_dict['slits_list'][0].slitord_id) - - + + def compute_offsets(self): """ @@ -1941,7 +1942,7 @@ def snr_report(self, snr_bar): Args: snr_bar (:obj:`numpy.ndarray`): Array of average S/N ratios for the brightest object in each exposure. Shape = (nexp,) - + Returns: """ diff --git a/pypeit/coadd3d.py b/pypeit/coadd3d.py index cef8d2bab5..8de83832b2 100644 --- a/pypeit/coadd3d.py +++ b/pypeit/coadd3d.py @@ -847,7 +847,7 @@ def add_grating_corr(self, flatfile, waveimg, slits, spat_flexure=None): scale_model = flatfield.illum_profile_spectral(flatframe, waveimg, slits, slit_illum_ref_idx=self.flatpar['slit_illum_ref_idx'], model=None, trim=self.flatpar['slit_trim'], - flexure=spat_flexure, + spat_flexure=spat_flexure, smooth_npix=self.flatpar['slit_illum_smooth_npix']) else: msgs.info("Using relative spectral illumination from FlatImages") @@ -951,7 +951,7 @@ def get_alignments(self, spec2DObj, slits, spat_flexure=None): msgs.info("Using slit edges for astrometric transform") # If nothing better was provided, use the slit edges if alignments is None: - left, right, _ = slits.select_edges(flexure=spat_flexure) + left, right, _ = slits.select_edges(spat_flexure=spat_flexure) locations = [0.0, 1.0] traces = np.append(left[:, None, :], right[:, None, :], axis=1) else: @@ -1016,8 +1016,8 @@ def load(self): # Initialise the slit edges msgs.info("Constructing slit image") slits = spec2DObj.slits - slitid_img = slits.slit_img(pad=0, flexure=spat_flexure) - slits_left, slits_right, _ = slits.select_edges(flexure=spat_flexure) + slitid_img = slits.slit_img(pad=0, spat_flexure=spat_flexure) + slits_left, slits_right, _ = slits.select_edges(spat_flexure=spat_flexure) # The order of operations below proceeds as follows: # (1) Get science image @@ -1110,7 +1110,8 @@ def load(self): crval_wv = self.cubepar['wave_min'] if self.cubepar['wave_min'] is not None else wave0 cd_wv = self.cubepar['wave_delta'] if self.cubepar['wave_delta'] is not None else dwv self.all_wcs.append(self.spec.get_wcs(spec2DObj.head0, slits, detector.platescale, crval_wv, cd_wv)) - ra_img, dec_img, minmax = slits.get_radec_image(self.all_wcs[ff], alignSplines, spec2DObj.tilts, flexure=spat_flexure) + ra_img, dec_img, minmax = slits.get_radec_image(self.all_wcs[ff], alignSplines, spec2DObj.tilts, + spat_flexure=spat_flexure) # Extract wavelength and delta wavelength arrays from the images wave_ext = waveimg[onslit_gpm] diff --git a/pypeit/core/flexure.py b/pypeit/core/flexure.py index 48e6253372..2d82cfc54b 100644 --- a/pypeit/core/flexure.py +++ b/pypeit/core/flexure.py @@ -18,7 +18,6 @@ from astropy import stats from astropy import units from astropy.io import ascii -from astropy.table import Table from astropy.stats import sigma_clipped_stats import scipy.signal import scipy.optimize as opt @@ -45,7 +44,7 @@ from IPython import embed -def spat_flexure_shift(sciimg, slits, bpm=None, maxlag=20, sigdetect=10., debug=False, qa_outfile=None, qa_vrange=None): +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. @@ -60,6 +59,14 @@ def spat_flexure_shift(sciimg, slits, bpm=None, maxlag=20, sigdetect=10., debug= Science image slits (:class:`pypeit.slittrace.SlitTraceSet`): Slits object + method (:obj:`str`, optional): + Method to use to calculate the spatial flexure shift. Options + are 'detector' (default), 'slit', and 'edge'. The 'detector' + method calculates the shift for all slits simultaneously, the + '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): @@ -79,6 +86,8 @@ def spat_flexure_shift(sciimg, slits, bpm=None, maxlag=20, sigdetect=10., debug= float: The spatial flexure shift relative to the initial slits """ + # 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) @@ -116,11 +125,11 @@ def spat_flexure_shift(sciimg, slits, bpm=None, maxlag=20, sigdetect=10., debug= 'consider either changing the maximum lag for the cross-correlation, ' 'or the "spat_flexure_sigdetect" parameter, or use the manual flexure correction.') - return 0. + return np.zeros((slits.nslits, 2), dtype=float) lag0_max = np.where(lags_max == 0)[0][0] shift = round(pix_max[0] - lag0_max, 3) - msgs.info('Spatial flexure measured: {}'.format(shift)) + msgs.info('Spatial flexure measured: {} pixels'.format(shift)) if debug: # 1D plot of the cross-correlation @@ -150,13 +159,12 @@ def spat_flexure_shift(sciimg, slits, bpm=None, maxlag=20, sigdetect=10., debug= 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 shift + 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 @@ -182,8 +190,8 @@ def spat_flexure_qa(img, slits, shift, gpm=None, vrange=None, outfile=None): vrange = None # 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) + left_slits, right_slits, mask_slits = slits.select_edges(initial=True, spat_flexure=None) + left_flex, right_flex, mask = slits.select_edges(initial=True, spat_flexure=shift) if debug: # where to start and end the plot in the spatial&spectral direction @@ -1513,7 +1521,7 @@ def sky_em_residuals(wave:np.ndarray, flux:np.ndarray, wline = [line-noff,line+noff] mw = (wave > wline[0]) & (wave < wline[1]) & good_ivar - # Reuire minimum number + # Require minimum number if np.sum(mw) <= nfit_min: continue diff --git a/pypeit/core/gui/skysub_regions.py b/pypeit/core/gui/skysub_regions.py index b8e268197c..16d2b45fca 100644 --- a/pypeit/core/gui/skysub_regions.py +++ b/pypeit/core/gui/skysub_regions.py @@ -41,7 +41,7 @@ class SkySubGUI: """ def __init__(self, canvas, image, frame, outname, det, slits, axes, pypeline, spectrograph, printout=False, - runtime=False, resolution=None, initial=False, flexure=None, overwrite=False): + runtime=False, resolution=None, initial=False, spat_flexure=None, overwrite=False): """Controls for the interactive sky regions definition tasks in PypeIt. The main goal of this routine is to interactively select sky background @@ -73,8 +73,8 @@ def __init__(self, canvas, image, frame, outname, det, slits, axes, pypeline, sp initial : bool, optional To use the initial edges regardless of the presence of the tweaked edges, set this to True. - flexure : float, optional - If provided, offset each slit by this amount + spat_flexure : float, optional + If provided, offset each slit in the spatial direction by this amount runtime : bool Is the GUI being launched during data reduction? resolution : int @@ -137,7 +137,7 @@ def __init__(self, canvas, image, frame, outname, det, slits, axes, pypeline, sp self._fitr = [] # Matplotlib shaded fit region self._fita = None - self.slits_left, self.slits_right, _ = slits.select_edges(initial=initial, flexure=flexure) + self.slits_left, self.slits_right, _ = slits.select_edges(initial=initial, spat_flexure=spat_flexure) self.initialize_menu() self.reset_regions() @@ -149,7 +149,7 @@ def __init__(self, canvas, image, frame, outname, det, slits, axes, pypeline, sp @classmethod def initialize(cls, det, frame, slits, pypeline, spectrograph, outname="skyregions.fits", overwrite=False, initial=False, - flexure=None, runtime=False, printout=False): + spat_flexure=None, runtime=False, printout=False): """ Initialize the 'ObjFindGUI' window for interactive object tracing @@ -167,6 +167,8 @@ def initialize(cls, det, frame, slits, pypeline, spectrograph, outname="skyregio Name of the spectrograph printout : bool Should the results be printed to screen + spat_flexure : float, optional + If provided, offset each slit in the spatial direction by this amount runtime : bool Is this GUI being launched during a data reduction? @@ -178,7 +180,7 @@ def initialize(cls, det, frame, slits, pypeline, spectrograph, outname="skyregio # NOTE: SlitTraceSet objects always store the left and right # traces as 2D arrays, even if there's only one slit. nslit = slits.nslits - lordloc, rordloc, _ = slits.select_edges(initial=initial, flexure=flexure) + lordloc, rordloc, _ = slits.select_edges(initial=initial, spat_flexure=spat_flexure) # Determine the scale of the image med = np.median(frame) @@ -210,7 +212,7 @@ def initialize(cls, det, frame, slits, pypeline, spectrograph, outname="skyregio # Initialise the object finding window and display to screen fig.canvas.manager.set_window_title('PypeIt - Sky regions') srgui = SkySubGUI(fig.canvas, image, frame, outname, det, slits, axes, pypeline, spectrograph, - printout=printout, runtime=runtime, initial=initial, flexure=flexure, overwrite=overwrite) + printout=printout, runtime=runtime, initial=initial, spat_flexure=spat_flexure, overwrite=overwrite) plt.show() return srgui diff --git a/pypeit/core/skysub.py b/pypeit/core/skysub.py index a6704253e7..87da5b81ad 100644 --- a/pypeit/core/skysub.py +++ b/pypeit/core/skysub.py @@ -1617,16 +1617,20 @@ def generate_mask(pypeline, skyreg, slits, slits_left, slits_right, spat_flexure A 2D array containing the pixel coordinates of the left slit edges slits_right : `numpy.ndarray`_ A 2D array containing the pixel coordinates of the right slit edges - resolution: int, optional - The percentage regions will be scaled to the specified resolution. The - resolution should probably correspond to the number of spatial pixels - on the slit. + spat_flexure (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, to apply to each slit. + This is used to correct for spatial flexure. The shape of the array should + be (nslits, 2), where the first column is the shift to apply to the left + edge of each slit and the second column is the shift to apply to the + right edge of each slit. Returns ------- mask : `numpy.ndarray`_ Boolean mask containing sky regions """ + # Check the input + _spat_flexure = np.zeros((slits.nslits, 2)) if spat_flexure is None else spat_flexure # Grab the resolution that was used to generate skyreg resolution = skyreg[0].size # Using the left/right slit edge traces, generate a series of traces that mark the @@ -1636,7 +1640,11 @@ def generate_mask(pypeline, skyreg, slits, slits_left, slits_right, spat_flexure # regardless of which slit the sky regions falls in) left_edg, righ_edg = np.zeros((slits.nspec, 0)), np.zeros((slits.nspec, 0)) spec_min, spec_max = np.array([]), np.array([]) + # Create a dummy spatial flexure that determines the spatial flexure at the edge of each sky region + skyreg_spat_flexure = np.zeros((0, 2)) for sl in range(slits.nslits): + this_spat_flexure_left = _spat_flexure[sl, 0] + this_spat_flexure_diff = _spat_flexure[sl, 1] - _spat_flexure[sl, 0] # Calculate the slit width diff = slits_right[:, sl] - slits_left[:, sl] # Break up the slit into `resolution` subpixels @@ -1654,6 +1662,10 @@ def generate_mask(pypeline, skyreg, slits, slits_left, slits_right, spat_flexure nreg += 1 spec_min = np.append(spec_min, slits.specmin[sl]) spec_max = np.append(spec_max, slits.specmax[sl]) + # Determine the spatial flexure at the edge of the sky regions -- linearly interpolate from the left edge + flex_left = this_spat_flexure_left + this_spat_flexure_diff * wl[rr]/(resolution-1.0) + flex_right = this_spat_flexure_left + this_spat_flexure_diff * wr[rr]/(resolution-1.0) + skyreg_spat_flexure = np.append(skyreg_spat_flexure, np.array([[flex_left, flex_right]]), axis=0) # Check if no regions were added if left_edg.shape[1] == 0: @@ -1666,7 +1678,7 @@ def generate_mask(pypeline, skyreg, slits, slits_left, slits_right, spat_flexure specmax=spec_max, binspec=slits.binspec, binspat=slits.binspat, pad=0) # Generate the mask, and return - return (slitreg.slit_img(use_spatial=False, flexure=spat_flexure) >= 0).astype(bool) + return (slitreg.slit_img(use_spatial=False, spat_flexure=skyreg_spat_flexure) >= 0).astype(bool) diff --git a/pypeit/core/tracewave.py b/pypeit/core/tracewave.py index 622094a9be..9f68ac2419 100644 --- a/pypeit/core/tracewave.py +++ b/pypeit/core/tracewave.py @@ -851,43 +851,118 @@ def fit_tilts(trc_tilt_dict, thismask, slit_cen, spat_order=3, spec_order=4, max # msgs.info("RMS/FWHM: {}".format(rms_real/fwhm)) -def fit2tilts(shape, coeff2, func2d, spat_shift=None): +def fit2tilts_prepareSlit(slit_left, slit_right, thismask_science, spat_flexure=None): """ - Evaluate the wavelength tilt model over the full image. + Prepare the slit for the fit2tilts function + + Parameters + ---------- + slit_left : `numpy.ndarray`_ + Left slit edge + slit_right : `numpy.ndarray`_ + Right slit edge + thismask_science : `numpy.ndarray`_ + Boolean mask for the science pixels in this slit (True = on slit) + spat_flexure : `numpy.ndarray`_, optional + Spatial flexure in pixels. This should be a two element array with the first element being the flexure at the + left edge of the slit and the second element being the flexure at the right edge of the slit. If None, no + flexure is applied. + + Returns + ------- + :obj:`tuple` of two `numpy.ndarray`_ + Tuple containing the spectral and spatial coordinates of the slit, including spatial flexure. + These variables are to be used in the fit2tilts function. + """ + # Get the image shape + img_shape = thismask_science.shape + # Check the spatial flexure input + if spat_flexure is not None and len(spat_flexure) != 2: + msgs.error('Spatial flexure must be a two element array') + _spat_flexure = np.zeros(2) if spat_flexure is None else spat_flexure + # Check dimensions + if len(slit_left) != len(slit_right): + msgs.error('Slit left and right edges must have the same length') + if len(slit_left) != img_shape[0]: + msgs.error('Slit edges must have the same length as the spectral dimension of the image') + + # Prepare the spatial and spectral coordinates + nspec, nspat = img_shape + spec_vec = np.arange(nspec) + spat_vec = np.arange(nspat) + spat_img, spec_img = np.meshgrid(spat_vec, spec_vec) + # Evaluate the slit coordinates + _spec_eval = spec_img[thismask_science] / (nspec - 1) + flex_coord = (spat_img - slit_left[:, None]) / (slit_right - slit_left)[:, None] + # Linearly interpolate the spatial flexure over the slit coordinates + spat_eval = spat_img + _spat_flexure[0] + flex_coord * (_spat_flexure[1] - _spat_flexure[0]) + _spat_eval = spat_eval[thismask_science] / (nspat - 1) + return _spec_eval, _spat_eval + + +def fit2tilts(coeff2, func2d, shape=None, spec_eval=None, spat_eval=None): + """ + Evaluate the wavelength tilt model over the full image. Note that this function + requires either shape or both spec_eval and spat_eval to be provided. If all three + are provided, spec_eval and spat_eval will be used. Parameters ---------- - shape : tuple of ints, - shape of image coeff2 : `numpy.ndarray`_, float result of griddata tilt fit func2d : str the 2d function used to fit the tilts - spat_shift : float, optional - Spatial shift to be added to image pixels before evaluation - If you are accounting for flexure, then you probably wish to - input -1*flexure_shift into this parameter. + shape : tuple of ints, optional + Shape of image. Only used if spat_eval and spec_eval are not provided. + spat_eval : `numpy.ndarray`_, optional + 1D array indicating how spatial pixel locations move across the + image. If spat_eval is provided, spec_eval must also be provided. + spat_shift is ignored when spat_eval and spec_eval are provided. + If you wish to account for spatial flexure, then you should include + the spatial flexure into this parameter. spat_eval should be given + as spatial_coordinates + spatial_flexure. + spec_eval : `numpy.ndarray`_, optional + 1D array indicating how spectral pixel locations move across the + image. If spec_eval is provided, spat_eval must also be provided. Returns ------- tilts : `numpy.ndarray`_, float - Image indicating how spectral pixel locations move across the - image. This output is used in the pipeline. - + Image indicating how spectral pixel locations move across the image. """ - # Init - _spat_shift = 0. if spat_shift is None else spat_shift + # Determine if coordinates have been supplied that include the spatial flexure. + if spec_eval is not None and spat_eval is not None: + _spec_eval = spec_eval + _spat_eval = spat_eval + else: + # Print a warning just in case only one was provided + if (spec_eval is None and spat_eval is not None) or (spec_eval is not None and spat_eval is None): + msgs.warn('Both spec_eval and spat_eval must be provided.' + msgs.newline() + + 'Only one variable provided, so a new (full) grid will be generated.') + # Print a warning if neither are provided + if spec_eval is None and spat_eval is None: + msgs.warn('No spatial and spectral coordinates provided.' + msgs.newline() + + 'A new (full) grid will be generated.') + # Print a warning is shape is not provided + if shape is None: + msgs.error('No shape provided for the image.' + msgs.newline() + + 'You must provide either `shape` or both `spat_eval` and `spec_eval`.') + msgs.warn("Assuming no spatial flexure.") + _spat_shift = 0.0 + # Setup the evaluation grid + nspec, nspat = shape + xnspecmin1 = float(nspec - 1) + xnspatmin1 = float(nspat - 1) + spec_vec = np.arange(nspec) + spat_vec = np.arange(nspat) - _spat_shift + spat_img, spec_img = np.meshgrid(spat_vec, spec_vec) + _spec_eval = spec_img / xnspecmin1 + _spat_eval = spat_img / xnspatmin1 + # Compute the tilts image - nspec, nspat = shape - xnspecmin1 = float(nspec - 1) - xnspatmin1 = float(nspat - 1) - spec_vec = np.arange(nspec) - spat_vec = np.arange(nspat) - _spat_shift - spat_img, spec_img = np.meshgrid(spat_vec, spec_vec) - # pypeitFit = fitting.PypeItFit(fitc=coeff2, minx=0.0, maxx=1.0, minx2=0.0, maxx2=1.0, func=func2d) - tilts = pypeitFit.eval(spec_img / xnspecmin1, x2=spat_img / xnspatmin1) + tilts = pypeitFit.eval(_spec_eval, x2=_spat_eval) # Added this to ensure that tilts are never crazy values due to extrapolation of fits which can break # wavelength solution fitting return np.fmax(np.fmin(tilts, 1.2), -0.2) diff --git a/pypeit/extraction.py b/pypeit/extraction.py index 68f31a3e3f..34d0e8a62e 100644 --- a/pypeit/extraction.py +++ b/pypeit/extraction.py @@ -160,8 +160,8 @@ def __init__(self, sciImg, slits, sobjs_obj, spectrograph, par, objtype, global_ # frames. Is that okay for this usage? # Flexure self.spat_flexure_shift = None - if (objtype == 'science' and self.par['scienceframe']['process']['spat_flexure_correct']) or \ - (objtype == 'standard' and self.par['calibrations']['standardframe']['process']['spat_flexure_correct']): + if (objtype == 'science' and self.par['scienceframe']['process']['spat_flexure_method'] != "skip") or \ + (objtype == 'standard' and self.par['calibrations']['standardframe']['process']['spat_flexure_method'] != "skip"): self.spat_flexure_shift = self.sciImg.spat_flexure elif objtype == 'science_coadd2d': self.spat_flexure_shift = None @@ -214,14 +214,15 @@ def __init__(self, sciImg, slits, sobjs_obj, spectrograph, par, objtype, global_ self.waveTilts = waveTilts self.waveTilts.is_synced(self.slits) # Deal with Flexure - if self.par['calibrations']['tiltframe']['process']['spat_flexure_correct']: + if self.par['calibrations']['tiltframe']['process']['spat_flexure_method'] != "skip": _spat_flexure = 0. if self.spat_flexure_shift is None else self.spat_flexure_shift # If they both shifted the same, there will be no reason to shift the tilts tilt_flexure_shift = _spat_flexure - self.waveTilts.spat_flexure else: tilt_flexure_shift = self.spat_flexure_shift msgs.info("Generating tilts image from fit in waveTilts") - self.tilts = self.waveTilts.fit2tiltimg(self.slitmask, flexure=tilt_flexure_shift) + self.tilts = self.waveTilts.fit2tiltimg(self.slitmask, self.slits_left, self.slits_right, + spat_flexure=tilt_flexure_shift) elif waveTilts is None and tilts is not None: msgs.info("Using user input tilts image") self.tilts = tilts @@ -311,7 +312,7 @@ def initialize_slits(self, slits, initial=False): # TODO JFH: his is an ugly hack for the present moment until we get the slits object sorted out self.slits_left, self.slits_right, _ \ - = self.slits.select_edges(initial=initial, flexure=self.spat_flexure_shift) + = self.slits.select_edges(initial=initial, spat_flexure=self.spat_flexure_shift) # This matches the logic below that is being applied to the slitmask. Better would be to clean up slits to # to return a new slits object with the desired selection criteria which would remove the ambiguity # about whether the slits and the slitmask are in sync. @@ -322,7 +323,7 @@ def initialize_slits(self, slits, initial=False): #self.slits_right = slits_right[:, gpm] # Slitmask - self.slitmask = self.slits.slit_img(initial=initial, flexure=self.spat_flexure_shift, + self.slitmask = self.slits.slit_img(initial=initial, spat_flexure=self.spat_flexure_shift, exclude_flag=self.slits.bitmask.exclude_for_reducing) # Now add the slitmask to the mask (i.e. post CR rejection in proc) # NOTE: this uses the par defined by EdgeTraceSet; this will diff --git a/pypeit/find_objects.py b/pypeit/find_objects.py index f2403a6d43..847533a12a 100644 --- a/pypeit/find_objects.py +++ b/pypeit/find_objects.py @@ -91,7 +91,13 @@ class FindObjects: slits (:class:`~pypeit.slittrace.SlitTraceSet`): sobjs_obj (:class:`pypeit.specobjs.SpecObjs`): Objects found - spat_flexure_shift (:obj:`float`): + spat_flexure_shift (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. tilts (`numpy.ndarray`_): WaveTilts images generated on-the-spot waveimg (`numpy.ndarray`_): @@ -152,8 +158,8 @@ def __init__(self, sciImg, slits, spectrograph, par, objtype, wv_calib=None, wav # frames. Is that okay for this usage? # Flexure self.spat_flexure_shift = None - if (objtype == 'science' and self.par['scienceframe']['process']['spat_flexure_correct']) or \ - (objtype == 'standard' and self.par['calibrations']['standardframe']['process']['spat_flexure_correct']): + if (objtype == 'science' and self.par['scienceframe']['process']['spat_flexure_method'] != "skip") or \ + (objtype == 'standard' and self.par['calibrations']['standardframe']['process']['spat_flexure_method'] != "skip"): self.spat_flexure_shift = self.sciImg.spat_flexure elif objtype == 'science_coadd2d': self.spat_flexure_shift = None @@ -221,14 +227,16 @@ def __init__(self, sciImg, slits, spectrograph, par, objtype, wv_calib=None, wav self.waveTilts = waveTilts self.waveTilts.is_synced(self.slits) # Deal with Flexure - if self.par['calibrations']['tiltframe']['process']['spat_flexure_correct']: - _spat_flexure = 0. if self.spat_flexure_shift is None else self.spat_flexure_shift + if self.par['calibrations']['tiltframe']['process']['spat_flexure_method'] != "skip": + _spat_flexure = np.zeros((slits.nslits, 2)) if self.spat_flexure_shift is None \ + else self.spat_flexure_shift # If they both shifted the same, there will be no reason to shift the tilts tilt_flexure_shift = _spat_flexure - self.waveTilts.spat_flexure else: tilt_flexure_shift = self.spat_flexure_shift msgs.info("Generating tilts image from fit in waveTilts") - self.tilts = self.waveTilts.fit2tiltimg(self.slitmask, flexure=tilt_flexure_shift) + self.tilts = self.waveTilts.fit2tiltimg(self.slitmask, self.slits_left, self.slits_right, + spat_flexure=tilt_flexure_shift) elif waveTilts is None and tilts is not None: msgs.info("Using user input tilts image") self.tilts = tilts @@ -296,7 +304,7 @@ def initialize_slits(self, slits, initial=False): # Select the edges to use # TODO JFH: his is an ugly hack for the present moment until we get the slits object sorted out self.slits_left, self.slits_right, _ \ - = self.slits.select_edges(initial=initial, flexure=self.spat_flexure_shift) + = self.slits.select_edges(initial=initial, spat_flexure=self.spat_flexure_shift) # This matches the logic below that is being applied to the slitmask. Better would be to clean up slits to # to return a new slits object with the desired selection criteria which would remove the ambiguity # about whether the slits and the slitmask are in sync. @@ -308,7 +316,7 @@ def initialize_slits(self, slits, initial=False): # Slitmask - self.slitmask = self.slits.slit_img(initial=initial, flexure=self.spat_flexure_shift, + self.slitmask = self.slits.slit_img(initial=initial, spat_flexure=self.spat_flexure_shift, exclude_flag=self.slits.bitmask.exclude_for_reducing+['BOXSLIT']) # Now add the slitmask to the mask (i.e. post CR rejection in proc) # NOTE: this uses the par defined by EdgeTraceSet; this will @@ -1083,8 +1091,8 @@ def joint_skysub(self, skymask=None, update_crmask=True, trim_edg=(0,0), model_ivar = self.sciImg.ivar sl_ref = self.par['calibrations']['flatfield']['slit_illum_ref_idx'] # Prepare the slitmasks for the relative spectral illumination - slitmask = self.slits.slit_img(pad=0, flexure=self.spat_flexure_shift) - slitmask_trim = self.slits.slit_img(pad=-3, flexure=self.spat_flexure_shift) + slitmask = self.slits.slit_img(pad=0, spat_flexure=self.spat_flexure_shift) + slitmask_trim = self.slits.slit_img(pad=-3, spat_flexure=self.spat_flexure_shift) for nn in range(numiter): msgs.info("Performing iterative joint sky subtraction - ITERATION {0:d}/{1:d}".format(nn+1, numiter)) # TODO trim_edg is in the parset so it should be passed in here via trim_edg=tuple(self.par['reduce']['trim_edge']), @@ -1230,8 +1238,8 @@ def calculate_flexure(self, global_sky): new_slitshift = self.slitshift + this_slitshift # Now report the flexure values for slit_idx, slit_spat in enumerate(self.slits.spat_id): - msgs.info("Flexure correction, slit {0:d} (spat id={1:d}): {2:.3f} pixels".format(1+slit_idx, slit_spat, - self.slitshift[slit_idx])) + msgs.info("Spectral flexure correction, slit {0:d} (spat id={1:d}): {2:.3f} pixels".format(1+slit_idx, slit_spat, + new_slitshift[slit_idx])) # Save QA # TODO :: Need to implement QA once the flexure code has been tidied up, and this routine has been moved # out of the find_objects() class. diff --git a/pypeit/flatfield.py b/pypeit/flatfield.py index e1f5d0623b..fb01e5b98b 100644 --- a/pypeit/flatfield.py +++ b/pypeit/flatfield.py @@ -360,15 +360,20 @@ def fit2illumflat(self, slits, frametype='illum', finecorr=False, initial=False, zeroth order correction (finecorr=False) initial (bool, optional): If True, the initial slit edges will be used - spat_flexure (float, optional): - Spatial flexure in pixels + spat_flexure (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. Returns: `numpy.ndarray`_: An image of the spatial illumination profile for all slits. """ # Check spatial flexure type - if spat_flexure is not None and not isinstance(spat_flexure, float): - msgs.error('Spatial flexure must be None or float.') + if spat_flexure is not None and not isinstance(spat_flexure, np.ndarray): + msgs.error('Spatial flexure must be None or a 2D numpy array.') # Initialise the returned array illumflat = np.ones(self.shape, dtype=float) # Load spatial bsplines @@ -386,12 +391,12 @@ def fit2illumflat(self, slits, frametype='illum', finecorr=False, initial=False, continue # Skip those without a bspline # DO it - _slitid_img = slits.slit_img(slitidx=slit_idx, initial=initial, flexure=spat_flexure) + _slitid_img = slits.slit_img(slitidx=slit_idx, initial=initial, spat_flexure=spat_flexure) onslit = _slitid_img == slits.spat_id[slit_idx] spat_coo = slits.spatial_coordinate_image(slitidx=slit_idx, initial=initial, slitid_img=_slitid_img, - flexure_shift=spat_flexure) + spat_flexure=spat_flexure) if finecorr: spec_coo = np.where(onslit)[0] / (slits.nspec - 1) illumflat[onslit] = spat_bsplines[slit_idx].eval(spat_coo[onslit], spec_coo) @@ -705,11 +710,12 @@ def build_waveimg(self): if self.wavetilts is None or self.wv_calib is None: msgs.error("Wavelength calib or tilts are not available. Cannot generate wavelength image.") else: - flex = self.wavetilts.spat_flexure - slitmask = self.slits.slit_img(initial=True, flexure=flex) - tilts = self.wavetilts.fit2tiltimg(slitmask, flexure=flex) + spat_flexure = self.wavetilts.spat_flexure + left, right, msk = self.slits.select_edges(initial=True, spat_flexure=spat_flexure) + slitmask = self.slits.slit_img(initial=True, spat_flexure=spat_flexure) + tilts = self.wavetilts.fit2tiltimg(slitmask, left, right, spat_flexure=spat_flexure) # Save to class attribute for inclusion in the Flat calibration frame - self.waveimg = self.wv_calib.build_waveimg(tilts, self.slits, spat_flexure=flex) + self.waveimg = self.wv_calib.build_waveimg(tilts, self.slits, spat_flexure=spat_flexure) def show(self, wcs_match=True): """ @@ -991,12 +997,16 @@ def fit(self, spat_illum_only=False, doqa=True, debug=False): # Create the tilts image for this slit if self.slitless: - tilts = np.tile(np.arange(rawflat.shape[0]) / rawflat.shape[0], (rawflat.shape[1], 1)).T + tilts = np.tile(np.arange(rawflat.shape[0]) / (rawflat.shape[0]-1), (rawflat.shape[1], 1)).T else: - # TODO -- JFH Confirm the sign of this shift is correct! - _flexure = 0. if self.wavetilts.spat_flexure is None else self.wavetilts.spat_flexure - tilts = tracewave.fit2tilts(rawflat.shape, self.wavetilts['coeffs'][:,:,slit_idx], - self.wavetilts['func2d'], spat_shift=-1*_flexure) + _flexure = np.zeros(2) if self.wavetilts.spat_flexure is None else self.wavetilts.spat_flexure[slit_idx,:] + _spec_eval, _spat_eval = tracewave.fit2tilts_prepareSlit(self.slits.left_init[:, slit_idx], + self.slits.right_init[:, slit_idx], + onslit_init, _flexure) + tilts = np.zeros(rawflat.shape, dtype=float) + tilts[onslit_init] = tracewave.fit2tilts(self.wavetilts['coeffs'][:,:,slit_idx], + self.wavetilts['func2d'], + spec_eval=_spec_eval, spat_eval=_spat_eval) # Convert the tilt image to an image with the spectral pixel index spec_coo = tilts * (nspec-1) @@ -1465,7 +1475,7 @@ def spatial_fit_finecorr(self, normed, onslit_tweak, slit_idx, slit_spat, gpm, # check id self.waveimg is available if self.waveimg is None: msgs.warn("Cannot perform the fine correction to the spatial illumination without the wavelength image.") - return + return np.ones_like(self.rawflatimg.image) # TODO :: Include fit_order in the parset?? fit_order = np.array([3, 6]) slit_txt = self.slits.slitord_txt @@ -1477,8 +1487,8 @@ def spatial_fit_finecorr(self, normed, onslit_tweak, slit_idx, slit_spat, gpm, onslit_tweak_trim = self.slits.slit_img(pad=-slit_trim, slitidx=slit_idx, initial=False) == slit_spat # Setup slitimg = (slit_spat + 1) * onslit_tweak.astype(int) - 1 # Need to +1 and -1 so that slitimg=-1 when off the slit - - left, right, msk = self.slits.select_edges(flexure=self.wavetilts.spat_flexure if self.wavetilts is not None else 0.0) + spat_flexure = self.wavetilts.spat_flexure if self.wavetilts is not None else np.zeros((self.slits.nslits, 2)) + left, right, msk = self.slits.select_edges(spat_flexure=spat_flexure) this_left = left[:, slit_idx] this_right = right[:, slit_idx] slitlen = int(np.median(this_right - this_left)) @@ -1488,7 +1498,7 @@ def spatial_fit_finecorr(self, normed, onslit_tweak, slit_idx, slit_spat, gpm, this_wave = self.waveimg[this_slit] xpos_img = self.slits.spatial_coordinate_image(slitidx=slit_idx, slitid_img=slitimg, - flexure_shift=self.wavetilts.spat_flexure if self.wavetilts is not None else 0.0) + spat_flexure=spat_flexure) # Generate the trimmed versions for fitting this_slit_trim = np.where(onslit_tweak_trim & self.rawflatimg.select_flag(invert=True)) this_wave_trim = self.waveimg[this_slit_trim] @@ -1585,10 +1595,11 @@ def extract_structure(self, rawflat_orig, slit_trim=3): # Now fit the spectral profile # TODO: Should this be *any* flag, or just BPM? gpm = self.rawflatimg.select_flag(flag='BPM', invert=True) + spat_flexure = self.wavetilts.spat_flexure if self.wavetilts is not None else np.zeros((self.slits.nslits, 2)) scale_model = illum_profile_spectral(rawflat, self.waveimg, self.slits, slit_illum_ref_idx=self.flatpar['slit_illum_ref_idx'], model=None, gpmask=gpm, skymask=None, trim=self.flatpar['slit_trim'], - flexure=self.wavetilts.spat_flexure if self.wavetilts is not None else 0.0, + spat_flexure=spat_flexure, smooth_npix=self.flatpar['slit_illum_smooth_npix']) # Trim the edges by a few pixels to avoid edge effects onslits_trim = gpm & (self.slits.slit_img(pad=-slit_trim, initial=False) != -1) @@ -1637,7 +1648,7 @@ def spectral_illumination(self, gpm=None, debug=False): # check if the waveimg is available if self.waveimg is None: msgs.warn("Cannot perform the spectral illumination without the wavelength image.") - return None + return np.ones_like(self.rawflatimg.image) msgs.info('Performing a joint fit to the flat-field response') # Grab some parameters trim = self.flatpar['slit_trim'] @@ -1648,12 +1659,11 @@ def spectral_illumination(self, gpm=None, debug=False): gpm = self.rawflatimg.select_flag(flag='BPM', invert=True) # Obtain relative spectral illumination + spat_flexure = self.wavetilts.spat_flexure if self.wavetilts is not None else np.zeros((self.slits.nslits, 2)) return illum_profile_spectral(rawflat, self.waveimg, self.slits, slit_illum_ref_idx=self.flatpar['slit_illum_ref_idx'], - model=None, gpmask=gpm, skymask=None, trim=trim, - flexure=self.wavetilts.spat_flexure if self.wavetilts is not None else 0.0, - smooth_npix=self.flatpar['slit_illum_smooth_npix'], - debug=debug) + model=None, gpmask=gpm, skymask=None, trim=trim, spat_flexure=spat_flexure, + smooth_npix=self.flatpar['slit_illum_smooth_npix'], debug=debug) def tweak_slit_edges(self, left, right, spat_coo, norm_flat, method='threshold', thresh=0.93, maxfrac=0.1, debug=False): @@ -2154,7 +2164,7 @@ def show_flats(image_list, wcs_match=True, slits=None, waveimg=None): # TODO :: This could possibly be moved to core.flat def illum_profile_spectral(rawimg, waveimg, slits, slit_illum_ref_idx=0, smooth_npix=None, polydeg=None, - model=None, gpmask=None, skymask=None, trim=3, flexure=None, maxiter=5, debug=False): + model=None, gpmask=None, skymask=None, trim=3, spat_flexure=None, maxiter=5, debug=False): """ Determine the relative spectral illumination of all slits. Currently only used for image slicer IFUs. @@ -2174,17 +2184,22 @@ def illum_profile_spectral(rawimg, waveimg, slits, slit_illum_ref_idx=0, smooth_ polydeg : int, optional Degree of polynomial to be used for determining relative spectral sensitivity. If None, coadd.smooth_weights will be used, with the smoothing length set to smooth_npix. - model : `numpy.ndarray`_, None + model : `numpy.ndarray`_, optional A model of the rawimg data. If None, rawimg will be used. - gpmask : `numpy.ndarray`_, None + gpmask : `numpy.ndarray`_, optional Good pixel mask - skymask : `numpy.ndarray`_, None + skymask : `numpy.ndarray`_, optional Sky mask trim : int Number of pixels to trim from the edges of the slit when deriving the spectral illumination - flexure : float, None - Spatial flexure + spat_flexure : `numpy.ndarray`_, optional + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. maxiter : :obj:`int` Maximum number of iterations to perform debug : :obj:`bool` @@ -2205,8 +2220,8 @@ def illum_profile_spectral(rawimg, waveimg, slits, slit_illum_ref_idx=0, smooth_ gpm = gpmask if (gpmask is not None) else np.ones_like(rawimg, dtype=bool) modelimg = model if (model is not None) else rawimg.copy() # Setup the slits - slitid_img = slits.slit_img(pad=0, flexure=flexure) - slitid_img_trim = slits.slit_img(pad=-trim, flexure=flexure) + slitid_img = slits.slit_img(pad=0, spat_flexure=spat_flexure) + slitid_img_trim = slits.slit_img(pad=-trim, spat_flexure=spat_flexure) scaleImg = np.ones_like(rawimg) modelimg_copy = modelimg.copy() # Obtain the minimum and maximum wavelength of all slits diff --git a/pypeit/images/buildimage.py b/pypeit/images/buildimage.py index b09581ed56..1642bb4c6f 100644 --- a/pypeit/images/buildimage.py +++ b/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,6 +250,7 @@ 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: @@ -254,13 +258,14 @@ def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm= 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']) diff --git a/pypeit/images/pypeitimage.py b/pypeit/images/pypeitimage.py index dc8e07b591..d69f174b22 100644 --- a/pypeit/images/pypeitimage.py +++ b/pypeit/images/pypeitimage.py @@ -94,7 +94,7 @@ class PypeItImage(datamodel.DataContainer): # TODO These docs are confusing. The __init__ method needs to be documented just as it is for # every other class that we have written in PypeIt, i.e. the arguments all need to be documented. They are not # documented here and instead we have the odd Args documentation above. - version = '1.3.0' + version = '1.3.1' """Datamodel version number""" datamodel = {'PYP_SPEC': dict(otype=str, descr='PypeIt spectrograph name'), @@ -128,9 +128,12 @@ class PypeItImage(datamodel.DataContainer): 'exptime': dict(otype=(int, float), descr='Effective exposure time (s)'), 'noise_floor': dict(otype=float, descr='Noise floor included in variance'), 'shot_noise': dict(otype=bool, descr='Shot-noise included in variance'), - 'spat_flexure': dict(otype=float, + 'spat_flexure': dict(otype=np.ndarray, atype=np.floating, descr='Shift, in spatial pixels, between this image ' - 'and SlitTrace'), + 'and SlitTrace. Shape is (nslits, 2), where' + 'spat_flexure[i,0] is the spatial shift of the left ' + 'edge of slit i and spat_flexure[i,1] is the spatial ' + 'shift of the right edge of slit i.'), 'filename': dict(otype=str, descr='Filename for the image'),} """Data model components.""" @@ -782,9 +785,16 @@ def sub(self, other): # Spatial flexure spat_flexure = self.spat_flexure if other.spat_flexure is not None and spat_flexure is not None \ - and other.spat_flexure != spat_flexure: - msgs.warn(f'Spatial flexure different for images being subtracted ({spat_flexure} ' - f'vs. {other.spat_flexure}). Adopting {np.max(np.abs([spat_flexure, other.spat_flexure]))}.') + and not np.array_equal(other.spat_flexure, spat_flexure): + msgs.warn(f'Spatial flexure different for images being subtracted. Adopting ' + f'the maximum spatial flexure of each individual edge.') + # Loop through all slit edges and find the largest flexure + for ii in range(spat_flexure.shape[0]): + # Assign the largest flexure (irrespective of sign) for each edge + if np.abs(other.spat_flexure[ii,0]) > np.abs(spat_flexure[ii,0]): + spat_flexure[ii,0] = other.spat_flexure[ii,0] + if np.abs(other.spat_flexure[ii,1]) > np.abs(spat_flexure[ii,1]): + spat_flexure[ii,1] = other.spat_flexure[ii,1] # Create a copy of the detector, if it is defined, to be used when # creating the new pypeit image below diff --git a/pypeit/images/rawimage.py b/pypeit/images/rawimage.py index d7e46f152d..3d8908db02 100644 --- a/pypeit/images/rawimage.py +++ b/pypeit/images/rawimage.py @@ -113,8 +113,12 @@ class RawImage: datasec_img (`numpy.ndarray`_): Image identifying which amplifier was used to read each section of the *processed* image. - spat_flexure_shift (:obj:`float`): - The spatial flexure shift in pixels, if calculated + spat_flexure_shift (`numpy.ndarray`_): + The spatial flexure shift in pixels, if calculated. This is a 2D array + of shape (nslits, 2) where spat_flexure_shift[i,0] is the shift in the + spatial direction for the left edge of slit i and spat_flexure_shift[i,1] + is the shift in the spatial direction for the right edge of slit i. + """ def __init__(self, ifile, spectrograph, det): @@ -239,7 +243,7 @@ def use_slits(self): """ if self.par is None: return False - return self.par['spat_flexure_correct'] or (self.use_flat and self.par['use_illumflat']) + return (self.par['spat_flexure_method'] != "skip") or (self.use_flat and self.par['use_illumflat']) def apply_gain(self, force=False): """ @@ -399,7 +403,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. @@ -504,7 +508,7 @@ def process(self, par, bpm=None, scattlight=None, flatimages=None, bias=None, sl Bias image for bias subtraction. slits (:class:`~pypeit.slittrace.SlitTraceSet`, optional): Used to calculate spatial flexure between the image and the - slits, if requested via the ``spat_flexure_correct`` parameter + slits, if requested via the ``spat_flexure_method`` parameter in :attr:`par`; see :func:`~pypeit.core.flexure.spat_flexure_shift`. Also used to construct the slit illumination profile, if requested via the @@ -517,6 +521,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. @@ -544,7 +551,7 @@ def process(self, par, bpm=None, scattlight=None, flatimages=None, bias=None, sl msgs.error('No dark available for dark subtraction!') if self.par['subtract_scattlight'] and scattlight is None: msgs.error('Scattered light subtraction requested, but scattered light model not provided.') - if self.par['spat_flexure_correct'] and slits is None: + if (self.par['spat_flexure_method'] != "skip") and slits is None: msgs.error('Spatial flexure correction requested but no slits provided.') if self.use_flat and flatimages is None: msgs.error('Flat-field corrections requested but no flat-field images generated ' @@ -669,8 +676,10 @@ 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, debug=debug) \ - if self.par['spat_flexure_correct'] 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, manual_spat_flexure=manual_spat_flexure, + debug=debug) # - Subtract scattered light... this needs to be done before flatfielding. if self.par['subtract_scattlight']: @@ -763,7 +772,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, debug=False): + def spatial_flexure_shift(self, slits, force=False, manual_spat_flexure=np.ma.masked, debug=False): """ Calculate a spatial shift in the edge traces due to flexure. @@ -771,16 +780,23 @@ def spatial_flexure_shift(self, slits, force=False, debug=False): :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. debug (:obj:`bool`, optional): Run in debug mode. Return: - float: The calculated flexure correction + `numpy.ndarray`_: The calculated flexure correction for the edge of each slit shape is (nslits, 2) """ step = inspect.stack()[0][3] @@ -792,19 +808,43 @@ def spatial_flexure_shift(self, slits, force=False, debug=False): msgs.error('CODING ERROR: Must use a single image (single detector or detector ' 'mosaic) to determine spatial flexure.') - # 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) + # 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: + # 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, bpm=self._bpm[0], + method=self.par['spat_flexure_method'], + 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]): + 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.spat_flexure_shift = flexure.spat_flexure_shift(self.image[0], slits, 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']) self.steps[step] = True # Return - return self.spat_flexure_shift + return spat_flexure def flatfield(self, flatimages, slits=None, force=False, debug=False): """ @@ -867,7 +907,7 @@ def flatfield(self, flatimages, slits=None, force=False, debug=False): illum_flat = flatimages.fit2illumflat(slits, spat_flexure=self.spat_flexure_shift, finecorr=False) illum_flat *= flatimages.fit2illumflat(slits, spat_flexure=self.spat_flexure_shift, finecorr=True) if debug: - left, right = slits.select_edges(flexure=self.spat_flexure_shift) + left, right = slits.select_edges(spat_flexure=self.spat_flexure_shift) viewer, ch = display.show_image(illum_flat, chname='illum_flat') display.show_slits(viewer, ch, left, right) # , slits.id) # @@ -1231,7 +1271,7 @@ def subtract_scattlight(self, msscattlight, slits, debug=False): f" {tmp[13]}, {tmp[14]}, {tmp[15]}]) # Polynomial terms (coefficients of spec**index)\n" print(strprint) pad = msscattlight.pad // spatbin - offslitmask = slits.slit_img(pad=pad, flexure=None) == -1 + offslitmask = slits.slit_img(pad=pad, spat_flexure=None) == -1 from matplotlib import pyplot as plt _frame = self.image[ii, ...] vmin, vmax = 0, np.max(scatt_img) @@ -1256,7 +1296,7 @@ def subtract_scattlight(self, msscattlight, slits, debug=False): elif self.par["scattlight"]["method"] == "frame": # Calculate a model specific for this frame pad = msscattlight.pad // spatbin - offslitmask = slits.slit_img(pad=pad, flexure=None) == -1 + offslitmask = slits.slit_img(pad=pad, spat_flexure=None) == -1 # Get starting parameters for the scattered light model x0, bounds = self.spectrograph.scattered_light_archive(binning, dispname) # Perform a fit to the scattered light @@ -1280,11 +1320,11 @@ def subtract_scattlight(self, msscattlight, slits, debug=False): # Check if a fine correction to the scattered light should be applied if do_finecorr: pad = self.par['scattlight']['finecorr_pad'] // spatbin - offslitmask = slits.slit_img(pad=pad, flexure=None) == -1 + offslitmask = slits.slit_img(pad=pad, spat_flexure=None) == -1 # Check if the user wishes to mask some inter-slit regions if self.par['scattlight']['finecorr_mask'] is not None: # Get the central trace of each slit - left, right, _ = slits.select_edges(flexure=None) + left, right, _ = slits.select_edges(spat_flexure=None) centrace = 0.5*(left+right) # Now mask user-defined inter-slit regions offslitmask = scattlight.mask_slit_regions(offslitmask, centrace, diff --git a/pypeit/inputfiles.py b/pypeit/inputfiles.py index ce8379cad0..63988dc392 100644 --- a/pypeit/inputfiles.py +++ b/pypeit/inputfiles.py @@ -371,14 +371,14 @@ def _read_data_file_table(lines, preserve_comments): ## Recast each as "object" in case the user has mucked with the Table ## e.g. a mix of floats and None - ## Also handle Masked columns -- fill with '' for key in tbl.keys(): # Object tbl[key] = tbl[key].data.astype(object) - if isinstance(tbl[key], column.MaskedColumn): + # RJC -- Not sure why we need to fill masked columns with empty data. Let's just retain the masked values. + # if isinstance(tbl[key], column.MaskedColumn): # Fill with empty string - tbl[key].fill_value = '' - tbl[key] = tbl[key].filled() + # tbl[key].fill_value = np.ma.masked + # tbl[key] = tbl[key].filled() # Build the table -- Old code # Because we allow (even encourage!) the users to modify entries by hand, @@ -472,20 +472,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: @@ -494,7 +495,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): diff --git a/pypeit/metadata.py b/pypeit/metadata.py index 0bbf644346..d43b74716d 100644 --- a/pypeit/metadata.py +++ b/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. @@ -1470,7 +1492,7 @@ def get_frame_types(self, flag_unknown=False, user=None, merge=True): msgs.info("Typing completed!") return self.set_frame_types(type_bits, merge=merge) - def set_pypeit_cols(self, write_bkg_pairs=False, write_manual=False, write_shift = False): + def set_pypeit_cols(self, write_bkg_pairs=False, write_manual=False, write_shift=False): """ Generate the list of columns to be included in the fitstbl (nearly the complete list). @@ -1540,7 +1562,7 @@ def set_combination_groups(self, assign_objects=True): if 'bkg_id' not in self.keys(): self['bkg_id'] = -1 if 'shift' not in self.keys(): - self['shift'] = 0 + self['shift'] = np.ma.masked # NOTE: Importantly, this if statement means that, if the user has # defined any non-negative combination IDs in their pypeit file, none of @@ -1574,9 +1596,10 @@ 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(): - self['shift'] = 0 + # Instantiate a masked array + self['shift'] = np.ma.array(np.zeros(len(self)), mask=np.ones(len(self), dtype=bool)) def write_sorted(self, ofile, overwrite=True, ignore=None, write_bkg_pairs=False, write_manual=False): @@ -1661,7 +1684,7 @@ def write_sorted(self, ofile, overwrite=True, ignore=None, def write_pypeit(self, output_path=None, cfg_lines=None, write_bkg_pairs=False, write_manual=False, - write_shift = False, + write_shift=False, configs=None, config_subdir=True, version_override=None, date_override=None): """ @@ -1741,7 +1764,7 @@ def write_pypeit(self, output_path=None, cfg_lines=None, # Grab output columns output_cols = self.set_pypeit_cols(write_bkg_pairs=write_bkg_pairs, write_manual=write_manual, - write_shift = write_shift) + write_shift=write_shift) # Write the pypeit files ofiles = [None]*len(cfg_keys) @@ -1792,7 +1815,7 @@ def write_pypeit(self, output_path=None, cfg_lines=None, pypeItFile = inputfiles.PypeItFile(cfg_lines, paths, subtbl, setup_dict) # Write pypeItFile.write(ofiles[j], version_override=version_override, - date_override=date_override) + date_override=date_override) # Return return ofiles diff --git a/pypeit/par/pypeitpar.py b/pypeit/par/pypeitpar.py index a88c46d4cb..09bf3f06b5 100644 --- a/pypeit/par/pypeitpar.py +++ b/pypeit/par/pypeitpar.py @@ -224,7 +224,7 @@ def __init__(self, trim=None, apply_gain=None, orient=None, use_pixelflat=None, use_illumflat=None, use_specillum=None, skip_write_2d=None, use_pattern=None, subtract_scattlight=None, scattlight=None, subtract_continuum=None, - spat_flexure_correct=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 @@ -366,14 +366,20 @@ def __init__(self, trim=None, apply_gain=None, orient=None, 'after ensuring the quality of the resulting reductions.' # Flexure - defaults['spat_flexure_correct'] = False - dtypes['spat_flexure_correct'] = bool - descr['spat_flexure_correct'] = 'Correct slits, illumination flat, etc. for flexure' - + defaults['spat_flexure_method'] = 'skip' + options['spat_flexure_method'] = ProcessImagesPar.valid_spatial_flexure() + dtypes['spat_flexure_method'] = str + descr['spat_flexure_method'] = 'Correct slits, illumination flat, etc. for spatial flexure. ' \ + 'Options are: {0}'.format(', '.join(options['spat_flexure_method'])) + \ + '"skip" means no correction is performed. ' \ + '"detector" means that a single shift is applied to all slits. ' \ + '"slit" means that each slit is shifted independently.' \ + '"edge" means that each slit edge is shifted independently.' + defaults['spat_flexure_maxlag'] = 20 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 ' \ @@ -484,7 +490,7 @@ 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_correct', 'spat_flexure_maxlag', 'spat_flexure_sigdetect', + 'spat_flexure_method', 'spat_flexure_maxlag', 'spat_flexure_sigdetect', 'spat_flexure_vrange', 'use_illumflat', 'use_specillum', 'skip_write_2d', 'empirical_rn', 'shot_noise', 'noise_floor', 'use_pixelflat', 'combine', 'scale_to_mean', 'correct_nonlinear', 'satpix', #'calib_setup_and_bit', @@ -517,6 +523,13 @@ def valid_combine_methods(): """ return ['median', 'mean' ] + @staticmethod + def valid_spatial_flexure(): + """ + Return the valid methods for combining frames. + """ + return ['skip', 'detector', 'slit', 'edge'] + @staticmethod def valid_saturation_handling(): """ @@ -3604,7 +3617,7 @@ def __init__(self, filt_iter=None, sobel_mode=None, edge_thresh=None, sobel_enha descr['niter_gaussian'] = 'The number of iterations of ' \ ':func:`~pypeit.core.trace.fit_trace` to use when using ' \ 'Gaussian weighting.' - + defaults['min_edge_side_sep'] = 5.0 dtypes['min_edge_side_sep'] = [int, float] descr['min_edge_side_sep'] = 'Minimum separation between same-side edges (e.g., the ' \ diff --git a/pypeit/pypeit.py b/pypeit/pypeit.py index 64d36eb7a7..6bb18cda5e 100644 --- a/pypeit/pypeit.py +++ b/pypeit/pypeit.py @@ -83,6 +83,7 @@ class PypeIt: fitstbl (:obj:`pypeit.metadata.PypeItMetaData`): holds the meta info """ + def __init__(self, pypeit_file, verbosity=2, overwrite=True, reuse_calibs=False, logname=None, show=False, redux_path=None, calib_only=False): @@ -117,14 +118,14 @@ def __init__(self, pypeit_file, verbosity=2, overwrite=True, reuse_calibs=False, # Build the meta data # - Re-initilize based on the file data msgs.info('Compiling metadata') - self.fitstbl = PypeItMetaData(self.spectrograph, self.par, + self.fitstbl = PypeItMetaData(self.spectrograph, self.par, files=self.pypeItFile.filenames, - usrdata=self.pypeItFile.data, + usrdata=self.pypeItFile.data, strict=True) # - Interpret automated or user-provided data from the PypeIt # file self.fitstbl.finalize_usr_build( - self.pypeItFile.frametypes, + self.pypeItFile.frametypes, self.pypeItFile.setup_name) # Other Internals @@ -395,7 +396,7 @@ def reduce_all(self): science_basename = [None]*len(grp_science) # Loop on unique comb_id u_combid = np.unique(self.fitstbl['comb_id'][grp_science]) - + for j, comb_id in enumerate(u_combid): # TODO: This was causing problems when multiple science frames # were provided to quicklook and the user chose *not* to stack @@ -509,11 +510,11 @@ def reduce_exposure(self, frames, bg_frames=None, std_outfile=None): # something other than the default of None. self.find_negative = (('science' in self.fitstbl['frametype'][bg_frames[0]]) | ('standard' in self.fitstbl['frametype'][bg_frames[0]])) \ - if self.par['reduce']['findobj']['find_negative'] is None else \ - self.par['reduce']['findobj']['find_negative'] + if self.par['reduce']['findobj']['find_negative'] is None else \ + self.par['reduce']['findobj']['find_negative'] else: self.bkg_redux = False - self.find_negative= False + self.find_negative = False # Container for all the Spec2DObj all_spec2d = spec2dobj.AllSpec2DObj() @@ -579,7 +580,7 @@ def reduce_exposure(self, frames, bg_frames=None, std_outfile=None): # global_sky, skymask and sciImg are needed in the extract loop initial_sky, sobjs_obj, sciImg, bkg_redux_sciimg, objFind = self.objfind_one( frames, self.det, bg_frames=bg_frames, std_outfile=std_outfile) - if len(sobjs_obj)>0: + if len(sobjs_obj) > 0: all_specobjs_objfind.add_sobj(sobjs_obj) initial_sky_list.append(initial_sky) sciImg_list.append(sciImg) @@ -589,10 +590,10 @@ def reduce_exposure(self, frames, bg_frames=None, std_outfile=None): # slitmask stuff if len(calibrated_det) > 0 and self.par['reduce']['slitmask']['assign_obj']: # get object positions from slitmask design and slitmask offsets for all the detectors - spat_flexure = np.array([ss.spat_flexure for ss in sciImg_list]) + spat_flexure = [ss.spat_flexure for ss in sciImg_list] # Grab platescale with binning bin_spec, bin_spat = parse.parse_binning(self.binning) - platescale = np.array([ss.detector.platescale*bin_spat for ss in sciImg_list]) + platescale = np.array([ss.detector.platescale * bin_spat for ss in sciImg_list]) # get the dither offset if available and if desired dither_off = None if self.par['reduce']['slitmask']['use_dither_offset']: @@ -629,18 +630,18 @@ def reduce_exposure(self, frames, bg_frames=None, std_outfile=None): # Extract all_spec2d[detname], tmp_sobjs \ - = self.extract_one(frames, self.det, sciImg_list[i], bkg_redux_sciimg_list[i], objFind_list[i], - initial_sky_list[i], all_specobjs_on_det) + = self.extract_one(frames, self.det, sciImg_list[i], bkg_redux_sciimg_list[i], objFind_list[i], + initial_sky_list[i], all_specobjs_on_det) # Hold em if tmp_sobjs.nobj > 0: all_specobjs_extract.add_sobj(tmp_sobjs) # Add calibration associations to the SpecObjs object all_specobjs_extract.calibs = calibrations.Calibrations.get_association( - self.fitstbl, self.spectrograph, self.calibrations_path, - self.fitstbl[frames[0]]['setup'], - self.fitstbl.find_frame_calib_groups(frames[0])[0], self.det, - must_exist=True, proc_only=True) + self.fitstbl, self.spectrograph, self.calibrations_path, + self.fitstbl[frames[0]]['setup'], + self.fitstbl.find_frame_calib_groups(frames[0])[0], self.det, + must_exist=True, proc_only=True) # JFH TODO write out the background frame? @@ -669,9 +670,9 @@ def get_sci_metadata(self, frame, det): # Set binning, obstime, basename, and objtype binning = self.fitstbl['binning'][frame] - obstime = self.fitstbl.construct_obstime(frame) + obstime = self.fitstbl.construct_obstime(frame) basename = self.fitstbl.construct_basename(frame, obstime=obstime) - types = self.fitstbl['frametype'][frame].split(',') + types = self.fitstbl['frametype'][frame].split(',') if 'science' in types: objtype_out = 'science' elif 'standard' in types: @@ -761,14 +762,14 @@ def objfind_one(self, frames, det, bg_frames=None, std_outfile=None): """ # Grab some meta-data needed for the reduction from the fitstbl self.objtype, self.setup, self.obstime, self.basename, self.binning \ - = self.get_sci_metadata(frames[0], det) + = self.get_sci_metadata(frames[0], det) msgs.info("Object finding begins for {} on det={}".format(self.basename, det)) # Is this a standard star? self.std_redux = self.objtype == 'standard' frame_par = self.par['calibrations']['standardframe'] \ - if self.std_redux else self.par['scienceframe'] + if self.std_redux else self.par['scienceframe'] # Get the standard trace if need be if self.std_redux is False and std_outfile is not None: @@ -778,6 +779,7 @@ 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, @@ -785,7 +787,8 @@ def objfind_one(self, frames, det, bg_frames=None, std_outfile=None): 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 @@ -797,6 +800,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. @@ -807,36 +811,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 = None - # 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_correct']) or \ - (self.objtype == 'standard' and self.par['calibrations']['standardframe']['process']['spat_flexure_correct']) or \ - manual_flexure: - if (manual_flexure or manual_flexure == 0) and not (np.issubdtype(self.fitstbl[frames[0]]["shift"], np.integer)): - msgs.info(f'Implementing manual flexure of {manual_flexure}') - spat_flexure = np.float64(manual_flexure) - sciImg.spat_flexure = spat_flexure - else: - msgs.info(f'Using auto-computed flexure') - spat_flexure = sciImg.spat_flexure - msgs.info(f'Flexure being used is: {spat_flexure}') # 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 @@ -896,8 +886,12 @@ def load_skyregions(self, initial_slits=False, scifile=None, frame=None, spat_fl frame : :obj:`int`, optional The index of the frame used to construct the calibration key. Only used if ``user_regions = user``. - spat_flexure : :obj:`float`, None, optional - The spatial flexure (measured in pixels) of the science frame relative to the trace frame. + spat_flexure (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, to apply to each slit. + This is used to correct for spatial flexure. The shape of the array should + be (nslits, 2), where the first column is the shift to apply to the left + edge of each slit and the second column is the shift to apply to the + right edge of each slit. Returns ------- @@ -913,9 +907,9 @@ def load_skyregions(self, initial_slits=False, scifile=None, frame=None, spat_fl if self.par['reduce']['skysub']['user_regions'] == 'user': # Build the file name calib_key = CalibFrame.construct_calib_key( - self.fitstbl['setup'][frame], - CalibFrame.ingest_calib_id(self.fitstbl['calib'][frame]), - self.spectrograph.get_det_name(self.det)) + self.fitstbl['setup'][frame], + CalibFrame.ingest_calib_id(self.fitstbl['calib'][frame]), + self.spectrograph.get_det_name(self.det)) regfile = buildimage.SkyRegions.construct_file_name(calib_key, calib_dir=self.calibrations_path, basename=io.remove_suffix(scifile)) @@ -934,9 +928,9 @@ def load_skyregions(self, initial_slits=False, scifile=None, frame=None, spat_fl # NOTE : Do not include spatial flexure here! # It is included when generating the mask in the return statement below slits_left, slits_right, _ \ - = self.caliBrate.slits.select_edges(initial=initial_slits, flexure=None) + = self.caliBrate.slits.select_edges(initial=initial_slits, spat_flexure=None) - maxslitlength = np.max(slits_right-slits_left) + maxslitlength = np.max(slits_right - slits_left) # Get the regions status, regions = skysub.read_userregions(skyregtxt, self.caliBrate.slits.nslits, maxslitlength) if status == 1: @@ -986,7 +980,7 @@ def extract_one(self, frames, det, sciImg, bkg_redux_sciimg, objFind, initial_sk """ # Grab some meta-data needed for the reduction from the fitstbl self.objtype, self.setup, self.obstime, self.basename, self.binning \ - = self.get_sci_metadata(frames[0], det) + = self.get_sci_metadata(frames[0], det) # Is this a standard star? self.std_redux = 'standard' in self.objtype @@ -1014,7 +1008,7 @@ def extract_one(self, frames, det, sciImg, bkg_redux_sciimg, objFind, initial_sk skymask = objFind.create_skymask(sobjs_obj) if skymask is None else skymask # DO NOT reinit_bpm, nor update_crmask bkg_redux_global_sky = objFind.global_skysub(skymask=skymask, bkg_redux_sciimg=bkg_redux_sciimg, - reinit_bpm=False, update_crmask=False, show=self.show) + reinit_bpm=False, update_crmask=False, show=self.show) scaleImg = objFind.scaleimg @@ -1073,7 +1067,7 @@ def extract_one(self, frames, det, sciImg, bkg_redux_sciimg, objFind, initial_sk # Tack on wavelength RMS for sobj in sobjs: iwv = np.where(self.caliBrate.wv_calib.spat_ids == sobj.SLITID)[0][0] - sobj.WAVE_RMS =self.caliBrate.wv_calib.wv_fits[iwv].rms + sobj.WAVE_RMS = self.caliBrate.wv_calib.wv_fits[iwv].rms # Construct table of spectral flexure spec_flex_table = Table() @@ -1102,10 +1096,10 @@ def extract_one(self, frames, det, sciImg, bkg_redux_sciimg, objFind, initial_sk spec2DObj.process_steps = sciImg.process_steps spec2DObj.calibs = calibrations.Calibrations.get_association( - self.fitstbl, self.spectrograph, self.calibrations_path, - self.fitstbl[frames[0]]['setup'], - self.fitstbl.find_frame_calib_groups(frames[0])[0], det, - must_exist=True, proc_only=True) + self.fitstbl, self.spectrograph, self.calibrations_path, + self.fitstbl[frames[0]]['setup'], + self.fitstbl.find_frame_calib_groups(frames[0])[0], det, + must_exist=True, proc_only=True) # QA spec2DObj.gen_qa() @@ -1212,8 +1206,8 @@ def save_exposure(self, frame:int, all_spec2d:spec2dobj.AllSpec2DObj, if self.par['rdx']['detnum'] is None: update_det = None elif isinstance(self.par['rdx']['detnum'], list): - update_det = [self.spectrograph.allowed_mosaics.index(d)+1 - if isinstance(d, tuple) else d for d in self.par['rdx']['detnum']] + update_det = [self.spectrograph.allowed_mosaics.index(d) + 1 + if isinstance(d, tuple) else d for d in self.par['rdx']['detnum']] else: update_det = self.par['rdx']['detnum'] @@ -1277,10 +1271,8 @@ def show_science(self): Simple print of science frames """ indx = self.fitstbl.find_frames('science') - print(self.fitstbl[['target','ra','dec','exptime','dispname']][indx]) + print(self.fitstbl[['target', 'ra', 'dec', 'exptime', 'dispname']][indx]) def __repr__(self): # Generate sets string return '<{:s}: pypeit_file={}>'.format(self.__class__.__name__, self.pypeit_file) - - diff --git a/pypeit/scattlight.py b/pypeit/scattlight.py index 675da0d07b..e4bb76909b 100644 --- a/pypeit/scattlight.py +++ b/pypeit/scattlight.py @@ -122,7 +122,7 @@ def show(self, image=None, slits=None, mask=False, wcs_match=True): wcs_match : :obj:`bool`, optional Use a reference image for the WCS and match all image in other channels to it. """ - offslitmask = slits.slit_img(pad=0, flexure=None) == -1 if mask else 1 + offslitmask = slits.slit_img(pad=0, spat_flexure=None) == -1 if mask else 1 # Prepare the frames _data = self.scattlight_raw if image is None else image diff --git a/pypeit/scripts/chk_flexure.py b/pypeit/scripts/chk_flexure.py index 74faf77b01..4dc93a14ee 100644 --- a/pypeit/scripts/chk_flexure.py +++ b/pypeit/scripts/chk_flexure.py @@ -29,7 +29,6 @@ def get_parser(cls, width=None): @staticmethod def main(args): - from IPython import embed from astropy.io import fits from pypeit import msgs from pypeit import specobjs @@ -41,7 +40,7 @@ def main(args): # Loop over the input files for in_file in args.input_file: - msgs.info(f'Checking fluxure for file: {in_file}') + msgs.info(f'Checking flexure for file: {in_file}') # What kind of file are we?? hdul = fits.open(in_file) diff --git a/pypeit/scripts/setup.py b/pypeit/scripts/setup.py index aee936337e..3d3a5e2777 100644 --- a/pypeit/scripts/setup.py +++ b/pypeit/scripts/setup.py @@ -146,7 +146,7 @@ def main(args): pypeit_files = ps.fitstbl.write_pypeit(output_path=output_path, cfg_lines=ps.user_cfg, write_bkg_pairs=args.background, write_manual=args.manual_extraction, - write_shift = args.flexure, + write_shift=args.flexure, configs=configs, version_override=args.version_override, date_override=args.date_override) diff --git a/pypeit/scripts/show_2dspec.py b/pypeit/scripts/show_2dspec.py index fd16e874fd..7aa4512406 100644 --- a/pypeit/scripts/show_2dspec.py +++ b/pypeit/scripts/show_2dspec.py @@ -271,7 +271,7 @@ def main(args): if spec2DObj.sci_spat_flexure is not None: msgs.info(f'Offseting slits by {spec2DObj.sci_spat_flexure}') slit_left, slit_right, slit_mask \ - = spec2DObj.slits.select_edges(flexure=spec2DObj.sci_spat_flexure) + = spec2DObj.slits.select_edges(spat_flexure=spec2DObj.sci_spat_flexure) slit_spat_id = spec2DObj.slits.spat_id slit_mask_id = spec2DObj.slits.maskdef_id slit_slid_IDs = spec2DObj.slits.slitord_id diff --git a/pypeit/scripts/skysub_regions.py b/pypeit/scripts/skysub_regions.py index 516d24a691..bb175af31d 100644 --- a/pypeit/scripts/skysub_regions.py +++ b/pypeit/scripts/skysub_regions.py @@ -96,7 +96,7 @@ def main(args): # Finally, initialise the GUI skyreg = SkySubGUI.initialize(det, frame, slits, pypeline, specname, outname=regfile, overwrite=args.overwrite, runtime=False, printout=True, - initial=args.initial, flexure=spat_flexure) + initial=args.initial, spat_flexure=spat_flexure) # Get the results msskyreg = skyreg.get_result() diff --git a/pypeit/slittrace.py b/pypeit/slittrace.py index ada5c686e3..9137d2d6fd 100644 --- a/pypeit/slittrace.py +++ b/pypeit/slittrace.py @@ -479,7 +479,7 @@ def get_slitlengths(self, initial=False, median=False): slitlen = right - left return np.median(slitlen, axis=1) if median else slitlen - def get_radec_image(self, wcs, alignSplines, tilts, slit_compute=None, slice_offset=None, initial=False, flexure=None): + def get_radec_image(self, wcs, alignSplines, tilts, slit_compute=None, slice_offset=None, initial=False, spat_flexure=None): """Generate an RA and DEC image for every pixel in the frame NOTE: This function is currently only used for SlicerIFU reductions. @@ -504,8 +504,8 @@ def get_radec_image(self, wcs, alignSplines, tilts, slit_compute=None, slice_off is set to 0.0. initial : bool Select the initial slit edges? - flexure : float, optional - If provided, offset each slit by this amount. + spat_flexure : float, optional + If provided, offset each slit in the spatial direction by this amount. Returns ------- @@ -542,7 +542,7 @@ def get_radec_image(self, wcs, alignSplines, tilts, slit_compute=None, slice_off decimg = np.zeros((self.nspec, self.nspat)) minmax = np.zeros((self.nslits, 2)) # Get the slit information - slitid_img_init = self.slit_img(pad=0, initial=initial, flexure=flexure) + slitid_img_init = self.slit_img(pad=0, initial=initial, spat_flexure=spat_flexure) for slit_idx, spatid in enumerate(self.spat_id): if slit_idx not in slit_compute: continue @@ -563,7 +563,7 @@ def get_radec_image(self, wcs, alignSplines, tilts, slit_compute=None, slice_off decimg[onslit] = world_dec.copy() return raimg, decimg, minmax - def select_edges(self, initial=False, flexure=None): + def select_edges(self, initial=False, spat_flexure=None): """ Select between the initial or tweaked slit edges and allow for flexure correction. @@ -578,8 +578,13 @@ def select_edges(self, initial=False, flexure=None): initial (:obj:`bool`, optional): To use the initial edges regardless of the presence of the tweaked edges, set this to True. - flexure (:obj:`float`, optional): - If provided, offset each slit by this amount + spat_flexure (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. Returns: tuple: Returns the full arrays containing the left and right @@ -593,15 +598,18 @@ def select_edges(self, initial=False, flexure=None): left, right = self.left_init, self.right_init # Add in spatial flexure? - if flexure: - self.left_flexure = left + flexure - self.right_flexure = right + flexure + self.left_flexure = left.copy() + self.right_flexure = right.copy() + if spat_flexure is not None: + for slit in range(self.nslits): + self.left_flexure[:,slit] += spat_flexure[slit, 0] + self.right_flexure[:,slit] += spat_flexure[slit, 1] left, right = self.left_flexure, self.right_flexure # Return return left.copy(), right.copy(), self.mask.copy() - def slit_img(self, pad=None, slitidx=None, initial=False, flexure=None, exclude_flag=None, + def slit_img(self, pad=None, slitidx=None, initial=False, spat_flexure=None, exclude_flag=None, use_spatial=True): r""" Construct an image identifying each pixel with its associated @@ -651,9 +659,13 @@ def slit_img(self, pad=None, slitidx=None, initial=False, flexure=None, exclude_ Warning -- This could conflict with input slitids, i.e. avoid using both use_spatial (bool, optional): If True, use self.spat_id value instead of 0-based indices - flexure (:obj:`float`, optional): - If provided, offset each slit by this amount - Done in select_edges() + spat_flexure (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. Returns: `numpy.ndarray`_: The image with the slit index @@ -673,7 +685,8 @@ def slit_img(self, pad=None, slitidx=None, initial=False, flexure=None, exclude_ spat = np.arange(self.nspat) spec = np.arange(self.nspec) - left, right, _ = self.select_edges(initial=initial, flexure=flexure) + # Get the edges + left, right, _ = self.select_edges(initial=initial, spat_flexure=spat_flexure) # Choose the slits to use if slitidx is not None: @@ -701,7 +714,7 @@ def slit_img(self, pad=None, slitidx=None, initial=False, flexure=None, exclude_ return slitid_img def spatial_coordinate_image(self, slitidx=None, full=False, slitid_img=None, - pad=None, initial=False, flexure_shift=None): + pad=None, initial=False, spat_flexure=None): r""" Generate an image with the normalized spatial coordinate within each slit. @@ -738,6 +751,13 @@ def spatial_coordinate_image(self, slitidx=None, full=False, slitid_img=None, :attr:`right`) are used. To use the nominal edges regardless of the presence of the tweaked edges, set this to True. See :func:`select_edges`. + spat_flexure (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. Returns: `numpy.ndarray`_: Array specifying the spatial coordinate @@ -759,7 +779,7 @@ def spatial_coordinate_image(self, slitidx=None, full=False, slitid_img=None, msgs.error('Provided slit ID image does not have the correct shape!') # Choose the slit edges to use - left, right, _ = self.select_edges(initial=initial, flexure=flexure_shift) + left, right, _ = self.select_edges(initial=initial, spat_flexure=spat_flexure) # Slit width slitwidth = right - left @@ -785,7 +805,7 @@ def spatial_coordinate_image(self, slitidx=None, full=False, slitid_img=None, coo_img = coo return coo_img - def spatial_coordinates(self, initial=False, flexure=None): + def spatial_coordinates(self, initial=False, spat_flexure=None): """ Return a fiducial coordinate for each slit. @@ -793,26 +813,33 @@ def spatial_coordinates(self, initial=False, flexure=None): :func:`slit_spat_pos`. Args: - original (:obj:`bool`, optional): + initial (:obj:`bool`, optional): By default, the method will use the tweaked slit edges if they have been defined. If they haven't - been, the nominal edges (:attr:`left` and - :attr:`right`) are used. To use the nominal edges + been, the initial edges (:attr:`left` and + :attr:`right`) are used. To use the initial edges regardless of the presence of the tweaked edges, set this to True. See :func:`select_edges`. + spat_flexure (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. Returns: `numpy.ndarray`_: Vector with the list of floating point spatial coordinates. """ - # TODO -- Confirm it makes sense to pass in flexure - left, right, _ = self.select_edges(initial=initial, flexure=flexure) + # TODO -- Confirm it makes sense to pass in spatial flexure + left, right, _ = self.select_edges(initial=initial, spat_flexure=spat_flexure) return SlitTraceSet.slit_spat_pos(left, right, self.nspat) @staticmethod def slit_spat_pos(left, right, nspat): r""" - Return a fidicial, normalized spatial coordinate for each slit. + Return a fiducial, normalized spatial coordinate for each slit. The fiducial coordinates are given by:: @@ -847,8 +874,13 @@ def mask_add_missing_obj(self, sobjs, spat_flexure, fwhm, boxcar_rad): Args: sobjs (:class:`~pypeit.specobjs.SpecObjs`): List of SpecObj that have been found and traced - spat_flexure (:obj:`float`): - Shifts, in spatial pixels, between this image and SlitTrace + spat_flexure (`numpy.ndarray`_): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. fwhm (:obj:`float`): FWHM in pixels to be used in the optimal extraction boxcar_rad (:obj:`float`): @@ -880,9 +912,9 @@ def mask_add_missing_obj(self, sobjs, spat_flexure, fwhm, boxcar_rad): cut_sobjs = sobjs # get slits edges init; includes flexure - left_init, _, _ = self.select_edges(initial=True, flexure=spat_flexure) + left_init, _, _ = self.select_edges(initial=True, spat_flexure=spat_flexure) # get slits edges tweaked; includes flexure - left_tweak, right_tweak, _ = self.select_edges(initial=False, flexure=spat_flexure) + left_tweak, right_tweak, _ = self.select_edges(initial=False, spat_flexure=spat_flexure) # midpoint in the spectral direction specmid = left_init[:,0].size//2 @@ -981,17 +1013,27 @@ def mask_add_missing_obj(self, sobjs, spat_flexure, fwhm, boxcar_rad): # Return return sobjs - def assign_maskinfo(self, sobjs, plate_scale, spat_flexure, TOLER=1.): + def assign_maskinfo(self, sobjs, plate_scale, spat_flexure, tolerance=1.): """ Assign RA, DEC, Name to objects. Modified in place. Args: - sobjs (:class:`pypeit.specobjs.SpecObjs`): List of SpecObj that have been found and traced. - plate_scale (:obj:`float`): platescale for the current detector. - spat_flexure (:obj:`float`): Shifts, in spatial pixels, between this image and SlitTrace. - det_buffer (:obj:`int`): Minimum separation between detector edges and a slit edge. - TOLER (:obj:`float`, optional): Matching tolerance in arcsec. + sobjs (:class:`pypeit.specobjs.SpecObjs`): + List of SpecObj that have been found and traced. + plate_scale (:obj:`float`): + platescale for the current detector. + spat_flexure (`numpy.ndarray`_): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. + det_buffer (:obj:`int`): + Minimum separation between detector edges and a slit edge. + tolerance (:obj:`float`, optional): + Matching tolerance in arcsec. Returns: :class:`pypeit.specobjs.SpecObjs`: Updated list of SpecObj that have been found and traced. @@ -1031,7 +1073,7 @@ def assign_maskinfo(self, sobjs, plate_scale, spat_flexure, TOLER=1.): 'Matching tolerance includes user-provided tolerance, slit tracing uncertainties and object size.') # get slits edges init - left_init, right_init, _ = self.select_edges(initial=True, flexure=spat_flexure) # includes flexure + left_init, right_init, _ = self.select_edges(initial=True, spat_flexure=spat_flexure) # includes flexure # midpoint in the spectral direction specmid = left_init[:, 0].size // 2 @@ -1084,7 +1126,7 @@ def assign_maskinfo(self, sobjs, plate_scale, spat_flexure, TOLER=1.): obj_fwhm = cut_sobjs[ipeak].FWHM*plate_scale else: obj_fwhm = 0. - in_toler = np.abs(separ*plate_scale) < (TOLER + cc_rms + obj_fwhm/2) + in_toler = np.abs(separ*plate_scale) < (tolerance + cc_rms + obj_fwhm / 2) if np.any(in_toler): # Find positive peakflux peak_flux = cut_sobjs[idx].smash_peakflux[in_toler] @@ -1240,17 +1282,28 @@ def get_maskdef_offset(self, sobjs, platescale, spat_flexure, slitmask_off, brig Determine the Slitmask offset (pixels) from position expected by the slitmask design Args: - sobjs (:class:`pypeit.specobjs.SpecObjs`): List of SpecObj that have been found and traced - platescale (:obj:`float`): Platescale - spat_flexure (:obj:`float`): Shifts, in spatial pixels, between this image and SlitTrace - slitmask_off (:obj:`float`): User provided slitmask offset in pixels - bright_maskdefid (:obj:`str`): User provided maskdef_id of a bright object to be used to measure offset - snr_thrshd (:obj:`float`): Objects detected above this S/N ratio threshold will be use to - compute the slitmask offset - use_alignbox (:obj:`bool`): Flag that determines if the alignment boxes are used to measure the offset - dither_off (:obj:`float`, optional): dither offset recorded in the header of the observations - - + sobjs (:class:`pypeit.specobjs.SpecObjs`): + List of SpecObj that have been found and traced + platescale (:obj:`float`): + Platescale of this detector + spat_flexure (`numpy.ndarray`_): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. + slitmask_off (:obj:`float`): + User provided slitmask offset in pixels + bright_maskdefid (:obj:`str`): + User provided maskdef_id of a bright object to be used to measure offset + snr_thrshd (:obj:`float`): + Objects detected above this S/N ratio threshold will be use to + compute the slitmask offset + use_alignbox (:obj:`bool`): + Flag that determines if the alignment boxes are used to measure the offset + dither_off (:obj:`float`, optional): + dither offset recorded in the header of the observations """ if self.maskdef_objpos is None: msgs.error('An array of object positions predicted by the slitmask design must be provided.') @@ -1292,7 +1345,7 @@ def get_maskdef_offset(self, sobjs, platescale, spat_flexure, slitmask_off, brig align_maskdef_ids = obj_maskdef_id[flag_align == 1] # get slits edges init - left_init, _, _ = self.select_edges(initial=True, flexure=spat_flexure) # includes flexure + left_init, _, _ = self.select_edges(initial=True, spat_flexure=spat_flexure) # includes flexure # midpoint in the spectral direction specmid = left_init[:, 0].size // 2 @@ -1558,13 +1611,25 @@ def get_maskdef_objpos_offset_alldets(sobjs, calib_slits, spat_flexure, platesca This info is recorded in the `SlitTraceSet` datamodel. Args: - sobjs (:class:`pypeit.specobjs.SpecObjs`): List of SpecObj that have been found and traced - calib_slits (:obj:`list`): List of `SlitTraceSet` with information on the traced slit edges - spat_flexure (:obj:`list`): List of shifts, in spatial pixels, between this image and SlitTrace - platescale (:obj:`list`): List of platescale for every detector - det_buffer (:obj:`int`): Minimum separation between detector edges and a slit edge - slitmask_par (:class:`pypeit.par.pypeitpar.PypeItPar`): slitmask PypeIt parameters - dither_off (:obj:`float`, optional): dither offset recorded in the header of the observations + sobjs (:class:`pypeit.specobjs.SpecObjs`): + List of SpecObj that have been found and traced + calib_slits (:obj:`list`): + List of `SlitTraceSet` with information on the traced slit edges + spat_flexure (:obj:`list`, optional): + If provided, this is a list of the shifts, in spatial pixels, + to apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. + platescale (:obj:`list`): + List of platescale for every detector + det_buffer (:obj:`int`): + Minimum separation between detector edges and a slit edge + slitmask_par (:class:`pypeit.par.pypeitpar.PypeItPar`): + slitmask PypeIt parameters + dither_off (:obj:`float`, optional): + dither offset recorded in the header of the observations Returns: List of `SlitTraceSet` with updated information on the traced slit edges @@ -1693,7 +1758,10 @@ def assign_addobjs_alldets(sobjs, calib_slits, spat_flexure, platescale, slitmas calib_slits (`numpy.ndarray`_): Array of `SlitTraceSet` with information on the traced slit edges. spat_flexure (:obj:`list`): - List of shifts, in spatial pixels, between this image and SlitTrace. + List of spatial flexure shifts, in spatial pixels, between this image and SlitTrace. + Each element of this list should be a `numpy.ndarray`_ of shape (nslits, 2), + where the first column is the shift to apply to the left edge of each slit + and the second column is the shift to apply to the right edge of each slit. platescale (:obj:`list`): List of platescale for every detector. slitmask_par (:class:`~pypeit.par.pypeitpar.PypeItPar`): @@ -1714,7 +1782,7 @@ def assign_addobjs_alldets(sobjs, calib_slits, spat_flexure, platescale, slitmas if calib_slits[i].maskdef_designtab is not None: # Assign slitmask design information to detected objects sobjs = calib_slits[i].assign_maskinfo(sobjs, platescale[i], spat_flexure[i], - TOLER=slitmask_par['obj_toler']) + tolerance=slitmask_par['obj_toler']) if slitmask_par['extract_missing_objs']: # Set the FWHM for the extraction of missing objects diff --git a/pypeit/spec2dobj.py b/pypeit/spec2dobj.py index 5e693e44c0..c065ade575 100644 --- a/pypeit/spec2dobj.py +++ b/pypeit/spec2dobj.py @@ -48,7 +48,7 @@ class Spec2DObj(datamodel.DataContainer): Primary header if instantiated from a FITS file """ - version = '1.1.1' + version = '1.1.2' # TODO 2d data model should be expanded to include: # waveimage -- flexure and heliocentric corrections should be applied to the final waveimage and since this is unique to @@ -72,8 +72,8 @@ class Spec2DObj(datamodel.DataContainer): 'tilts': dict(otype=np.ndarray, atype=np.floating, descr='2D tilts image (float64)'), 'scaleimg': dict(otype=np.ndarray, atype=np.floating, - descr='2D multiplicative scale image [or a single scalar as an array] that has been applied to ' - 'the science image (float32)'), + descr='2D multiplicative scale image [or a single scalar as an array] ' + 'that has been applied to the science image (float32)'), 'waveimg': dict(otype=np.ndarray, atype=np.floating, descr='2D wavelength image in vacuum (float64)'), 'bpmmask': dict(otype=imagebitmask.ImageBitMaskArray, @@ -85,12 +85,15 @@ class Spec2DObj(datamodel.DataContainer): descr='Table with WaveCalib diagnostic info'), 'maskdef_designtab': dict(otype=table.Table, descr='Table with slitmask design and object info'), - 'sci_spat_flexure': dict(otype=float, - descr='Shift, in spatial pixels, between this image ' - 'and SlitTrace'), + 'sci_spat_flexure': dict(otype=np.ndarray, atype=np.floating, + descr='Shift, in spatial pixels, between this image ' + 'and SlitTrace. Shape is (nslits, 2), where ' + 'spat_flexure[i,0] is the spatial shift of the left ' + 'edge of slit i and spat_flexure[i,1] is the spatial ' + 'shift of the right edge of slit i.'), 'sci_spec_flexure': dict(otype=table.Table, - descr='Global shift of the spectrum to correct for spectral' - 'flexure (pixels). This is based on the sky spectrum at' + descr='Global shift of the spectrum to correct for spectral ' + 'flexure (pixels). This is based on the sky spectrum at ' 'the center of each slit'), 'vel_type': dict(otype=str, descr='Type of reference frame correction (if any). ' 'Options are listed in the routine: ' @@ -244,6 +247,9 @@ def _bundle(self): # maskdef_designtab elif key == 'maskdef_designtab': d.append(dict(maskdef_designtab=self.maskdef_designtab)) + # Spatial flexure + elif key == 'sci_spat_flexure': + d.append(dict(sci_spat_flexure=self.sci_spat_flexure)) # Spectral flexure elif key == 'sci_spec_flexure': d.append(dict(sci_spec_flexure=self.sci_spec_flexure)) @@ -331,7 +337,7 @@ def update_slits(self, spec2DObj): self.slits.mask[gpm] = spec2DObj.slits.mask[gpm] # Slitmask - slitmask = spec2DObj.slits.slit_img(flexure=spec2DObj.sci_spat_flexure, + slitmask = spec2DObj.slits.slit_img(spat_flexure=spec2DObj.sci_spat_flexure, exclude_flag=spec2DObj.slits.bitmask.exclude_for_reducing) # Fill in the image for slit_idx, spat_id in enumerate(spec2DObj.slits.spat_id[gpm]): @@ -574,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? @@ -661,7 +675,7 @@ def write_to_fits(self, outfile, pri_hdr=None, update_det=None, msgs.error("Original spec2D object has a different version. Too risky to continue. Rerun both") # Generate the slit "mask" slitmask = _allspecobj[det].slits.slit_img( - flexure=_allspecobj[det].sci_spat_flexure) + spat_flexure=_allspecobj[det].sci_spat_flexure) # Save the new one in a copy new_Spec2DObj = deepcopy(self[det]) # Replace with the old @@ -762,11 +776,21 @@ def flexure_diagnostics(self, flexure_type='spat'): return_flex[det] = spec_flex # get and print the spatial flexure if flexure_type == 'spat': - spat_flex = self[det].sci_spat_flexure - # print the value - print(f'Spat shift: {spat_flex}') + spat_flexure = self[det].sci_spat_flexure + if np.all(spat_flexure == spat_flexure[0, 0]): + # print the value + print(f'Spatial shift: {spat_flexure}') + elif np.array_equal(spat_flexure[:,0],spat_flexure[:,1]): + # print the value of each slit + for ii in range(spat_flexure.shape[0]): + print(f' Slit {ii+1} spatial shift: {spat_flexure[ii,0]}') + else: + # print the value for the edge of each slit + for ii in range(spat_flexure.shape[0]): + print(' Slit {0:2d} -- left edge spatial shift: {1:f}'.format(ii+1, spat_flexure[ii,0])) + print(' -- right edge spatial shift: {0:f}'.format(spat_flexure[ii,1])) # return the value - return_flex[det] = spat_flex + return_flex[det] = spat_flexure return return_flex diff --git a/pypeit/specobjs.py b/pypeit/specobjs.py index 624168addc..3e21d2e811 100644 --- a/pypeit/specobjs.py +++ b/pypeit/specobjs.py @@ -833,7 +833,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)): diff --git a/pypeit/spectrographs/gemini_gnirs.py b/pypeit/spectrographs/gemini_gnirs.py index c878a8393b..8b8301d632 100644 --- a/pypeit/spectrographs/gemini_gnirs.py +++ b/pypeit/spectrographs/gemini_gnirs.py @@ -640,7 +640,7 @@ def default_pypeit_par(cls): par['scienceframe']['process']['objlim'] = 1.5 par['scienceframe']['process']['use_illumflat'] = False # illumflat is applied when building the relative scale image in reduce.py, so should be applied to scienceframe too. par['scienceframe']['process']['use_specillum'] = False # apply relative spectral illumination - par['scienceframe']['process']['spat_flexure_correct'] = False # don't correct for spatial flexure - varying spatial illumination profile could throw this correction off. Also, there's no way to do astrometric correction if we can't correct for spatial flexure of the contbars frames + par['scienceframe']['process']['spat_flexure_method'] = "skip" # don't correct for spatial flexure - varying spatial illumination profile could throw this correction off. Also, there's no way to do astrometric correction if we can't correct for spatial flexure of the contbars frames par['scienceframe']['process']['use_biasimage'] = False par['scienceframe']['process']['use_darkimage'] = False par['calibrations']['flatfield']['slit_illum_finecorr'] = False diff --git a/pypeit/spectrographs/gtc_osiris.py b/pypeit/spectrographs/gtc_osiris.py index 1135d6b604..862ae7a209 100644 --- a/pypeit/spectrographs/gtc_osiris.py +++ b/pypeit/spectrographs/gtc_osiris.py @@ -602,7 +602,7 @@ def default_pypeit_par(cls): par['scienceframe']['process']['objlim'] = 1.5 par['scienceframe']['process']['use_illumflat'] = False # illumflat is applied when building the relative scale image in reduce.py, so should be applied to scienceframe too. par['scienceframe']['process']['use_specillum'] = False # apply relative spectral illumination - par['scienceframe']['process']['spat_flexure_correct'] = False # don't correct for spatial flexure - varying spatial illumination profile could throw this correction off. Also, there's no way to do astrometric correction if we can't correct for spatial flexure of the contbars frames + par['scienceframe']['process']['spat_flexure_method'] = "skip" # don't correct for spatial flexure - varying spatial illumination profile could throw this correction off. Also, there's no way to do astrometric correction if we can't correct for spatial flexure of the contbars frames par['scienceframe']['process']['use_biasimage'] = False par['scienceframe']['process']['use_darkimage'] = False par['calibrations']['flatfield']['slit_illum_finecorr'] = False diff --git a/pypeit/spectrographs/keck_kcwi.py b/pypeit/spectrographs/keck_kcwi.py index e377c31fec..7c91183e6c 100644 --- a/pypeit/spectrographs/keck_kcwi.py +++ b/pypeit/spectrographs/keck_kcwi.py @@ -326,7 +326,7 @@ def default_pypeit_par(cls): # Illumination corrections par['scienceframe']['process']['use_illumflat'] = True # illumflat is applied when building the relative scale image in reduce.py, so should be applied to scienceframe too. par['scienceframe']['process']['use_specillum'] = True # apply relative spectral illumination - par['scienceframe']['process']['spat_flexure_correct'] = False # don't correct for spatial flexure - varying spatial illumination profile could throw this correction off. Also, there's no way to do astrometric correction if we can't correct for spatial flexure of the contbars frames + par['scienceframe']['process']['spat_flexure_method'] = "skip" # don't correct for spatial flexure - varying spatial illumination profile could throw this correction off. Also, there's no way to do astrometric correction if we can't correct for spatial flexure of the contbars frames par['scienceframe']['process']['use_biasimage'] = True # Need to use bias frames for KCWI, because the bias level varies monotonically with spatial and spectral direction par['scienceframe']['process']['use_darkimage'] = False diff --git a/pypeit/spectrographs/keck_lris.py b/pypeit/spectrographs/keck_lris.py index fe16205599..1d37fa4abd 100644 --- a/pypeit/spectrographs/keck_lris.py +++ b/pypeit/spectrographs/keck_lris.py @@ -104,8 +104,8 @@ def default_pypeit_par(cls): # Always correct for spatial flexure on science images # TODO -- Decide whether to make the following defaults # May not want to do them for LongSlit - par['scienceframe']['process']['spat_flexure_correct'] = True - par['calibrations']['standardframe']['process']['spat_flexure_correct'] = True + par['scienceframe']['process']['spat_flexure_method'] = "detector" + par['calibrations']['standardframe']['process']['spat_flexure_method'] = "detector" par['scienceframe']['exprng'] = [61, None] diff --git a/pypeit/tests/test_spec2dobj.py b/pypeit/tests/test_spec2dobj.py index 2a639153fe..ae30d78b49 100644 --- a/pypeit/tests/test_spec2dobj.py +++ b/pypeit/tests/test_spec2dobj.py @@ -25,8 +25,9 @@ def init_dict(): sciimg = np.ones((1000,1000)).astype(float) # Slits - left = np.full((1000, 3), 2, dtype=float) - right = np.full((1000, 3), 8, dtype=float) + nslits = 3 + left = np.full((1000, nslits), 2, dtype=float) + right = np.full((1000, nslits), 8, dtype=float) left[:,1] = 15. right[:,1] = 21. left[:,2] = 25. @@ -37,6 +38,7 @@ def init_dict(): spec_flex_table = Table() spec_flex_table['spat_id'] = slits.spat_id spec_flex_table['sci_spec_flexure'] = np.zeros(left.shape[1]) + spat_flexure = np.full((nslits, 2), 3.5) # return dict(sciimg = sciimg, ivarraw = 0.1 * np.ones_like(sciimg), @@ -52,7 +54,7 @@ def init_dict(): maskdef_designtab=None, tilts=np.ones_like(sciimg).astype(float), #tilts=wavetilts.WaveTilts(**test_wavetilts.instant_dict), - sci_spat_flexure=3.5, + sci_spat_flexure=spat_flexure, sci_spec_flexure=spec_flex_table, vel_type='HELIOCENTRIC', vel_corr=1.0+1.0e-5) diff --git a/pypeit/wavecalib.py b/pypeit/wavecalib.py index 18600f8a38..b4e69d9581 100644 --- a/pypeit/wavecalib.py +++ b/pypeit/wavecalib.py @@ -270,18 +270,23 @@ def build_fwhmimg(self, tilts, slits, initial=False, spat_flexure=None): Properties of the slits initial (bool, optional): If True, the initial slit locations will be used. Otherwise, the tweaked edges will be used. - spat_flexure (float, optional): - Spatial flexure correction in pixels. + spat_flexure (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. Returns: `numpy.ndarray`_: The spectral FWHM image. """ # Check spatial flexure type - if (spat_flexure is not None) and (not isinstance(spat_flexure, float)): - msgs.error("Spatial flexure must be None or float") + if (spat_flexure is not None) and (not isinstance(spat_flexure, np.ndarray)): + msgs.error("Spatial flexure must be None or numpy.ndarray.") # Generate the slit mask and slit edges - pad slitmask by 1 for edge effects - slitmask = slits.slit_img(pad=1, initial=initial, flexure=spat_flexure) - slits_left, slits_right, _ = slits.select_edges(initial=initial, flexure=spat_flexure) + slitmask = slits.slit_img(pad=1, initial=initial, spat_flexure=spat_flexure) + slits_left, slits_right, _ = slits.select_edges(initial=initial, spat_flexure=spat_flexure) # need to exclude slits that are masked (are bad) bad_slits = slits.bitmask.flagged(slits.mask, and_not=slits.bitmask.exclude_for_reducing) ok_spat_ids = slits.spat_id[np.logical_not(bad_slits)] @@ -308,8 +313,13 @@ def build_waveimg(self, tilts, slits, spat_flexure=None, spec_flexure=None): Image holding tilts slits (:class:`pypeit.slittrace.SlitTraceSet`): Properties of the slits - spat_flexure (float, optional): - Spatial flexure correction in pixels. + spat_flexure (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, to + apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. spec_flexure (float, `numpy.ndarray`_, optional): Spectral flexure correction in pixels. If a float, the same spectral flexure correction will be applied @@ -322,8 +332,8 @@ def build_waveimg(self, tilts, slits, spat_flexure=None, spec_flexure=None): `numpy.ndarray`_: The wavelength image. """ # Check spatial flexure type - if (spat_flexure is not None) and (not isinstance(spat_flexure, float)): - msgs.error("Spatial flexure must be None or float") + if (spat_flexure is not None) and (not isinstance(spat_flexure, np.ndarray)): + msgs.error("Spatial flexure must be None or numpy.ndarray") # Check spectral flexure type if spec_flexure is None: spec_flex = np.zeros(slits.nslits) elif isinstance(spec_flexure, float): spec_flex = spec_flexure*np.ones(slits.nslits) @@ -338,7 +348,7 @@ def build_waveimg(self, tilts, slits, spat_flexure=None, spec_flexure=None): # image = np.zeros_like(tilts) # Grab slit_img - slitmask = slits.slit_img(flexure=spat_flexure, exclude_flag=slits.bitmask.exclude_for_reducing) + slitmask = slits.slit_img(spat_flexure=spat_flexure, exclude_flag=slits.bitmask.exclude_for_reducing) # Separate detectors for the 2D solutions? if self.par['ech_separate_2d']: @@ -554,7 +564,7 @@ def __init__(self, msarc, slits, spectrograph, par, lamps, # Load up slits # TODO -- Allow for flexure - slits_left, slits_right, mask = self.slits.select_edges(initial=True, flexure=None) # Grabs all, init slits + flexure + slits_left, slits_right, mask = self.slits.select_edges(initial=True, spat_flexure=None) # Grabs all, init slits + flexure self.orders = self.slits.ech_order # Can be None # self.spat_coo = self.slits.spatial_coordinates() # All slits, even masked # Internal mask for failed wv_calib analysis @@ -566,7 +576,7 @@ def __init__(self, msarc, slits, spectrograph, par, lamps, self.wvc_bpm_init = self.wvc_bpm.copy() # Slitmask -- Grabs only unmasked, initial slits #self.slitmask_science = self.slits.slit_img(initial=True, flexure=None, exclude_flag=['BOXSLIT']) - self.slitmask_science = self.slits.slit_img(initial=True, flexure=None) + self.slitmask_science = self.slits.slit_img(initial=True, spat_flexure=None) # Resize self.shape_science = self.slitmask_science.shape self.shape_arc = self.msarc.image.shape diff --git a/pypeit/wavetilts.py b/pypeit/wavetilts.py index b638cfc4c1..b75d7a5320 100644 --- a/pypeit/wavetilts.py +++ b/pypeit/wavetilts.py @@ -46,7 +46,7 @@ class WaveTilts(calibframe.CalibFrame): included in the output. """ - version = '1.2.0' + version = '1.2.1' # Calibration frame attributes calib_type = 'Tilts' @@ -72,7 +72,12 @@ class WaveTilts(calibframe.CalibFrame): 'spec_order': dict(otype=np.ndarray, atype=np.integer, descr='Order for spectral fit (nslit)'), 'func2d': dict(otype=str, descr='Function used for the 2D fit'), - 'spat_flexure': dict(otype=float, descr='Flexure shift from the input TiltImage'), + 'spat_flexure': dict(otype=np.ndarray, atype=np.floating, + descr='Spatial flexure shift, in spatial pixels, between TiltImage ' + 'and SlitTrace. Shape is (nslits, 2), where ' + 'spat_flexure[i,0] is the spatial shift of the left ' + 'edge of slit i and spat_flexure[i,1] is the spatial ' + 'shift of the right edge of slit i.'), 'slits_filename': dict(otype=str, descr='Path to SlitTraceSet file. This helps to ' 'find the Slits calibration file when running ' 'pypeit_chk_tilts()'), @@ -124,38 +129,52 @@ def is_synced(self, slits): msgs.error('Your tilt solutions are out of sync with your slits. Remove calibrations ' 'and restart from scratch.') - def fit2tiltimg(self, slitmask, flexure=None): + def fit2tiltimg(self, slitmask, slits_left, slits_right, spat_flexure=None): """ Generate a tilt image from the fit parameters - - Mainly to allow for flexure + Mainly to allow for spatial flexure Args: slitmask (`numpy.ndarray`_): - ?? - flexure (float, optional): - Spatial shift of the tilt image onto the desired frame - (typically a science image) + An image identifying the slit/order at each pixel. Pixels + without a slit are marked with -1. + slits_left (`numpy.ndarray`_): + Left slit edges + slits_right (`numpy.ndarray`_): + Right slit edges + spat_flexure (`numpy.ndarray`_, optional): + If provided, this is the shift, in spatial pixels, of the tilt + image onto the desired frame (typically a science image). The + shifts apply to each slit. This is used to correct for spatial + flexure. The shape of the array should be (nslits, 2), + where the first column is the shift to apply to the + left edge of each slit and the second column is the + shift to apply to the right edge of each slit. Returns: `numpy.ndarray`_: New tilt image - """ msgs.info("Generating a tilts image from the fit parameters") - _flexure = 0. if flexure is None else flexure + # Check the optional inputs + _spat_flexure = np.zeros((slits_left.shape[1], 2)) if spat_flexure is None else spat_flexure + # Setup the output image final_tilts = np.zeros_like(slitmask).astype(float) gdslit_spat = np.unique(slitmask[slitmask >= 0]).astype(int) - # Loop + # Loop through all good slits for slit_spat in gdslit_spat: + # Get the slit index slit_idx = self.spatid_to_zero(slit_spat) - # Calculate - coeff_out = self.coeffs[:self.spec_order[slit_idx]+1,:self.spat_order[slit_idx]+1,slit_idx] - _tilts = tracewave.fit2tilts(final_tilts.shape, coeff_out, self.func2d, spat_shift=-1*_flexure) - # Fill - thismask_science = slitmask == slit_spat - final_tilts[thismask_science] = _tilts[thismask_science] + # Prepare the coefficients + coeff_out = self.coeffs[:self.spec_order[slit_idx]+1, :self.spat_order[slit_idx]+1, slit_idx] + # Extract the spectral and spatial coordinates for this slit + thismask_science = (slitmask == slit_spat) + _spec_eval, _spat_eval = tracewave.fit2tilts_prepareSlit(slits_left[:, slit_idx], slits_right[:, slit_idx], + thismask_science, _spat_flexure[slit_idx, :]) + # Calculate the tilts + final_tilts[thismask_science] = tracewave.fit2tilts(coeff_out, self.func2d, + spec_eval=_spec_eval, spat_eval=_spat_eval) # Return return final_tilts @@ -166,10 +185,11 @@ def spatid_to_zero(self, spat_id): Args: spat_id (int): + Slit spat_id Returns: - int: - + int: index of slit corresponding to spat_id + TODO :: This code is a direct copy of the slits method (and only used in the function above. Should be consolidated... """ mtch = self.spat_id == spat_id return np.where(mtch)[0][0] @@ -217,8 +237,8 @@ def show(self, waveimg=None, wcs_match=True, in_ginga=True, show_traces=False, cal_file = Path(_calib_dir).absolute() / self.slits_filename if cal_file.exists(): slits = slittrace.SlitTraceSet.from_file(cal_file, chk_version=chk_version) - _slitmask = slits.slit_img(initial=True, flexure=self.spat_flexure) - _left, _right, _mask = slits.select_edges(flexure=self.spat_flexure) + _slitmask = slits.slit_img(initial=True, spat_flexure=self.spat_flexure) + _left, _right, _mask = slits.select_edges(spat_flexure=self.spat_flexure) gpm = _mask == 0 # resize slitmask = arc.resize_mask2arc(tilt_img_dict.image.shape, _slitmask) @@ -234,7 +254,7 @@ def show(self, waveimg=None, wcs_match=True, in_ginga=True, show_traces=False, wv_calib_name = wavecalib.WaveCalib.construct_file_name(self.calib_key, calib_dir=_calib_dir) if Path(wv_calib_name).absolute().exists(): wv_calib = wavecalib.WaveCalib.from_file(wv_calib_name, chk_version=chk_version) - tilts = self.fit2tiltimg(slitmask, flexure=self.spat_flexure) + tilts = self.fit2tiltimg(slitmask, _left, _right, spat_flexure=self.spat_flexure) waveimg = wv_calib.build_waveimg(tilts, slits, spat_flexure=self.spat_flexure) else: msgs.warn('Could not load Wave image to show with tilts image.') @@ -334,7 +354,7 @@ def __init__(self, mstilt, slits, spectrograph, par, wavepar, det=1, qa_path=Non self.slits = slits self.det = det self.qa_path = qa_path - self.spat_flexure = spat_flexure + self.spat_flexure = spat_flexure if spat_flexure is not None else np.zeros((slits.nslits, 2), dtype=float) self.measured_fwhms = measured_fwhms if measured_fwhms is not None else np.array([None] * slits.nslits) # Get the non-linear count level @@ -352,7 +372,7 @@ def __init__(self, mstilt, slits, spectrograph, par, wavepar, det=1, qa_path=Non # TODO -- Tidy this up into one or two methods? # Load up all slits # TODO -- Discuss further with JFH - all_left, all_right, mask = self.slits.select_edges(initial=True, flexure=self.spat_flexure) # Grabs all, initial slits + all_left, all_right, mask = self.slits.select_edges(initial=True, spat_flexure=self.spat_flexure) # Grabs all, initial slits # self.tilt_bpm = np.invert(mask == 0) # At this point of the reduction the only bitmask flags that may have been generated are 'USERIGNORE', # 'SHORTSLIT', 'BOXSLIT' and 'BADWVCALIB'. Here we use only 'USERIGNORE' and 'SHORTSLIT' to create the bpm mask @@ -360,7 +380,7 @@ def __init__(self, mstilt, slits, spectrograph, par, wavepar, det=1, qa_path=Non self.tilt_bpm_init = self.tilt_bpm.copy() # Slitmask # TODO -- Discuss further with JFH - self.slitmask_science = self.slits.slit_img(initial=True, flexure=self.spat_flexure, exclude_flag=['BOXSLIT']) # All unmasked slits + self.slitmask_science = self.slits.slit_img(initial=True, spat_flexure=self.spat_flexure, exclude_flag=['BOXSLIT']) # All unmasked slits # Resize # TODO: Should this be the bpm or *any* flag? gpm = self.mstilt.select_flag(flag='BPM', invert=True) if self.mstilt is not None \ @@ -375,7 +395,6 @@ def __init__(self, mstilt, slits, spectrograph, par, wavepar, det=1, qa_path=Non # Key Internals self.mask = None self.all_trace_dict = [None]*self.slits.nslits - self.tilts = None # 2D fits are stored as a dictionary rather than list because we will jsonify the dict self.all_fit_dict = [None]*self.slits.nslits self.steps = [] @@ -738,6 +757,9 @@ def run(self, doqa=True, debug=False, show=False): self.spat_order = np.zeros(self.slits.nslits, dtype=int) self.spec_order = np.zeros(self.slits.nslits, dtype=int) + # Grab the slit edges + slits_left, slits_right, _ = self.slits.select_edges(initial=True, spat_flexure=self.spat_flexure) + # Loop on all slits for slit_idx, slit_spat in enumerate(self.slits.spat_id): if self.tilt_bpm[slit_idx]: @@ -812,17 +834,18 @@ def run(self, doqa=True, debug=False, show=False): # Tilts are created with the size of the original slitmask, # which corresonds to the same binning as the science # images, trace images, and pixelflats etc. - self.tilts = tracewave.fit2tilts(self.slitmask_science.shape, coeff_out, - self.par['func2d']) + thismask_science = self.slitmask_science == slit_spat + _spec_eval, _spat_eval = tracewave.fit2tilts_prepareSlit(slits_left[:, slit_idx], slits_right[:, slit_idx], + thismask_science, self.spat_flexure[slit_idx, :]) + tilts = tracewave.fit2tilts(coeff_out, self.par['func2d'], spec_eval=_spec_eval, spat_eval=_spat_eval) # Check that the tilts image has values that span a reasonable range # TODO: Is this the right threshold? - if np.nanmax(self.tilts) - np.nanmin(self.tilts) < 0.8: + if np.nanmax(tilts) - np.nanmin(tilts) < 0.8: msgs.warn('Tilts image fit not good. This slit/order will not be reduced!') self.slits.mask[slit_idx] = self.slits.bitmask.turn_on(self.slits.mask[slit_idx], 'BADTILTCALIB') continue # Save to final image - thismask_science = self.slitmask_science == slit_spat - self.final_tilts[thismask_science] = self.tilts[thismask_science] + self.final_tilts[thismask_science] = tilts if show: viewer, ch = display.show_image(self.mstilt.image * (self.slitmask > -1), chname='tilts')