0

I've been using matplotlib.pyplot to plot a spectrum from fits files in Python and getting the intensity versus pixel, but what I actually need is to convert the pixels to wavelength. I've seen similar questions that got me in the right path of what I need to do (e.g. similar question,RGB example) but I still feel lost in the process.

I have FITS files with wavelengths around 3500 and 6000 (A), in float32 format and dimensions (53165,).

So as I understand, I need to calibrate the pixel positions to wavelength. I have my rest wavelength header (RESW), my "step" wavelength header (STW) and I would need to get:

x = RESW + (number of pixels * STW)

and the plot it. Here is what I got in my code so far.

import os, glob
from glob import glob
from pylab import *
from astropy.io import ascii
import scipy.constants as constants
import matplotlib.pylab as plt
from astropy.io import fits

#add directory you have your files in
dir = '' 


#OPEN ALL FILES AND STORE THEM INTO A LIST
files= glob(dir + '*.fits')
for fi in (files):
    print(fi)
    name = fi[:-len('.fits')] #Remove '.fits' from the file name

    with fits.open(dir + fi) as hdu:
        hdu.info()
        data = hdu[0].data
        hdr = hdu[0].header #added to try2
        step = hdr['CDELT1'] #added to try2
        restw = hdr['CRVAL1'] #added to try2
        #step = fits.getheader('STW') #added to try
        #restw = fits.getheader('RESW') #added to try
        spectra = restw + (data * step) #added to try

plt.clf()
plt.plot(spectra)
plt.savefig(name + '.pdf')

I've tried using fits.getheader('') but I don't know where or how to put it because this way is not working right.

Could someone please help? Thanks in advance!

Ale
  • 1
  • 4
  • isn't it, "hdr = hdu[0].header" and then "step=hdr('STW')"? See: http://docs.astropy.org/en/stable/io/fits/usage/headers.html – Gus Jul 04 '18 at 13:49
  • Thank you for your help! So I just tried it with that (edit above) but its still not working. – Ale Jul 04 '18 at 14:33
  • what isn't working? from what I can tell you want to plot "x" (wavelength) versus "y" flux. your matplotlib call really should be `plt.plot(wave, flux)` where `flux=hdu[0].data` and `wave = restw + (numpy.arange(len(flux)) * step)`. `numpy.arange()` will return the array index of the pixel values for the spectrum flux, which you then scale with the step and zeropoint. i'm not sure I can help without more information on the spectrum. – Gus Jul 05 '18 at 12:08
  • It doesn't work as in I don't get anything back, even by just putting `plt.show()` instead of just saving the figure. I've tried adding `numpy.arange()` as you said, but still nothing :/ – Ale Jul 06 '18 at 14:20
  • I'm still not sure what you mean by "nothing" but I'll point out that in your loop you reassign the `spectra` value in `spectra = restw + (data * step)` on each loop, so your `plot()` call is only plotting the data from the last file. – Iguananaut Jul 07 '18 at 22:30
  • I don't get my pdfs ("nothing") like I did when I first run the script to plot the flux versus pixels. And yes, each loop needs to plot the data only from the last file, so that is okay. – Ale Jul 10 '18 at 13:58

0 Answers0