Fix for Tabular-fits reader/writer fix (#1113)

* stand alone old tabular fits tests

* modified tests to desired header behavior

* write header to HDU0

* test header data roundtrips in hdu0

* put tests back in big file for easier comparison

* restore old legacy tests.

* codestyle fix

* modify the reader to get the header from HDU0

* tests pass now, removed print statements.

* specify format while reading file

* Debugging compressed tabular fits

* Remove format spec in read in test

* Fix rebase

* Remove commented prints, restore NAXIS check

* Update changelog

* Fix comment

* Update CHANGES.rst

* Add option to store data header instead

* Move header storage inside open file context

* Update CHANGES.rst

Co-authored-by: Derek Homeier <709020+dhomeier@users.noreply.github.com>

* Add `store_data_header` option to tabular-fits writer

* Fix `store_data_header` tests and add roundtrip

---------

Co-authored-by: Simon Torres <simon.torres@noirlab.edu>
Co-authored-by: Ricky O'Steen <rosteen@stsci.edu>
Co-authored-by: Ricky O'Steen <39831871+rosteen@users.noreply.github.com>
Co-authored-by: Derek Homeier <709020+dhomeier@users.noreply.github.com>
Co-authored-by: Derek Homeier <dhomeier@users.noreply.github.com>

CHANGES.rst

@@ -10,6 +10,11 @@ Bug Fixes
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
+- Changed the ``tabular-fits`` reader/writer to round-trip the primary header.
+   The reader now by default reads the primary header into ``meta``; the old behaviour
+   can be restored by setting the option ``store_data_header=True``. The writer is writing
+   copies of ``meta['header']`` to both primary and data extension headers. [#1113]
+
 1.15.0 (2024-05-01)
 -------------------
 

specutils/io/default_loaders/tabular_fits.py

@@ -41,7 +41,7 @@ def identify_tabular_fits(origin, *args, **kwargs):
 
 @data_loader("tabular-fits", identifier=identify_tabular_fits,
              dtype=Spectrum1D, extensions=['fits', 'fit'], priority=6)
-def tabular_fits_loader(file_obj, column_mapping=None, hdu=1, **kwargs):
+def tabular_fits_loader(file_obj, column_mapping=None, hdu=1, store_data_header=False, **kwargs):
     """
     Load spectrum from a FITS file.
 
@@ -52,6 +52,9 @@ def tabular_fits_loader(file_obj, column_mapping=None, hdu=1, **kwargs):
             or HDUList (as resulting from astropy.io.fits.open()).
     hdu: int
         The HDU of the fits file (default: 1st extension) to read from
+    store_data_header: bool
+        Defaults to ``False``, which stores the primary header in ``Spectrum1D.meta['header']``.
+        Set to ``True`` to instead store the header from the specified data HDU.
     column_mapping : dict
         A dictionary describing the relation between the FITS file columns
         and the arguments of the `Spectrum1D` class, along with unit
@@ -71,8 +74,11 @@ def tabular_fits_loader(file_obj, column_mapping=None, hdu=1, **kwargs):
     # routines to search for spectral axis information in the file.
     with read_fileobj_or_hdulist(file_obj, **kwargs) as hdulist:
         wcs = WCS(hdulist[hdu].header)
-
-    tab = Table.read(file_obj, format='fits', hdu=hdu)
+        tab = Table.read(hdulist[hdu])
+        if store_data_header:
+            tab.meta = hdulist[hdu].header
+        else:
+            tab.meta = hdulist[0].header
 
     # Minimal checks for wcs consistency with table data -
     # assume 1D spectral axis (having shape (0, NAXIS1),
@@ -90,7 +96,7 @@ def tabular_fits_loader(file_obj, column_mapping=None, hdu=1, **kwargs):
 
 
 @custom_writer("tabular-fits")
-def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, **kwargs):
+def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, store_data_header=False, **kwargs):
     """
     Write spectrum to BINTABLE extension of a FITS file.
 
@@ -102,7 +108,11 @@ def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, **kwarg
     hdu: int
         Header Data Unit in FITS file to write to (currently only extension HDU 1)
     update_header: bool
-        Update FITS header with all compatible entries in `spectrum.meta`
+        Write all compatible items in ``Spectrum1D.meta`` directly to FITS header;
+        this will overwrite any identically named keys from ``Spectrum1D.meta['header']``.
+    store_data_header: bool
+        If ``True``, store ``Spectrum1D.meta['header']`` in the header of the target data HDU
+        instead of the primary header (default ``False``).
     wunit : str or `~astropy.units.Unit`
         Unit for the spectral axis (wavelength or frequency-like)
     funit : str or `~astropy.units.Unit`
@@ -117,6 +127,7 @@ def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, **kwarg
 
     header = spectrum.meta.get('header', fits.header.Header()).copy()
 
+    # Should this warn about incompatible items or overwriting existing ones?
     if update_header:
         hdr_types = (str, int, float, complex, bool,
                      np.floating, np.integer, np.complexfloating, np.bool_)
@@ -181,8 +192,18 @@ def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, **kwarg
         if columns[c].ndim > 1:
             columns[c] = columns[c].T
 
-    tab = Table(columns, names=colnames, meta=header)
+    tab = Table(columns, names=colnames)
+    if store_data_header:
+        hdu0 = fits.PrimaryHDU()
+        hdu1 = fits.BinTableHDU(data=tab, header=header)
+    else:
+        hdu0 = fits.PrimaryHDU(header=header)
+        hdu1 = fits.BinTableHDU(data=tab)
+
+    # This will overwrite any 'EXTNAME' previously read from a valid header; should it?
+    hdu1.header.update(EXTNAME='DATA')
+
+    hdulist = fits.HDUList([hdu0, hdu1])
 
-    # Todo: support writing to other HDUs than the default (1st)
-    #       and an 'update' mode so different HDUs can be written to separately
-    tab.write(file_name, format="fits", **kwargs)
+    # TODO: Use output_verify options to check for valid FITS
+    hdulist.writeto(file_name, **kwargs)

specutils/tests/test_loaders.py

@@ -657,38 +657,86 @@ def test_tabular_fits_mask(tmp_path, mask_type):
         assert sp1.mask.dtype == sp2.mask.dtype
 
 
-def test_tabular_fits_maskheader(tmp_path):
+@pytest.mark.parametrize("metadata_hdu", [0, 1])
+def test_tabular_fits_maskheader(tmp_path, metadata_hdu):
     # Create a small data set + header with reserved FITS keywords
     disp = np.linspace(1, 1.2, 21) * u.AA
     flux = np.random.normal(0., 1.0e-14, disp.shape[0]) * u.Jy
     hdr = fits.header.Header({'TELESCOP': 'Leviathan', 'APERTURE': 1.8,
-                              'OBSERVER': 'Parsons', 'NAXIS': 1, 'NAXIS1': 8})
+                              'OBSERVER': 'Parsons'})
 
     spectrum = Spectrum1D(flux=flux, spectral_axis=disp, meta={'header': hdr})
     tmpfile = str(tmp_path / '_tst.fits')
