From 85effbc1cb43a7c8b1d7aa2d0102043f6d42c643 Mon Sep 17 00:00:00 2001 From: "Timothy P. Ellsworth Bowers" Date: Fri, 23 Jun 2023 11:35:22 -0700 Subject: [PATCH 1/8] Add sensible error message When users try to ingest a `spec1d` file into a Spectrum1D object, the current scheme ends badly with an obscure error message like: ``` OSError: Could not identify column containing the wavelength, frequency or energy ``` Because PypeIt `spec1d` files may have multiple extensions, by convention we load these into SpectrumList objects in specutils, even if there is only a single spectrum. This commit catches the case when a user tries to load a `spec1d` file into a Spectrum1D object. Previously, such an attempt ended up with specutils loading in the `spec1d` file as a tabular FITS format, then didn't know what to do about the spectral axis (leading to the above obscure error). This commit specifically includes a loader for `spec1d` files into Spectrum1D objects, but emits a sensible error message in this case rather than trying to actually load the file. As it happens, in order to properly do this, I needed to make an adjustment in specutils itself. A PR has been submitted (astropy/specutils#1068), but this commit should not be considered as a PypeIt PR until the specutils change has incorporated into an official release. At that time, the specutils optional dependency will need to have a minimum version specified. So... this commit will quietly sit in a branch until specutils makes its move. modified: pypeit/specutils/pypeit_loaders.py --- pypeit/specutils/pypeit_loaders.py | 42 +++++++++++++++++++++++++++--- 1 file changed, 39 insertions(+), 3 deletions(-) diff --git a/pypeit/specutils/pypeit_loaders.py b/pypeit/specutils/pypeit_loaders.py index b7cd629c99..e051ef81b3 100644 --- a/pypeit/specutils/pypeit_loaders.py +++ b/pypeit/specutils/pypeit_loaders.py @@ -11,6 +11,7 @@ - 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) """ @@ -80,7 +81,8 @@ def identify_pypeit_onespec(origin, *args, **kwargs): identifier=identify_pypeit_spec1d, extensions=["fits"], priority=10, - dtype=SpectrumList) + dtype=SpectrumList, + autogenerate_spectrumlist=False) def pypeit_spec1d_loader(filename, extract=None, fluxed=True, **kwargs): """ Load spectra from a PypeIt spec1d file into a SpectrumList. @@ -109,7 +111,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 ' @@ -160,7 +162,7 @@ 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 ' @@ -188,5 +190,39 @@ def pypeit_onespec_loader(filename, grid=False, **kwargs): 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.') From a9cc26ac5d2dfc202ef3ff3104eba44f76b06592 Mon Sep 17 00:00:00 2001 From: "Timothy P. Ellsworth Bowers" Date: Mon, 28 Aug 2023 12:00:56 -0700 Subject: [PATCH 2/8] Update specutils version to accomodate this PR modified: setup.cfg --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index bbb6f42214..ca496fe4d9 100644 --- a/setup.cfg +++ b/setup.cfg @@ -59,7 +59,7 @@ scripts = scikit-image = scikit-image specutils = - specutils + specutils>=1.12 test = pytest>=6.0.0 pytest-astropy From 1b6f0acd3c9de3af0bfda6ec1521e6b141f5af1c Mon Sep 17 00:00:00 2001 From: Kyle Westfall Date: Thu, 19 Oct 2023 12:59:41 -0700 Subject: [PATCH 3/8] monotonic check --- pypeit/specutils/pypeit_loaders.py | 83 +++++++++++++++++++++++++++--- 1 file changed, 77 insertions(+), 6 deletions(-) diff --git a/pypeit/specutils/pypeit_loaders.py b/pypeit/specutils/pypeit_loaders.py index 3cd8cee4c8..0408d31185 100644 --- a/pypeit/specutils/pypeit_loaders.py +++ b/pypeit/specutils/pypeit_loaders.py @@ -16,6 +16,8 @@ from IPython import embed +import numpy as np + import astropy.io.fits import astropy.nddata import astropy.units @@ -35,6 +37,57 @@ 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 not np.any(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!') + indx = np.append([True], np.logical_not(indx)) + return wave[indx], flux[indx], None if ivar is None else ivar[indx] + + # Identifier Functions =======================================================# def _identify_pypeit(*args, **kwargs): """ @@ -81,7 +134,7 @@ def identify_pypeit_onespec(origin, *args, **kwargs): extensions=["fits"], priority=10, dtype=SpectrumList) -def pypeit_spec1d_loader(filename, extract=None, fluxed=True, **kwargs): +def pypeit_spec1d_loader(filename, extract=None, fluxed=True, strict=True, **kwargs): """ Load spectra from a PypeIt spec1d file into a SpectrumList. @@ -98,6 +151,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 ------- @@ -121,7 +182,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 +199,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 +210,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 ------- @@ -168,6 +237,8 @@ def pypeit_onespec_loader(filename, grid=False, **kwargs): 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,9 +250,9 @@ 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), From b02e23c76fa844ba31a299dd6526a0cf4a292a35 Mon Sep 17 00:00:00 2001 From: Kyle Westfall Date: Thu, 19 Oct 2023 13:25:15 -0700 Subject: [PATCH 4/8] need to iteratively remove points --- pypeit/specutils/pypeit_loaders.py | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/pypeit/specutils/pypeit_loaders.py b/pypeit/specutils/pypeit_loaders.py index 0408d31185..279f167104 100644 --- a/pypeit/specutils/pypeit_loaders.py +++ b/pypeit/specutils/pypeit_loaders.py @@ -67,8 +67,8 @@ def enforce_monotonic_wavelengths(wave, flux, ivar, strict=True): Edited inverse variance vector. This may be an unchanged reference to the original vector. """ - indx = np.diff(wave) < 0 - if not np.any(indx): + indx = np.diff(wave) > 0 + if np.all(indx): # Wavelengths are monotonic, so we're done return wave, flux, ivar @@ -84,8 +84,20 @@ def enforce_monotonic_wavelengths(wave, flux, ivar, strict=True): '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!') - indx = np.append([True], np.logical_not(indx)) - return wave[indx], flux[indx], None if ivar is None else ivar[indx] + + # 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. + _wave = wave.copy() + 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 =======================================================# From 33a152039de41919396d11fa2235e7096275ea4f Mon Sep 17 00:00:00 2001 From: Kyle Westfall Date: Thu, 19 Oct 2023 14:12:15 -0700 Subject: [PATCH 5/8] changes --- doc/releases/1.14.1dev.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/doc/releases/1.14.1dev.rst b/doc/releases/1.14.1dev.rst index 8c59aae6af..efce08cc5a 100644 --- a/doc/releases/1.14.1dev.rst +++ b/doc/releases/1.14.1dev.rst @@ -16,6 +16,10 @@ 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. Instrument-specific Updates --------------------------- From 4379d29a2fb40af15a633f240cafada5e2186774 Mon Sep 17 00:00:00 2001 From: Kyle Westfall Date: Thu, 19 Oct 2023 14:26:49 -0700 Subject: [PATCH 6/8] test --- pypeit/tests/test_specutils.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/pypeit/tests/test_specutils.py b/pypeit/tests/test_specutils.py index 87308f1da1..ae3b6f5647 100644 --- a/pypeit/tests/test_specutils.py +++ b/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() + + From 795230d6609e483241de8e3b81ca544ae77d8d60 Mon Sep 17 00:00:00 2001 From: "Timothy P. Ellsworth Bowers" Date: Thu, 19 Oct 2023 22:52:54 -0700 Subject: [PATCH 7/8] Update release notes modified: doc/releases/1.14.1dev.rst modified: setup.cfg --- doc/releases/1.14.1dev.rst | 4 ++++ setup.cfg | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/doc/releases/1.14.1dev.rst b/doc/releases/1.14.1dev.rst index efce08cc5a..8f3c998841 100644 --- a/doc/releases/1.14.1dev.rst +++ b/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 ---------------------------------------------------- @@ -20,6 +21,9 @@ Functionality/Performance Improvements and Additions 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 --------------------------- diff --git a/setup.cfg b/setup.cfg index d7e0b56625..e46e893c82 100644 --- a/setup.cfg +++ b/setup.cfg @@ -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 From e00758a6c6c344c7f0aa93785693782d54ba8a2d Mon Sep 17 00:00:00 2001 From: Kyle Westfall Date: Mon, 23 Oct 2023 09:39:19 -0700 Subject: [PATCH 8/8] PR comments --- pypeit/specutils/pypeit_loaders.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pypeit/specutils/pypeit_loaders.py b/pypeit/specutils/pypeit_loaders.py index aceea013c9..f5d8e6e3d5 100644 --- a/pypeit/specutils/pypeit_loaders.py +++ b/pypeit/specutils/pypeit_loaders.py @@ -38,7 +38,7 @@ from pypeit import onespec -def enforce_monotonic_wavelengths(wave, flux, ivar, strict=True): +def _enforce_monotonic_wavelengths(wave, flux, ivar, strict=True): """ Force the spectrum to have a monotonically increasing wavelength vector. @@ -90,7 +90,6 @@ def enforce_monotonic_wavelengths(wave, flux, ivar, strict=True): # want to be more acceptable, we should consider instead fitting a low-order # polynomial to the pixel vs. wavelength function and rejecting strong # outliers. - _wave = wave.copy() pix = np.arange(wave.size) indx = np.append([True], indx) while not np.all(indx): @@ -196,7 +195,7 @@ def pypeit_spec1d_loader(filename, extract=None, fluxed=True, strict=True, **kwa _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) + _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), @@ -251,7 +250,7 @@ def pypeit_onespec_loader(filename, grid=False, strict=True, **kwargs): 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) + 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