如何正确对.fits文件执行rebinning以匹配目标像素尺度？

我正在处理两张FITS图像，二者像素尺度（pixscale）不同：图像1的CDELT为±0.0013888888888889度，图像2的CDELT为±0.00277778000000度。需要将图像1的像素尺度匹配到图像2（即重采样（rebinning）），但尝试了网络及AI提供的代码后，结果数据混乱、绘图无意义，且输出文件远大于原文件（降低分辨率不应出现此情况）。我最后尝试的代码如下：

import numpy as np
from scipy.ndimage import zoom
from astropy.io import fits

# Load the data from the second image
hdul = fits.open('ha_image.fits')
data = hdul[0].data
header = hdul[0].header

# Calculate the rebinning factors
scale_factor_x = 0.00277778000000 / 0.0013888888888889
scale_factor_y = 0.00277778000000 / 0.0013888888888889

# Rebin the data using zoom
rebinned_data = zoom(data, (scale_factor_y, scale_factor_x), order=1)

# Update the WCS information
header['CDELT1'] = -0.00277778000000
header['CDELT2'] = 0.00277778000000

# Create a new PrimaryHDU object with the rebinned data and updated header
rebinned_hdu = fits.PrimaryHDU(rebinned_data, header=header)

# Save the rebinned image to a new FITS file
rebinned_hdu.writeto('rebinned_ha_image.fits', overwrite=True)

问题分析 #

你的代码存在两个核心错误：

1. 缩放方向完全搞反：要将高分辨率（小像素）的图像1缩放到低分辨率（大像素）的图像2尺度，缩放因子应为原像素尺度/目标像素尺度（即≈0.5），但你用了目标/原得到≈2的放大因子，导致图像被放大4倍，文件体积自然增大。
2. WCS头信息不完整：仅修改CDELT无法让头文件与新数据匹配，还必须同步更新图像尺寸（NAXIS1/NAXIS2），否则WCS坐标与实际像素错位，导致绘图混乱。

正确解决方案 #

根据你的像素尺度（刚好是约2倍整数关系），推荐两种可靠方法：

方法1：整数倍重采样（精度损失最小）

利用整数倍合并像素的方式，直接对原像素取平均或求和，保留数据统计特性：

import numpy as np
from astropy.io import fits

# 加载需要重采样的图像1（高分辨率）
hdul = fits.open('image1.fits')
data = hdul[0].data
header = hdul[0].header.copy()  # 复制头文件避免修改原数据

# 计算重采样倍数：原像素尺度/目标像素尺度，这里≈2
bin_factor = int(round(abs(header['CDELT1']) / 0.00277778000000))
# 保持CDELT的正负方向与原图像一致
target_cdelt1 = -0.00277778000000 if header['CDELT1'] < 0 else 0.00277778000000
target_cdelt2 = 0.00277778000000 if header['CDELT2'] > 0 else -0.00277778000000

# 裁剪图像到能被重采样倍数整除的尺寸
new_shape = (data.shape[0] // bin_factor, data.shape[1] // bin_factor)
cropped_data = data[:new_shape[0]*bin_factor, :new_shape[1]*bin_factor]

# 对N×N原像素取平均（若需保留总通量可改为sum）
rebinned_data = cropped_data.reshape(new_shape[0], bin_factor, new_shape[1], bin_factor).mean(axis=(1, 3))

# 更新头文件关键信息
header['NAXIS1'] = new_shape[1]
header['NAXIS2'] = new_shape[0]
header['CDELT1'] = target_cdelt1
header['CDELT2'] = target_cdelt2
# 同步调整参考像素位置，保持WCS坐标正确
if 'CRPIX1' in header:
    header['CRPIX1'] = (header['CRPIX1'] - 0.5) / bin_factor + 0.5
if 'CRPIX2' in header:
    header['CRPIX2'] = (header['CRPIX2'] - 0.5) / bin_factor + 0.5

# 保存结果
rebinned_hdu = fits.PrimaryHDU(rebinned_data, header=header)
rebinned_hdu.writeto('rebinned_image1.fits', overwrite=True)
hdul.close()

方法2：基于WCS的通用重采样（支持非整数倍）

如果像素尺度不是整数倍，使用astropy的reproject工具，直接匹配目标图像的WCS，确保坐标完全对齐：

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

# 加载源图像（图像1）和目标图像（图像2，用于获取目标WCS）
source_hdul = fits.open('image1.fits')
target_hdul = fits.open('image2.fits')

# 获取源和目标的WCS信息与图像尺寸
source_wcs = WCS(source_hdul[0].header)
target_wcs = WCS(target_hdul[0].header)
target_shape = (target_hdul[0].header['NAXIS2'], target_hdul[0].header['NAXIS1'])

# 执行重采样（线性插值，可选nearest/spline等方法）
rebinned_data, _ = reproject_interp(
    (source_hdul[0].data, source_wcs),
    target_wcs,
    shape_out=target_shape
)

# 使用目标图像的头文件保存结果
rebinned_header = target_hdul[0].header.copy()
rebinned_hdu = fits.PrimaryHDU(rebinned_data, header=rebinned_header)
rebinned_hdu.writeto('rebinned_image1_to_image2.fits', overwrite=True)

# 关闭文件
source_hdul.close()
target_hdul.close()

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