使用光谱立方体

问题描述

我正在尝试使用 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,wz/vunit)

但这对主要问题并不重要,我认为)

现在对于 THINGS data sets(例如 NGC 628)的情况,这会打印出我期望的准确速度值,范围从 588 到 735 公里/秒(用 kvis 和 miriad 验证,两者都是其中非常可靠)。如果我将速度约定更改为无线电,结果会按预期改变。

但是对于 AGES data sets(例如 VC2),我得到了截然不同的值。我预计范围为 -2277 - 20108 公里/秒;我实际得到的是 -2350 - 19370,在上限处超过 700 公里/秒!有趣的是,如果我不使用 .with_spectral_unit,即如果我只是这样做:

vcube.spectral_axis

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

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

解决方法

好吧,经过一个星期的困惑,我找到了解决方案!

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

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

使用 velsw 是一种快速简便的解决方法,因为它只转换标头值(它不会重新网格化数据)。缺点是首先必须转换为 miriad 自己的格式并返回到 fits(对于不熟悉 miriad 的任何人都使用 fits task)。我想应该可以直接使用光谱立方体转换标题值来跳过这一步,但如果我不知道如何做到这一点,我会发布一个单独的问题。