使用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 验证，两者都是超可靠的）。如果我将速度约定更改为无线电，结果会按预期更改。

但是对于AGES data sets（例如VC2），我得到的值大不相同。我预计范围 -2277 - 20108 km/s；我实际得到的是 -2350 - 19370，在上端时速超过 700 公里/秒！有趣的是，如果我不使用 .with_spectral_unit，即如果我只是这样做：

vcube.spectral_axis

...然后我得到正确的结果。所以这与单位转换有关，但我不知道是什么。我尝试将速度绘制为通道的函数。与正确速度的差异遵循抛物线，但最低差异不在参考通道处。

我唯一的怀疑是它可能与立方体的网格化方式有关。 AGES 网格使通道在频率上具有恒定大小，因此每个通道的速度宽度略有变化。我相信 THINGS 使用的是恒定的速度间隔。那么，spectral_cube 可以处理必要的转换，还是我找错了树？

Answer 1

好的，经过一个星期的困惑，我找到了解决方案！

问题确实是网格化。我尝试过的每个非 AGES 立方体都没有与天体或光谱立方体有关的光谱坐标问题。我不认为将数据网格化以具有恒定的频率通道宽度但不同的速度宽度会如此不寻常，但显然它是。真正让我感到困惑的是，如果没有对轴应用任何转换，那么这些值是正确的，但是如果为 with_spectral_unit 命令提供了任何关键字——即使只是为了将立方体保持在其本机单位中——那么这些值都是错误的。

在尝试了我能想到的每个标题的调整后，我发现了 miriad 任务velsw，它可以在不同的速度轴之间进行转换。直接设置我想要的速度约定（光学）不起作用，给光谱立方体转换带来了类似的 - 尽管不相同 - 错误。但是，任务说明会警告“非线性轴仅对参考点的一阶正确”。所以答案是转换成频率，在这个数据中是线性的。然后，spectral-cube 可以以近乎完美的精度处理转换回速度。

使用 velsw 是一种快速简便的解决方法，因为它只转换标题值（它不会重新网格化数据）。不利的一面是，首先必须转换为 miriad 自己的格式并恢复为适合（对于不熟悉 miriad 的任何人使用fits task）。我想应该可以直接使用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是我们要转换的像素坐标。