WFSS 光谱第三部分：发射线图谱

WFSS 光谱第三部分：发射线图谱#

用例： 最优提取光谱；红移测量；发射线图。简化版的 JDox 科学用例 # 33。

数据： JWST 模拟的 NIRISS 图像来自 MIRAGE，经过 JWST 校准管道处理；星系团。

工具： specutils, astropy, pandas, emcee, lmfit, corner, h5py。

跨仪器： NIRSpec

文档： 本笔记本是 STScI 更大后处理数据分析工具生态系统的一部分。

引言#

本笔记本是关于 NIRISS WFSS 数据的一组 4 个笔记本中的第 4 个：

1D 最优提取，因为 JWST 管道仅提供框提取。最优提取提高了微弱源光谱的信噪比（S/N）。
合并和归一化 1D 光谱。
将星系与模板进行交叉关联以获取红移。
空间分辨的发射线图。

本笔记本从 NIRISS WFSS 光谱的 2D 剪切图中获取发射线图。本笔记本将需要 NIRISS WFSS 光谱的 2D 剪切图（例如，l3_nis_f150w_G150C_s00004_cal.fits；管道 spec3 产品）。

注意： 在这个例子中，笔记本仅使用在一个抖动位置直接从当前版本的管道（构建 7.5）获得的 2D 矫正光谱，但也可以从所有抖动位置堆叠的光谱开始，以提高信噪比。

%matplotlib inline  # 在Jupyter Notebook中内联显示Matplotlib图形

import os  # 导入操作系统模块，用于文件和目录操作

import numpy as np  # 导入NumPy库，用于数值计算和数组操作

from astropy.io import fits  # 从Astropy库导入FITS模块，用于处理FITS文件

import astropy.units as u  # 导入Astropy单位模块，用于物理量的单位处理

from astropy.modeling import models  # 从Astropy建模模块导入模型类

from astropy.modeling.polynomial import Chebyshev1D  # 导入一维切比雪夫多项式模型

from specutils.spectra import Spectrum1D  # 从specutils库导入一维光谱类

from specutils.fitting import fit_lines  # 从specutils库导入拟合谱线的函数

from specutils.fitting import continuum  # 从specutils库导入连续谱拟合的函数

from astropy import __version__ as astropy_version  # 获取Astropy库的版本信息

print('astropy', astropy_version)  # 打印Astropy库的版本信息

import matplotlib.pyplot as plt  # 导入绘图库 matplotlib 的 pyplot 模块

import matplotlib as mpl  # 导入 matplotlib 库

# 设置保存图像时的分辨率为 80 dpi
mpl.rcParams['savefig.dpi'] = 80  

# 设置图像显示时的分辨率为 80 dpi
mpl.rcParams['figure.dpi'] = 80

0. 下载笔记本 01 产品#

这些也可以通过运行笔记本来获取。

import os  # 导入os模块，用于文件和目录操作

# 检查输出目录是否存在
if not os.path.exists('./output'):
    
    import zipfile  # 导入zipfile模块，用于处理zip文件
    import urllib.request  # 导入urllib.request模块，用于下载文件

    # 设置要下载的zip文件的链接
    boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/NIRISS_lensing_cluster/output.zip'
    
    # 设置下载后保存的文件名
    boxfile = './output.zip'
    
    # 下载zip文件
    urllib.request.urlretrieve(boxlink, boxfile)
    
    # 打开下载的zip文件
    zf = zipfile.ZipFile(boxfile, 'r')
    
    # 获取zip文件内所有文件的名称列表
    list_names = zf.namelist()
    
    # 遍历文件名称列表
    for item in list_names:
        # 提取zip文件中的每个文件，使用extract比extractall更安全
        zf.extract(member=item)  # 提取文件

else:
    # 如果输出目录已存在，打印提示信息
    print('Already exists')

# 下载文件，如果尚不存在