-    spectrum.write(tmpfile, format='tabular-fits')
+    spectrum.write(tmpfile, format='tabular-fits', store_data_header=bool(metadata_hdu))
 
     # Read it in and check against the original
     with fits.open(tmpfile) as hdulist:
+
+        # Test HDU0 header
         assert hdulist[0].header['NAXIS'] == 0
+
+        # keys relevant to datashape are in HDU1 header
         assert hdulist[1].header['NAXIS'] == 2
         assert hdulist[1].header['NAXIS2'] == disp.shape[0]
-        assert hdulist[1].header['OBSERVER'] == 'Parsons'
+
+        # Test storage of metadata in selected HDU header
+        assert hdulist[metadata_hdu].header['OBSERVER'] == 'Parsons'
+
+
+@pytest.mark.parametrize("metadata_hdu", [0, 1])
+def test_tabular_fits_roundtrip_header(tmp_path, metadata_hdu):
+    """Test roundtrip of full header."""
+    disp = np.linspace(1, 1.2, 21) * u.AA
+    flux = np.random.normal(0., 1.0e-14, disp.shape[0]) * u.erg / (u.s * u.cm**2 * u.AA)
+    hdr = fits.header.Header({'TELESCOP': 'Crystal', 'OBSERVER': 'Cruz'})
+    spec = Spectrum1D(flux=flux, spectral_axis=disp, meta={'header': hdr})
+    tmpfile = str(tmp_path / '_tst.fits')
+    spec.write(tmpfile, format='tabular-fits', store_data_header=bool(metadata_hdu))
+    spec = Spectrum1D.read(tmpfile, format='tabular-fits', store_data_header=bool(metadata_hdu))
+
+    # Confirm HDU-specific header cards are read back in.
+    assert spec.meta['header']['NAXIS'] == metadata_hdu * 2
+    assert spec.meta['header']['TELESCOP'] == 'Crystal'
+
+    # Write it out and read back again to compare full headers.
+    tmpfile = str(tmp_path / '_tst2.fits')
+    spec.write(tmpfile, format='tabular-fits', store_data_header=bool(metadata_hdu))
+    spectrum = Spectrum1D.read(tmpfile, format='tabular-fits', store_data_header=bool(metadata_hdu))
+
+    # 'EXTNAME' is rewritten by the writer and may pop up in a different location.
+    if metadata_hdu == 1:
+        assert spectrum.meta['header'].pop('EXTNAME') == spec.meta['header'].pop('EXTNAME')
+    assert spectrum.meta['header'] == spec.meta['header']
+
+
+@pytest.mark.parametrize("metadata_hdu", [0, 1])
+def test_tabular_fits_update_header(tmp_path, metadata_hdu):
+    """Test modification of header items using ``update_header`` option."""
+    disp = np.linspace(1, 1.2, 21) * u.AA
+    flux = np.random.normal(0., 1.0e-14, disp.shape[0]) * u.erg / (u.s * u.cm**2 * u.AA)
+    hdr = fits.header.Header({'TELESCOP': 'Crystal', 'OBSERVER': 'Cruz'})
+    spec = Spectrum1D(flux=flux, spectral_axis=disp, meta={'header': hdr})
+    tmpfile = str(tmp_path / '_tst.fits')
+    spec.write(tmpfile, format='tabular-fits', store_data_header=bool(metadata_hdu))
+    spectrum = Spectrum1D.read(tmpfile, format='tabular-fits', store_data_header=bool(metadata_hdu))
+
+    assert spectrum.meta['header']['OBSERVER'] == 'Cruz'
+    assert spectrum.meta['header']['TELESCOP'] == 'Crystal'
 
     # Now write with updated header information from spectrum.meta
