working with stellar spectra in fits format in python

Question

Hey I am new to working with astronomical data in Python. I wanted to start working with stellar spectra and I am having trouble with the data. To get a first look I just wanted to plot a spectra (flux over wavelength). I downloaded the Miles library (http://research.iac.es/proyecto/miles/pages/stellar-libraries/miles-library.php) and started working with a single fits file. My first steps where as follows :

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits

hdul = fits.open('s0013.fits')
data = hdul[0].data
h1 = hdul[0].header  

I read the header using print(repr(h1)) and know I wanted to plot the spectrum. The data has shape (1, 4367) and I am not sure how to proceed to achieve 2 arrays one with the wavelength and one with the flux to plot the data. I am sorry if this question is stupid but I can not figure it out.

Cheers

Answer 1

The information you need to recreate the wavelength array is in the World Coordinate System (WCS) of the header, specifically:

CRPIX1  =                 1.00
CRVAL1  =            3500.0000 / central wavelength of first pixel
CDELT1  =             0.900000 / linear dispersion (Angstrom/pixel)

which lists the starting/reference pixel of the wavelength array (1.0), the wavelength value at the start point (3500 angstroms (assumed)) and the step per pixel (0.9 Angstrom/pixel). To read this information, it is best to use a WCS library rather than trying to interpret them directly as they can be more complicated and there are many subformats of FITS WCS.

Fortunately astropy has a module to make this easy (starting from your code above and extending it):

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS

hdul = fits.open('s0013.fits')
data = hdul[0].data
h1 = hdul[0].header
obj_name = h1.get('OBJECT', 'Unknown')

flux = data[0]
w = WCS(h1, naxis=1, relax=False, fix=False)
lam = w.wcs_pix2world(np.arange(len(flux)), 0)[0]

plt.plot(lam, flux)
plt.ylim(0, )
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Normalized flux')
plt.savefig(obj_name + '.png')

This will produce the following plot:

If you want to do more extended manipulation of spectra, particularly with the bewildering variety of wavelength and flux units, it might be worth looking at synphot and specutils which build on Astropy and add more direct support for spectra beyond simple numpy arrays. For example, you could make a synphot SourceSpectrum from the above by doing:

from astropy import units as u
from synphot import units, SourceSpectrum
from synphot.spectrum import Empirical1D

source_spec = SourceSpectrum(Empirical1D, points=lam*u.AA, lookup_table=flux,
                   keep_neg=True, meta={'header': h1})