如何为Numpy数组生成的FITS图像添加WCS坐标以在DS9中显示

嗨，这个需求我太熟悉了！根本不用手动逐个处理像素坐标，只要给生成的FITS文件加上符合标准的WCS（世界坐标系统）头信息，ds9就会自动识别，鼠标悬停时直接显示你要的银道、赤道坐标。下面给你一步步讲怎么实现：

核心原理 #

ds9读取FITS文件时，会解析头文件里的WCS关键字，根据这些参数自动计算每个像素对应的天球坐标。我们只需要用astropy库（天文数据处理的必备工具）生成正确的WCS头，再和你的Numpy数组一起保存成FITS就行。

具体实现步骤 #

1. 准备依赖库 #

如果你还没装astropy，先通过pip安装：

pip install astropy

2. 创建WCS对象并配置参数 #

假设你的ZEA_N_sky是128x128的ZEA投影（从变量名看应该是天顶等面积投影），我们以银道坐标系为例配置WCS，ds9后续可以自动转换成赤道坐标：

import numpy as np
from astropy.io import fits
from astropy.wcs import WCS

# 你的原始数组，假设n=128
n = 128
ZEA_N_sky = np.zeros((n, n))
# （这里省略你给像素赋值的代码）

# 初始化WCS对象，指定2维图像
wcs = WCS(naxis=2)

# 关键配置：设置投影类型和坐标系
# GLON-ZEA表示银经（Galactic Longitude）用ZEA投影，GLAT-ZEA是银纬
wcs.wcs.ctype = ["GLON-ZEA", "GLAT-ZEA"]

# 设置参考像素：FITS采用1-based索引，128x128图像的中心像素是(64.5, 64.5)
wcs.wcs.crpix = [n/2 + 0.5, n/2 + 0.5]

# 参考像素对应的天球坐标：这里设为银心（银经0°，银纬0°），你可以改成自己需要的参考点
wcs.wcs.crval = [0.0, 0.0]

# 像素大小：单位是度，这里设为0.5度/像素，经度方向负号是因为天体图像通常左到右对应经度减小，可根据你的需求调整
wcs.wcs.cdelt = [-0.5, 0.5]

# 可选：添加历元信息，比如J2000标准历元
wcs.wcs.equinox = 2000.0

3. 保存带WCS的FITS文件 #

把WCS信息写入FITS头，再和数据一起保存：

# 创建PrimaryHDU（主数据单元），传入你的数组
hdu = fits.PrimaryHDU(data=ZEA_N_sky)

# 把WCS的关键字扩展到FITS头里
hdu.header.extend(wcs.to_header())

# 保存文件，overwrite=True允许覆盖已存在的文件
hdul = fits.HDUList([hdu])
hdul.writeto('ZEA_N_sky_with_wcs.fits', overwrite=True)

验证效果 #

用ds9打开生成的ZEA_N_sky_with_wcs.fits：

• 鼠标悬停在图像上，状态栏会自动显示当前像素对应的银经、银纬；
• 如果你想看赤经、赤纬，只需要在ds9顶部菜单选Frame -> Coordinate System -> Equatorial，鼠标悬停就会切换显示赤道坐标了。

额外提示 #

• 如果你的投影不是ZEA，把ctype里的-ZEA改成对应的投影代码即可，比如-CAR是平板投影、-TAN是切向投影，具体可以参考FITS WCS标准；
• 像素大小cdelt如果是弧秒单位，记得转换成度（1弧秒 = 1/3600度）；
• 如果你需要更复杂的WCS配置（比如包含自行、视差等），可以用astropy.wcs.WCS的其他属性进一步设置。

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