<p>The information you need to recreate the wavelength array is in the World Coordinate System (WCS) of the header, specifically:</p>
<pre><code>CRPIX1  =                 1.00
CRVAL1  =            3500.0000 / central wavelength of first pixel
CDELT1  =             0.900000 / linear dispersion (Angstrom/pixel)
</code></pre>
<p>which lists the starting/reference pixel of the wavelength array (<code>1.0</code>), the wavelength value at the start point (<code>3500</code> angstroms (assumed)) and the step per pixel (<code>0.9</code> 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.</p>
<p>Fortunately <code>astropy</code> has a module to make this easy (starting from your code above and extending it):</p>
<pre><code>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')
</code></pre>
<p>This will produce the following plot: <a href="https://i.sstatic.net/s3Lps.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/s3Lps.png" alt="Plotted spectrum of HD001461 from the MILES library" /></a></p>
<p>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 <a href="https://synphot.readthedocs.io/en/latest/" rel="nofollow noreferrer">synphot</a> and <a href="https://specutils.readthedocs.io/en/stable/index.html" rel="nofollow noreferrer">specutils</a> which build on Astropy and add more direct support for spectra beyond simple numpy arrays. For example, you could make a <em>synphot</em> <code>SourceSpectrum</code> from the above by doing:</p>
<pre><code>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})
</code></pre>

CC BY-SA 4.0
2021-04-26T13:19:02.107Z
2021-04-26T13:19:02.107Z
astrosnapper
<p>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 (<a href="http://research.iac.es/proyecto/miles/pages/stellar-libraries/miles-library.php" rel="nofollow noreferrer">http://research.iac.es/proyecto/miles/pages/stellar-libraries/miles-library.php</a>) and started working with a single fits file. My first steps where as follows :</p>
<pre><code>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  
</code></pre>
<p>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.</p>
<p>Cheers</p>

<data-analysis><spectra><python><fits>
working with stellar spectra in fits format in python
trynerror
# working with stellar spectra in fits format in python

(Asked by: username_0 on 2021-04-26)

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

**Tags:** data-analysis, spectra, python, fits


---

## 1 answers

### Answer by username_1 on 2021-04-26. Score: 2 | Accepted answer

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


---

