MIRI LRS 窄缝光谱学：使用 JWST 管道进行光谱提取

MIRI LRS 窄缝光谱学：使用 JWST 管道进行光谱提取#

2023年7月

用例： 使用 JWST 校准管道提取窄缝光谱的光谱数据。

数据： 公开可用的科学数据

工具： jwst, matplotlib, astropy.

跨仪器： NIRSpec, MIRI.

文档： 本笔记本是 STScI 更大后处理数据分析工具生态系统的一部分，可以直接从 JDAT Notebook Github 目录下载。

引言：JWST 校准管道中的光谱提取#

JWST 校准管道对所有光谱数据执行光谱提取，使用基本的默认假设，这些假设经过调整以生成大多数科学案例所需的准确校准光谱。此默认方法是简单的固定宽度箱体提取，其中光谱在有效波长范围内沿交叉色散轴的多个像素上进行求和。在光谱的每个像素上应用孔径修正，以考虑由于有限宽度孔径而损失的通量。

extract_1d 步骤使用以下输入进行其算法：

光谱提取参考文件：这是一个 JSON 格式的文件，可以作为参考文件从 JWST CRDS 系统获取。
边界框：assign_wcs 步骤将边界框定义附加到数据上，定义了可用有效校准的区域。我们将在下面演示如何可视化该区域。

然而，extract_1d 步骤具有执行更复杂光谱提取的能力，这需要手动编辑参数并重新运行管道步骤。

目标#

本笔记本将演示如何使用不同设置重新运行光谱提取步骤，以展示 JWST 校准管道的能力。

假设#

我们将在重新采样的校准光谱图像上演示光谱提取方法。基本演示和两个示例运行在 Level 3 数据上，其中节点曝光已合并为单个光谱图像。两个示例将使用 Level 2b 数据 - 其中一个是节点曝光。

测试数据#

本笔记本使用的数据是 Jha 等人在 PID 2072（Obs 1）中观察的 Ia 型超新星 SN2021aefx。这些数据在零独占访问期内获取，并在 Kwok et al 2023 中发布。您可以从这个 Box 文件夹获取数据，我们建议您将文件放置在本存储库的 data/ 文件夹中，或在运行之前更改笔记本中的目录设置。

当然，您也可以使用自己的数据，而不是演示数据。

JWST 管道版本和 CRDS 上下文#

本笔记本是使用校准管道版本 1.10.2 编写的。我们将 CRDS 上下文明确设置为 1089，以匹配 MAST 中当前最新版本。如果您使用不同的管道版本或 CRDS 上下文，请阅读相关的发行说明（管道的说明，CRDS 的说明）以了解可能相关的更改。

导入包#

astropy.io 用于访问FITS文件
os 用于管理系统路径
matplotlib 用于绘制数据
urllib 用于下载数据
tarfile 用于解压数据
numpy 用于基本数组操作
jwst 用于运行JWST管道和处理数据产品
json 用于处理json文件
crds 用于处理JWST参考文件

# 首先设置CRDS变量

import os  # 导入os模块以便进行操作系统相关的功能

# 设置CRDS上下文为jwst_1089.pmap
os.environ['CRDS_CONTEXT'] = 'jwst_1089.pmap'

# 设置CRDS缓存路径为用户主目录下的crds_cache文件夹
os.environ['CRDS_PATH'] = os.environ['HOME'] + '/crds_cache'

# 设置CRDS服务器的URL
os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds.stsci.edu'

# 打印CRDS缓存位置
print(f'CRDS cache location: {os.environ["CRDS_PATH"]}')

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

import urllib.request  # 导入用于处理URL请求的库

import tarfile  # 导入用于处理tar文件的库

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

import matplotlib.pyplot as plt  # 导入Matplotlib库，用于绘图

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

import astropy.units as u  # 导入Astropy库中的单位模块

from astropy.modeling import models, fitting  # 导入Astropy建模和拟合模块

import jwst  # 导入JWST相关库

from jwst import datamodels  # 从JWST库中导入数据模型模块

from jwst.extract_1d import Extract1dStep  # 从JWST库中导入一维提取步骤模块

from matplotlib.patches import Rectangle  # 从Matplotlib库中导入矩形补丁类

import json  # 导入JSON库，用于处理JSON数据

import crds  # 导入CRDS库，用于处理数据校准

print(f'Using JWST calibration pipeline version {jwst.__version__}')  # 打印当前使用的JWST校准管道版本

