使用astropy的fit_wcs_from_points为FITS文件添加WCS的问题及解决

操作背景 #

通过pixel_to_world工具获取了5颗恒星的赤经(RA)、赤纬(DEC)，同时已掌握这些恒星在目标图像上的xy坐标，因此选用fit_wcs_from_points作为目标图像的WCS计算方案。

测试拟合代码 #

import numpy as np
from astropy.wcs.utils import fit_wcs_from_points
from astropy.coordinates import SkyCoord

stars = np.array([[1246.63, 1372.83, 1455.12, 1611.95, 1644.85],
               [1588.42, 1502.92, 1677.24, 1610.39, 1325.12]])

known_coords = SkyCoord([(125.66419083, -42.96809252), (125.67730695, -42.98209958),
                         (125.65082259, -42.9914015), (125.6611325, -43.01438513), (125.70471982, -43.01228167)],
                        frame="icrs", unit="deg")

fit_wcs_from_points(xy = stars, world_coords = known_coords, projection='TAN')

拟合输出结果 #

WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN'  'DEC--TAN'  
CRVAL : 125.67776105611648  -42.99124198597975  
CRPIX : 1451.3389930261899  1501.555748491856  
CD1_1 CD1_2  : 6.547003719386885e-07  -0.00011159088984525699  
CD2_1 CD2_2  : -0.00012238492146305314  -1.1457490550730293e-05  
NAXIS : 1645  1678

参考图像WCS参数 #

手头参考图像与当前图像为180度旋转关系，其WCS参数如下，拟合得到的CRVAL数值符合预期：

CTYPE1  = 'RA---TAN'           / Pixel coordinate system                        
CTYPE2  = 'DEC--TAN'           / Pixel coordinate system                        
CRPIX1  =    -6457.18566867618 / Ref. pixel of center of rotation               
CRPIX2  =    -4673.96071610528 / Ref. pixel of center of rotation               
CRVAL1  =     125.116590442276 / [deg] [deg] 08:24:31.5 Value of ref pixel      
CRVAL2  =    -43.3619270216124 / [deg] [deg] -42:34:50.2 Value of ref pixel     
CD1_1   =  5.92477204932042E-5 / WCS transform matrix element                   
CD2_1   =  3.73014240770773E-8 / WCS transform matrix element                   
CD1_2   =  -3.4695999305039E-8 / WCS transform matrix element                   
CD2_2   =  5.92440992599392E-5 / WCS transform matrix element    

初始错误写入代码 #

hdulist = fits.open('Swirl06p_1.fits')
header = hdulist[0].header

header.keys
header.set('CTYPE1', 'RA---TAN')
header.set('CTYPE2', 'DEC--TAN')
header.set('CRVAL1', 125.67776105611648)
header.set('CRVAL2', -42.99124198597975)
header.set('CRPIX1', 1451.3389930261899)
header.set('CRPIX2', 1501.555748491856)
header.set('CD1_1', 6.547003719386885e-07)
header.set('CD1_2', -0.00011159088984525699)
header.set('CD2_1', -0.00012238492146305314)
header.set('CD2_2', -1.1457490550730293e-05)
header.set('NAXIS1', 1645, after=3)
header.set('NAXIS2', 1678, after=4)

hdulist.writeto('Swirl06p_1WCS.fits')
hdulist.close()

问题表现 #

写入后生成的图像WCS结果异常：赤经约为180、赤纬约为-3，与预期的126左右、-43左右完全不符。

运行环境 #

• Anaconda 4.10.3
• astropy 4.3.1
• numpy 1.20.3
• 操作系统：Ubuntu 20.04.3 LTS

最终解决方案 #

错误原因是手动逐个写入WCS关键字时，没有清理FITS头中原有的旧WCS相关字段，导致新老字段冲突。正确做法是直接将拟合得到的WCS对象转换为标准头字段，批量更新到FITS头中，避免字段冲突或遗漏。

可运行的最终代码 #

hdulist = fits.open('Swirl06p_1.fits')
header = hdulist[0].header
stars = np.array([[1246.63, 1372.83, 1455.12, 1611.95, 1644.85],
               [1588.42, 1502.92, 1677.24, 1610.39, 1325.12]])

known_coords = SkyCoord([(125.66419083, -42.96809252), (125.67730695, -42.98209958),
                         (125.65082259, -42.9914015), (125.6611325, -43.01438513), (125.70471982, -43.01228167)],
                        frame="icrs", unit="deg")

w = fit_wcs_from_points(xy = stars, world_coords = known_coords, projection='TAN')
hdulist[0].header.update(w.to_header())
hdulist.writeto('Swirl06p_1WCS.fits')
hdulist.close()

内容的提问来源于stack exchange，提问作者Jazamatazz