python - Write 3d Numpy array to FITS file with Astropy

Question

I have a 3D NumPy array (i.e. (10, 256, 256)) representing 256x256 images. I would like to write this array to a FITS file using astropy.io.fits so that I can open the file using ds9 -mecube and move through the frames. My attempt is shown below

export_array = numpy.array(images) #Create an array from a list of images
print export_array.shape ## (10, 256, 256)

hdu = fits.PrimaryHDU(export_array)
hdulist = fits.HDUList([hdu])
hdulist.writeto(out_file_name)
hdulist.close()

This will give me a FITS file which does in fact contain the 3D array. However if I open with ds9 -mecube I can only see the first image. Is there anyway to create the FITS file with this functionality using astropy.io.fits? Or is there perhaps some functionality with ds9 that I am missing?

Answer 1

I don't use ds9, but apparently the -mecube option means "multi-extension cube". The docs say "Load a multi-extension FITS file as a data cube. You're just writing a single array as a data cube. To write it as a multi-extension FITS you might do something like:

hdul = fits.HDUList()
hdul.append(fits.PrimaryHDU())

for img in export_array:
    hdul.append(fits.ImageHDU(data=img))

hdul.writeto('output.fits')

(You don't need to call hdul.close()--that only does anything if the HDUList was loaded from an existing file and you want to close the underlying file object; it has no effect for an HDUList created from scratch in memory).

I don't know exactly what ds9 is expecting for loading a multi-extension FITS file as a data cube--this isn't any specific FITS convention and the docs arent't clear. But it's probably something like that.

All that said, according to the ds9 docs you don't need to use this at all. If you don't use the -mecube option it will just read the 3D array in the primary HDU as a data cube.