if not os.path.exists('./pipeline_products'):  # 检查目录是否存在

    import zipfile  # 导入zipfile模块，用于处理zip文件
    import urllib.request  # 导入urllib.request模块，用于下载文件

    # 定义要下载的zip文件的链接
    boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/NIRISS_lensing_cluster/pipeline_products.zip'

    # 定义下载后保存的文件名
    boxfile = './pipeline_products.zip'

    # 下载zip文件
    urllib.request.urlretrieve(boxlink, boxfile)

    # 打开下载的zip文件
    zf = zipfile.ZipFile(boxfile, 'r')

    # 获取zip文件内的文件名列表
    list_names = zf.namelist()

    # 遍历文件名列表，逐个提取文件
    for item in list_names:
        zf.extract(member=item)  # 使用extract提取文件，避免使用extractall以防止文件路径问题

# 选择目标

DIR_OUT = './output/'  # 输出目录

filt = 'f200w'  # 选择的滤光片

grism = 'G150C'  # 选择的光栅

# grism = 'G150R'  # 可选的另一种光栅

id = '00004'  # 目标ID

1. 获取Hα发射线图；#

# 打开二维文件;

# 可以在 No.01a 下载。

DIR_DATA = './pipeline_products/'  # 数据目录

# 零索引的dither编号 --- 测试数据有两个dither位置。

ndither = 0

# 文件名

file_2d = f'{DIR_DATA}l3_nis_{filt}_{grism}_s{id}_cal.fits'  # 构建文件路径

hdu_2d = fits.open(file_2d)  # 打开FITS文件

# 对齐光栅方向

#   - x方向 = 色散（波长）方向。

#   - y方向 = 交叉色散。

# 在这个笔记本中。
    
if grism == 'G150C':  # 如果光栅是G150C

    # 如果光谱是水平的;

    data_2d = hdu_2d[ndither*7+1].data  # 获取二维数据

    dq_2d = hdu_2d[ndither*7+2].data  # 获取数据质量数组

    err_2d = hdu_2d[ndither*7+3].data  # 获取误差数组

    wave_2d = hdu_2d[ndither*7+4].data  # 获取波长数组

else:  # 如果光栅不是G150C

    data_2d = rotate(hdu_2d[ndither*7+1].data, 90)  # 将数据旋转90度

    dq_2d = rotate(hdu_2d[ndither*7+2].data, 90)  # 将数据质量数组旋转90度

    err_2d = rotate(hdu_2d[ndither*7+3].data, 90)  # 将误差数组旋转90度

    wave_2d = rotate(hdu_2d[ndither*7+4].data, 90)  # 将波长数组旋转90度

plt.imshow(data_2d[:, :])  # 显示二维数据数组的图像

在每个y位置减去连续背景：#

根据源的亮度，有两种选择；

1. 在每个交叉色散列的位置减去估计的连续背景。这适用于亮源。
1. 减去一个单一的连续背景（“主连续背景”），假设其形状在交叉色散方向上与一维提取的形状相同。这对于暗源是必需的。

由于我们清楚地看到源的连续背景，这里我们演示选项1。

# 我们希望在y位置查看光谱

# 作为示例；
yy = 10  # 设置y位置为10

spec_unit = u.dimensionless_unscaled  # 定义光谱单位为无量纲

# 创建掩膜，选择特定波长范围内的数据
mask_line = ((wave_2d[yy, :] > 1.75) & (wave_2d[yy, :] < 1.97)) | ((wave_2d[yy, :] > 2.08) & (wave_2d[yy, :] < 2.23))

# 创建光谱对象，包含波长和对应的flux数据
obs = Spectrum1D(spectral_axis=wave_2d[yy, :][mask_line]*u.um, flux=data_2d[yy, :][mask_line]*spec_unit)

# 拟合连续谱，使用7阶切比雪夫多项式模型
cont = continuum.fit_generic_continuum(obs, model=Chebyshev1D(7))

# 绘制观察到的光谱数据
plt.plot(wave_2d[yy, :], data_2d[yy, :], color='r', label=f'Observed at y={yy}')

