How to read in IRAF multispec spectra?

Question

I have a spectrum in a fits file that I generated with Iraf. The wavelength axis is encoded in the header as:

WAT0_001= 'system=multispec'
WAT1_001= 'wtype=multispec label=Wavelength units=angstroms'
WAT2_001= 'wtype=multispec spec1 = "1 1 2 1. 2.1919422441886 4200 0. 452.53 471'
WAT3_001= 'wtype=linear'
WAT2_002= '.60 1. 0. 3 3 1. 4200.00000000001 1313.88904209266 1365.65012876239 '
WAT2_003= '1422.67911152069 1479.0560707956 1535.24082980747 1584.94609332243'

Is there an easy way to load this into python?

Ary

Do you have Pyfits or astropy installed? – Felske

@Tom. We discussed this on a Facebook group for astronomers and found the answer ourself. We wanted to open a post and give the answer for future references. – Overshine

Answer 1

I have been using this code, which was given to me by Rick White. However, the specutils package is probably the better way to do it:

from specutils.io import read_fits
spectra_list = read_fits.read_fits_spectrum1d('mymultispec.fits')
print spectra_list[0]
   Spectrum1D([ 338.06109619,  395.59234619,  326.0012207 , ...,
         660.0098877 ,  686.54498291,  689.58374023])

print spectra_list[1].dispersion
   <Quantity [ 8293.44875263, 8293.40459999, 8293.36044556,...,
        8166.53073537, 8166.48250242, 8166.43426803] Angstrom>

Ary

Can you mark this as the correct answer? This way, a newcomer will know that this works. – Uniaxial

I agree, specutils is the way to go. – Whoreson

specutils has now stopped supporting the multispec format as they restructured their code and haven't so far rewritten the multispec part. They may add it back in the future. – Educational

I would like to point out that specutils has been updated and the method explained here is not working anymore. – Argo

Answer 2

Follow-up on the answer from @kgully - and the related discussion: specutils supports reading multispec files again (has for several years, even though it wasn't working for a bit in between these answers). The relevant code snippet (for any IRAF-formatted file) is now:

>>> from specutils import Spectrum1D
>>> spec = Spectrum1D.read('mymultispec.fits')
>>> print(spec)
Spectrum1D (length=762078)
flux:             [ 0.4408 adu, ..., -1.0 adu ],  mean=0.74146 adu
spectral axis:    [ 3726.7 Angstrom, ..., 9300.0 Angstrom ],  
mean=6094.4 Angstrom

and you replace Spectrum1D with SpectrumCollection if your multispec file actually has multiple spectra.

Ozzy