Merge pull request #818 from havok2063/jwstc1d

Adds support for JWST c1d spectroscopic files

Author: eteq

specutils/io/default_loaders/jwst_reader.py

@@ -3,6 +3,7 @@
 from astropy.table import Table
 from astropy.io import fits
 from astropy.nddata import StdDevUncertainty
+import logging
 import numpy as np
 import asdf
 from gwcs.wcstools import grid_from_bounding_box
@@ -14,23 +15,51 @@
 
 __all__ = ["jwst_x1d_single_loader", "jwst_x1d_multi_loader"]
 
+log = logging.getLogger(__name__)
 
-def identify_jwst_x1d_fits(origin, *args, **kwargs):
+
+def _identify_spec1d_fits(origin, extname, *args, **kwargs):
+    """ Generic spec 1d identifier function """
+    is_jwst = _identify_jwst_fits(*args)
+    with read_fileobj_or_hdulist(*args, memmap=False, **kwargs) as hdulist:
+        return (is_jwst and extname in hdulist and (extname, 2) not in hdulist)
+
+
+def _identify_spec1d_multi_fits(origin, extname, *args, **kwargs):
     """
-    Check whether the given file is a JWST x1d spectral data product.
+    Check whether the given file is a JWST c1d/x1d spectral data product with many slits.
     """
     is_jwst = _identify_jwst_fits(*args)
     with read_fileobj_or_hdulist(*args, memmap=False, **kwargs) as hdulist:
-        return (is_jwst and 'EXTRACT1D' in hdulist and ('EXTRACT1D', 2) not in hdulist)
+        return is_jwst and (extname, 2) in hdulist
+
+
+def identify_jwst_c1d_fits(origin, *args, **kwargs):
+    """
+    Check whether the given file is a JWST c1d spectral data product.
+    """
+    return _identify_spec1d_fits(origin, 'COMBINE1D', *args, **kwargs)
+
+
+def identify_jwst_c1d_multi_fits(origin, *args, **kwargs):
+    """
+    Check whether the given file is a JWST c1d spectral data product with many slits.
+    """
+    return _identify_spec1d_multi_fits(origin, 'COMBINE1D', *args, **kwargs)
+
+
+def identify_jwst_x1d_fits(origin, *args, **kwargs):
+    """
+    Check whether the given file is a JWST x1d spectral data product.
+    """
+    return _identify_spec1d_fits(origin, 'EXTRACT1D', *args, **kwargs)
 
 
 def identify_jwst_x1d_multi_fits(origin, *args, **kwargs):
     """
     Check whether the given file is a JWST x1d spectral data product with many slits.
     """
-    is_jwst = _identify_jwst_fits(*args)
-    with read_fileobj_or_hdulist(*args, memmap=False, **kwargs) as hdulist:
-        return is_jwst and ('EXTRACT1D', 2) in hdulist
+    return _identify_spec1d_multi_fits(origin, 'EXTRACT1D', *args, **kwargs)
 
 
 def identify_jwst_s2d_fits(origin, *args, **kwargs):
@@ -75,6 +104,50 @@ def _identify_jwst_fits(*args):
         return False
 
 
+@data_loader("JWST c1d", identifier=identify_jwst_c1d_fits, dtype=Spectrum1D,
+             extensions=['fits'])
+def jwst_c1d_single_loader(file_obj, **kwargs):
+    """
+    Loader for JWST c1d 1-D spectral data in FITS format
+
+    Parameters
+    ----------
+    filename : str
+        The path to the FITS file
+
+    Returns
+    -------
+    Spectrum1D
+        The spectrum contained in the file.
+    """
+    spectrum_list = _jwst_spec1d_loader(file_obj, extname='COMBINE1D', **kwargs)
+    if len(spectrum_list) == 1:
+        return spectrum_list[0]
+    else:
+        raise RuntimeError(f"Input data has {len(spectrum_list)} spectra. "
+                           "Use SpectrumList.read() instead.")
+
+
+@data_loader("JWST c1d multi", identifier=identify_jwst_c1d_multi_fits,
+             dtype=SpectrumList, extensions=['fits'])
+def jwst_c1d_multi_loader(file_obj, **kwargs):
+    """
+    Loader for JWST x1d 1-D spectral data in FITS format
+
+    Parameters
+    ----------
+    file_obj: str, file-like, or HDUList
+          FITS file name, object (provided from name by Astropy I/O Registry),
+          or HDUList (as resulting from astropy.io.fits.open()).
+
+    Returns
+    -------
+    SpectrumList
+        A list of the spectra that are contained in the file.
+    """
+    return _jwst_spec1d_loader(file_obj, extname='COMBINE1D', **kwargs)
+
+
 @data_loader("JWST x1d", identifier=identify_jwst_x1d_fits, dtype=Spectrum1D,
              extensions=['fits'])
 def jwst_x1d_single_loader(file_obj, **kwargs):