import os  # 导入操作系统模块
import urllib.request  # 导入urllib请求模块
import tarfile  # 导入tarfile模块，用于处理tar文件

data_tar_url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MIRI_LRS_notebook/data.tar.gz'  # 数据压缩包的URL地址

# 下载并解压数据（如果需要）

if not os.path.exists("data.tar.gz"):  # 如果当前目录下不存在data.tar.gz文件
    print("Downloading Data")  # 打印下载数据的提示
    urllib.request.urlretrieve(data_tar_url, 'data.tar.gz')  # 从URL下载数据并保存为data.tar.gz

if not os.path.exists("data/"):  # 如果当前目录下不存在data文件夹
    print("Unpacking Data")  # 打印解压数据的提示
    with tarfile.open('./data.tar.gz', "r:gz") as tar:  # 打开data.tar.gz文件
        tar.extractall(filter='data')  # 解压所有文件到data文件夹中

1. 第三级数据产品 #

让我们开始绘制主要的默认第三级输出产品：

s2d 文件：这是由共同加权的重采样单个点曝光构建的二维图像。
x1d 文件：这是从第三级 s2d 文件中提取的1维光谱。

s2d 图像显示了一个明亮的中央轨迹，旁边有两个负轨迹。这些轨迹是由于点曝光的组合而产生的，每个曝光也包含一个正轨迹和一个负轨迹，因为它们在进行背景减除时是相互减去的。

我们将x轴的短波长端限制在5微米，因为在此波长以下我们的校准非常差。第三级光谱是从重采样、抖动组合和校准的曝光中提取的。

l3_s2d_file = 'data/jw02072-o001_t010_miri_p750l_s2d_1089.fits'  # 定义JWST数据文件路径

l3_s2d = datamodels.open(l3_s2d_file)  # 打开JWST数据文件

fig, ax = plt.subplots(figsize=[2, 8])  # 创建一个2x8英寸的图形和坐标轴

im2d = ax.imshow(l3_s2d.data, origin='lower', aspect='auto', cmap='gist_gray')  # 显示2D光谱图像，设置原点在下方，自动调整纵横比，使用灰色调色板

ax.set_xlabel('column')  # 设置x轴标签为'column'

ax.set_ylabel('row')  # 设置y轴标签为'row'

ax.set_title('SN2021aefx - Level 3 resampled 2D spectral image')  # 设置图形标题

fig.colorbar(im2d)  # 添加颜色条以显示数据值的颜色映射

fig.show()  # 显示图形

l3_file = 'data/jw02072-o001_t010_miri_p750l_x1d_1089.fits'  # 定义JWST数据文件路径

l3_spec = datamodels.open(l3_file)  # 打开JWST数据文件并加载数据模型

fig2, ax2 = plt.subplots(figsize=[12, 4])  # 创建一个12x4英寸的图形和坐标轴

ax2.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['FLUX'])  # 绘制波长与通量的关系图

ax2.set_xlabel('wavelength (um)')  # 设置x轴标签为波长（微米）

ax2.set_ylabel('flux (Jy)')  # 设置y轴标签为通量（焦耳每秒每平方米）

ax2.set_title('SN2021aefx - Level 3 spectrum in MAST (pmap 1089)')  # 设置图形标题

ax2.set_xlim(5., 14.)  # 设置x轴范围为5到14微米

fig2.show()  # 显示图形

光谱提取参考文件 #

告诉 extract_1d 算法使用哪些参数的参考文件是一个使用 json 格式的文本文件，该文件可以在 CRDS 中找到。提取中使用的第二个参考文件是孔径校正；它用于修正因所用提取孔径大小而导致的波长函数下的光通量损失。您可以使用 x1d 文件的数模属性来检查算法调用了哪个提取参考文件。

我们下面展示如何以编程方式检查文件，以查看用于生成上述默认 Level 3 光谱的孔径。注意：这个 json 文件可以很容易地用简单的文本编辑器打开和编辑。

extract_1d 参考文件的完整文档可以在这里找到。我们建议您仔细阅读此页面及其中的任何链接，以了解文件中的参数是如何应用于数据的。

# 打印用于光谱提取的参考文件名称
print(f'Spectral extraction reference file used: {l3_spec.meta.ref_file.extract1d.name}')

file_path = 'data/jw02072-o001_t010_miri_p750l_x1d_1089.fits'  # 定义JWST数据文件的路径

