如何从FITS文件读取光谱数据？SpecUtils使用报错求助

问题背景 #

作为SpecUtils新用户，参照官方文档代码读取FITS光谱文件时，出现IndexError，无法通过loglam和flux索引获取数据，报错提示仅整数、切片等可作为有效索引。

1. 查看FITS文件的表头与数据结构 #

报错根源在于对FITS文件的HDU（Header Data Unit）结构和数据类型判断错误。先通过以下代码查看文件的完整结构：

from astropy.io import fits

# 打开FITS文件（用with语句自动关闭）
with fits.open('GIRAF.2021-02-14T01:00:34.723.fits') as fits_file:
    # 打印所有HDU的基本信息
    fits_file.info()
    # 查看目标HDU的数据类型
    target_hdu = fits_file[0]
    print(f"数据类型: {type(target_hdu.data)}")
    # 如果是表格型数据，打印所有列名
    if hasattr(target_hdu.data, 'columns'):
        print(f"表格列名: {[col.name for col in target_hdu.data.columns]}")
    # 打印表头元数据（包含波长校准等关键信息）
    print("\n表头信息:")
    print(repr(target_hdu.header))

通过输出可以确认：

• 光谱数据实际所在的HDU索引（原代码注释提到在第二个HDU，但用了f[0]，可能索引错误）
• 数据是普通numpy数组还是带列名的表格结构
• 表头中是否包含CRVAL1、CDELT1等波长校准关键字

2. 正确读取光谱数据 #

根据上述查看结果，调整读取逻辑：

from astropy.io import fits
from astropy import units as u
import numpy as np
from specutils import Spectrum1D

with fits.open('GIRAF.2021-02-14T01:00:34.723.fits') as fits_file:
    # 优先检查是否存在第二个HDU（索引1），若存在则使用它
    hdu = fits_file[1] if len(fits_file) > 1 else fits_file[0]
    
    # 情况1：数据是带列名的表格结构（如原文档示例）
    if hasattr(hdu.data, 'columns'):
        if 'loglam' in [col.name for col in hdu.data.columns] and 'flux' in [col.name for col in hdu.data.columns]:
            spectral_axis = 10**hdu.data['loglam'] * u.AA
            flux = hdu.data['flux'] * 10**-17 * u.Unit('erg cm-2 s-1 AA-1')
        else:
            raise ValueError("表格中未找到loglam或flux列，请检查列名")
    
    # 情况2：数据是普通数组，通过表头关键字计算波长
    else:
        # 从表头获取波长校准参数
        crval1 = hdu.header.get('CRVAL1')
        cdelt1 = hdu.header.get('CDELT1')
        naxis1 = hdu.header.get('NAXIS1')
        
        if all([crval1, cdelt1, naxis1]):
            # 生成光谱轴（波长）
            spectral_axis = (crval1 + cdelt1 * np.arange(naxis1)) * u.AA
            # 通量数据单位根据实际情况调整
            flux = hdu.data * u.Unit('erg cm-2 s-1 AA-1')
        else:
            raise ValueError("表头缺少CRVAL1、CDELT1或NAXIS1关键字，无法计算波长")

# 创建SpecUtils的Spectrum1D对象
spec = Spectrum1D(spectral_axis=spectral_axis, flux=flux)

3. 绘制光谱图 #

使用matplotlib结合SpecUtils完成绘图：

from matplotlib import pyplot as plt
from astropy.visualization import quantity_support

# 启用单位支持，自动在坐标轴显示单位
quantity_support()

# 基础绘图
plt.figure(figsize=(10, 6))
plt.plot(spec.spectral_axis, spec.flux, linewidth=0.5)
plt.xlabel(f"波长 ({spec.spectral_axis.unit})")
plt.ylabel(f"通量 ({spec.flux.unit})")
plt.title("GIRAF观测光谱")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# （可选）用SpecUtils内置方法绘图，支持更多光谱操作
from specutils.manipulation import FluxConservingResampler

# 重采样到均匀波长轴，优化绘图效果
new_wavelengths = np.linspace(spec.spectral_axis.value.min(), 
                              spec.spectral_axis.value.max(), 
                              1500) * spec.spectral_axis.unit
resampler = FluxConservingResampler()
resampled_spec = resampler(spec, new_wavelengths)

plt.figure(figsize=(10, 6))
resampled_spec.plot(linewidth=0.8)
plt.title("重采样后的GIRAF光谱")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

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