Collapsing / Flattening a FITS data cube in python

Question

I've looked all over the place and am not finding a solution to this issue. I feel like it should be fairly straightforward, but we'll see.

I have a .FITS format data cube and I need to collapse it into a 2D FITS image. The data cube has two spacial dimensions and one spectral/velocity dimension.

Just looking for a simple python routine to load in the cube and flatten all these layers (i.e. integrate them along the spectral/velocity axis). Thanks for any help.

Answer 1

OK, this seems to work:

import pyfits
import numpy as np

hdulist = pyfits.open(filename)  
header = hdulist[0].header  
data = hdulist[0].data   
data = np.nan_to_num(data)   
new_data = data[0]

for i in range(1,84):                #this depends on number of layers or pages
    new_data += data[i]

hdu = pyfits.PrimaryHDU(new_data)
hdu.writeto(new_filename)

One problem with this routine is that WCS coordinates (which are attached to the original data cube) are lost during this conversion.

Answer 2

This tutorial on pyfits is a little old, but still basically correct. The key is that the output of opening a FITS cube with pyfits (or astropy.io.fits) is that you have a 3 dimensional numpy array.

import pyfits
# if you are using astropy then for this example
# from astropy.io import fits as pyfits

data_cube, header_data_cube = pyfits.getdata("data_cube.fits", 0, header=True)
data_cube.shape
# (Z, X, Y) 

You then have to decided how to flatten/integrate cube along the Z axis, and there are plenty of resources out there to help you decide the right (hopefully based in some analysis framework) to do that.

Answer 3

This is a bit of an old question, but spectral-cube now provides a better solution for this.

Example, based on Teachey's answer:

from spectral_cube import SpectralCube
cube = SpectralCube.read(filename)

summed_image = cube.sum(axis=0)
summed_image.hdu.writeto(new_filename)