@@ -91,7 +164,7 @@ def jwst_x1d_single_loader(file_obj, **kwargs):
     Spectrum1D
         The spectrum contained in the file.
     """
-    spectrum_list = _jwst_x1d_loader(file_obj, **kwargs)
+    spectrum_list = _jwst_spec1d_loader(file_obj, extname='EXTRACT1D', **kwargs)
     if len(spectrum_list) == 1:
         return spectrum_list[0]
     else:
@@ -116,33 +189,38 @@ def jwst_x1d_multi_loader(file_obj, **kwargs):
     SpectrumList
         A list of the spectra that are contained in the file.
     """
-    return _jwst_x1d_loader(file_obj, **kwargs)
+    return _jwst_spec1d_loader(file_obj, extname='EXTRACT1D', **kwargs)
 
 
-def _jwst_x1d_loader(file_obj, **kwargs):
+def _jwst_spec1d_loader(file_obj, extname='EXTRACT1D', **kwargs):
     """Implementation of loader for JWST x1d 1-D spectral data in FITS format
 
     Parameters
     ----------
     file_obj: str, file-like, or HDUList
           FITS file name, object (provided from name by Astropy I/O Registry),
           or HDUList (as resulting from astropy.io.fits.open()).
+    extname: str
+        The name of the science extension. Either "COMBINE1D" or "EXTRACT1D".  By default "EXTRACT1D".
 
     Returns
     -------
     SpectrumList
         A list of the spectra that are contained in the file.
     """
 
+    if extname not in ['COMBINE1D', 'EXTRACT1D']:
+        raise ValueError('Incorrect extname given for 1d spectral data.')
+
     spectra = []
 
     with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist:
 
         primary_header = hdulist["PRIMARY"].header
 
         for hdu in hdulist:
-            # Read only the BinaryTableHDUs named EXTRACT1D and SCI
-            if hdu.name != 'EXTRACT1D':
+            # Read only the BinaryTableHDUs named COMBINE1D/EXTRACT1D and SCI
+            if hdu.name != extname:
                 continue
 
             data = Table.read(hdu)
@@ -155,18 +233,26 @@ def _jwst_x1d_loader(file_obj, **kwargs):
             # SRCTYPE used to be in primary header, but was moved around. As
             # per most recent pipeline definition, it should be in the
             # EXTRACT1D extension.
+            #
+            # SRCTYPE should either be POINT or EXTENDED.  In some cases, it is UNKNOWN
+            # or missing.  If that's the case, default to using POINT as the SRCTYPE.
+            # Error out only when SRCTYPE is a bad value.
             srctype = None
             if "srctype" in hdu.header:
-                srctype = hdu.header.get("srctype")
+                srctype = hdu.header.get("srctype", None)
+
+            # checking if SRCTPYE is missing or UNKNOWN
+            if not srctype or srctype == 'UNKNOWN':
+                log.warning('SRCTYPE is missing or UNKNOWN in JWST x1d loader. '
+                            'Defaulting to srctype="POINT".')
+                srctype = 'POINT'
 
             if srctype == "POINT":
                 flux = Quantity(data["FLUX"])
                 uncertainty = StdDevUncertainty(data["ERROR"])
-
             elif srctype == "EXTENDED":
                 flux = Quantity(data["SURF_BRIGHT"])
                 uncertainty = StdDevUncertainty(hdu.data["SB_ERROR"])
