Flux with electron/s cannot be read back in from file (roundtrip problem)

[pllim]

Consider this use case:

import numpy as np
from astropy import units as u
from specutils import Spectrum1D

flux = np.random.random(10) * (u.electron / u.s)
wave = (np.arange(10) + 1) * u.m
sp = Spectrum1D(flux=flux, spectral_axis=wave)  # Works fine

But then when writing it back out into a FITS file:

>>> sp.write('flux_es.fits')
WARNING: The unit 'electron / s' could not be saved in native FITS format
and cannot be recovered in reading. It can roundtrip within astropy by using
QTable both to write and read back, though one has to enable the unit before
reading. [astropy.io.fits.convenience]

Now, that file cannot be read back in, as the warning indicated. But why can't specutils write it out in a way that can roundtrip?

>>> sp2 = Spectrum1D.read('flux_es.fits')
.../astropy/nddata/mixins/ndio.py:59, in NDDataRead.__call__(self, *args, **kwargs)
     58 def __call__(self, *args, **kwargs):
---> 59     return self.registry.read(self._cls, *args, **kwargs)

.../astropy/io/registry/core.py:219, in UnifiedInputRegistry.read(self, cls, format, cache, *args, **kwargs)
    214     format = self._get_valid_format(
    215         "read", cls, path, fileobj, args, kwargs
    216     )
    218 reader = self.get_reader(format, cls)
--> 219 data = reader(*args, **kwargs)
    221 if not isinstance(data, cls):
    222     # User has read with a subclass where only the parent class is
    223     # registered.  This returns the parent class, so try coercing
    224     # to desired subclass.
    225     try:

.../specutils/io/default_loaders/tabular_fits.py:86, in tabular_fits_loader(file_obj, column_mapping, hdu, **kwargs)
     83 # If no column mapping is given, attempt to parse the file using
     84 # unit information
     85 if column_mapping is None:
---> 86     return generic_spectrum_from_table(tab, wcs=wcs, **kwargs)
     88 return spectrum_from_column_mapping(tab, column_mapping, wcs=wcs)

.../specutils/io/parsing_utils.py:258, in generic_spectrum_from_table(table, wcs, **kwargs)
    256 flux_column = _find_spectral_column(table, colnames, spectral_axis)
    257 if flux_column is None:
--> 258     raise IOError("Could not identify column containing the flux")
    259 flux = table[flux_column].to(table[flux_column].unit)
    260 colnames.remove(flux_column)

OSError: Could not identify column containing the flux

Desired behavior: specutils should know how to workaround this limitation for proper I/O roundtripping of its own product.

🐱