Having problems using specutils to write Spectrum1D object to FITS file

[kakirastern]

I have tried using the approach as suggested on https://specutils.readthedocs.io/en/stable/custom_loading.html#creating-a-custom-writer to save a Spectrum1D object to the FITS file format, but have been encountered various problems.

The Spectrum1D object concerned is an integrated spectrum obtained from some data cube of a well-known planetary nebula after splicing a two spectra (one on the red side and one on the blue side) into one in order to derive some optical spectrum for the scientific target. The given Spectrum1D object can be specified by its flux and spectral_axis as follows:

>>> type(final_spec)
specutils.spectra.spectrum1d.Spectrum1D

>>> final_spec.flux.value
array([1.37377823e-15, 1.42879234e-15, 1.51984429e-15, ...,
       2.73236522e-17, 3.10544382e-17, 2.67768878e-42])

>>> final_spec.spectral_axis.value
array([3500.        , 3500.77431644, 3501.54863288, ..., 6999.12060362,
       6999.55980181, 6999.999     ])

At first I tried the plain old write method without specifying the supposedly optional format parameter. But I got the following error with traceback:

---------------------------------------------------------------------------
IORegistryError                           Traceback (most recent call last)
<ipython-input-120-2f4909935542> in <module>
----> 1 final_spec.write('/Users/krisstern/Desktop/junk.fits')

~/astropy/astropy/nddata/mixins/ndio.py in __call__(self, *args, **kwargs)
     97 
     98     def __call__(self, *args, **kwargs):
---> 99         registry.write(self._instance, *args, **kwargs)
    100 
    101 

~/astropy/astropy/io/registry.py in write(data, format, *args, **kwargs)
    561 
    562         format = _get_valid_format(
--> 563             'write', data.__class__, path, fileobj, args, kwargs)
    564 
    565     writer = get_writer(format, data.__class__)

~/astropy/astropy/io/registry.py in _get_valid_format(mode, cls, path, fileobj, args, kwargs)
    601                               " 'format' argument.\n"
    602                               "The available formats are:\n"