with fits.open(file_path) as hdul:  # 打开FITS文件

    header = hdul[0].header  # 获取FITS文件的头部信息

    json_ref_default = crds.getreferences(header)['extract1d']  # 获取与头部信息相关的参考数据

    with open(json_ref_default) as json_ref:  # 打开参考数据的JSON文件

        x1dref_default = json.load(json_ref)  # 读取JSON文件内容

        print('Settings for SLIT data: {}'.format(x1dref_default['apertures'][0]))  # 打印SLIT数据的设置

        print('  ')  # 打印空行以便于阅读

        print('Settings for SLITLESS data: {}'.format(x1dref_default['apertures'][1]))  # 打印SLITLESS数据的设置

让我们来看看这个文件中的内容。

id: 识别标签，在本例中指定将应用于的曝光类型的参数。
region_type: 可选项，如果包含，必须设置为 ‘target’。
disp_axis: 定义色散方向（1 表示 x 轴，2 表示 y 轴）。对于 MIRI LRS，这应该始终设置为 2。
xstart (int): 水平方向上的第一个像素（x 轴；0 索引）
xstop (int): 水平方向上的最后一个像素（x 轴；0 索引；限制为包含）
bkg_order:
use_source_posn (True/False): 如果为 True，这将使用目标坐标在场中定位目标，并将提取孔径偏移到该位置。我们建议将其设置为 False。
bkg_order: 用于背景拟合的多项式阶数。如果未提供伴随参数 bkg_coeff，则不会执行背景拟合。对于 MIRI LRS 裂缝数据，默认的背景减法是在 Spec2Pipeline 中通过相互减去 nod 曝光来实现的。

对于 MIRI LRS，色散方向是垂直的（即 disp_axis = 2），提取孔径的宽度通过坐标 xstart 和 xstop 指定。如果未提供坐标 ystart 和 ystop，则光谱将提取整个 s2d 裁剪区域的高度。我们可以在 Level 3 s2d 文件上说明默认提取参数。

# 从x1dref_default字典中获取第一个光圈的x起始位置
xstart = x1dref_default['apertures'][0]['xstart']

# 从x1dref_default字典中获取第一个光圈的x结束位置
xstop = x1dref_default['apertures'][0]['xstop']

# 获取l3_s2d数据的高度（行数）
ap_height = np.shape(l3_s2d.data)[0]

# 计算光圈的宽度
ap_width = xstop - xstart + 1

# 创建一个矩形对象，表示光圈区域
x1d_rect = Rectangle(xy=(xstart, 0), width=ap_width, height=ap_height, angle=0., edgecolor='red',
                     facecolor='None', ls='-', lw=1.5)

# 创建一个新的图形和坐标轴，设置图形大小
fig, ax = plt.subplots(figsize=[2, 8])

# 在坐标轴上显示l3_s2d数据，设置原点在下方，纵横比自动，使用灰色调色板
im2d = ax.imshow(l3_s2d.data, origin='lower', aspect='auto', cmap='gist_gray')

# 将光圈矩形添加到坐标轴上
ax.add_patch(x1d_rect)

# 设置x轴标签
ax.set_xlabel('column')

# 设置y轴标签
ax.set_ylabel('row')

# 设置图形标题
ax.set_title('SN2021aefx - Level 3 resampled 2D spectral image')

# 添加颜色条以表示数据值
fig.colorbar(im2d)

# 显示图形
fig.show()

示例 1：改变提取宽度 #

在这个示例中，我们演示如何改变提取宽度，超出默认值。我们将提取12个像素，而不是8个，同时保持光谱孔径居中于轨迹上。

我们将在这个笔记本中使用Python修改json文件中的值，但该文件也可以直接在文本编辑器中进行编辑。

xstart2 = xstart - 2  # 将xstart减去2，得到新的xstart值

xstop2 = xstop + 2    # 将xstop加上2，得到新的xstop值

print('New xstart, xstop values = {0},{1}'.format(xstart2, xstop2))  # 打印新的xstart和xstop值

with open(json_ref_default) as json_ref:  # 打开默认的JSON参考文件

    x1dref_default = json.load(json_ref)  # 读取JSON文件内容并加载到x1dref_default变量中

    x1dref_ex1 = x1dref_default.copy()  # 复制默认的参考数据到x1dref_ex1变量中

    x1dref_ex1['apertures'][0]['xstart'] = xstart2  # 更新x1dref_ex1中的xstart值为新的xstart2

    x1dref_ex1['apertures'][0]['xstop'] = xstop2  # 更新x1dref_ex1中的xstop值为新的xstop2

