使用 spectral_cube 转换 FITS 文件的光谱轴
Converting the spectral axis of a FITS file using spectral_cube
我正在尝试使用 astropy/spectral_cube 提取 FITS 数据集的光谱轴。具体来说，我想将通道值转换为速度值，考虑到不同的 radio/optical 约定和正在检查的谱线的静止频率。
这适用于一个 FITS 文件，但不适用于另一个文件。我认为这是由于 header 的差异造成的，但我无法弄清楚关键差异是什么。
我的代码：
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' 是从 header 获得的。两个立方体都是 HI 数据。
然后我打印出光谱轴：
vcube.spectral_axis
（当愤怒使用它时，我会做额外的步骤，因为我需要在 non-integer 通道值和速度之间来回转换，例如：
cubewcs = vcube.wcs
wx, wy, wz = cubewcs.all_pix2world(150.0,125,0.0,0)
print(wx,wy,wz/vunit)
但我认为这不是主要问题的关键）
现在对于 THINGS data sets (e.g. NGC 628 的情况，它打印出我期望的确切速度值，范围从 588 到 735 km/s（已使用 kvis 和 miriad 验证，两者都是超级可靠）。如果我将速度约定更改为无线电，结果会按预期更改。
但是对于 AGES data sets (e.g. VC2)，我得到了截然不同的值。我希望范围为 -2277 - 20108 km/s；我实际得到的是 -2350 - 19370，在上端偏离 700 km/s！有趣的是，如果我不使用 .with_spectral_unit，即如果我只是使用 :
vcube.spectral_axis
...然后我得到了正确的结果。所以这与单位转换有关，但我不知道是什么。我尝试将速度绘制为通道的函数。与正确速度的差异遵循抛物线，但最低差异不在参考通道。
我唯一的怀疑是它可能与立方体的网格化方式有关。 AGES 网格使通道在频率上具有恒定大小，因此每个通道的速度宽度略有变化。我相信 THINGS 使用的是等速间隔。那么，spectral_cube 可以处理必要的转换，还是我找错树了？
好吧，在百思不得其解的一周之后，我找到了解决办法！
确实是网格问题。我尝试过的每个 non-AGES 立方体都没有天体或 spectral-cube 的光谱坐标问题。我不认为将数据网格化以具有恒定的频率通道宽度但变化的速度宽度会如此不寻常，但显然它是。真正让我感到困惑的是，如果没有对轴应用任何变换，那么值是正确的，但是如果向 with_spectral_unit 命令提供任何关键字——即使只是为了将立方体保持在其原始单位——那么值都错了。
在尝试了所有我能想到的header调整后，我找到了可以在不同速度轴之间转换的miriad任务velsw。直接设置我想要的速度约定（光学）不起作用，给出了与 spectral-cube 转换类似但不完全相同的错误。但是，任务说明给出警告“non-linear 轴仅在参考点处对 first-order 正确”。所以答案是转换为频率，在这个数据中频率是线性的。然后 spectral-cube 可以以 near-as-dammit 完美的精度处理回速度的转换。
使用 velsw 是一种快速简便的解决方法，因为它只转换 header 值（它不会重新网格化数据）。不利之处在于，首先必须转换为 miriad 自己的格式，然后再恢复为适合的格式（对不熟悉 miriad 的人使用 fits task）。我想应该可以直接使用 spectral-cube 来转换 header 值来跳过这一步，但是如果我不知道该怎么做，我会 post 一个单独的问题.
编辑：使用 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 文件的 header 值：
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]
我们现在再次读入 spectral-cube :
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是我们要转换的像素坐标。
我正在尝试使用 astropy/spectral_cube 提取 FITS 数据集的光谱轴。具体来说，我想将通道值转换为速度值，考虑到不同的 radio/optical 约定和正在检查的谱线的静止频率。 这适用于一个 FITS 文件，但不适用于另一个文件。我认为这是由于 header 的差异造成的，但我无法弄清楚关键差异是什么。
我的代码：
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' 是从 header 获得的。两个立方体都是 HI 数据。
然后我打印出光谱轴：
vcube.spectral_axis
（当愤怒使用它时，我会做额外的步骤，因为我需要在 non-integer 通道值和速度之间来回转换，例如：
cubewcs = vcube.wcs
wx, wy, wz = cubewcs.all_pix2world(150.0,125,0.0,0)
print(wx,wy,wz/vunit)
但我认为这不是主要问题的关键）
现在对于 THINGS data sets (e.g. NGC 628 的情况，它打印出我期望的确切速度值，范围从 588 到 735 km/s（已使用 kvis 和 miriad 验证，两者都是超级可靠）。如果我将速度约定更改为无线电，结果会按预期更改。
但是对于 AGES data sets (e.g. VC2)，我得到了截然不同的值。我希望范围为 -2277 - 20108 km/s；我实际得到的是 -2350 - 19370，在上端偏离 700 km/s！有趣的是，如果我不使用 .with_spectral_unit，即如果我只是使用 :
vcube.spectral_axis
...然后我得到了正确的结果。所以这与单位转换有关，但我不知道是什么。我尝试将速度绘制为通道的函数。与正确速度的差异遵循抛物线，但最低差异不在参考通道。
我唯一的怀疑是它可能与立方体的网格化方式有关。 AGES 网格使通道在频率上具有恒定大小，因此每个通道的速度宽度略有变化。我相信 THINGS 使用的是等速间隔。那么，spectral_cube 可以处理必要的转换，还是我找错树了？
好吧，在百思不得其解的一周之后，我找到了解决办法！
确实是网格问题。我尝试过的每个 non-AGES 立方体都没有天体或 spectral-cube 的光谱坐标问题。我不认为将数据网格化以具有恒定的频率通道宽度但变化的速度宽度会如此不寻常，但显然它是。真正让我感到困惑的是，如果没有对轴应用任何变换，那么值是正确的，但是如果向 with_spectral_unit 命令提供任何关键字——即使只是为了将立方体保持在其原始单位——那么值都错了。
在尝试了所有我能想到的header调整后，我找到了可以在不同速度轴之间转换的miriad任务velsw。直接设置我想要的速度约定（光学）不起作用，给出了与 spectral-cube 转换类似但不完全相同的错误。但是，任务说明给出警告“non-linear 轴仅在参考点处对 first-order 正确”。所以答案是转换为频率，在这个数据中频率是线性的。然后 spectral-cube 可以以 near-as-dammit 完美的精度处理回速度的转换。
使用 velsw 是一种快速简便的解决方法，因为它只转换 header 值（它不会重新网格化数据）。不利之处在于，首先必须转换为 miriad 自己的格式，然后再恢复为适合的格式（对不熟悉 miriad 的人使用 fits task）。我想应该可以直接使用 spectral-cube 来转换 header 值来跳过这一步，但是如果我不知道该怎么做，我会 post 一个单独的问题.
编辑：使用 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 文件的 header 值：
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]
我们现在再次读入 spectral-cube :
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是我们要转换的像素坐标。