如何从fits格式天文星图中移除饱和像素并输出处理后图像？

核心思路 #

将二维数组中等于最大像素值的饱和点替换为np.nan（空值），matplotlib渲染时会自动将空值区域显示为透明孔洞，完全匹配需求效果。由于FITS原始数据多为无符号16位整型（uint16），不支持空值存储，需要先将数组转换为浮点型再做替换。

完整可运行代码 #

import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt

fits_image_filename = "Acturus_V_2s.fit"
hdul = fits.open(fits_image_filename)

data = hdul[0].data
hdul.close() # 读取完成后关闭FITS文件句柄

datacut = data[610:710,755:855].astype(float) # 转换为浮点型支持NaN赋值
MaxPixelValue = np.amax(datacut)
print(MaxPixelValue) 

# 替换饱和像素为NaN
datacut[datacut == MaxPixelValue] = np.nan

# 展示处理后的图像
plt.figure(figsize=(6,6))
plt.imshow(datacut, origin="lower")
plt.colorbar(label="像素值")
plt.title("移除饱和像素后的星点子图")
plt.show()

关键逻辑说明 #

• 布尔索引datacut == MaxPixelValue会自动筛选出所有值为65535的饱和像素位置，直接对该索引区域赋值即可完成批量替换，不需要循环遍历
• 替换为np.nan的区域在渲染时不会显示颜色，会直接透出画布底色，形成清晰的中心孔洞
• 如果你不需要透明孔洞，也可以将替换值改为子图的背景值（比如非饱和区域的最小值np.amin(datacut[datacut != MaxPixelValue])），效果同样符合要求

可选操作：保存处理后的FITS文件 #

如果需要导出处理后的子图为新的FITS文件，可以追加以下代码：

# 构造新的FITS HDU
new_hdu = fits.PrimaryHDU(datacut)
new_hdul = fits.HDUList([new_hdu])
new_hdul.writeto("Acturus_V_2s_desaturated.fits", overwrite=True)
new_hdul.close()

内容的提问来源于stack exchange，提问作者Johanna Björk