# 绘制拟合的连续谱
plt.plot(wave_2d[yy, :][mask_line]*u.um, cont(wave_2d[yy, :][mask_line]*u.um), color='b', label='Fitted continuum')

# 添加图例
plt.legend(loc=0)

# 设置x轴范围
plt.xlim(1.6, 2.3)

# 计算并返回经过掩膜处理的二维波长数据的值，单位为微米
cont(wave_2d[yy, :][mask_line]*u.um).value

# 沿y轴重复此操作;

flux_cont2d = data_2d[:, :] * 0  # 初始化一个与data_2d相同形状的数组，所有元素为0

for yy in range(len(data_2d[:, 0])):  # 遍历data_2d的每一行（每个y坐标）

    # 创建掩膜，选择波长在1.75到1.97和2.08到2.23之间的元素
    mask_line = ((wave_2d[yy, :] > 1.75) & (wave_2d[yy, :] < 1.97)) | ((wave_2d[yy, :] > 2.08) & (wave_2d[yy, :] < 2.23))

    # 创建一个Spectrum1D对象，包含选定波长范围内的光谱数据
    obs = Spectrum1D(spectral_axis=wave_2d[yy, :][mask_line]*u.um, flux=data_2d[yy, :][mask_line]*spec_unit)

    # 使用Chebyshev1D模型拟合连续谱
    cont = continuum.fit_generic_continuum(obs, model=Chebyshev1D(7))

    # 将拟合的连续谱值存储到flux_cont2d的相应行
    flux_cont2d[yy, :] = cont(wave_2d[yy, :]*u.um).value

plt.imshow(flux_cont2d[:, :])  # 显示二维光谱数据的图像

plt.title('Fitted continuum')  # 设置图像标题为“拟合的连续谱”

plt.imshow(data_2d[:, :] - flux_cont2d[:, :])  # 显示二维数据图像，减去连续谱数据

plt.title('Continuum subtracted spectrum')  # 设置图像标题为“减去连续谱的光谱”

从去除连续谱的光谱中提取Hα发射；#

# 切割Ha线图

rsq = data_2d.shape[0]  # 获取数据的行数，假设为正方形

cut_ha = np.zeros((rsq, rsq), 'float32')  # 创建一个全零的数组，用于存储切割后的Ha图

zin = 2.1  # 从Notebook No.3中获取的红移值
lamcen = 0.6564 * (1. + zin)  # 计算中心波长，单位为微米

for yy in range(len(data_2d[:, 0])):  # 遍历每一行（每个y像素）

    # 由于波长数组可能倾斜，因此必须在每个y像素上执行此操作
    index_lamcen = np.argmin(np.abs(lamcen - wave_2d[yy, :]))  # 找到与中心波长最接近的波长索引

    # 从数据中减去连续光通量，并切割出相应的Ha线数据
    cut_ha[yy, :] = (data_2d - flux_cont2d)[yy, int(index_lamcen - rsq / 2.):int(index_lamcen + rsq / 2.)]

plt.imshow(cut_ha)  # 显示切割后的Ha图
plt.title('H$\\alpha$ map')  # 设置图表标题为Hα图

2. 获取Hβ和OIII图#

这更具挑战性，因为这些谱线彼此靠得很近。理想情况下，应该优先考虑迭代过程，但在这里我们使用Specutils的双高斯成分拟合，方法与Hα相似。

filt = 'f150w'  # 定义滤波器类型为'f150w'

file_2d = f'{DIR_DATA}l3_nis_{filt}_{grism}_s{id}_cal.fits'  # 构建2D数据文件路径

hdu_2d = fits.open(file_2d)  # 打开FITS文件并读取数据

# 对于光栅方向进行对齐
#   - x方向 = 色散（波长）方向
#   - y方向 = 横向色散
# 在这个笔记本中进行处理

