如何将FITS格式太阳图像质心移至指定坐标以完成对齐？

不用转PNG格式，直接基于numpy数组处理FITS原始数据即可完成对齐，以下是针对你需求的具体实现步骤：

核心思路 #

计算当前图像质心与模板圆心的偏移量，对FITS数据数组执行平移裁剪操作，全程保留FITS格式，无需转换为图片格式。

准备工作 #

确保已安装numpy和astropy（天文领域处理FITS文件的标准库），未安装的话执行以下命令：

pip install numpy astropy

具体代码实现 #

1. 定义基础参数 #

替换为你自己的模板参数和文件路径：

import numpy as np
from astropy.io import fits

# 你的模板参数：圆心坐标(x,y)、太阳盘面半径
target_x = 500  # 替换为模板中心x坐标
target_y = 500  # 替换为模板中心y坐标
template_r = 200  # 替换为模板半径

# 读取单张FITS示例文件
hdul = fits.open('your_sun_image.fits')
data = hdul[0].data
original_header = hdul[0].header
hdul.close()

2. 计算当前图像质心（复用你的代码） #

# 提取太阳盘面像素并计算质心
y_indices, x_indices = np.where(data >= 1e4)
centroid_y = np.mean(y_indices)
centroid_x = np.mean(x_indices)

3. 执行平移裁剪对齐 #

利用偏移量调整图像索引范围，裁剪出以模板圆心为中心的太阳盘面区域：

# 计算质心到目标圆心的偏移量
dx = target_x - centroid_x
dy = target_y - centroid_y

# 确定目标裁剪区域的边界（以模板圆心为中心，覆盖整个太阳盘面）
crop_start_y = int(target_y - template_r)
crop_end_y = int(target_y + template_r)
crop_start_x = int(target_x - template_r)
crop_end_x = int(target_x + template_r)

# 根据偏移量计算当前图像中对应的裁剪范围
shifted_start_y = int(crop_start_y - dy)
shifted_end_y = int(crop_end_y - dy)
shifted_start_x = int(crop_start_x - dx)
shifted_end_x = int(crop_end_x - dx)

# 确保索引不超出图像边界，避免报错
shifted_start_y = max(0, shifted_start_y)
shifted_end_y = min(data.shape[0], shifted_end_y)
shifted_start_x = max(0, shifted_start_x)
shifted_end_x = min(data.shape[1], shifted_end_x)

# 裁剪得到对齐后的图像数据
aligned_data = data[shifted_start_y:shifted_end_y, shifted_start_x:shifted_end_x]

4. 保存对齐后的FITS文件 #

# 复用原文件的头信息，保存对齐后的图像
new_hdul = fits.PrimaryHDU(aligned_data, header=original_header)
new_hdul.writeto('aligned_sun_image.fits', overwrite=True)

5. 批量处理200+图像的循环示例 #

import os

# 输入输出文件夹路径
input_dir = 'your_fits_input_folder'
output_dir = 'aligned_fits_output'

# 创建输出文件夹
os.makedirs(output_dir, exist_ok=True)

# 遍历所有FITS文件
for filename in os.listdir(input_dir):
    if not filename.endswith('.fits'):
        continue
    
    file_path = os.path.join(input_dir, filename)
    # 读取文件
    hdul = fits.open(file_path)
    data = hdul[0].data
    header = hdul[0].header
    hdul.close()
    
    # 计算质心
    y_indices, x_indices = np.where(data >= 1e4)
    centroid_y = np.mean(y_indices)
    centroid_x = np.mean(x_indices)
    
    # 计算偏移并裁剪
    dx = target_x - centroid_x
    dy = target_y - centroid_y
    
    shifted_start_y = int(target_y - template_r - dy)
    shifted_end_y = int(target_y + template_r - dy)
    shifted_start_x = int(target_x - template_r - dx)
    shifted_end_x = int(target_x + template_r - dx)
    
    # 边界检查
    shifted_start_y = max(0, shifted_start_y)
    shifted_end_y = min(data.shape[0], shifted_end_y)
    shifted_start_x = max(0, shifted_start_x)
    shifted_end_x = min(data.shape[1], shifted_end_x)
    
    aligned_data = data[shifted_start_y:shifted_end_y, shifted_start_x:shifted_end_x]
    
    # 保存结果
    output_path = os.path.join(output_dir, f'aligned_{filename}')
    new_hdul = fits.PrimaryHDU(aligned_data, header=header)
    new_hdul.writeto(output_path, overwrite=True)

注意事项 #

• 若不想裁剪，可使用numpy.pad函数用背景值（比如1e3）填充平移后的空白区域，但裁剪方案更高效且符合你的需求。
• 调整template_r时确保能完整覆盖太阳盘面，避免裁剪丢失有效数据。

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