如何基于坐标对3D FITS数据立方的NumPy数组进行掩膜

实现步骤 #

1. 读取FITS数据并提取WCS坐标信息
2. 生成像素网格对应的银道坐标
3. 将银道坐标转换为赤道坐标以获取赤纬δ
4. 构建满足所有条件的二维掩膜
5. 扩展掩膜维度并应用到3D数据立方

代码实现 #

import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
import astropy.units as u

# 读取FITS数据立方
hdul = fits.open("目标FITS文件路径.fits")
data_cube = hdul[0].data  # 形状为[v, y, x]
wcs_info = WCS(hdul[0].header)

# 生成像素坐标网格
y_dim, x_dim = data_cube.shape[1], data_cube.shape[2]
y_pix, x_pix = np.mgrid[0:y_dim, 0:x_dim]

# 转换像素坐标为银道坐标(l, b)
gal_coords = wcs_info.pixel_to_world(x_pix, y_pix)
l_deg = gal_coords.l.degree
b_deg = gal_coords.b.degree

# 转换银道坐标为赤道坐标，提取赤纬δ
eq_coords = gal_coords.transform_to('icrs')
delta_deg = eq_coords.dec.degree

# 构建二维掩膜
mask_2d = (l_deg > 229) & (l_deg < 234) & (b_deg > 4) & (b_deg < 14) & (delta_deg > -13)

# 扩展掩膜至3D以匹配数据立方维度
mask_3d = np.expand_dims(mask_2d, axis=0)
# 应用掩膜，用NaN填充非目标区域（也可替换为0或其他值）
masked_cube = np.where(mask_3d, data_cube, np.nan)

# 保存掩膜后的数据（可选）
hdul[0].data = masked_cube
hdul.writeto("掩膜后的数据立方.fits", overwrite=True)
hdul.close()

关键说明 #

• WCS坐标转换: 利用astropy的WCS工具自动处理像素与天球坐标的映射，无需手动计算像素尺度，避免人为误差。
• 跨坐标系转换: 通过transform_to('icrs')直接将银道坐标转为赤道坐标系，轻松获取每个像素的赤纬值。
• 掩膜扩展: 将二维空间掩膜扩展为三维，确保所有速度切片都应用相同的空间筛选条件。
• 数据标记: 使用NaN填充非目标区域，后续分析工具（如astropy、matplotlib）可自动忽略这些无效值。

注意事项 #

• 确认FITS文件的WCS头信息完整，否则坐标转换会失败。
• 若数据立方的轴顺序不是[v,y,x]，需调整代码中data_cube.shape的索引，确保y、x对应正确的空间轴。
• 处理超大数据立方时，可考虑分块处理以降低内存占用。

内容的提问来源于stack exchange，提问作者Firestar-Reimu