with open('x1d_reffile_example1.json', 'w') as jsrefout:  # 打开一个新的JSON文件以写入更新后的数据

    json.dump(x1dref_ex1, jsrefout, indent=4)  # 将更新后的数据以JSON格式写入文件，缩进为4个空格

ap_width2 = xstop2 - xstart2 + 1  # 计算第二个光圈的宽度

# 创建第一个光圈的矩形对象，设置位置、宽度、高度、颜色等属性
x1d_rect1 = Rectangle(xy=(xstart, 0), width=ap_width, height=ap_height, angle=0., edgecolor='red',
                      facecolor='None', ls='-', lw=1, label='8-px aperture (default)')  # 默认8像素光圈

# 创建第二个光圈的矩形对象，设置位置、宽度、高度、颜色等属性
x1d_rect2 = Rectangle(xy=(xstart2, 0), width=ap_width2, height=ap_height, angle=0., edgecolor='cyan',
                      facecolor='None', ls='-', lw=1, label='12-px aperture')  # 修改后的12像素光圈

# 创建一个新的图形和坐标轴，设置图形大小
fig4, ax4 = plt.subplots(figsize=[2, 8])

# 在坐标轴上显示二维图像，设置原点位置、纵横比和颜色映射
im2d = ax4.imshow(l3_s2d.data, origin='lower', aspect='auto', cmap='gist_gray')

# ax4.add_collection(aps_collection)  # 如果需要，可以添加光圈集合（注释掉）

# 在坐标轴上添加第一个光圈矩形
ax4.add_patch(x1d_rect1)

# 在坐标轴上添加第二个光圈矩形
ax4.add_patch(x1d_rect2)

# 设置坐标轴的标签
ax4.set_xlabel('column')  # 设置x轴标签为'column'
ax4.set_ylabel('row')     # 设置y轴标签为'row'
ax4.set_title('Example 1: Default vs modified extraction aperture')  # 设置图形标题

# 添加图例，位置为3
ax4.legend(loc=3)

# 为图形添加颜色条
fig.colorbar(im2d)

# 显示图形
fig.show()

接下来，我们运行光谱提取步骤，使用这个修改过的参考文件。注意：当单独运行某个步骤时，文件名后缀与我们完整运行Spec3Pipeline时不同。提取的光谱现在将包含extract1dstep.fits在文件名中。我们传递给步骤调用的自定义参数如下：

output_file：我们为这个示例提供了一个自定义的输出文件名（包括输出文件名使得save_results参数变得不再必要）
override_extract1d：在这里，我们提供了上面创建的自定义参考文件的名称

我们将输出与默认提取的产品进行比较。我们预计光谱几乎是相同的；在较长波长处可能会出现差异，因为我们的路径损失校正在这个低信噪比区域的校准较差。

# 调用Extract1dStep类的call方法，进行一维提取
sp3_ex1 = Extract1dStep.call(
    l3_s2d,  # 输入的二级科学数据
    output_dir='data/',  # 输出目录设置为'data/'
    output_file='lrs_slit_extract_example1',  # 输出文件名设置为'lrs_slit_extract_example1'
    override_extract1d='x1d_reffile_example1.json'  # 覆盖提取配置文件，使用'x1d_reffile_example1.json'
)

# 打印变量 sp3_ex1 的值
print(sp3_ex1)

fig5, ax5 = plt.subplots(figsize=[12, 4])  # 创建一个12x4英寸的图形和坐标轴

ax5.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['FLUX'], label='8-px aperture')  # 绘制8像素光圈的光谱数据

ax5.plot(sp3_ex1.spec[0].spec_table['WAVELENGTH'], sp3_ex1.spec[0].spec_table['FLUX'], label='12-px aperture')  # 绘制12像素光圈的光谱数据

ax5.set_xlabel('wavelength (um)')  # 设置x轴标签为波长（微米）

ax5.set_ylabel('flux (Jy)')  # 设置y轴标签为通量（焦耳/秒）

ax5.set_title('Example 1: Difference aperture sizes')  # 设置图表标题

ax5.set_xlim(5., 14.)  # 设置x轴的范围为5到14微米

ax5.legend()  # 显示图例