if grism == 'G150C':  # 如果光栅类型为'G150C'
    # 如果光谱是水平的；
    data_2d = hdu_2d[ndither*7+1].data  # 读取数据部分
    dq_2d = hdu_2d[ndither*7+2].data  # 读取数据质量部分
    err_2d = hdu_2d[ndither*7+3].data  # 读取误差部分
    wave_2d = hdu_2d[ndither*7+4].data  # 读取波长部分
else:  # 如果光栅类型不是'G150C'
    data_2d = rotate(hdu_2d[ndither*7+1].data, 90)  # 将数据旋转90度
    dq_2d = rotate(hdu_2d[ndither*7+2].data, 90)  # 将数据质量部分旋转90度
    err_2d = rotate(hdu_2d[ndither*7+3].data, 90)  # 将误差部分旋转90度
    wave_2d = rotate(hdu_2d[ndither*7+4].data, 90)  # 将波长部分旋转90度

# !! 注意提取的光谱波长方向已翻转 !!
plt.imshow(data_2d[:, ::-1])  # 显示数据图像，反转列以修正波长方向
plt.title(f'{filt}')  # 设置图像标题为滤波器类型

在上面的图中，您可以看到Oiii双线和Hbeta。#

2a. 获取连续谱；#

yy = 10  # 示例行，选择y轴的索引为10

spec_unit = u.dimensionless_unscaled  # 设置光谱单位为无量纲

# 创建掩膜，选择特定波长范围内的线
mask_line = ((wave_2d[yy, :] > 1.35) & (wave_2d[yy, :] < 1.48)) | ((wave_2d[yy, :] > 1.6) & (wave_2d[yy, :] < 1.65))

# 使用掩膜选择的波长和数据创建光谱对象
obs = Spectrum1D(spectral_axis=wave_2d[yy, :][mask_line]*u.um, flux=data_2d[yy, :][mask_line]*spec_unit)

# 拟合连续谱，使用7阶切比雪夫多项式模型
cont = continuum.fit_generic_continuum(obs, model=Chebyshev1D(7))

# 绘制观测数据，y轴为yy的索引
plt.plot(wave_2d[yy, :], data_2d[yy, :], color='r', label=f'Observed at y={yy}')

# 绘制拟合的连续谱
plt.plot(wave_2d[yy, :][mask_line]*u.um, cont(wave_2d[yy, :][mask_line]*u.um), color='b', label='Fitted continuum')

# 添加图例
plt.legend(loc=0)

# 设置x轴范围
plt.xlim(1.2, 1.8)

# 沿y轴重复此操作;

flux_cont2d_150 = data_2d[:, :] * 0  # 初始化一个与data_2d相同形状的数组，所有元素为0

for yy in range(len(data_2d[:, 0])):  # 遍历data_2d的每一行（即每个y坐标）

    # 创建一个掩膜，选择波长在1.35到1.48微米和1.6到1.65微米之间的元素
    mask_line = ((wave_2d[yy, :] > 1.35) & (wave_2d[yy, :] < 1.48)) | ((wave_2d[yy, :] > 1.6) & (wave_2d[yy, :] < 1.65))

    # 创建一个Spectrum1D对象，包含选定波长范围内的光谱数据
    obs = Spectrum1D(spectral_axis=wave_2d[yy, :][mask_line]*u.um, flux=data_2d[yy, :][mask_line]*spec_unit)

    # 使用Chebyshev1D模型拟合光谱的连续背景
    cont = continuum.fit_generic_continuum(obs, model=Chebyshev1D(7))

    # 将拟合的连续背景值存储到flux_cont2d_150的相应行中
    flux_cont2d_150[yy, :] = cont(wave_2d[yy, :]*u.um).value

# 使用imshow函数显示flux_cont2d_150的图像，y轴翻转
plt.imshow(flux_cont2d_150[:, ::-1])

line_2d = data_2d[:, :] - flux_cont2d_150[:, :]  # 从二维数据中减去150波段的连续光通量

