Merge pull request #1708 from pypeit/specutil_wave

Allow a way to circumvent monotonic wavelength requirement in Spectrum1D

doc/releases/1.14.1dev.rst

@@ -6,6 +6,7 @@ Dependency Changes
 ------------------
 
 - Removes use of deprecated ``pkg_resources``.
+- Require version ``>=1.12`` for ``specutils``.
 
 Functionality/Performance Improvements and Additions
 ----------------------------------------------------
@@ -16,6 +17,13 @@ Functionality/Performance Improvements and Additions
 - Added support for Keck/KCRM RL data reduction.
 - Allow one to turn off reinit of reduce_bpm in global_skysub and made this
   the default for the final pass
+- Allow pypeit Spectrum1D loaders to circumvent the requirement that reduced
+  spectrum have monotonically increasing wavelength data.  Non-monotonic
+  wavelength data usually indicate a problem with the data reduction, but this
+  at least lets users ingest the spectrum.
+- Add a sensible error message to the pypeit Spectrum1D loaders in the event a
+  user inadvertently tries to use Spectrum1D instead of SpectrumList for a
+  ``spec1d`` file.
 
 Instrument-specific Updates
 ---------------------------

pypeit/specutils/pypeit_loaders.py

@@ -11,11 +11,14 @@
 - 2022-02-04: Initial Version (tbowers)
 - 2022-09-16: Correct an import error and add module docstring (tbowers)
 - 2023-03-09: Moved into the main pypeit repo and refactored (KBW)
+- 2023-06-23: Added sensible error message for incorrect spec1d loading (tbowers)
 
 """
 
 from IPython import embed
 
+import numpy as np
+
 import astropy.io.fits
 import astropy.nddata
 import astropy.units
@@ -35,6 +38,68 @@
 from pypeit import onespec
 
 
+def _enforce_monotonic_wavelengths(wave, flux, ivar, strict=True):
+    """
+    Force the spectrum to have a monotonically increasing wavelength vector.
+
+    Parameters
+    ----------
+    wave : `numpy.ndarray`_
+        Spectrum wavelengths
+    flux : `numpy.ndarray`_
+        Spectrum flux
+    ivar : `numpy.ndarray`_
+        Spectrum inverse variance. Can be None.
+    strict : bool, optional
+        Check that the wavelength vector is monotonically increasing.  If not,
+        raise an error (as would be done by the `specutils.SpectrumList`_ class).
+        If False, wavelengths that are *not* monotonically increasing are masked
+        in the construction of the returned `specutils.SpectrumList`_ object.
+
+    Returns
+    -------
+    wave : `numpy.ndarray`_
+        Edited wavelength vector.  This may be an unchanged reference to the
+        original vector.
+    flux : `numpy.ndarray`_
+        Edited flux vector.  This may be an unchanged reference to the original
+        vector.
+    ivar : `numpy.ndarray`_
+        Edited inverse variance vector.  This may be an unchanged reference to
+        the original vector.
+    """
+    indx = np.diff(wave) > 0
+    if np.all(indx):
+        # Wavelengths are monotonic, so we're done
+        return wave, flux, ivar
+
+    if strict:
+        # Wavelengths are not monotonic, but the user expects them to be, so
+        # fault.
+        msgs.error('Wavelengths are not monotonically increasing!  Circumvent this fault by '
+                   'setting strict=False, but BEWARE that this is likely the result of an '
+                   'error in the data reduction!')
+
+    # Wavelengths are not monotonic, but the user wants to keep going.
+    msgs.warn('Wavelengths are not monotonically increasing!  Strict was set to False, so '
+              'measurements after a negative step in wavelength are removed from the constructed '
+              'spectrum.  BEWARE that this is likely the result of an error in the data '
+              'reduction!')
+
+    # NOTE: This is the brute force approach.  If this becomes something that we
+    # want to be more acceptable, we should consider instead fitting a low-order
+    # polynomial to the pixel vs. wavelength function and rejecting strong
+    # outliers.
+    pix = np.arange(wave.size)
+    indx = np.append([True], indx)
+    while not np.all(indx):
+        pix = pix[indx]
+        wave = wave[indx]
+        indx = np.append([True], np.diff(wave) > 0)
+
+    return wave, flux[pix], None if ivar is None else ivar[pix]
+
+
 # Identifier Functions =======================================================#
 def _identify_pypeit(*args, **kwargs):
     """
@@ -80,8 +145,9 @@ def identify_pypeit_onespec(origin, *args, **kwargs):
              identifier=identify_pypeit_spec1d,
              extensions=["fits"],
              priority=10,