fig5.show()  # 显示图形

示例 2：改变光圈位置 #

在这个示例中，我们将演示如何在狭缝的不同位置进行光谱提取。一个好的应用场景是从其中一个偏移曝光中提取光谱，这在将偏移合并到Spec3Pipeline之前进行。我们将使用来自Spec2Pipeline的s2d输出，并提取光谱。在偏移1的曝光中，我们看到光谱峰位于第13列（0索引），并提取一个默认的8像素固定宽度光圈。

l2_s2d_file = 'data/jw02072001001_06101_00001_mirimage_s2d.fits'  # 定义JWST二级数据文件的路径

l2_s2d = datamodels.open(l2_s2d_file)  # 使用datamodels库打开该文件并加载数据

xstart3 = 9  # 定义x轴起始值

xstop3 = 17  # 定义x轴结束值

with open(json_ref_default) as json_ref:  # 打开默认的JSON参考文件
    x1dref_default = json.load(json_ref)  # 加载JSON数据到字典中
    x1dref_ex2 = x1dref_default.copy()  # 复制默认数据以便修改
    x1dref_ex2['apertures'][0]['xstart'] = xstart3  # 更新xstart值
    x1dref_ex2['apertures'][0]['xstop'] = xstop3  # 更新xstop值

with open('x1d_reffile_example2.json', 'w') as jsrefout:  # 打开输出文件以写入
    json.dump(x1dref_ex2, jsrefout, indent=4)  # 将修改后的数据写入JSON文件，格式化为4个空格缩进

ap_width3 = xstop3 - xstart3 + 1  # 计算光圈的宽度

# 创建一个矩形对象，表示光圈的位置和大小
x1d_rect3 = Rectangle(xy=(xstart3, 0), width=ap_width3, height=ap_height, angle=0., edgecolor='red',
                      facecolor='None', ls='-', lw=1, label='8-px aperture at nod 1 location')

# 创建一个新的图形和坐标轴
fig6, ax6 = plt.subplots(figsize=[2, 8])

# 在坐标轴上显示二维图像，设置原点在左下角，保持宽高比，使用灰色调色板
im2d = ax6.imshow(l2_s2d.data, origin='lower', aspect='auto', cmap='gist_gray')

# 在坐标轴上添加光圈矩形
ax6.add_patch(x1d_rect3)

# 设置x轴和y轴的标签
ax6.set_xlabel('column')  # x轴标签为'column'
ax6.set_ylabel('row')     # y轴标签为'row'

# 设置图形的标题
ax6.set_title('Example 2: Different aperture location')  # 图形标题

# 添加图例，位置设置为3
ax6.legend(loc=3)

# 为图形添加颜色条
fig6.colorbar(im2d)

# 显示图形
fig6.show()

# 调用Extract1dStep类的call方法进行一维提取
sp2_ex2 = Extract1dStep.call(
    l2_s2d_file,  # 输入的二级空间数据文件
    output_dir='data/',  # 输出目录设置为'data/'
    output_file='lrs_slit_extract_example2',  # 输出文件名设置为'lrs_slit_extract_example2'
    override_extract1d='x1d_reffile_example2.json'  # 指定覆盖的一维提取参考文件
)

让我们再次将输出与默认提取的产品进行对比绘图。我们预期这个1-nod光谱会更嘈杂，但与组合产品相比不会有显著差异。光谱中可能会有更多的坏像素，这些坏像素会表现为光谱中的尖峰或凹陷。

fig7, ax7 = plt.subplots(figsize=[12, 4])  # 创建一个12x4英寸的图形和坐标轴

ax7.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['FLUX'], label='default location (nods combined)')  # 绘制默认位置（合并nods）的光谱数据

ax7.plot(sp2_ex2.spec[0].spec_table['WAVELENGTH'], sp2_ex2.spec[0].spec_table['FLUX'], label='nod 1 location (single nod)')  # 绘制单个nod位置的光谱数据

ax7.set_xlabel('wavelength (um)')  # 设置x轴标签为波长（微米）

ax7.set_ylabel('flux (Jy)')  # 设置y轴标签为通量（Jy）

ax7.set_title('Example 2: Different aperture locations')  # 设置图形标题

ax7.set_xlim(5., 14.)  # 设置x轴范围为5到14微米

ax7.legend()  # 显示图例

fig7.show()  # 显示图形

示例 3：带背景减除的提取 #