plt.imshow(line_2d[:, ::-1])  # 显示反转后的二维数据图像

2b. 拟合发射线，使用多组分高斯函数；#

yy = 10  # 设置要处理的光谱行索引

# 拟合光谱

con = (1.4 < wave_2d[yy, :]) & (wave_2d[yy, :] < 1.65)  # 选择波长范围在1.4到1.65微米之间的像素

spectrum_cut = Spectrum1D(flux=line_2d[yy, :][con]*spec_unit,  # 创建一个Spectrum1D对象，包含选定波长范围内的光谱数据
                          spectral_axis=wave_2d[yy, :][con]*u.um)  # 设置光谱轴为选定波长范围

# !!! 可能需要对初始值进行一些调整，以成功运行拟合；

# 对于Hb（氢β线）

g1_init = models.Gaussian1D(amplitude=100*spec_unit, mean=1.51*u.um, stddev=0.009*u.um)  # 初始化Hb的高斯模型

# 对于O3蓝（氧III蓝线）

g2_init = models.Gaussian1D(amplitude=80.*spec_unit, mean=1.53*u.um, stddev=0.006*u.um)  # 初始化O3蓝的高斯模型

# 对于O3红（氧III红线）

g3_init = models.Gaussian1D(amplitude=200.*spec_unit, mean=1.55*u.um, stddev=0.006*u.um)  # 初始化O3红的高斯模型

g123_fit = fit_lines(spectrum_cut, g1_init+g2_init+g3_init)  # 拟合组合的高斯模型到光谱数据
y_fit = g123_fit(wave_2d[yy, :]*u.um)  # 计算拟合后的光谱值

print(g123_fit)  # 输出拟合结果

# 单独绘制图形？

plt.plot(wave_2d[yy, :], line_2d[yy, :], marker='.', ls='', color='r', label=f'Observed at y={yy}')  # 绘制观测数据，y轴位置为yy，红色点标记

plt.plot(wave_2d[yy, :], y_fit, color='b', label='Fit', zorder=-2, alpha=0.4, lw=5)  # 绘制拟合曲线，蓝色，透明度0.4，线宽5

y_fit1 = g123_fit[0](wave_2d[yy, :]*u.um)  # 计算第一个拟合函数的值

plt.plot(wave_2d[yy, :], y_fit1, color='g', label='1')  # 绘制第一个拟合曲线，绿色

y_fit1 = g123_fit[1](wave_2d[yy, :]*u.um)  # 计算第二个拟合函数的值

plt.plot(wave_2d[yy, :], y_fit1, color='orange', label='2')  # 绘制第二个拟合曲线，橙色

y_fit1 = g123_fit[2](wave_2d[yy, :]*u.um)  # 计算第三个拟合函数的值

plt.plot(wave_2d[yy, :], y_fit1, color='purple', label='3')  # 绘制第三个拟合曲线，紫色

plt.xlim(1.4, 1.6)  # 设置x轴范围为1.4到1.6

plt.title('Single fit peak')  # 设置图表标题

plt.grid(True)  # 显示网格

plt.legend(loc=2)  # 显示图例，位置在左上角

# 打印g123_fit列表中的第一个元素
print(g123_fit[0])  

# 打印g123_fit列表中的第二个元素
print(g123_fit[1])  

# 打印g123_fit列表中的第三个元素
print(g123_fit[2])

拟合中央阵列效果良好。沿y轴重复此操作，并获取发射线图，方式与Hα相同。#

# 切割Hb和Oiii线图

rsq = data_2d.shape[0]  # 获取数据的维度

cut_hb = np.zeros((rsq, rsq), 'float32')  # 初始化Hb切割图像数组

cut_o3b = np.zeros((rsq, rsq), 'float32')  # 初始化O3蓝切割图像数组

cut_o3r = np.zeros((rsq, rsq), 'float32')  # 初始化O3红切割图像数组

zin = 2.1  # 从Notebook No.2获得的红移估计值，通过交叉相关得到

