如何将含天区坐标的多维数组写入FITS文件并设置正确轴格式？

嘿，这个需求我熟！用Python的astropy.io.fits库就能完美解决，它是天文数据处理的标配工具。我给你一步步拆解实现，保证轴的配置完全符合你的要求：

核心思路 #

你的array2是[100,100,50]的三维数组，前两轴对应天区坐标，第三轴对应array1的50个值。我们需要：

• 将array2作为FITS的图像数据写入
• 在FITS表头中配置前两轴的天区坐标元数据（用标准WCS格式）
• 把array1的50个值关联到第三轴，让后续读取时能明确对应关系

完整代码实现 #

1. 导入依赖库 #

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

2. 准备数据（模拟你的输入，实际使用时替换成你的真实数据） #

# 模拟array1：50个自定义值（比如能量、波段等）
array1 = np.linspace(0.5, 10.0, 50)
# 模拟array2：100x100x50的天区数据
array2 = np.random.rand(100, 100, 50)

3. 配置WCS与FITS表头 #

WCS（World Coordinate System）是天文数据中描述像素与真实天区坐标映射的标准，我们用它来定义前两轴的天区信息，同时关联第三轴与array1：

# 创建3轴WCS对象
wcs = WCS(naxis=3)

# 配置前两轴（天区坐标：RA/赤经、DEC/赤纬）
wcs.wcs.ctype = ["RA---TAN", "DEC--TAN", "CHANNEL"]  # 第三轴自定义为"通道"
wcs.wcs.crval = [120.0, 30.0, array1[0]]  # 参考点坐标：RA=120°，DEC=30°，第三轴起始值
wcs.wcs.crpix = [50.5, 50.5, 1.0]  # 参考点的像素位置（FITS像素从1开始计数，这里取中心）
wcs.wcs.cdelt = [-0.1, 0.1, array1[1]-array1[0]]  # 像素增量：RA每个像素减0.1°，DEC加0.1°，第三轴步长
wcs.wcs.cunit = ["deg", "deg", "keV"]  # 单位：角度用度，第三轴假设为keV（可根据你的array1含义修改）

# 将WCS信息转为FITS表头格式
wcs_header = wcs.to_header()

4. 创建FITS文件并写入数据 #

# 创建ImageHDU（图像数据单元），传入array2和配置好的表头
hdu = fits.ImageHDU(data=array2, header=wcs_header)

# 额外把array1的所有值存入表头，方便直接读取
hdu.header['CHAN_VALS'] = array1.tolist()

# 创建HDU列表（必须包含PrimaryHDU），写入文件
hdulist = fits.HDUList([fits.PrimaryHDU(), hdu])
hdulist.writeto('sky_data.fits', overwrite=True)  # overwrite=True覆盖已有文件
hdulist.close()

关键说明 #

• 前两轴的WCS配置：CTYPE1和CTYPE2用RA---TAN和DEC--TAN表示切向投影的赤经赤纬，你可以根据实际天区投影方式修改（比如GAL-TAN代表银道坐标）。
• 第三轴关联：通过CTYPE3定义轴类型，同时把array1的完整值存入CHAN_VALS关键字，这样读取时可以直接获取每个通道对应的原始值。
• 读取验证：后续读取时，你可以用fits.getdata('sky_data.fits')获取array2，用fits.getheader('sky_data.fits')['CHAN_VALS']获取array1，也可以用WCS(fits.getheader('sky_data.fits'))解析轴的坐标信息。

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