对于 LRS（低分辨率光谱）狭缝观测，默认的背景减除策略是在 Spec2Pipeline 的 background 步骤中执行的；两个点位曝光相互减去，结果返回每个 2D 光谱图像，包含一个正向和一个负向的轨迹，并且背景已被减除。

然而，对于非标准情况或无狭缝的 LRS 数据，可以在 extract_1d 中作为光谱提取的一部分减去背景。在 extract_1d 参考文件中，我们可以传递特定的背景参数：

bkg_coeff（列表或浮点数列表）：用于作为背景的区域。这是进行背景减除所需的主要参数
bkg_fit（字符串）：背景计算的类型或方法。（例如：None、’poly’、’mean’ 或 ‘median’）
bkg_order（整数）：拟合背景区域的多项式阶数。如果 bkg_fit 未设置为 ‘poly’，则此参数将被忽略。
smoothing_length（奇数整数；可选）：用于平滑背景信号的箱形滤波器的宽度，沿色散方向。这在数据噪声较大的情况下可以提供更好的质量。

对于 bkg_fit 参数的 ‘poly’ 选项，将取背景区域中给定行的所有像素值，并对它们拟合一个阶数为 bkg_order 的多项式。此选项在背景存在梯度的情况下可能会很有用。

我们在这里使用的数据已经减去了背景，因此我们预计其影响将是最小的，但我们提供一个使用点位 1、级别 2b 光谱图像的演示。在此示例中，我们将从两个 4 列窗口计算背景，将 bkg_fit 设置为 ‘median’。

rows = [140, 200, 325]  # 定义要绘制的行号

fig8, ax8 = plt.subplots(figsize=[8, 4])  # 创建一个8x4英寸的图形和坐标轴

ncols = np.shape(l2_s2d.data)[1]  # 获取数据的列数

pltx = np.arange(ncols)  # 创建一个从0到ncols-1的数组，用于x轴

for rr in rows:  # 遍历每一行
    label = 'row {}'.format(rr)  # 创建当前行的标签

    ax8.plot(pltx, l2_s2d.data[rr, :], label=label)  # 绘制当前行的数据

# 添加垂直线，表示背景区域
ax8.axvline(x=1, ymin=0, ymax=1, ls='--', lw=1., color='coral', label='background regions')  # 在x=1处添加虚线

ax8.axvline(x=5, ymin=0, ymax=1, ls='--', lw=1., color='coral')  # 在x=5处添加虚线

ax8.axvline(x=39, ymin=0, ymax=1, ls='--', lw=1., color='coral')  # 在x=39处添加虚线

ax8.axvline(x=43, ymin=0, ymax=1, ls='--', lw=1., color='coral')  # 在x=43处添加虚线

ax8.legend()  # 显示图例

fig8.show()  # 显示图形

with open(json_ref_default) as json_ref:  # 打开默认的JSON参考文件

    x1dref_default = json.load(json_ref)  # 加载JSON数据到变量x1dref_default

    x1dref_ex3 = x1dref_default.copy()  # 复制默认的参考数据到新的变量x1dref_ex3

    x1dref_ex3['apertures'][0]['xstart'] = xstart3  # 设置第一个光圈的起始x坐标

    x1dref_ex3['apertures'][0]['xstop'] = xstop3  # 设置第一个光圈的结束x坐标

    x1dref_ex3['apertures'][0]['bkg_coeff'] = [[0.5], [4.5], [38.5], [43.5]]  # 设置背景系数

    x1dref_ex3['apertures'][0]['bkg_fit'] = 'median'  # 设置背景拟合方法为中位数

with open('x1d_reffile_example3.json', 'w') as jsrefout:  # 打开一个新的JSON文件以写入

    json.dump(x1dref_ex3, jsrefout, indent=4)  # 将修改后的数据写入新文件，格式化为4个空格缩进

# 调用Extract1dStep类的方法，进行一维提取
sp2_ex3 = Extract1dStep.call(
    l2_s2d_file,  # 输入的二级空间数据文件
    output_dir='data/',  # 输出目录
    output_file='lrs_slit_extract_example3',  # 输出文件名
    override_extract1d='x1d_reffile_example3.json'  # 覆盖提取配置文件
)

当 extract_1d 步骤执行背景减法时，背景光谱是输出产品的一部分，因此您可以检查被减去的内容。在下面的图中，我们可以看到，正如预期的那样，这个特定曝光的背景接近于零（除了噪声较大的长波长端），因为减法已经执行过。