-    spectrum.meta.update({'OBSERVER': 'Rosse', 'EXPTIME': 32.1, 'NAXIS2': 12})
+    spectrum.meta['header'].update({'OBSERVER': 'Rosse', 'EXPTIME': 32.1, 'NAXIS2': 12})
     spectrum.write(tmpfile, format='tabular-fits', overwrite=True,
-                   update_header=True)
+                   update_header=True, store_data_header=bool(metadata_hdu))
 
     with fits.open(tmpfile) as hdulist:
         assert hdulist[1].header['NAXIS2'] == disp.shape[0]
-        assert hdulist[1].header['OBSERVER'] == 'Rosse'
-        assert_allclose(hdulist[1].header['EXPTIME'], 3.21e1)
+        assert hdulist[metadata_hdu].header['OBSERVER'] == 'Rosse'
+        assert_allclose(hdulist[metadata_hdu].header['EXPTIME'], 3.21e1)
 
     # Test that unsupported types (dict) are not added to written header
     spectrum.meta['MYHEADER'] = {'OBSDATE': '1848-02-26', 'TARGET': 'M51'}
     spectrum.write(tmpfile, format='tabular-fits', overwrite=True,
-                   update_header=True)
+                   update_header=True, store_data_header=bool(metadata_hdu))
 
     with fits.open(tmpfile) as hdulist:
         assert 'MYHEADER' not in hdulist[0].header
@@ -703,7 +751,7 @@ def test_tabular_fits_autowrite(tmp_path):
     disp = np.linspace(1, 1.2, 21) * u.AA
     flux = np.random.normal(0., 1.0e-14, disp.shape[0]) * u.W / (u.m**2 * u.AA)
     hdr = fits.header.Header({'TELESCOP': 'Leviathan', 'APERTURE': 1.8,
-                              'OBSERVER': 'Parsons', 'NAXIS': 1, 'NAXIS1': 8})
+                              'OBSERVER': 'Parsons'})
 
     spectrum = Spectrum1D(flux=flux, spectral_axis=disp, meta={'header': hdr})
     tmpfile = str(tmp_path / '_tst.fits')