--> 603                               "{}".format(format_table_str))
    604     elif len(valid_formats) > 1:
    605         raise IORegistryError(

IORegistryError: Format could not be identified based on the file name or contents, please provide a 'format' argument.
The available formats are:
    Format    Read Write Auto-identify
------------- ---- ----- -------------
  fits-writer   No   Yes            No
 tabular-fits  Yes   Yes           Yes

However, if I first define the following:

from astropy.table import Table
from specutils.io.registers import custom_writer

@custom_writer("fits-writer")
def generic_fits(spectrum, file_name, **kwargs):
    flux = spectrum.flux.value
    disp = spectrum.spectral_axis.value
    meta = spectrum.meta

    tab = Table([disp, flux], names=("spectral_axis", "flux"), meta=meta)

    tab.write(file_name, format="fits")

then do

>>> final_spec.write("/Users/krisstern/Desktop/junk.fits", format="fits-writer")

even though no error is thrown at this stage, later on when I tried to open the "junk.fits" file using SAOImageDS9 I got the error "Unable to load fits /Users/krisstern/Desktop/junk.fits". So that is not very helpful.

I did later on resorted to using the astropy.io.fits's writeto method in a roundabout way to tackle the same task of saving the Spectrum1D object without major issues, but I would love to see some specutils writing solution for Spectrum1D into FITS file. That would be a good feature to have.

[dhomeier]

Could you list the values of spectrum.spectral_axis.unit, specrum.flux.unit and report what happens if you directly try to write it as table with

tab = Table([final_spec.spectral_axis.value, final_spec.flux.value], names=("spectral_axis", "flux"), meta=final_spec.meta)
tab.write("junk.fits", format="fits")

[kakirastern]

After some discussion on Astropy's Slack it was agreed that it would be better to clarity the file format and the reason why SAOImageDS9 cannot load the perfectly good file, maybe some documentation can be worked on and some write feature(s) could be added to the codebase. Will follow up with some PR's soon.

[kakirastern]

@dhomeier Sure, I will follow up once I have tested the issue out more shortly.

[dhomeier]

So fits.open() and Spectrum1D.read() are working just fine on the file? Agree that it sounds more like a DS9 issue then. Trying to confirm this with my probably very old (8.0.1) SAOImageDS9 installation, I found that it does not only fail with the same nondescriptive error on all test output files from Spectrum1D.write(), but on every single reference file from the tests (SDSS, MUSCLES etc.) except for those storing the spectra in IMAGE_HDU format. Are you sure that DS9 even supports reading from BINARY_TBL?

[dhomeier]

Just checked against 8.1, and indeed the only tables DS9 seems to support are Binary Events or HEALPIX Tables. So to be able to read the files with SAOImageDS9 you would need to create a writer saving the spectrum as a 1D-IMAGE - there was at one point a wcs1d-fits template that had been removed in #530, apparently because it was also just writing in BINTABLE format (i.e. identical to tabular-fits).

[kakirastern]

I am still experimenting with writing a FITS file with Spectrum1D. At one point I was able to save the spectrum according to the instructions in the official documentations, but am still having problems opening the file written this way with Spectrum1D.read() as described at https://specutils.readthedocs.io/en/stable/spectrum1d.html#reading-from-a-file. Seems to be some issue arising from the format parameter setting. Will continue experimenting and will open some PR(s) to fix this unexpected behavior accordingly.

[kakirastern]

For example, I am getting the following error when trying to read in/load a Spectrum1D object saved using the @custom_writer("fits-writer") approach as shown at https://specutils.readthedocs.io/en/stable/custom_loading.html#creating-a-custom-writer if I approach the task naively:

---------------------------------------------------------------------------
IORegistryError                           Traceback (most recent call last)
<ipython-input-184-ce5565b615ac> in <module>
----> 1 Spectrum1D.read('/Users/krisstern/Desktop/ic418_edited_integrated_spectrum_cube1.fits', format='fits')

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/nddata/mixins/ndio.py in __call__(self, *args, **kwargs)
     54 
     55     def __call__(self, *args, **kwargs):
---> 56         return registry.read(self._cls, *args, **kwargs)
     57 
     58 

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/io/registry.py in read(cls, format, *args, **kwargs)
    520                 'read', cls, path, fileobj, args, kwargs)
    521 
--> 522         reader = get_reader(format, cls)
    523         data = reader(*args, **kwargs)
    524 

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/io/registry.py in get_reader(data_format, data_class)
    453             "No reader defined for format '{}' and class '{}'.\n\nThe "
    454             "available formats are:\n\n{}".format(
--> 455                 data_format, data_class.__name__, format_table_str))
    456 
    457 

IORegistryError: No reader defined for format 'fits' and class 'Spectrum1D'.

The available formats are:

      Format      Read Write Auto-identify
----------------- ---- ----- -------------
    APOGEE apStar  Yes    No           Yes
   APOGEE apVisit  Yes    No           Yes
APOGEE aspcapStar  Yes    No           Yes
            ASCII  Yes    No           Yes
             ECSV  Yes    No           Yes
          HST/COS  Yes    No           Yes
         HST/STIS  Yes    No           Yes
             IPAC  Yes    No           Yes
 SDSS-I/II spSpec  Yes    No           Yes
 SDSS-III/IV spec  Yes    No           Yes
 Subaru-pfsObject  Yes    No           Yes
             iraf  Yes    No           Yes
      muscles-sed  Yes    No           Yes
     tabular-fits  Yes   Yes           Yes
       wcs1d-fits  Yes    No           Yes

[kakirastern]

I could open the saved Spectrum1D object saved in the "fits-writer" format with astropy's astropy.io.fits machinery though, however the loaded object looks a bit unnatural at first and appears as something as a astropy.io.fits.fitsrec.FITS_rec object like the following:

FITS_rec([(3500.        , 1.3737783e-15), (3500.77431644, 1.4287923e-15),
          (3501.54863288, 1.5198443e-15), ...,
          (6999.12060362, 2.7232715e-17), (6999.55980181, 3.1819597e-17),
          (6999.999     , 3.3631163e-42)],
         dtype=(numpy.record, [('spectral_axis', '>f8'), ('flux', '>f4')]))

[dhomeier]

I thought reading back as Spectrum1D was already working.

I could open the saved Spectrum1D object saved in the "fits-writer" format with astropy's astropy.io.fits machinery though, however the loaded object looks a bit unnatural at first and appears as something as a astropy.io.fits.fitsrec.FITS_rec object like the following:

That is just what is expected, basically a recarray of wavelength (spectral_axis) and flux.

If you want to read it in as a format specified by your custom_loader, that loader also needs to define a @custom_writer("fits") (or however you want to call it) as described in the docs. What happens when you try to read it without specifying a format, or with format='tabular-fits'?

[kakirastern]

I thought reading back as Spectrum1D was already working.

I could open the saved Spectrum1D object saved in the "fits-writer" format with astropy's astropy.io.fits machinery though, however the loaded object looks a bit unnatural at first and appears as something as a astropy.io.fits.fitsrec.FITS_rec object like the following:

That is just what is expected, basically a recarray of wavelength (spectral_axis) and flux.

Thanks! That's good to know.

If you want to read it in as a format specified by your custom_loader, that loader also needs to define a @custom_writer("fits") (or however you want to call it) as described in the docs. What happens when you try to read it without specifying a format, or with format='tabular-fits'?

I did try the custom_loader approach, but got some error originally. But I was able to get it working after adapting the sample code to the following:

from astropy.io import fits
from astropy.nddata import StdDevUncertainty
from astropy.table import Table
from astropy.units import Unit
from astropy.wcs import WCS

from specutils.io.registers import data_loader
from specutils import Spectrum1D

# Define an optional identifier. If made specific enough, this circumvents the
# need to add ``format="fits-format"`` in the ``Spectrum1D.read`` call.

def identify_generic_fits(origin, *args, **kwargs):
    return (isinstance(args[0], str) and
            os.path.splitext(args[0].lower())[1] == '.fits')

@data_loader("fits-format", identifier=identify_generic_fits,
             extensions=['fits'])
def generic_fits(file_name, **kwargs):
    with fits.open(file_name, **kwargs) as hdulist:
        header = hdulist[0].header

        tab = Table.read(file_name)

        meta = {'header': header}
        wcs = WCS(hdulist[0].header)
        data = tab["flux"] * Unit("ST")
        lamb = tab["spectral_axis"] * Unit("Angstrom")

    return Spectrum1D(flux=data, spectral_axis=lamb, wcs=wcs, meta=meta)

Otherwise if I blindly follow the example like a novice would I would be getting some error when I call spec as defined in

spec = Spectrum1D.read("/Users/krisstern/Desktop/ic418_edited_integrated_spectrum_cube1.fits", format="my-format")

[kakirastern]

BTW, using the format='tabular-fits' as a required parameter just gave me some flat-out mistake like

OSError: Could not identify column containing the flux

[kakirastern]

The corresponding traceback for the format='tabular-fits' ΟSError is as follows, for your information:

---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
<ipython-input-248-2c82bd8c0f2c> in <module>
----> 1 test_spec = Spectrum1D.read("/Users/krisstern/Desktop/ic418_edited_integrated_spectrum_cube1.fits", format="tabular-fits")

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/nddata/mixins/ndio.py in __call__(self, *args, **kwargs)
     54 
     55     def __call__(self, *args, **kwargs):
---> 56         return registry.read(self._cls, *args, **kwargs)
     57 
     58 

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/io/registry.py in read(cls, format, *args, **kwargs)
    521 
    522         reader = get_reader(format, cls)
--> 523         data = reader(*args, **kwargs)
    524 
    525         if not isinstance(data, cls):

~/specutils/specutils/io/default_loaders/tabular_fits.py in tabular_fits_loader(file_name, column_mapping, hdu, **kwargs)
     73     # unit information
     74     if column_mapping is None:
---> 75         return generic_spectrum_from_table(tab, wcs=wcs, **kwargs)
     76 
     77     return spectrum_from_column_mapping(tab, column_mapping, wcs=wcs)

~/specutils/specutils/io/parsing_utils.py in generic_spectrum_from_table(table, wcs, **kwargs)
    183     flux_column = _find_spectral_column(table,colnames,spectral_axis)
    184     if flux_column is None:
--> 185         raise IOError("Could not identify column containing the flux")
    186     flux = table[flux_column].to(table[flux_column].unit)
    187     colnames.remove(flux_column)

OSError: Could not identify column containing the flux

[dhomeier]

Thanks for the detailed output! I was a bit confused initially by the IOError turned into an OSError, but it looks like this was done within Astropy registry, so it is a bit beyond our control. Good to hear that it is working with your loader, though. I think the docs implicitly assume that you are substituting my-format with the respective name you've chosen, but they are unfortunately not very up-to-date on a number of points.

To track the tabular-fits problem further down, could you try it as

test_spec = Spectrum1D.read("/Users/krisstern/Desktop/ic418_edited_integrated_spectrum_cube1.fits", format="tabular-fits,
                        column_mapping={'flux': ('flux', 'ST'), 'spectral_axis': ('spectral_axis', 'AA')})

I suspect the problem is the flux unit you used in writing the spectrum to file – what unit is 'ST'?The loader expects some flux unit that can be converted to Jy (but can also be flux per wavelength interval); if your spectrum is using an unrecognised unit, that might also cause hiccups in other parts of specutils.

[nmearl]

For example, I am getting the following error when trying to read in/load a Spectrum1D object saved using the @custom_writer("fits-writer") approach as shown at https://specutils.readthedocs.io/en/stable/custom_loading.html#creating-a-custom-writer if I approach the task naively:

I'd like to point out that you're attempting to read in the written file using Spectrum1D.read('/Users/krisstern/Desktop/ic418_edited_integrated_spectrum_cube1.fits', format='fits') -- fits is not a defined reader format, which is why you're seeing the error. You'd have to use either tabular-fits or wcs1d-fits.

[dhomeier]

I'd like to point out that you're attempting to read in the written file using Spectrum1D.read('/Users/krisstern/Desktop/ic418_edited_integrated_spectrum_cube1.fits', format='fits') -- fits is not a defined reader format, which is why you're seeing the error. You'd have to use either tabular-fits or wcs1d-fits.

But the successful attempt I assume was using fits-format, which is defined in the custom_loader as advised in the documentation. Anyhow those docs should urgently be revised to first pointing the user to one of the builtin loaders above instead of leading them directly to step two with writing their own custom readers and writers.

[nmearl]

I suspect the problem is the flux unit you used in writing the spectrum to file – what unit is 'ST'?The loader expects some flux unit that can be converted to Jy (but can also be flux per wavelength interval); if your spectrum is using an unrecognised unit, that might also cause hiccups in other parts of specutils.

ST is a recognized astropy unit. This may very well be a bug if it's not being identified as such. I.e. the following works fine:

In [15]: wavelengths = np.arange(10) * u.AA
    ...: flux = np.random.sample(10) * u.Unit("ST")
    ...: flux.to("Jy", equivalencies=u.spectral_density(wavelengths))
Out[15]:
<Quantity [0.00000000e+00, 4.83452598e-05, 1.07028980e-04, 8.17246872e-04,
           1.13296289e-03, 2.13623159e-03, 1.80384816e-03, 5.82781596e-03,
           6.39521064e-03, 9.72076522e-03] Jy>

[dhomeier]

    flux.to("Jy", equivalencies=u.spectral_density(wavelengths))

Right, that is almost exactly what _find_spectral_column is doing – I got confused by attempting this with a spectrum with an incompatible uncertainty unit.