lamcen_hb = 0.4862680 * (1. + zin)  # 计算Hb线的中心波长

lamcen_o3b = 0.4960295 * (1. + zin)  # 计算O3蓝线的中心波长

lamcen_o3r = 0.5008240 * (1. + zin)  # 计算O3红线的中心波长

for yy in range(len(data_2d[:, 0])):  # 遍历每一行数据

    # 拟合光谱

    con = (1.4 < wave_2d[yy, :]) & (wave_2d[yy, :] < 1.65)  # 选择波长范围

    spectrum_cut = Spectrum1D(flux=line_2d[yy, :][con]*spec_unit,  # 创建光谱对象
                              spectral_axis=wave_2d[yy, :][con]*u.um)  # 设置光谱轴

    # !!! 可能需要对初始值进行一些调整，以成功运行拟合

    # 对于Hb线

    g1_init = models.Gaussian1D(amplitude=100*spec_unit, mean=1.51*u.um, stddev=0.009*u.um)  # 初始化Hb高斯模型

    # 对于O3蓝线

    g2_init = models.Gaussian1D(amplitude=80.*spec_unit, mean=1.53*u.um, stddev=0.006*u.um)  # 初始化O3蓝高斯模型

    # 对于O3红线

    g3_init = models.Gaussian1D(amplitude=200.*spec_unit, mean=1.55*u.um, stddev=0.006*u.um)  # 初始化O3红高斯模型

    g123_fit = fit_lines(spectrum_cut, g1_init+g2_init+g3_init)  # 拟合所有高斯模型

    y_fit = g123_fit(wave_2d[yy, :]*u.um)  # 计算拟合结果

    # 这必须在每个y像素上进行，因为波长数组可能会倾斜

    index_lamcen_hb = np.argmin(np.abs(lamcen_hb - wave_2d[yy, :]))  # 找到Hb中心波长的索引

    cut_hb[yy, :] = g123_fit[0](wave_2d[yy, :]*u.um)[int(index_lamcen_hb-rsq/2.):int(index_lamcen_hb+rsq/2.)]  # 提取Hb切割图像

    index_lamcen_o3b = np.argmin(np.abs(lamcen_o3b - wave_2d[yy, :]))  # 找到O3蓝中心波长的索引

    cut_o3b[yy, :] = g123_fit[1](wave_2d[yy, :]*u.um)[int(index_lamcen_o3b-rsq/2.):int(index_lamcen_o3b+rsq/2.)]  # 提取O3蓝切割图像

    index_lamcen_o3r = np.argmin(np.abs(lamcen_o3r - wave_2d[yy, :]))  # 找到O3红中心波长的索引

    cut_o3r[yy, :] = g123_fit[2](wave_2d[yy, :]*u.um)[int(index_lamcen_o3r-rsq/2.):int(index_lamcen_o3r+rsq/2.)]  # 提取O3红切割图像

plt.imshow(cut_hb)  # 显示切割后的Hβ图像

plt.title('H$\\beta$ map')  # 设置图像标题为Hβ图

plt.imshow(cut_o3b)  # 显示Oiii 4960的图像数据

plt.title('Oiii 4960')  # 设置图像标题为'Oiii 4960'

plt.imshow(cut_o3r)  # 显示切割后的Oiii 5008图像

plt.title('Oiii 5008')  # 设置图像标题为'Oiii 5008'

摘要；#

如上所示，除了中心区域以外的区域看起来不太对。这主要是由于多组分拟合失败，特别是对于Oiii双线。

为了改善拟合，可以采取以下措施：

在每个y轴上反复拟合并检查拟合，直到收敛为止，
使用MCMC进行更深入的拟合，这也可以通过设置先验来固定两个Oiii线的比率。
或者减少组分的数量，特别是在源位置边缘的Oiii双线区域，那里线条是混合的。

WFSS 光谱 第三部分：发射线图谱

Contents

WFSS 光谱 第三部分：发射线图谱#

引言#