2

I am new to python and FITS image files, as such I am running into issues. I have two FITS files; the first FITS file is pixels/counts and the second FITS file (calibration file) is pixels/wavelength. I need to convert pixels/counts into wavelength/counts. Once this is done, I need to output wavelength/counts as a new FITS file for further analysis. So far I have managed to array the required data as shown in the code below.

import numpy as np
from astropy.io import fits

# read the images 

image_file = ("run_1.fits")
image_calibration = ("cali_1.fits")

hdr = fits.getheader(image_file)
hdr_c = fits.getheader(image_calibration)

#   print headers
sp = fits.open(image_file)
print('\n\nHeader of the spectrum :\n\n', sp[0].header, '\n\n')

sp_c = fits.open(image_calibration)
print('\n\nHeader of the spectrum :\n\n', sp_c[0].header, '\n\n')

# generation of arrays with the wavelengths and counts
count = np.array(sp[0].data)
wave = np.array(sp_c[0].data)

I do not understand how to save two separate arrays into one FITS file. I tried an alternative approach by creating list as shown in this code

file_list = fits.open(image_file)
calibration_list = fits.open(image_calibration)


image_data = file_list[0].data
calibration_data = calibration_list[0].data


# make a list to hold images
img_list = []
img_list.append(image_data)
img_list.append(calibration_data)

# list to numpy array
img_array = np.array(img_list)

# save the array as fits - image cube
fits.writeto('mycube.fits', img_array)

However I could only save as a cube, which is not correct because I just need wavelength and counts data. Also, I lost all the headers in the newly created FITS file. To say I am lost is an understatement! Could someone point me in the right direction please? Thank you.

I am still working on this problem. I have now managed (I think) to produce a FITS file containing the wavelength and counts using this website:

https://www.mubdirahman.com/assets/lecture-3---numerical-manipulation-ii.pdf

This is my code:

# Making a Primary HDU (required):
primaryhdu = fits.PrimaryHDU(flux) # Makes a header # or if you have a header that you’ve created: primaryhdu = fits.PrimaryHDU(arr1, header=head1)
# If you have additional extensions:
secondhdu = fits.ImageHDU(wave)
# Making a new HDU List:
hdulist1 = fits.HDUList([primaryhdu, secondhdu])
# Writing the file:
hdulist1.writeto("filename.fits", overwrite=True)

image = ("filename.fits")
hdr = fits.open(image)
image_data = hdr[0].data
wave_data = hdr[1].data

I am sure this is not the correct format for wavelength/counts. I need both wavelength and counts to be contained in hdr[0].data

pixel_image
  • 23
  • 4

1 Answers1

0

If you are working with spectral data, it might be useful to look into specutils which is designed for common tasks associated with reading/writing/manipulating spectra.

It's common to store spectral data in FITS files using tables, rather than images. For example you can create a table containing wavelength, flux, and counts columns, and include the associated units in the column metadata.

The docs include an example on how to create a generic "FITS table" writer with wavelength and flux columns. You could start from this example and modify it to suit your exact needs (which can vary quite a bit from case to case, which is probably why a "generic" FITS writer is not built-in).

You might also be able to use the fits-wcs1d format.

If you prefer not to use specutils, that example still might be useful as it demonstrates how to create an Astropy Table from your data and output it to a well-formatted FITS file.

Iguananaut
  • 15,675
  • 4
  • 43
  • 50