@kakirastern in your original post your custom writer just created the Table from spectrum.flux.value and spectrum.spectral_axis.value, without units. Could you post the FITS header of the saved file to check which spectral and flux units are defined there?

[kakirastern]

@dhomeier Essentially I have something like the following for the WCS (sorry for seemingly going on a tangent) and the meta (with the header info):

>>> spec.wcs
<WCS(output_frame=SpectralFrame, input_frame=CoordinateFrame, forward_transform=Model: SpectralTabular1D
N_inputs: 1
N_outputs: 1
Parameters: 
  points: (array([   0,    1,    2, ..., 6483, 6484, 6485]),)
  lookup_table: [3500.         3500.77431644 3501.54863288 ... 6999.12060362 6999.55980181
 6999.999     ] Angstrom
  method: linear
  fill_value: nan
  bounds_error: True)>

>>> spec.meta
{'header': SIMPLE  =                    T / conforms to FITS standard                      
 BITPIX  =                    8 / array data type                                
 NAXIS   =                    0 / number of array dimensions                     
 EXTEND  =                    T                                                  }

[dhomeier]

Thanks, I specifically meant the full header from fits.open()[1].header - looks like your reader already stripped the relevant TUNIT* cards.

[kakirastern]

The header I got with the command fits.open('some_file_name.fits')[1].header is as follows (which may not be so important as the TUNIT* cards are indeed already stripped like you have suggested):

XTENSION= 'BINTABLE'           / binary table extension                         
BITPIX  =                    8 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                   12 / length of dimension 1                          
NAXIS2  =                 6486 / length of dimension 2                          
PCOUNT  =                    0 / number of group parameters                     
GCOUNT  =                    1 / number of groups                               
TFIELDS =                    2 / number of table fields                         
TTYPE1  = 'spectral_axis'                                                       
TFORM1  = 'D       '                                                            
TTYPE2  = 'flux    '                                                            
TFORM2  = 'E       '

[dhomeier]

@kakirastern no, fits.open does not remove any header cards! So you wrote the spectrum (presumably following the outdated docs) without any unit information, and this would cause generic_spectrum_from_table to fail as it needs to know how to transform between F_nu and F_lambda flux scales. tabular-fits should still work for reading with the column_mapping={'flux': ('flux', 'ST'), 'spectral_axis': ('spectral_axis', 'AA')} option, but I would rather recommend to update your custom writer to write the full quantities – it should suffice to modify it to

    flux = spectrum.flux
    disp = spectrum.spectral_axis

[kakirastern]

Sorry @dhomeier. Somewhere along the line the flux and disp were stripped of the units while I was following the outdated docs. So I propose me opening a PR to update the docs to clear the confusion.

But when I modified the code into

>>> flux = spectrum.flux
>>> disp = spectrum.spectral_axis

before I attempt to save the Spectrum1D object to file, I got the following error with traceback:

---------------------------------------------------------------------------
UnitScaleError                            Traceback (most recent call last)
/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/io/fits/convenience.py in table_to_hdu(table, character_as_bytes)
    517             try:
--> 518                 col.unit = unit.to_string(format='fits')
    519             except UnitScaleError:

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/units/core.py in to_string(self, format)
    604         f = unit_format.get_format(format)
--> 605         return f.to_string(self)
    606 

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/units/format/fits.py in to_string(cls, unit)
    123                     "that are not powers of 10.  Multiply your data by "
--> 124                     "{:e}.".format(unit.scale))
    125             elif unit.scale != 1.0:

UnitScaleError: The FITS unit format is not able to represent scales that are not powers of 10.  Multiply your data by 3.630781e-09.

During handling of the above exception, another exception occurred:

UnitScaleError                            Traceback (most recent call last)
<ipython-input-147-4eb030df73d7> in <module>
----> 1 new_spec.write('/Users/krisstern/Desktop/ic418_edited_integrated_spectrum_junk.fits', format='fits-writer-1')

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/nddata/mixins/ndio.py in __call__(self, *args, **kwargs)
     97 
     98     def __call__(self, *args, **kwargs):
