Revamp Spectrum1D tabular-fits writer

[kelle]

I'm having a lot of trouble figuring out how to write a Spectrum1D object to a file. A few of my questions:

• What file types are supported? FITS, ASDF?
• How do I also make sure the header data is also saved?

I'm happy to work on this once I understand it. :)

[kelle]

As a result of digging into this, I discovered that header data does not currently roundtrip using Spectrum1D.write. I propose that this be implemented using something like this:

# Make a dedicated HDU for the header and spectrum
hdu0 = fits.PrimaryHDU(header=header)
hdu1 = fits.BinTableHDU(data=NIR_table)
hdu1.header['EXTNAME'] = 'SPECTRUM'
hdu1.header.set('OBJECT', NIRSpec_dict['object_name'], 'Object Name')

# Combine the two HDUs into a multi-extension FITS (MEF)
spectrum_mef = fits.HDUList([hdu0, hdu1])  # hdu0 is header and hdu1 is data

# Write out the MEF to a file
spectrum_mef.writeto(fits_filename, output_verify="exception")

# Demonstrate how to read the file back in using specutils
spec1d_NIR = Spectrum1D.read(fits_filename, format='tabular-fits')
header_NIR = fits.getheader(fits_filename)

If folks agree, I would be happy to implement this given some guidance.

[pllim]

@dhomeier might also be interested in this issue since he is working on:

• Fix and extend wcs1d-fits loader for multi-D WCS #1009

[kelle]

Summarizing relevant discussions from Slack:

• If someone modifies the spectrum, the keywords might become inaccurate. My proposal is to add a warning that says something like, "Header keywords are the responsibility of the user and are not guaranteed to be accurate." We can open issues and what not on specutils routines which modify header keywords and gradually make this less true. But I think that header round tripping is important enough to tolerate this shortcoming.

• I am not sure how to implement this for Spectrum2D, 3D or SpectrumList. I am happy to just do it for Spectrum1D -- this would be immensely useful -- OR to figure out how to make it extendable to these other classes.

[dhomeier]

First I am not sure this issue really is a duplicate or, or completes #908 – the latter is primarily on better documentation of the existing options for writing spectra, while this one at least with the recent discussion is heading more into proposing new functionality.

I think the current "generic" read/write loaders (in fact the only ones implementing write functionality) have already been constructed with some round-tripping support in mind, though incomplete yet.
E.g. #1009 tries to improve write and read-back of the mask and uncertainty attributes.

One important point here is that Spectrum1D has no header property as such; that is rather tied to the specific FITS file layout.
The location to hold general information of that type is meta, which could be constructed directly from header information of a FITS HDU, or be a superset of the same, or could also contain objects that are very hard or impossible to serialise to a FITS header. But round-trip of all information that can be represented as a FITS header would certainly be a reasonable target. I think the main problems here are

1. Existing FITS loaders, if they save header information, variously do so by extracting the entire header directly into spectrum.meta or creating a spectrum.meta['header'] dictionary that contains only the FITS header information. The latter is more robust as other content could be added to meta later, but within that scheme the FITS writers would probably be restricted to writing only meta['header'] and no other metadata information (otherwise I foresee running into an infinite recursion...).
2. With the above convention there is still the question whether to load the header data from the HDU containing the spectrum, the primary HDU (which at least for tabular-fits will be different, but still may contain relevant information), or all of them; and then where to write them back.
3. Settling onto one of those schemes has already been discussed at various time, but I think with no clear conclusion yet. Also there has been some concern about always saving all header information, potentially creating inflated meta attributes with content most people might never look at.
4. As noted from the Slack discussions, many FITS keywords really are connected directly to the data format and representation and thus are already present in Spectrum1D e.g. as flux.dtype, flux.unit. Still storing them in meta in a way is redundant and may get out of sync if the data parts are modified. I think to some extent the header parsers already try to take care of that and are stripping cards like BUNIT (which would be overwritten when writing back to a HDU anyway).

A good starting point could perhaps be to review what info is already preserved in round-trip tests like

specutils/specutils/tests/test_loaders.py

Lines 608 to 611 in ca3153f

	# Now write with updated header information from spectrum.meta
	spectrum.meta.update({'OBSERVER': 'Rosse', 'EXPTIME': 32.1, 'NAXIS2': 12})
	spectrum.write(tmpfile, format='tabular-fits', overwrite=True,
	update_header=True)

and what is missing that you would like to be saved. As sketched in your proposed implementation, getting header data back from tabular-fits would require bypassing the Table.read interface in the loader, as it provides to my knowledge no option for saving header data. wcs1d-fits would already do this, but then again is inconsistent with even its own writer, since it is using the meta['header'] scheme, while both writers are trying to write back all of meta with update_header=True set. Again I think there has been discussion why that would be problematic to change, and break existing behaviour.

[kelle]

We've tested the spectrum.meta.update method and it's pretty close! Problems include

• extra END keywords. This makes many header keywords inaccessible by other tools. For example,

header=fits.getheader(tmpfile)

results in a header which just has the minimal keywords and not any of the things we added. (I assume this is because of the no-loneger-necessary END keyword)