fig9, ax9 = plt.subplots(nrows=2, ncols=1, figsize=[12, 4])  # 创建一个包含2行1列的子图，图形大小为12x4英寸

# ax9.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['FLUX'], label='default location (nods combined)')  # 注释掉的代码，绘制默认位置的光谱

ax9[0].plot(sp2_ex2.spec[0].spec_table['WAVELENGTH'], sp2_ex2.spec[0].spec_table['FLUX'], label='nod 1 spectrum - no bkg sub')  # 绘制未进行背景减除的nod 1光谱

ax9[0].plot(sp2_ex3.spec[0].spec_table['WAVELENGTH'], sp2_ex3.spec[0].spec_table['FLUX'], label='nod 1 spectrum - with bkg sub')  # 绘制已进行背景减除的nod 1光谱

ax9[1].plot(sp2_ex3.spec[0].spec_table['WAVELENGTH'], sp2_ex3.spec[0].spec_table['BACKGROUND'], label='background')  # 绘制背景光谱

ax9[1].set_xlabel('wavelength (um)')  # 设置第二个子图的x轴标签为“波长（微米）”

ax9[0].set_ylabel('flux (Jy)')  # 设置第一个子图的y轴标签为“通量（Jy）”

ax9[0].set_title('Example 3: Extraction with background subtraction')  # 设置第一个子图的标题

ax9[0].set_xlim(5., 14.)  # 设置第一个子图的x轴范围为5到14微米

ax9[1].set_xlim(5., 14.)  # 设置第二个子图的x轴范围为5到14微米

ax9[0].legend()  # 在第一个子图中显示图例

ax9[1].legend()  # 在第二个子图中显示图例

fig9.show()  # 显示图形

示例 4：锥形列提取 #

在这个示例中，我们将使用JWST校准管道在锥形列孔径中进行光谱提取。在提取参考文件中指定这一点的方法是使用src_soeff参数，而不是更简单的xstart和xstop设置。src_coeff参数可以采用多项式系数，而不是固定的像素值。在这个示例中，我们将定义一个对应于空间轮廓的3倍全宽半最大值（FWHM）的锥形列孔径。

提取孔径的多项式定义可以作为像素或波长的函数进行指定，这在independent_var参数中定义。

我们将使用预先测量的FWHM值作为波长的函数，拟合FWHM(\(\lambda\))轮廓的直线，并根据这个拟合设置提取参数。当然，FWHM也可以直接从数据中测量。

def calc_xap_fit():

    # 这些是从调试数据中测得的值。FWHM以角秒为单位。

    lam = [5.0, 7.5, 10.0, 12.0]  # 波长列表

    fwhm = [0.29, 0.3, 0.36, 0.42]  # FWHM值列表

    # 使用MIRI像素缩放0.11角秒/像素将角秒转换为像素

    fwhm_px = fwhm / (0.11*u.arcsec/u.pixel)  # 将FWHM转换为像素

    # 我们想提取3 * FWHM，这意味着在轨迹的两侧各1.5 * FWHM

    xap_pix = 1.5 * fwhm_px  # 计算提取的像素范围

    # 现在我们想对这些点拟合一条直线

    line_init = models.Linear1D()  # 初始化线性模型

    fit = fitting.LinearLSQFitter()  # 初始化线性最小二乘拟合器

    fitted_line = fit(line_init, lam, xap_pix.value)  # 拟合直线

    print(fitted_line)  # 打印拟合结果

    fig, ax = plt.subplots(figsize=[8, 4])  # 创建图形和坐标轴

    xplt = np.linspace(4.0, 14., num=50)  # 创建用于绘图的x值

    ax.plot(lam, xap_pix.value, 'rx', label='1.5 * FWHM(px)')  # 绘制1.5 * FWHM的点

    ax.plot(xplt, fitted_line(xplt), 'b-', label='best-fit line')  # 绘制拟合直线

    ax.set_xlabel('wavelength')  # 设置x轴标签

    ax.set_ylabel('px')  # 设置y轴标签

    ax.legend()  # 显示图例

    return fitted_line  # 返回拟合的直线

poly_pos = calc_xap_fit()  # 调用函数计算XAP拟合，返回结果赋值给poly_pos

print(poly_pos.slope, poly_pos.intercept)  # 打印poly_pos的斜率和截距