-             dtype=SpectrumList)
-def pypeit_spec1d_loader(filename, extract=None, fluxed=True, **kwargs):
+             dtype=SpectrumList,
+             autogenerate_spectrumlist=False)
+def pypeit_spec1d_loader(filename, extract=None, fluxed=True, strict=True, **kwargs):
     """
     Load spectra from a PypeIt spec1d file into a SpectrumList.
 
@@ -98,6 +164,14 @@ def pypeit_spec1d_loader(filename, extract=None, fluxed=True, **kwargs):
         If True, return the flux-calibrated spectrum, if it exists.  If the flux
         calibration hasn't been performed or ``fluxed=False``, the spectrum is
         returned in counts.
+    strict : bool, optional
+        Check that the wavelength vector is monotonically increasing.  If not,
+        raise an error (as would be done by the `specutils.SpectrumList`_ class).
+        If False, wavelengths that are *not* monotonically increasing are masked
+        in the construction of the returned `specutils.SpectrumList`_ object.
+    kwargs : dict, optional
+        **Ignored!**  Used to catch spurious arguments passed to the base class
+        *that are ignored by this function.
 
     Returns
     -------
@@ -109,7 +183,7 @@ def pypeit_spec1d_loader(filename, extract=None, fluxed=True, **kwargs):
         sobjs = specobjs.SpecObjs.from_fitsfile(filename, chk_version=False)
     except PypeItError:
         file_pypeit_version = astropy.io.fits.getval(filename, 'VERSPYP', 'PRIMARY')
-        msgs.error(f'Unable to ingest {filename} using pypeit.specobjs module from your version '
+        msgs.error(f'Unable to ingest {filename.name} using pypeit.specobjs module from your version '
                    f'of PypeIt ({__version__}).  The version used to write the file is '
                    f'{file_pypeit_version}.  If these are different, you may need to re-reduce '
                    'your data using your current PypeIt version or install the matching version '
@@ -121,7 +195,7 @@ def pypeit_spec1d_loader(filename, extract=None, fluxed=True, **kwargs):
         _ext, _cal = sobj.best_ext_match(extract=extract, fluxed=fluxed)
         _wave, _flux, _ivar, _gpm = sobj.get_box_ext(fluxed=_cal) if _ext == 'BOX' \
                                         else sobj.get_opt_ext(fluxed=_cal)
-
+        _wave, _flux, _ivar = _enforce_monotonic_wavelengths(_wave, _flux, _ivar, strict=strict)
         flux_unit = astropy.units.Unit("1e-17 erg/(s cm^2 Angstrom)" if _cal else "electron")
         spec += \
             [Spectrum1D(flux=astropy.units.Quantity(_flux * flux_unit),
@@ -138,7 +212,7 @@ def pypeit_spec1d_loader(filename, extract=None, fluxed=True, **kwargs):
              extensions=["fits"],
              priority=10,
              dtype=Spectrum1D)
-def pypeit_onespec_loader(filename, grid=False, **kwargs):
+def pypeit_onespec_loader(filename, grid=False, strict=True, **kwargs):
     """
     Load a spectrum from a PypeIt OneSpec file into a Spectrum1D object.
 
@@ -149,6 +223,14 @@ def pypeit_onespec_loader(filename, grid=False, **kwargs):
     grid : bool, optional
         Use the uniform grid wavelengths, instead of the contribution-weighted
         center.
+    strict : bool, optional
+        Check that the wavelength vector is monotonically increasing.  If not,
+        raise an error (as would be done by the `specutils.Spectrum1D`_ class).
+        If False, wavelengths that are *not* monotonically increasing are masked
+        in the construction of the returned `specutils.Spectrum1D`_ object.
+    kwargs : dict, optional
+        **Ignored!**  Used to catch spurious arguments passed to the base class
+        *that are ignored by this function.
 
     Returns
     -------
@@ -160,14 +242,16 @@ def pypeit_onespec_loader(filename, grid=False, **kwargs):
         spec = onespec.OneSpec.from_file(filename)
     except PypeItError:
         file_pypeit_version = astropy.io.fits.getval(filename, 'VERSPYP', 'PRIMARY')