• COMMENT and HISTORY are stripped. I assume this is on purpose. Is there a way to keep them?

[dhomeier]

That's what you get with fits.getheader()? That would be an io.fits issue, but I cannot reproduce any stripping of COMMENT cards.
Looking back at

specutils/specutils/io/default_loaders/tabular_fits.py

Lines 117 to 127 in ca3153f

	header = spectrum.meta.get('header', fits.header.Header()).copy()
	if update_header:
	hdr_types = (str, int, float, complex, bool,
	np.floating, np.integer, np.complexfloating, np.bool_)
	header.update([keyword for keyword in spectrum.meta.items() if
	isinstance(keyword[1], hdr_types)])
	# Strip header of FITS reserved keywords
	for keyword in ['NAXIS', 'NAXIS1', 'NAXIS2']:
	header.remove(keyword, ignore_missing=True)

for anything you put into a meta['header'] directory you should actually not even need the update_header option; but from either source it should keep everything that has types it can handle, except for NAXIS*.
That's also something perhaps to be matched in wcs1d-fits, which by default only initialises from the WCS header.

[kelle]

• In my case, the COMMENT and HISTORY cards are stripped because they are <class 'astropy.io.fits.header._HeaderCommentaryCards'>
• Even though I use spec.meta.update(header), the meta doesn't roundtrip unless I use update_header = True. Not sure why. (This is not a big deal to me at the moment.)
• The meta dictionary doesn't keep the comment that goes along with each keyword/value pair. A FITS header a list of three elements, not just two. The keyword, the value, AND the comment. So as soon as I do spec1d.meta.update(header), I've lost information.
• The extra ENDs must come from io.fits or from Table. I've posted on Slack to track this one down.

I think the bottom line is that meta doesn't maintain the functionality that one expects of FITS headers.

We could continue to go down this road or maybe I should just write a custom reader/writer that is based on tabular-fits but expects HDU0 to have the header data.

[kelle]

Just to show y'all what I'm talking about, here's the degradation of the header information using the current tooling.

The header data I have compiled and added using spec1d.meta.update(header):

OBJECT  = 'VHS1256b'           / Name of observed object                        
RA      = '194.007637'         / [deg] Pointing position                        
DEC     = '-12.957692'         / [deg] Pointing position                        
TELESCOP= 'JWST    '           / name of telescope                              
INSTRUME= 'NIRSpec '           / name of instrument                             
SPECBAND= 'NIR     '           / SED.bandpass                                   
SPEC_VAL=                 2.25 / [um] Characteristic spec coord                 
SPEC_BW =                  4.5 / [um] Width of spectrum                         
TDMIN1  =                  0.0 / [um] Starting wavelength                       
TDMAX1  =                  4.5 / [um] Ending wavelength                         
DATE    = '2023-02-24'         / Date of file creation                          
COMMENT To read this file in with specutils use Spectrum1D.read() with format = 
COMMENT 'tabular-fits'                                                          
HISTORY Data Reduced by JWST pipeline version 1.7.2, CRDS version 11.16.12, CRDS
HISTORY  context jwst_0977.pmap  

The header obtained via meta in the read back in file:

OrderedDict([('OBJECT', 'VHS1256b'),
             ('RA', '194.007637'),
             ('DEC', '-12.957692'),
             ('TELESCOP', 'JWST'),
             ('INSTRUME', 'NIRSpec'),
             ('SPECBAND', 'NIR'),
             ('SPEC_VAL', 2.25),
             ('SPEC_BW', 4.5),
             ('TDMIN1', 0.0),
             ('TDMAX1', 4.5),
             ('DATE', '2023-02-24')])

The header obtained using fits.getheader:

SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                    8 / array data type                                
NAXIS   =                    0 / number of array dimensions                     
EXTEND  =                    T   

And the header directly from the FITS file

SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                    8 / array data type                                
NAXIS   =                    0 / number of array dimensions                     
EXTEND  =                    T                                                  
END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             XTENSION= 'BINTABLE'           / binary table extension                         
BITPIX  =                    8 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                   24 / length of dimension 1                          
NAXIS2  =                   10 / length of dimension 2                          
PCOUNT  =                    0 / number of group parameters                     
GCOUNT  =                    1 / number of groups                               
TFIELDS =                    3 / number of table fields                         
TTYPE1  = 'wavelength'                                                          
TFORM1  = 'D       '                                                            
TUNIT1  = 'um      '                                                            
TTYPE2  = 'flux    '                                                            
TFORM2  = 'D       '                                                            
TUNIT2  = 'Jy      '                                                            
TTYPE3  = 'uncertainty'                                                         
TFORM3  = 'D       '                                                            
TUNIT3  = 'Jy      '                                                            
OBJECT  = 'VHS1256b'                                                            
RA      = '194.007637'                                                          
DEC     = '-12.957692'                                                          
TELESCOP= 'JWST    '                                                            
INSTRUME= 'NIRSpec '                                                            
SPECBAND= 'NIR     '                                                            
SPEC_VAL=                 2.25                                                  
SPEC_BW =                  4.5                                                  
TDMIN1  =                  0.0                                                  
TDMAX1  =                  4.5                                                  
DATE    = '2023-02-24'                                                          
END