From a6fd06f975306b1c7a15337461e1e3873fe178a7 Mon Sep 17 00:00:00 2001 From: Kelle Cruz Date: Thu, 11 Jul 2024 11:40:58 -0400 Subject: [PATCH] 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 Co-authored-by: Ricky O'Steen 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 --- CHANGES.rst | 5 ++ specutils/io/default_loaders/tabular_fits.py | 39 ++++++++--- specutils/tests/test_loaders.py | 68 +++++++++++++++++--- 3 files changed, 93 insertions(+), 19 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index bb5bf08a8..7c7621907 100644 --- a/CHANGES.rst +++ b/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) ------------------- diff --git a/specutils/io/default_loaders/tabular_fits.py b/specutils/io/default_loaders/tabular_fits.py index ad3c6b8aa..c11c072c9 100644 --- a/specutils/io/default_loaders/tabular_fits.py +++ b/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) diff --git a/specutils/tests/test_loaders.py b/specutils/tests/test_loaders.py index c033c4097..e3c42be2d 100644 --- a/specutils/tests/test_loaders.py +++ b/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')