-        msgs.error(f'Unable to ingest {filename} using pypeit.specobjs module from your version '
+        msgs.error(f'Unable to ingest {filename.name} using pypeit.specobjs module from your version '
                    f'of PypeIt ({__version__}).  The version used to write the file is '
                    f'{file_pypeit_version}.  If these are different, you may need to re-reduce '
                    'your data using your current PypeIt version or install the matching version '
                    'of PypeIt (e.g., pip install pypeit==1.11.0).')
 
     flux_unit = astropy.units.Unit("1e-17 erg/(s cm^2 Angstrom)" if spec.fluxed else "ct/s")
     wave = spec.wave_grid_mid if grid else spec.wave
+    wave, flux, ivar = _enforce_monotonic_wavelengths(wave, spec.flux, spec.ivar, strict=strict)
+
     # If the input filename is actually a string, assign it as the spectrum
     # name.  Otherwise, try assuming it's a _io.FileIO object, and if that
     # doesn't work assign an empty string as the name.
@@ -179,14 +263,48 @@ def pypeit_onespec_loader(filename, grid=False, **kwargs):
         except AttributeError:
             name = ''
 
-    return Spectrum1D(flux=astropy.units.Quantity(spec.flux * flux_unit),
-                      uncertainty=None if spec.ivar is None 
-                        else astropy.nddata.InverseVariance(spec.ivar / flux_unit**2),
+    return Spectrum1D(flux=astropy.units.Quantity(flux * flux_unit),
+                      uncertainty=None if ivar is None 
+                        else astropy.nddata.InverseVariance(ivar / flux_unit**2),
                       meta={'name': name, 'extract': spec.ext_mode, 'fluxed': spec.fluxed,
                             'grid': grid},
                       spectral_axis=astropy.units.Quantity(wave * astropy.units.angstrom),
                       velocity_convention="doppler_optical",
                       bin_specification="centers")
 
+# Warning Function ===========================================================#
+@data_loader('PypeIt spec1d nolist',
+             identifier=identify_pypeit_spec1d,
+             extensions=["fits"],
+             priority=10,
+             dtype=Spectrum1D,
+             autogenerate_spectrumlist=False)
+def pypeit_spec1d_loader_nolist(filename, extract=None, fluxed=True, **kwargs):
+    """
+    Sensible error message if a user tries to load spectra from a PypeIt spec1d
+    file into a Spectrum1D.
 
+    This is not allowed because spec1d files may contain mutliple spectra.  This
+    function accepts all arguments as the SpectrumList version, but only outputs
+    a PypeIt Error with a sensible message.
 
+    This avoids receiving unhelpful error messages such as::
+
+        OSError: Could not identify column containing the wavelength, frequency or energy
+
+    Parameters
+    ----------
+    filename : str
+        The path to the FITS file
+    extract : str, optional
+        The extraction used to produce the spectrum.  Must be either None,
+        ``'BOX'`` (for a boxcar extraction), or ``'OPT'`` for optimal
+        extraction.  If None, the optimal extraction will be returned, if it
+        exists, otherwise the boxcar extraction will be returned.
+    fluxed : bool, optional
+        If True, return the flux-calibrated spectrum, if it exists.  If the flux
+        calibration hasn't been performed or ``fluxed=False``, the spectrum is
+        returned in counts.
+    """
+    msgs.error(f'The spec1d file {filename.name} cannot be ingested into a Spectrum1D object.'
+               f'{msgs.newline()}Please use the SpectrumList object for spec1d files.')

pypeit/tests/test_specutils.py

@@ -14,11 +14,13 @@
 from pypeit import specobjs
 from pypeit.specutils import Spectrum1D, SpectrumList
 from pypeit.tests import tstutils
+from pypeit.pypmsgs import PypeItError
 
 import pytest
 specutils_required = pytest.mark.skipif(Spectrum1D is None or SpectrumList is None,
                                         reason='specutils not installed')
 
+
 @specutils_required
 def test_onespec_io():
     rng = np.random.default_rng()
@@ -134,3 +136,32 @@ def test_spec1d_io():
 
     ofile.unlink()
 
+
+@specutils_required
+def test_onespec_monotonic():
+    rng = np.random.default_rng(999)
+    grid_wave = np.linspace(3500,10000,1000)
+    wave = grid_wave + (10*rng.uniform(size=grid_wave.size) - 1)
+    flux = np.ones(1000, dtype=float)
+    # TODO: PYP_SPEC is required if we want to be able to read the file!
+    spec = onespec.OneSpec(wave, grid_wave, flux, PYP_SPEC='shane_kast_blue')
+    ofile = Path(tstutils.data_path('tmp.fits')).resolve()
+    spec.to_file(str(ofile), overwrite=True)
+
+    with pytest.raises(PypeItError):
+        # Should fault because the wavelength vector is not monotonically
+        # increasing
+        _spec = Spectrum1D.read(ofile)
+
+    # This will be fine because the grid *is* monotonically increasing
+    _spec = Spectrum1D.read(ofile, grid=True)
+
+    # This should be fine because reader will remove non-monotonically
+    # increasing wavelength measurements.
+    __spec = Spectrum1D.read(ofile, strict=False)
+
+    assert _spec.shape[0] > __spec.shape[0], 'Strict should remove data'
+
+    ofile.unlink()
+
+

setup.cfg

@@ -61,7 +61,7 @@ scripts =
 scikit-image =
     scikit-image
 specutils =
-    specutils
+    specutils>=1.12
 test =
     pytest>=6.0.0
     pytest-astropy
@@ -75,7 +75,7 @@ docs =
     sphinx_rtd_theme==1.2.2
 dev =
     scikit-image
-    specutils
+    specutils>=1.12
     pytest>=6.0.0
     pytest-astropy
     tox