2

I have a bunch of fits files that can be read using the below script

from astropy.io import fits
hdu = fits.open('file.fits')
data = hdu[0].data

I am trying to make an image cube using the data read from multiple fits files. (An image cube is a 3D image that contains data from multiple fits file where the x and y axis is the 2D image dimension and the 3rd axis is time or frequency)

I believe it can be done usin spectral _cube module, however most of the documentation only talks about how to read an image cube and not how to make one using individual fits files.

So far I have tried the following script.

#In the below script data is a 3D numpy array
from spectral_cube import SpectralCube
cube = SpectralCube(data=data)
cube.write('new_cube.fits', format='fits')

However, the above script gives an error saying 3 arguments required while only 2 given.

DavidG
  • 24,279
  • 14
  • 89
  • 82
Mc Missile
  • 715
  • 11
  • 25
  • Which part gives that error? When posting problems with Python code please always post the exact exception message and as much of the traceback you think is relevant. – Iguananaut Feb 04 '19 at 18:53
  • According to the docs you also need to provide a [WCS object](https://spectral-cube.readthedocs.io/en/latest/creating_reading.html#direct-initialization) . Normally when reading a cube from an existing FITS file it will read the full WCS out of the FITS header. In this case you'll have to construct it yourself, providing WCS coordinates for each of your axes. The details of that depend largely on your application. Spectral WCS are outside my wheelhouse so maybe someone else can comment there. Looking more at the spectral_cube sources might help as well. – Iguananaut Feb 04 '19 at 19:03
  • @Iguananaut Thank you for looking into my problem. I figured out an alternative way to make cubes. Your intention to help is much appreciated :) – Mc Missile Feb 05 '19 at 07:06

2 Answers2

3

The simplest way to do this is to simply put the images you want to have in your cube into a numpy array and then save that array as a fits file. You can also save them into a numpy array directly, but appending lists is easier if you do it in a for loop, instead of doing it explicitly for each image, like I do it here.

import numpy as np
from astropy import fits

# Read the images you want to concatenate into a cube
img1 = fits.getdata('img1.fits')
img2 = fits.getdata('img2.fits')

# Make a list that will hold all your images
img_list = []
img_list.append(img1)
img_list.append(img2)

# Cast the list into a numpy array
img_array = np.array(img_list)

# Save the array as fits - it will save it as an image cube
fits.writeto('mycube.fits', img_array)
Iva Laginja
  • 202
  • 1
  • 10
  • Great! Very simple answer. I almost done with my app. Now I just have to deal with the header. Thanks a lot ! – edison Nov 03 '20 at 18:30
2

Apparently, one does not need to use the spectral_cube module to make image cubes. It can be much easily done using the AstroPy python module. Below is the script.

from astropy.io import fits
import numpy as np

cube = np.zeros((50,1000,1000)) #Here 1000x1000 is the dimension of the individual fits images and 50 is the third perpendicular axis(time/freq)

for i in range(50):
    hdu = fits.open('image' + str(i) + '.fits') #The fits images that you want to combine have the name string 'image' + str(i) + '.fits'
    data = hud[0].data[:,:]
    cube[i,:,:] = data

hdu_new = fits.PrimaryHDU(cube)
hdu_new.writeto('cube.fits')
ascendants
  • 2,123
  • 3
  • 11
  • 22
Mc Missile
  • 715
  • 11
  • 25