---> 99         registry.write(self._instance, *args, **kwargs)
    100 
    101 

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/io/registry.py in write(data, format, *args, **kwargs)
    564 
    565     writer = get_writer(format, data.__class__)
--> 566     writer(data, *args, **kwargs)
    567 
    568 

<ipython-input-146-64a1ec58dfd6> in generic_fits(spectrum, file_name, **kwargs)
     11                )
     12 
---> 13     tab.write(file_name, format="fits")

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/table/connect.py in __call__(self, serialize_method, *args, **kwargs)
    111         instance = self._instance
    112         with serialize_method_as(instance, serialize_method):
--> 113             registry.write(instance, *args, **kwargs)

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/io/registry.py in write(data, format, *args, **kwargs)
    564 
    565     writer = get_writer(format, data.__class__)
--> 566     writer(data, *args, **kwargs)
    567 
    568 

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/io/fits/connect.py in write_table_fits(input, output, overwrite)
    402     input = _encode_mixins(input)
    403 
--> 404     table_hdu = table_to_hdu(input, character_as_bytes=True)
    405 
    406     # Check if output file already exists

/opt/anaconda3/envs/astro/lib/python3.7/site-packages/astropy/io/fits/convenience.py in table_to_hdu(table, character_as_bytes)
    520                 scale = unit.scale
    521                 raise UnitScaleError(
--> 522                     f"The column '{col.name}' could not be stored in FITS "
    523                     f"format because it has a scale '({str(scale)})' that "
    524                     f"is not recognized by the FITS standard. Either scale "

UnitScaleError: The column 'flux' could not be stored in FITS format because it has a scale '(1.0)' that is not recognized by the FITS standard. Either scale the data or change the units.

I suspect this is particular to my case since I am using a non-typical flux scale.

[dhomeier]

Thanks for following up on this @kakirastern – now Astropy io.fits really does not let you get any further. Although ST is an accepted Astropy unit, it is not recognised as a valid unit by the FITS standard (see also https://github.com/astropy/astropy/issues/2933), and is thus rejected:

>>> u.Unit('ST').get_format_name('fits')
'ST'
>>> u.format.fits.Fits._validate_unit(u.Unit('ST').get_format_name('fits'))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/derek/lib/python3.8/site-packages/astropy/units/format/fits.py", line 86, in _validate_unit
    raise ValueError(
ValueError: Unit 'ST' not supported by the FITS standard.

although it tries its best with

u.format.utils.decompose_to_known_units(u.Unit('ST'), u.format.fits.Fits._get_unit_name)
Unit("3.63078e-09 erg / (Angstrom cm2 s)")

but since this can not be represented with a powers-of-ten prefix, it gives you the slightly more obscure error above.

Seems valid as ST is not listed in Table 4 of the standard paper. Interestingly io.fits did accept it as BUNIT for writing in IMAGE format, though.

So you would have to resort in your custom loader to rescaling the actual flux data with something like flux = spectrum.flux.to('W/(m2 nm)') (or any other F_lambda unit you fancy).

Note to self: tabular_fits_writer might check this earlier on to provide a more helpful error message.

[kakirastern]

Hi @dhomeier Thank you for the info! I finally figured out a solution and replaced the astropy flux unit u.ST with the synphot unit unit.FLAM (or Unit("flam") in the context of specutils). And it worked both for the custom writer and custom loader in specutils for writing and reading the concerned FITS files. As promised I will follow up with some PR(s) to update the documentation to clarify things a bit, probably over the coming weekend.

[dhomeier]

unit u.ST with the synphot unit unit.FLAM (or Unit("flam") in the context of specutils). And it worked both for the custom writer and custom loader in specutils for writing and reading the concerned FITS files. As promised I will follow up with some PR(s) to update the documentation to clarify things a bit, probably over the coming weekend.

Good to hear, although I was not able to reproduce the solution using u.Unit("flam") - did you define it manually or via one of the format options? I am a bit surprised that this is accepted by io.fits, as I cannot find it in the lists of bases or simple_units, nor is it listed in the standard. The default loaders for tabular and wcs1d reject it as well.

Anyway, does "FLAM" have a defined physical scale? Otherwise at least for the docs update I'd vote for recommending conversion to a standard unit, so one could always convert back to the correct values in u.ST.

[kakirastern]

Oh sorry, I meant I have used Unit("flam") in the custom data loader as follows:

from astropy.io import fits
from astropy.nddata import StdDevUncertainty
from astropy.table import Table
from astropy.units import Unit
from astropy.wcs import WCS

from specutils.io.registers import data_loader
from specutils import Spectrum1D

# Define an optional identifier. If made specific enough, this circumvents the
# need to add ``format="fits-format"`` in the ``Spectrum1D.read`` call.

def identify_generic_fits(origin, *args, **kwargs):
    return (isinstance(args[0], str) and
            os.path.splitext(args[0].lower())[1] == '.fits')

@data_loader("fits-format", identifier=identify_generic_fits,
             extensions=['fits'])
def generic_fits(file_name, **kwargs):
    with fits.open(file_name, **kwargs) as hdulist:
        header = hdulist[0].header

        tab = Table.read(file_name)

        meta = {'header': header}
        wcs = WCS(hdulist[0].header)
        data = tab["flux"] * Unit("flam")
        lamb = tab["spectral_axis"] * Unit("Angstrom")

    return Spectrum1D(flux=data, spectral_axis=lamb, wcs=wcs, meta=meta)

And, for the corresponding custom data writer I have used the following code snippet:

from astropy.table import Table
from specutils.io.registers import custom_writer

@custom_writer("fits-writer")
def generic_fits(spectrum, file_name, **kwargs):
    flux = spectrum.flux
    disp = spectrum.spectral_axis
    meta = spectrum.meta

    tab = Table([disp, flux], names=("spectral_axis", "flux"), meta=meta)

    tab.write(file_name, format="fits")

Also, the synphot.units.FLAM unit variable corresponds toerg * s^{−1} * cm^{−2} * Angstrom^{−1}.

[dhomeier]

Also, the synphot.units.FLAM unit variable corresponds toerg * s^{−1} * cm^{−2} * Angstrom^{−1}.

OK, so basically it's sort of an alias for u.Unit('erg / (AA cm2 s)'). Somehow I still cannot reproduce your usage, though, I get

>>> from astropy.units import Unit
>>> Unit("flam")
Traceback (most recent call last):
  File "/Users/derek/lib/python3.8/site-packages/astropy/units/format/generic.py", line 566, in _do_parse
    return cls._parse_unit(s, detailed_exception=False)
  File "/Users/derek/lib/python3.8/site-packages/astropy/units/format/generic.py", line 487, in _parse_unit
    raise ValueError()
ValueError

[loads of further traceback]

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/derek/lib/python3.8/site-packages/astropy/units/core.py", line 1878, in __call__
    raise ValueError(msg)
ValueError: 'flam' did not parse as unit: At col 0, flam is not a valid unit. Did you mean flm? If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see http://docs.astropy.org/en/latest/units/combining_and_defining.html

and I also cannot find it defined in any of the special unit modules like cds or imperial either. Are you sure you did not run some u.add_enabled_units(..) code prior to this?

[kakirastern]

I checked and it works for me when I do the following:

>>> from astropy.units import Unit
>>> from synphot import units
>>> Unit("flam")
flam

So you will also need to import synphot for my trick to work.

[kakirastern]

PR #641 opened to clarify the custom loader example confusion raised in this issue previously.

[kakirastern]

I think this issue is related and it may be useful to record it here: https://github.com/spacetelescope/synphot_refactor/issues/156. It is about separating out the units module of synphot, which is supported by the HST Exposure Time Calculator, such that all HST units defined in the module are guaranteed to work. It explains why the synphot "FLAM" unit worked in the case given while astropy's equivalent "ST" unit failed. Even if we do not end up to include this detail it would be good to have some reference here.

[kelle]

I think this issue can be closed.