如何用Astropy将JWST FITS图像重投影至指定笛卡尔坐标系？

核心方案 #

使用天文领域标准库astropy和reproject实现图像重投影，这两个工具比你之前尝试的FITS_tools更稳定且功能完整，核心步骤分为读取图像与WCS、定义目标坐标系、执行重投影三个环节。

1. 安装依赖 #

先确保安装必要的库：

pip install astropy reproject

2. 读取FITS文件与WCS信息 #

从FITS文件中提取图像数据和对应的世界坐标系（WCS）信息：

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

# 读取MIRI图像及WCS
miri_hdu = fits.open("jw02107-o039_t018_miri_f1130w_i2d.fits")[0]
miri_data = miri_hdu.data
miri_wcs = WCS(miri_hdu.header)

# 读取NIRcam图像及WCS
nircam_hdu = fits.open("jw02107-o040_t018_nircam_clear-f335m_i2d.fits")[0]
nircam_data = nircam_hdu.data
nircam_wcs = WCS(nircam_hdu.header)

3. 定义目标坐标系与尺寸 #

要统一到x=5000、y=2000的笛卡尔坐标系，需基于两张图像的共同天区构建目标WCS（采用天文常用的CAR笛卡尔投影，ICRS参考系）：

# 计算两张图像的共同天区中心坐标
# 获取MIRI图像四角的天球坐标
miri_corners = miri_wcs.pixel_to_world(
    [0, 0, miri_data.shape[1], miri_data.shape[1]],
    [0, miri_data.shape[0], 0, miri_data.shape[0]]
)
# 获取NIRcam图像四角的天球坐标
nircam_corners = nircam_wcs.pixel_to_world(
    [0, 0, nircam_data.shape[1], nircam_data.shape[1]],
    [0, nircam_data.shape[0], 0, nircam_data.shape[0]]
)

# 计算共同天区的RA/Dec中心
all_ra = [c.ra.degree for c in miri_corners + nircam_corners]
all_dec = [c.dec.degree for c in miri_corners + nircam_corners]
center_ra = (min(all_ra) + max(all_ra)) / 2
center_dec = (min(all_dec) + max(all_dec)) / 2

# 构建目标WCS
target_wcs = WCS(naxis=2)
target_wcs.wcs.crpix = [2500, 1000]  # 目标图像中心像素（x=2500, y=1000）
target_wcs.wcs.crval = [center_ra, center_dec]  # 对应天球中心坐标
target_wcs.wcs.cdelt = [-miri_wcs.wcs.cdelt[0], miri_wcs.wcs.cdelt[1]]  # 保持原像素尺度（负号匹配RA方向）
target_wcs.wcs.ctype = ["RA---CAR", "DEC--CAR"]  # 笛卡尔投影
target_wcs.wcs.radesys = "ICRS"  # JWST默认参考系

# 目标图像尺寸：numpy数组维度为(y, x)，对应(2000, 5000)
target_shape = (2000, 5000)

4. 执行重投影 #

使用reproject_exact（高精度）或reproject_interp（快速插值）完成图像转换：

from reproject import reproject_exact

# 重投影MIRI图像
miri_reproj, miri_footprint = reproject_exact(
    (miri_data, miri_wcs), target_wcs, shape_out=target_shape
)

# 重投影NIRcam图像
nircam_reproj, nircam_footprint = reproject_exact(
    (nircam_data, nircam_wcs), target_wcs, shape_out=target_shape
)

5. 验证结果 #

用matplotlib可视化重投影后的图像，确认尺寸和对齐情况：

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

ax1.imshow(miri_reproj, origin='lower', cmap='inferno')
ax1.set_title('MIRI Reprojected (5000x2000)')
ax1.set_xlabel('X Pixel')
ax1.set_ylabel('Y Pixel')

ax2.imshow(nircam_reproj, origin='lower', cmap='inferno')
ax2.set_title('NIRcam Reprojected (5000x2000)')
ax2.set_xlabel('X Pixel')
ax2.set_ylabel('Y Pixel')

plt.tight_layout()
plt.show()

问题排查 #

• WCS模块无合适方法：你之前使用的WCS工具可能不是astropy.wcs，这是当前天文领域的标准WCS实现，功能覆盖所有重投影需求。
• FITS_tools变量未初始化报错：该错误通常是函数参数传递错误或文件读取失败导致，替换为astropy+reproject组合可避免这类问题。

注意事项 #

• 若目标是纯粹的笛卡尔坐标系（无需天球映射），可直接定义线性WCS，跳过天区中心计算环节。
• 大图像建议用reproject_interp替代reproject_exact以提升速度。
• 确保两张图像的天球参考系一致（JWST图像默认均为ICRS，无需额外转换）。

内容的提问来源于stack exchange，提问作者Kutluhan Aygüzel