上述多项式定义了波长与提取像素数量之间的关系。为了确保提取位置位于光谱的中心位置，我们在截距值上加上轨迹的中心位置，即第30.5列。

在下一个单元格中，我们将这些参数提供给提取参考文件中的src_coeff参数。注意：如果同时存在src_coeff、xstart和xstop这三个参数，则src_coeff参数具有优先权；为清晰起见，我们将后两者从我们的参考文件中移除。

trace_cen = 30.5  # 设置中心轨迹值为30.5

with open(json_ref_default) as json_ref:  # 打开默认的JSON参考文件

    x1dref_default = json.load(json_ref)  # 读取JSON数据并加载到x1dref_default变量中

    x1dref_ex4 = x1dref_default.copy()  # 复制默认的参考数据到x1dref_ex4

    x1dref_ex4['apertures'][0]['xstart'] = None  # 将第一个光圈的起始x值设置为None

    x1dref_ex4['apertures'][0]['xstop'] = None  # 将第一个光圈的结束x值设置为None

    x1dref_ex4['apertures'][0]['independent_var'] = 'wavelength'  # 设置独立变量为波长

    # 设置源系数，包含两个点的坐标
    x1dref_ex4['apertures'][0]['src_coeff'] = [[-1*poly_pos.intercept.value + trace_cen, -1*poly_pos.slope.value], 
                                                [poly_pos.intercept.value + trace_cen, poly_pos.slope.value]]

with open('x1d_reffile_example4.json', 'w') as jsrefout:  # 打开一个新的JSON文件以写入

    json.dump(x1dref_ex4, jsrefout, indent=4)  # 将修改后的数据写入JSON文件，格式化为4个空格缩进

# 调用Extract1dStep类的call方法进行一维提取
sp3_ex4 = Extract1dStep.call(
    l3_s2d,  # 输入的二级数据
    output_dir='data/',  # 输出目录
    output_file='lrs_slit_extract_example4',  # 输出文件名
    override_extract1d='x1d_reffile_example4.json'  # 覆盖提取配置文件
)

fig10, ax10 = plt.subplots(figsize=[12, 4])  # 创建一个12x4英寸的图形和坐标轴

ax10.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['FLUX'], label='default fixed-width aperture')  # 绘制默认固定宽度光圈的光谱数据

ax10.plot(sp3_ex4.spec[0].spec_table['WAVELENGTH'], sp3_ex4.spec[0].spec_table['FLUX'], label='tapered column aperture')  # 绘制锥形光柱光圈的光谱数据

ax10.set_xlabel('wavelength (um)')  # 设置x轴标签为波长（微米）

ax10.set_ylabel('flux (Jy)')  # 设置y轴标签为通量（焦耳每秒每平方米）

ax10.set_title('Example 4: Tapered column vs. fixed-width extraction aperture')  # 设置图表标题

ax10.set_xlim(5., 14.)  # 设置x轴范围为5到14微米

ax10.legend()  # 显示图例

fig10.show()  # 显示图形

输出光谱还包含了一个关于用于提取的像素数量与波长的关系。让我们也将其可视化。

fig11, ax11 = plt.subplots(figsize=[12, 4])  # 创建一个12x4英寸的图形和坐标轴

ax11.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['NPIXELS'], label='default fixed-width aperture')  # 绘制默认固定宽度光圈的像素数量与波长的关系

ax11.plot(sp3_ex4.spec[0].spec_table['WAVELENGTH'], sp3_ex4.spec[0].spec_table['NPIXELS'], label='tapered column aperture')  # 绘制锥形光圈的像素数量与波长的关系

ax11.set_xlabel('wavelength (um)')  # 设置x轴标签为波长（微米）

ax11.set_ylabel('number of pixels')  # 设置y轴标签为像素数量

ax11.set_title('Example 4: Number of pixels extracted')  # 设置图表标题

ax11.set_xlim(5., 14.)  # 设置x轴的范围为5到14微米

ax11.legend()  # 显示图例

fig11.show()  # 显示图形

摘要#

我们希望这个笔记本能帮助您理解JWST（詹姆斯·韦伯太空望远镜）在光谱提取方面的校准能力。上述示例并不是所有可能性的详尽列表：不同的源和背景提取方法可以结合起来进行更复杂的提取操作。

如果您有任何问题、意见或希望进一步演示这些功能的请求，请联系 JWST 帮助中心。