使用spectral_cube转换FITS文件的光谱轴

Question

我正在尝试使用astropy/ spectral _cube提取FITS数据集的谱轴。具体地说，我希望将通道值转换为速度值，以考虑正在检查的谱线的不同无线电/光学约定和静止频率。这对于一个FITS文件工作得很好，但对于另一个文件就不行了。我认为这是由于头部的不同，但我不知道关键的不同是什么。

我的代码：

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) 

其中'fitsfilename‘显然设置了要加载的文件的名称。“rest_value”是从报头中获取的。两个多维数据集都是HI数据。

然后，我使用以下命令打印光谱轴：

vcube.spectral_axis

(在愤怒中使用时，我会执行额外的步骤，因为我需要在非整数通道值和速度之间来回转换，例如：

cubewcs = vcube.wcs
wx, wy, wz = cubewcs.all_pix2world(150.0,125,0.0,0)
print(wx,wy,wz/vunit)

但我认为这并不是主要问题的关键。)

现在，对于THINGS data sets (例如NGC 628)，这将打印出我期望的准确速度值，范围从588到735 km/s (通过kvis和miriad验证，这两个都是超级可靠的)。如果我将速度约定更改为radio，结果将按预期更改。

但是对于AGES data sets (例如VC2)，我得到的值有很大的不同。我期望的范围-2277 - 20108公里/秒；我实际得到的是-2350 - 19370，在上端超过700公里/秒！有趣的是，如果我不使用.with_spectral_unit，也就是说，我只使用：

vcube.spectral_axis

..。然后我就得到了正确的结果。所以这跟单位转换有关，但我不知道是什么。我试着把速度画成通道的函数。与正确速度的差异遵循抛物线，但最小的差异不在参考通道上。

我唯一的怀疑是，这可能与立方体的网格化方式有关。调整栅格的年龄，使通道在频率上具有恒定的大小，以便每个通道的速度宽度略有变化。我相信事物使用的是一个恒定的速度间隔。那么，spectral_cube可以处理必要的转换吗，或者我找错了人？

Answer 1

好了，在完全困惑了一周之后，我找到了一个解决方案！

问题确实出在网格上。我尝试过的每个非AGES立方体都没有星座或光谱立方体的光谱坐标问题。我不认为栅格数据具有恒定的频率通道宽度但变化的速度宽度是不寻常的，但很明显它是不寻常的。真正让我困惑的是，如果没有对轴应用任何转换，那么值就是正确的，但是如果给with_spectral_unit命令提供了任何关键字-即使只是为了保持立方体的原生单位-那么值都是错误的。

在尝试了所有我能想到的对头的调整后，我找到了miriad任务velsw，它可以在不同的速度轴之间转换。直接设置我想要的速度约定(光学)不起作用，给出了一个类似的-尽管不是完全相同的-光谱立方体转换的误差。然而，任务指令给出了一个警告，“非线性轴只对参考点的一阶正确”。所以答案是转换成频率，在这个数据中它是线性的。然后，spectral cube可以以近乎完美的精度处理转换回速度的过程。

使用velsw是一种快速而简单的变通方法，因为它只转换头文件的值(它不会重新格式化数据)。缺点是，人们必须首先转换成miriad自己的格式，然后再转换回fits (对于任何不熟悉miriad的人来说，都可以使用fits task )。我想应该可以直接使用spectral cube转换头文件的值来跳过这一步，但是如果我不知道怎么做，我会发一个单独的问题。

编辑:使用spectral cube执行此操作的代码如下所示。首先，我们像往常一样加载立方体，并将其转换为频率：

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)

然后，我们使用为频率立方体生成的值，在内存中转换FITS文件的头值：

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]

我们现在再一次读入光谱立方体：

scube = SpectralCube.read(fitsfile)

现在我们可以像预期的那样进行转换，例如：

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)

其中cx，cy和cz是我们想要转换的像素坐标。