-
             else:
                 raise RuntimeError(f"Keyword SRCTYPE is {srctype}.  It should "
                                    "be 'POINT' or 'EXTENDED'. Can't decide between `flux` and "

specutils/io/default_loaders/tests/test_jwst_reader.py

@@ -13,15 +13,15 @@
 from specutils import Spectrum1D, SpectrumList
 
 
-# The x1d reader tests --------------------------
+# The c1d/x1d reader tests --------------------------
 
-def create_spectrum_hdu(data_len, srctype=None, ver=1):
+def create_spectrum_hdu(data_len, srctype=None, ver=1, name='EXTRACT1D'):
     """Mock a JWST x1d BinTableHDU"""
     data = np.random.random((data_len, 5))
     table = Table(data=data, names=['WAVELENGTH', 'FLUX', 'ERROR', 'SURF_BRIGHT',
         'SB_ERROR'])
 
-    hdu = fits.BinTableHDU(table, name='EXTRACT1D')
+    hdu = fits.BinTableHDU(table, name=name)
     hdu.header['TUNIT1'] = 'um'
     hdu.header['TUNIT2'] = 'Jy'
     hdu.header['TUNIT3'] = 'Jy'
@@ -48,27 +48,46 @@ def x1d_single():
 
 
 @pytest.fixture(scope="function")
-def x1d_multi():
-    """Mock a JWST x1d multispec HDUList with 3 spectra"""
+def spec_single(request):
+    """Mock a JWST c1d/x1d HDUList with a single spectrum"""
+    name = request.param
+    hdulist = fits.HDUList()
+    hdulist.append(fits.PrimaryHDU())
+    hdulist["PRIMARY"].header["TELESCOP"] = ("JWST", "comment")
+    # Add a BinTableHDU that contains spectral data
+    hdulist.append(create_spectrum_hdu(100, 'POINT', ver=1, name=name))
+    # Mock the ASDF extension
+    hdulist.append(fits.BinTableHDU(name='ASDF'))
+
+    return hdulist
+
+
+@pytest.fixture(scope="function")
+def spec_multi(request):
+    """Mock a JWST c1d/x1d multispec HDUList with 3 spectra"""
+    name = request.param
     hdulist = fits.HDUList()
     hdulist.append(fits.PrimaryHDU())
     hdulist["PRIMARY"].header["TELESCOP"] = "JWST"
     # Add a few BinTableHDUs that contain spectral data
-    hdulist.append(create_spectrum_hdu(100, 'POINT', ver=1))
-    hdulist.append(create_spectrum_hdu(120, 'EXTENDED', ver=2))
-    hdulist.append(create_spectrum_hdu(110, 'POINT', ver=3))
+    hdulist.append(create_spectrum_hdu(100, 'POINT', ver=1, name=name))
+    hdulist.append(create_spectrum_hdu(120, 'EXTENDED', ver=2, name=name))
+    hdulist.append(create_spectrum_hdu(110, 'POINT', ver=3, name=name))
     # Mock the ASDF extension
     hdulist.append(fits.BinTableHDU(name='ASDF'))
 
     return hdulist
 
 
-def test_jwst_x1d_multi_reader(tmpdir, x1d_multi):
-    """Test SpectrumList.read for JWST x1d multi data"""
+@pytest.mark.parametrize('spec_multi, format',
+                         [('EXTRACT1D', 'JWST x1d multi'),
+                          ('COMBINE1D', 'JWST c1d multi')], indirect=['spec_multi'])
+def test_jwst_1d_multi_reader(tmpdir, spec_multi, format):
+    """Test SpectrumList.read for JWST c1d/x1d multi data"""
     tmpfile = str(tmpdir.join('jwst.fits'))
-    x1d_multi.writeto(tmpfile)
+    spec_multi.writeto(tmpfile)
 
-    data = SpectrumList.read(tmpfile, format='JWST x1d multi')
+    data = SpectrumList.read(tmpfile, format=format)
     assert type(data) is SpectrumList
     assert len(data) == 3
 
@@ -80,20 +99,39 @@ def test_jwst_x1d_multi_reader(tmpdir, x1d_multi):
     assert data[2].shape == (110,)
 
 
-def test_jwst_x1d_single_reader(tmpdir, x1d_single):
+@pytest.mark.parametrize('spec_single, format',
+                         [('EXTRACT1D', 'JWST x1d'),
+                          ('COMBINE1D', 'JWST c1d')], indirect=['spec_single'])
+def test_jwst_1d_single_reader(tmpdir, spec_single, format):
     """Test Spectrum1D.read for JWST x1d data"""
     tmpfile = str(tmpdir.join('jwst.fits'))
+    spec_single.writeto(tmpfile)
+
+    data = Spectrum1D.read(tmpfile, format=format)
+    assert type(data) is Spectrum1D
+    assert data.shape == (100,)
+
+
+@pytest.mark.parametrize("srctype", [None, "UNKNOWN"])
+def test_jwst_srctpye_defaults(tmpdir, x1d_single, srctype):
+    """ Test """
+    tmpfile = str(tmpdir.join('jwst.fits'))
+
+    # Add a spectrum with missing or UNKNOWN SRCTYPE (mutate the fixture)
+    x1d_single['EXTRACT1D'].header['SRCTYPE'] == srctype
     x1d_single.writeto(tmpfile)
 
     data = Spectrum1D.read(tmpfile, format='JWST x1d')
     assert type(data) is Spectrum1D
     assert data.shape == (100,)
+    assert x1d_single['EXTRACT1D'].header['SRCTYPE'] == "POINT"
 
 
-def test_jwst_x1d_single_reader_no_format(tmpdir, x1d_single):
-    """Test Spectrum1D.read for JWST x1d data without format arg"""
+@pytest.mark.parametrize('spec_single', ['EXTRACT1D', 'COMBINE1D'], indirect=['spec_single'])
+def test_jwst_1d_single_reader_no_format(tmpdir, spec_single):
+    """Test Spectrum1D.read for JWST c1d/x1d data without format arg"""
     tmpfile = str(tmpdir.join('jwst.fits'))
-    x1d_single.writeto(tmpfile)
+    spec_single.writeto(tmpfile)
 
     data = Spectrum1D.read(tmpfile)
     assert type(data) is Spectrum1D
@@ -102,10 +140,11 @@ def test_jwst_x1d_single_reader_no_format(tmpdir, x1d_single):
     assert data.spectral_axis.unit == u.um
 
 
-def test_jwst_x1d_multi_reader_no_format(tmpdir, x1d_multi):
-    """Test Spectrum1D.read for JWST x1d data without format arg"""
+@pytest.mark.parametrize('spec_multi', ['EXTRACT1D', 'COMBINE1D'], indirect=['spec_multi'])
+def test_jwst_1d_multi_reader_no_format(tmpdir, spec_multi):
+    """Test Spectrum1D.read for JWST c1d/x1d data without format arg"""
     tmpfile = str(tmpdir.join('jwst.fits'))
-    x1d_multi.writeto(tmpfile)
+    spec_multi.writeto(tmpfile)
 
     data = SpectrumList.read(tmpfile)
     assert type(data) is SpectrumList
@@ -115,39 +154,42 @@ def test_jwst_x1d_multi_reader_no_format(tmpdir, x1d_multi):
         assert isinstance(item, Spectrum1D)
 
 
-def test_jwst_x1d_multi_reader_check_units(tmpdir, x1d_multi):
-    """Test units for Spectrum1D.read for JWST x1d data"""
+@pytest.mark.parametrize('spec_multi', ['EXTRACT1D', 'COMBINE1D'], indirect=['spec_multi'])
+def test_jwst_1d_multi_reader_check_units(tmpdir, spec_multi):
+    """Test units for Spectrum1D.read for JWST c1d/x1d data"""
     tmpfile = str(tmpdir.join('jwst.fits'))
-    x1d_multi.writeto(tmpfile)
+    spec_multi.writeto(tmpfile)
 
     data = SpectrumList.read(tmpfile)
     assert data[0].unit == u.Jy
     assert data[1].unit == u.MJy / u.sr
     assert data[2].unit == u.Jy
 
 
-def test_jwst_x1d_reader_meta(tmpdir, x1d_single):
-    """Test that the Primary and EXTRACT1D extension headers are merged in meta"""
+@pytest.mark.parametrize('spec_single', ['EXTRACT1D', 'COMBINE1D'], indirect=['spec_single'])
+def test_jwst_1d_reader_meta(tmpdir, spec_single):
+    """Test that the Primary and COMBINE1D/EXTRACT1D extension headers are merged in meta"""
     tmpfile = str(tmpdir.join('jwst.fits'))
-    x1d_single.writeto(tmpfile)
+    spec_single.writeto(tmpfile)
 
     data = Spectrum1D.read(tmpfile)
     assert ('TELESCOP', 'JWST') in data.meta['header'].items()
     assert ('SRCTYPE', 'POINT') in data.meta['header'].items()
 
 
-def test_jwst_x1d_single_reader_fail_on_multi(tmpdir, x1d_multi):
-    """Make sure Spectrum1D.read on JWST x1d with many spectra errors out"""
+@pytest.mark.parametrize('spec_multi', ['EXTRACT1D', 'COMBINE1D'], indirect=['spec_multi'])
+def test_jwst_1d_single_reader_fail_on_multi(tmpdir, spec_multi):
+    """Make sure Spectrum1D.read on JWST c1d/x1d with many spectra errors out"""
     tmpfile = str(tmpdir.join('jwst.fits'))
-    x1d_multi.writeto(tmpfile)
+    spec_multi.writeto(tmpfile)
 
     with pytest.raises(IORegistryError):
         Spectrum1D.read(tmpfile)
 
 
-@pytest.mark.parametrize("srctype", [None, "UNKNOWN"])
+@pytest.mark.parametrize("srctype", ["BADVAL"])
 def test_jwst_reader_fail(tmpdir, x1d_single, srctype):
-    """Check that the reader fails when SRCTYPE is not set or is UNKNOWN"""
+    """Check that the reader fails when SRCTYPE is a BADVAL"""
     tmpfile = str(tmpdir.join('jwst.fits'))
     hdulist = x1d_single
     # Add a spectrum with bad SRCTYPE (mutate the fixture)