1

I've some problems and I could not find any answer to my problem.

I'm trying to create a datacube in python, where the three axis are (RA,DEC,z), that is 2 sky position and red shift. I think my code for generating the cube works, I define the cube as:

cube     = np.zeros([int(size_x),int(size_y),int(Nchannel)])

where x and y are pixel coordinates and the redshift is sliced in channels. Having this cube I'm filling it with intensity of some lines. At the end I define my .fits header as follows:

hdr = fits.Header()

hdr['EQUINOX'] = 2000
hdr['CRPIX1']  = round(size_ra*3600./pix_size/2.) 
hdr['CRPIX2']  = round(size_dec*3600./pix_size/2.)
hdr['CRPIX3']  = 0
hdr['CRVAL1']  = ra0
hdr['CRVAL2']  = dec0
hdr['CRVAL3']  = z_min
hdr['CD1_1']   =  pix_size/3600.
hdr['CD1_2']   =  0.
hdr['CD2_1']   =  0.
hdr['CD2_2']   =  pix_size/3600.
hdr['CTYPE1']  = "RA---TAN"
hdr['CTYPE2']  = "DEC--TAN"
hdr['CTYPE3']  = "Z"
hdr['BUNIT']   = "Jy/pixel"

fits.writeto('cube.fits',cube,hdr,overwrite=True)

And here is the problem, my cube.fits is in the "bad" direction. When I open it using ds9 the z-axis is not the redshift z... I'm suspecting a bad header, but where can I specify the axis in the fits header? Cheers

SF..MJ
  • 862
  • 8
  • 19
GreatMonark
  • 33
  • 3

1 Answers1

2

The axes are indeed inverted, FITS uses the Fortran convention (column-major order) whereas Python/Numpy uses the C convention (row-major order). http://docs.astropy.org/en/latest/io/fits/appendix/faq.html#what-convention-does-astropy-use-for-indexing-such-as-of-image-coordinates

So for your cube you need to define the axes as (z, y, x):

In [1]: import numpy as np

In [2]: from astropy.io import fits

In [3]: fits.ImageHDU(data=np.zeros((5,4,3))).header
Out[3]: 
XTENSION= 'IMAGE   '           / Image extension                                
BITPIX  =                  -64 / array data type                                
NAXIS   =                    3 / number of array dimensions                     
NAXIS1  =                    3                                                  
NAXIS2  =                    4                                                  
NAXIS3  =                    5                                                  
PCOUNT  =                    0 / number of parameters                           
GCOUNT  =                    1 / number of groups
saimn
  • 443
  • 2
  • 10
  • Thank you for the answer, I've thought it was possible to modify this convention in the header but it seems that it's not possible. – GreatMonark Mar 20 '18 at 13:17
  • The axis ordering convention is what it is--there's way to change that. But do please take note of the discussion in https://github.com/astropy/astropy/issues/3836 In general when storing arrays in FITS one doesn't not have to (and in fact almost always *must not*) manually set the `NAXIS*` keywords. Those will be generated appropriately from the array object. The array in Python land uses Python/C-like axis ordering, but that will be mapped to FITS ordering when storing the file. So if the axes in Python are `(x, y, z)` they will still be stored as `z -> NAXIS1` and so on. – Iguananaut Mar 20 '18 at 17:34
  • Likewise you must invert the axis numbers when writing WCS headers. The WCS interface will take care of this correctly as well. – Iguananaut Mar 20 '18 at 17:35