如何对FITS文件进行圆形裁剪并保留WCS坐标信息

你当前操作仅替换了HDU中的数据数组为裁剪后的小尺寸数组，未同步更新对应原始全图的WCS头信息，是WCS失效的核心原因，可参考以下两种解决思路：

思路1（推荐）：结合Cutout2D和圆形掩膜实现 #

先通过astropy的Cutout2D生成刚好容纳圆形区域的矩形裁剪片，Cutout2D会自动计算并生成对应裁剪后区域的正确WCS信息，再在该裁剪片上应用圆形掩膜，即可同时满足圆形区域提取和WCS保留的需求。
示例代码如下：

from astropy.io import fits
from astropy.wcs import WCS
from astropy.nddata import Cutout2D
from regions import PixCoord, CirclePixelRegion
import matplotlib.pyplot as plt

# 用with上下文管理器自动处理文件关闭，避免读写异常
with fits.open('/Users/Wesley/OneDrive/Documents/Honours/Data and Code/M101-10/2-Luminosity/M101-LMIPS24.fits') as hdulist:
    hdu = hdulist[0]
    orig_wcs = WCS(hdu.header)
    orig_data = hdu.data

    # 定义圆形区域参数
    center = PixCoord(144, 144)
    radius = 85
    # 先做矩形cutout，尺寸刚好容纳圆形，自动生成对应WCS
    cutout_size = 2 * radius # 也可根据需求留少量余量，比如2*radius+2
    cutout = Cutout2D(orig_data, center, cutout_size, wcs=orig_wcs)
    
    # 在裁剪后的小图上生成圆形掩膜
    new_center = PixCoord(radius, radius) # 裁剪后圆心移到新图中心位置
    aperture = CirclePixelRegion(new_center, radius)
    mask = aperture.to_mask(mode='exact')
    weighted_data = mask.multiply(cutout.data)

    # 生成保存用HDU，写入裁剪后的正确WCS头和掩膜后的数据
    new_hdu = fits.PrimaryHDU(weighted_data, header=cutout.wcs.to_header())
    new_hdu.writeto('/Users/Wesley/OneDrive/Documents/Honours/Data and Code/M101-10/2-Luminosity/M101-Test.fits', overwrite=True)

# 可视化验证
plt.title("Data cutout multiplied by mask", size=9)
plt.imshow(weighted_data, cmap=plt.cm.viridis, origin='lower')
plt.show()

思路2：手动更新原始WCS头 #

如果你不想使用Cutout2D，也可以手动调整原始WCS的参考像素坐标：

• 先确定裁剪区域在原始大图上的左下角像素坐标(x_min, y_min)，本次裁剪的x范围是59_{229、y范围是59}229，对应x_min=59, y_min=59（注意astropy像素默认从0开始计数，如果你的WCS头符合FITS标准从1开始计数，偏移量要额外加1）
• 把原始WCS头里的CRPIX1减去x_min，CRPIX2减去y_min，其余WCS参数（CRVAL、CD矩阵等）保持不变
• 将更新后的WCS头写入新的HDU，替换原始头即可

额外注意事项 #

原有代码中在hdulist.close()之后仍修改hdulist[0].data并写入文件的操作存在隐患，close后访问的仅为内存缓存数据，建议使用with上下文管理器管理文件打开关闭流程，避免读写异常。

效果参考：

内容的提问来源于stack exchange，提问作者Wesley Van Kempen