I'm trying to extract the spectral axis of a FITS data set using astropy/spectral_cube. Specifically I want to convert the channel values into velocity values, accounting for the different radio/optical conventions and rest frequency of the spectral line being examined. This is working just fine for one FITS file but not for another. I presume this is due to a difference in the headers, but I can't figure out what the critical difference is.
My code :
from astropy.io import fits as pyfits
from astropy.wcs import WCS
from astropy import units
import spectral_cube
from spectral_cube import SpectralCube
fitsfile = pyfits.open(fitsfilename)
scube = SpectralCube.read(fitsfile)
vcube = scube.with_spectral_unit(units.km / units.s, velocity_convention='optical', rest_value=1420405750.0 * units.Hz)
Where 'fitsfilename' obviously sets the name of the file to load. The 'rest_value' is obtained from the header. Both cubes are HI data.
I then print out the spectral axis with :
vcube.spectral_axis
(When using this in anger, I would do additional steps, as I need to convert back and forth between non-integer channel values and velocities, e.g. :
cubewcs = vcube.wcs
wx, wy, wz = cubewcs.all_pix2world(150.0,125,0.0,0)
print(wx,wy,wz/vunit)
But this is not crucial to the main problem, I think)
Now for the case of the THINGS data sets (e.g. NGC 628) this prints out the exact velocity values I'm expecting, ranging from 588 to 735 km/s (verified with kvis and miriad, both of which are ultra reliable). If I change the velocity convention to radio, the results change as expected.
But for the AGES data sets (e.g. VC2), I get substantially different values. I expect the range -2277 - 20108 km/s; what I actually get is -2350 - 19370, off by over 700 km/s at the upper end ! Interestingly if I don't use .with_spectral_unit, i.e. if I just do :
vcube.spectral_axis
... then I get the correct result. So it's something to do with the unit conversion, but I don't know what. I tried plotting the velocities as a function of channel. The differences from the correct velocity follows a parabola, but the lowest difference is NOT at the reference channel.
My only suspicion is that it may relate to how the cubes are gridded. AGES grids so that the channels are of constant size in frequency, so that the velocity width of each channel varies slightly. I believe THINGS uses a constant velocity interval instead. So, can spectral_cube handle the necessary conversion, or am I barking up the wrong tree ?
Okay, after a week of being utterly perplexed, I found a solution !
The problem was indeed the gridding. Every non-AGES cube I tried had no spectral coordinate problems with astropy or spectral-cube at all. I didn't think it would be so unusual to grid the data to have constant frequency channel width but varying velocity width, but apparently it is. What really baffled me was that if no transformation is applied to the axis then the values are correct, but if any keywords are provided to the with_spectral_unit command - even if only to keep the cube in its native units - then the values are all wrong.
After trying every adjustment to the header I could think of, I found the miriad task velsw which can convert between different velocity axes. Directly setting the velocity convention that I wanted (optical) did not work, giving a similar - though not identical - error to the spectral-cube conversion. However, the task instructions give a warning that "non-linear axes are only correct to first-order at the reference point". So the answer is to convert to frequency, which in this data is linear. Then spectral-cube can handle the conversion back to velocity with near-as-dammit perfect accuracy.
Using velsw is a fast and easy workaround as it only transforms the header values (it doesn't regrid the data). The downside is that one first has to convert to miriad's own format and back to fits (using the fits task for anyone unfamiliar with miriad). I imagine it should be possible to convert the header values directly using spectral-cube to skip this step, but I'll post a separate question if I can't figure out how to do that.
EDIT : The code to do this with spectral-cube is as follows. First we load the cube as usual and convert it to frequency :
from astropy.io import fits as pyfits
from astropy.wcs import WCS
fitsfile = pyfits.open(fitsfilename)
scube = SpectralCube.read(fitsfile)
fcube = scube.with_spectral_unit(units.Hz)
Then we convert the header values of the FITS file in memory, using the values generated for the frequency cube :
fitsfile[0].header['CRVAL3'] = fcube.wcs.wcs.crval[2]
fitsfile[0].header['CDELT3'] = fcube.wcs.wcs.cdelt[2]
fitsfile[0].header['CTYPE3'] = fcube.wcs.wcs.ctype[2]
We now again read into spectral-cube :
scube = SpectralCube.read(fitsfile)
And now we can convert as expected, e.g. :
vcube = scube.with_spectral_unit(units.MHz, velocity_convention='optical', rest_value=rfq_value * units.Hz)
cubewcs = vcube.wcs
wx, wy, wz = cubewcs.all_pix2world(cx,cy,cz,0)
Where cx, cy and cz are the pixel coordinates we want to convert.