<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: [![Plotted spectrum of HD001461 from the MILES library](https://i.sstatic.net/s3Lps.png)](https://i.sstatic.net/s3Lps.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](https://synphot.readthedocs.io/en/latest/) and [specutils](https://specutils.readthedocs.io/en/stable/index.html) 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})

```


---

