Working with stellar spectra in fits format in python

Working with stellar spectra in fits format in python

We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

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 ( and started working with a single fits file. My first steps where as follows :

import numpy as np import matplotlib.pyplot as plt from import fits hdul ='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.


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 (3500angstroms (assumed)) and the step per pixel (0.9Angstrom/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.

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

import numpy as np import matplotlib.pyplot as plt from import fits from astropy.wcs import WCS hdul ='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">

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 synphotSourceSpectrumfrom 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})