使用 NSClean 清理 NIRSpec BOTS 产品中的残余 1/f 噪声#
笔记本目标#
本笔记本的目标是通过去除残余的 1/f 噪声来生成清理后的 BOTS (_rateint.fits) 文件。这些清理后的文件将作为 3 (TSO3Pipeline) 管道的输入。
目录#
1. 介绍 #
JWST NIRSpec 仪器有许多特性和特点，观察者在规划观测和解释数据时应当注意。其中一个显著特征是在提取的 1-D 光谱中出现负值和/或过剩通量，通常伴有不规则的波长依赖性起伏。这种伪影的原因是相关噪声，称为 1/f 噪声，源于低水平的探测器热不稳定性，在 2-D 计数率图像中表现为垂直条带，尤其是在 NRS2 探测器的曝光中。虽然 IRS2 读出模式减少了这种效应，但并未完全消除。
为了解决这个问题，JWST 科学校准管道集成了由 Bernard Rauscher 开发的外部包，称为 NSClean，并在 Spec2Pipeline 中通过 NSCleanStep 实现。该算法利用探测器的暗区在傅里叶空间中拟合背景模型。它需要一个输入掩模来识别探测器的所有暗区。掩模越全面、越完整，背景拟合效果越好。
在本笔记本中，我们将使用集成到管道中的 NSClean 算法，利用默认参数动态生成的掩模来去除 1/f 噪声。在某些情况下，这个掩模可能不够完整或过于严格，无法实现最佳的噪声去除。为了解决这个问题，我们演示了如何手动修改默认掩模，以及如何通过调整 NSCleanStep 参数 创建替代掩模。如有需要，请参阅 NSClean 文档，获取有关手动创建自定义掩模的一些建议。
本笔记本利用了在 BOTS G395H 光栅中观察到的 WASP-39b 的过境数据，作为 JWST 早期发布科学计划 ERS-1366 观测 3 的示例。
2. 导入库 #
# ---------- 设置CRDS环境变量 ----------
import os # 导入os模块，用于操作系统相关功能
import jwst # 导入jwst模块，用于JWST相关功能
# 设置CRDS上下文为指定的.pmap文件
os.environ['CRDS_CONTEXT'] = 'jwst_1215.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"]}')
# 打印JWST校准管道的版本信息
print("JWST Calibration Pipeline Version={}".format(jwst.__version__))
# 打印当前操作的CRDS上下文（注释掉的代码）
# print("Current Operational CRDS Context = {}".format(crds.get_default_context()))
# ------ General Imports ------
import numpy as np # 导入NumPy库，用于数值计算
import time as tt # 导入time库，用于时间相关操作
import logging # 导入logging库，用于记录日志
import warnings # 导入warnings库，用于警告管理
import json # 导入json库，用于处理JSON数据
# ------ JWST Calibration Pipeline Imports ------
from crds.client import api # 从CRDS客户端导入API模块，用于数据检索
from stpipe import crds_client # 从stpipe导入CRDS客户端，用于数据处理
from jwst.pipeline.calwebb_spec2 import Spec2Pipeline # 导入JWST的Spec2Pipeline，用于光谱数据处理
# ------ Plotting/Stats Imports ------
from matplotlib import pyplot as plt # 从matplotlib导入pyplot，用于绘图
from astropy.io import fits # 从Astropy导入fits模块，用于处理FITS文件
from utils import get_jwst_file, plot_dark_data, plot_cleaned_data, plot_spectra # 导入自定义工具函数
# Hide all log and warning messages. # 隐藏所有日志和警告消息
logging.disable(logging.ERROR) # 禁用错误级别的日志
warnings.simplefilter("ignore", RuntimeWarning) # 忽略运行时警告
3. 下载数据 #
本笔记本的输入数据包含通过NIRSpec BOTS模式使用G395H光栅观测的WASP-39b的过境。该数据集是JWST早期发布科学计划ERS-1366的一部分，具体为观测3。它由456次积分组成，每次积分有70组，使用SUB2048子阵列（每次积分每组2048 x 32像素的阵列大小）。整个观测持续了10.56小时，并分为三个部分。本笔记本将以第二部分（seg002）为例。然而，需要注意的是，在进入TSO3Pipeline之前，所有部分必须首先通过Spec2Pipeline进行处理。
有关当前和计划中的JWST时间序列观测（TSO）程序的全面概述，特别是涉及使用NIRSpec BOTS模式的过境系外行星的程序，请参考TrExoLiSTS。
# 定义下载目录
mast_products_dir = "./mast_products/"
# 检查目录是否存在
if not os.path.exists(mast_products_dir):
# 如果目录不存在，则创建该目录
os.makedirs(mast_products_dir)
# 该笔记本专注于第二个分段。
obs_ids = ["jw01366003001_04101_00001-seg002"] # 观测ID列表
detectors = [1, 2] # 包含探测器NRS1和NRS2的列表
# 指定计数率产品的名称。
rateints_names = [] # 存储计数率产品名称的列表
for obs_id in obs_ids: # 遍历每个观测ID
for detector in detectors: # 遍历每个探测器
rateints_names.append(f"{obs_id}_nrs{detector}_rateints.fits") # 生成计数率文件名并添加到列表中
# 下载所有的FITS文件。
for name in rateints_names: # 遍历每个文件名
print(f"Downloading {name}") # 打印正在下载的文件名
get_jwst_file( # 调用函数下载文件
name, # 文件名
mast_api_token=None, # API令牌，默认为None
save_directory=mast_products_dir, # 保存目录
)
4. 运行 Spec2Pipeline 而不进行 NSClean（原始数据） #
下面的单元格执行 Spec2Pipeline，在处理过程中明确跳过 NSClean 步骤。生成的二级产品将作为参考点，展示在未去除 1/f 噪声的情况下，计数率图像和最终提取光谱的外观。
# 设置用于运行管道的目录，跳过NSClean处理
stage2_nsclean_skipped_dir = "./stage2_nsclean_skipped/"
# 检查目录是否存在
if not os.path.exists(stage2_nsclean_skipped_dir):
os.makedirs(stage2_nsclean_skipped_dir) # 如果目录不存在，则创建该目录
4.1 修改 EXTRACT1D 参考文件 #
我们调整 extract_1d 步骤的默认参数参考文件，特别是修改提取宽度孔径。这一修改是必要的，因为当前的处理流程并未对 2D 光谱进行重采样以校正光谱（去除曲率）。因此，我们依赖于较宽的提取孔径，尤其是对于 NRS2。
# EXTRACT1D 参考文件，用于 BOTS 数据 NRS1。
refs = api.dump_references(crds_client.get_context_used('jwst'), # 获取JWST上下文的参考文件
['jwst_nirspec_extract1d_0006.json']) # 指定需要的参考文件
extract_1d_ref = refs['jwst_nirspec_extract1d_0006.json'] # 获取提取的参考文件路径
# 以读取模式打开 EXTRACT1D 参考文件。
with open(extract_1d_ref, "r") as ref_file: # 打开参考文件
params = json.load(ref_file) # 读取 JSON 数据
params["apertures"][0][
"extract_width"
] = 27 # 提取孔径半径（单位：像素），可以修改。
params["apertures"][0]["xstart"] = 100 # 设置起始 x 坐标
params["apertures"][0].pop("nod2_offset") # 移除 nod2_offset。
params["apertures"][0].pop("nod3_offset") # 移除 nod3_offset。
params["apertures"][0].pop("nod5_offset") # 移除 nod5_offset。
params["apertures"][1][
"extract_width"
] = 27 # 提取孔径半径（单位：像素），可以修改。
params["apertures"][1]["xstart"] = 100 # 设置起始 x 坐标
params["apertures"][1].pop("nod2_offset") # 移除 nod2_offset。
params["apertures"][1].pop("nod3_offset") # 移除 nod3_offset。
params["apertures"][1].pop("nod5_offset") # 移除 nod5_offset。
# 将更改写入新文件。
newData = json.dumps(params, indent=4) # 将参数转换为格式化的 JSON 字符串
# 添加后缀 '_bots' 以区分文件与默认版本。
with open(os.path.basename(extract_1d_ref)[:-5] + "_nrs1_bots.json", "w") as file: # 创建新文件
file.write(newData) # 写入新数据
# EXTRACT1D参考文件用于BOTS数据NRS2。
refs = api.dump_references(crds_client.get_context_used('jwst'), # 获取JWST上下文使用的参考文件
['jwst_nirspec_extract1d_0006.json']) # 指定要获取的参考文件名
extract_1d_ref = refs['jwst_nirspec_extract1d_0006.json'] # 提取指定的EXTRACT1D参考文件
# 以读模式打开EXTRACT1D参考文件。
with open(extract_1d_ref, "r") as ref_file: # 打开文件进行读取
params = json.load(ref_file) # 将文件内容加载为JSON格式的字典
params["apertures"][0][ # 访问第一个光圈的参数
"extract_width"
] = 27 # 提取光圈半径，单位为像素，可以修改。
params["apertures"][0].pop("nod2_offset") # 移除nod2_offset参数。
params["apertures"][0].pop("nod3_offset") # 移除nod3_offset参数。
params["apertures"][0].pop("nod5_offset") # 移除nod5_offset参数。
params["apertures"][1][ # 访问第二个光圈的参数
"extract_width"
] = 27 # 提取光圈半径，单位为像素，可以修改。
params["apertures"][1].pop("nod2_offset") # 移除nod2_offset参数。
params["apertures"][1].pop("nod3_offset") # 移除nod3_offset参数。
params["apertures"][1].pop("nod5_offset") # 移除nod5_offset参数。
# 将更改写入新文件。
newData = json.dumps(params, indent=4) # 将修改后的参数转换为JSON格式字符串，并格式化为4个空格缩进
# 添加后缀'_bots'以区分该文件与默认版本。
with open(os.path.basename(extract_1d_ref)[:-5] + "_nrs2_bots.json", "w") as file: # 创建新文件并写入
file.write(newData) # 将新数据写入文件
# 原始数据（未应用NSClean）。
# 估计运行时间：10分钟（可能会有所不同）。
start = tt.time() # 记录开始时间
for i in rateints_names: # 遍历每个速率积分文件名
print(f"Processing {i}...") # 打印正在处理的文件名
if "nrs1" in i: # 检查文件名中是否包含"nrs1"
extract1dref = os.path.basename(extract_1d_ref)[:-5] + "_nrs1_bots.json" # 设置提取1D参考文件名
else: # 如果不是"nrs1"
extract1dref = os.path.basename(extract_1d_ref)[:-5] + "_nrs2_bots.json" # 设置提取1D参考文件名
Spec2Pipeline.call( # 调用Spec2Pipeline处理数据
mast_products_dir + i, # 输入文件路径
save_results=True, # 保存结果
steps={ # 处理步骤配置
# 为光谱数据分配波长解决方案。
"assign_wcs": {"skip": False}, # 不跳过此步骤
# 从NIRSpec图像中去除相关读噪声（1/f噪声）。
"nsclean": {"skip": True}, # 跳过此步骤
"extract_2d": {"skip": False}, # 获取2D切片。
"srctype": {"skip": False}, # 分配源类型：'point'或'extended'。
# 我们选择跳过此步骤，因为它可能引入虚假像素
# 这会影响光曲线的散射。
"flat_field": {"skip": True}, # 跳过此步骤
# 我们也跳过此步骤，因为BOTS系外行星观测是相对的。
"photom": {"skip": True}, # 跳过此步骤
# 估计在2D光谱中被标记为DO_NOT_USE的像素的通量值。
"pixel_replace": { # 像素替换配置
"skip": False, # 不跳过此步骤
"n_adjacent_cols": 5, # 邻近列数
"algorithm": "fit_profile", # 使用的算法
},
"extract_1d": { # 1D提取配置
"skip": False, # 不跳过此步骤
"override_extract1d": extract1dref, # 覆盖提取1D参考文件
"apply_apcorr": False, # 不应用光圈校正
"subtract_background": False, # 不减去背景
},
# 'dispaxis': 1, 'xstart': 100, 'xstop': 150, 'extract_width': 10},
},
output_dir=stage2_nsclean_skipped_dir, # 输出目录
)
print(f"Saved {i[:-13]}" + "calints.fits") # 打印保存的calints文件名
print(f"Saved {i[:-13]}" + "x1dints.fits") # 打印保存的x1dints文件名
end = tt.time() # 记录结束时间
print("Run time: ", round(end - start, 1) / 60.0, " min") # 打印运行时间
5. 使用 NSClean 清理 1/f 噪声（默认管道掩模） #
如果在 Spec2Pipeline 的 NSClean 步骤 中未提供用户自定义的掩模文件，管道将根据默认参数生成一个掩模。该掩模将识别任何未被照亮的像素。也就是说，掩模必须包含 True 和 False 值，其中 True 表示该像素是黑暗的，False 表示该像素是照亮的（不是黑暗的）。
默认情况下，管道将以下探测器区域标记为照亮的非黑暗区域（False）：
被指定为 FS 数据的科学区域的像素。
5-σ 异常值（默认值）
在速率数据中设置为 NaN 的任何像素。
要调整掩模中的异常值检测，可以尝试修改 n_sigma 参数（将在下一节中探讨）。较高的值将识别较少的异常值，而较低的值将识别更多的异常值。
默认生成的掩模将在下面保存和分析。
# 设置用于运行NSClean的目录，使用默认参数。
stage2_nsclean_default_dir = "./stage2_nsclean_default/" # 定义默认目录路径
if not os.path.exists(stage2_nsclean_default_dir): # 检查目录是否存在
os.makedirs(stage2_nsclean_default_dir) # 如果目录不存在，则创建该目录。
# 1/f 噪声清理数据（默认 NSClean 管道掩膜）。
# 估计运行时间：36 分钟（可能会有所不同）。
start = tt.time() # 记录开始时间
for i in rateints_names: # 遍历每个速率积分名称
print(f"Processing {i}...") # 打印正在处理的文件名
if "nrs1" in i: # 如果文件名中包含 "nrs1"
extract1dref = os.path.basename(extract_1d_ref)[:-5] + "_nrs1_bots.json" # 设置提取1D参考文件名
else: # 否则
extract1dref = os.path.basename(extract_1d_ref)[:-5] + "_nrs2_bots.json" # 设置提取1D参考文件名
Spec2Pipeline.call( # 调用 Spec2Pipeline 处理数据
mast_products_dir + i, # 输入文件路径
save_results=True, # 保存结果
steps={ # 定义处理步骤
# 为光谱数据分配波长解决方案。
"assign_wcs": {"skip": False, "save_results": True}, # 赋值 WCS
# 从 NIRSpec 图像中去除相关读噪声（1/f 噪声）。
"nsclean": {"skip": False, "save_mask": True, "save_results": True}, # 噪声清理
"extract_2d": {"skip": False}, # 获取 2D 剪切图。
"srctype": {"skip": False}, # 分配源类型：'point' 或 'extended'。
# 我们选择跳过此步骤，因为它可能会引入虚假像素
# 这会影响光曲线的散射。
"flat_field": {"skip": True}, # 跳过平场校正
# 我们也跳过此步骤，因为 BOTS 系外行星观测是相对的。
"photom": {"skip": True}, # 跳过光度校正
# 估计在 2D 光谱中标记为 DO_NOT_USE 的像素的通量值。
"pixel_replace": { # 像素替换步骤
"skip": False, # 不跳过
"n_adjacent_cols": 5, # 邻近列数
"algorithm": "fit_profile", # 使用的算法
},
"extract_1d": { # 提取 1D 数据
"skip": False, # 不跳过
"override_extract1d": extract1dref, # 覆盖提取1D参考
"apply_apcorr": False, # 不应用光圈校正
"subtract_background": False, # 不减去背景
},
# 'dispaxis': 1, 'xstart': 100, 'xstop': 150, 'extract_width': 10},
},
output_dir=stage2_nsclean_default_dir, # 输出目录
)
print(f"Saved {i[:-13]}" + "assign_wcs.fits") # 打印保存的文件名
print(f"Saved {i[:-13]}" + "mask.fits") # 打印保存的掩膜文件名
print(f"Saved {i[:-13]}" + "nsclean.fits") # 打印保存的噪声清理文件名
print(f"Saved {i[:-13]}" + "calints.fits") # 打印保存的校准积分文件名
print(f"Saved {i[:-13]}" + "x1dints.fits") # 打印保存的1D积分文件名
end = tt.time() # 记录结束时间
print("Run time: ", round(end - start, 1) / 60.0, " min") # 打印运行时间
警告：
在某些情况下，NSClean步骤可能无法找到背景噪声的拟合。如果掩模中没有足够的暗数据（标记为True），则可能会发生这种失败。特别是，掩模中的每一列（除了前4列和最后4列）必须包含一些标记为True的像素。背景拟合过程会逐列考虑，因此如果某一列没有数据可供拟合，则会崩溃。如果发生失败，请检查下面图像中的掩模，确保每一列至少有一些True值。
5.1 验证掩模（默认管道掩模） #
检查掩模与速率数据，以确保它仅保留探测器的暗区。
请注意，仍然存在一些剩余的照明区域，主要是由于瞬态伪影，如宇宙射线和雪球。
此外，对于某些数据集，探测器上可能存在多个未被WCS边界框掩盖的照明区域。在此数据中，请注意掩模未覆盖光谱轨迹的蓝端区域。该区域是被照明的，但尚未校准。
# 绘制带有屏蔽区域的速率数据。
# 从管道中构建的即时屏蔽列表。
nsclean_default_masks = [
stage2_nsclean_default_dir + "jw01366003001_04101_00001-seg002_nrs1_mask.fits", # NRS1的屏蔽文件路径
stage2_nsclean_default_dir + "jw01366003001_04101_00001-seg002_nrs2_mask.fits", # NRS2的屏蔽文件路径
]
# 绘制每组相关的速率数据和屏蔽文件。
for rateint_file, mask_file in zip(rateints_names, nsclean_default_masks):
plot_dark_data(
mast_products_dir + rateint_file, # 速率数据文件的完整路径
mask_file, # 对应的屏蔽文件
slice_index=100, # 切片索引
axis=0, # 绘图的轴
layout="rows", # 布局方式
)
默认的管道掩模（pipeline mask）对于NRS2掩模会遮盖整个计数率图像（countrate image）。
5.2 比较原始数据与清理后的数据（默认管道掩模） #
现在我们可以将清理后的数据（使用默认管道掩模）与原始速率文件进行比较，并验证1/f噪声是否已减少。
在许多情况下，清理过程会在速率文件中引入新的伪影。这些伪影应仔细检查，并与噪声减少的好处进行权衡。如果瞬态伪影（如雪球）干扰了清理过程，手动编辑掩模以将这些区域排除在背景拟合之外可能是有益的。为此，可以尝试调整异常值检测阈值或直接编辑掩模数组中的特定像素（将在接下来的几个部分中探讨）。否则，请参考NSClean文档以获取有关手动编辑的其他建议。
请注意，在下面的图像中，有一些值与原始速率文件存在较大的相对差异（在下面的相对差异图像中显示）。这些是清理过程的伪影。
此外，还有更广泛的低水平残余背景效应（在右侧的相对差异图像中显示，散布的异常值通过σ裁剪识别，并被掩模隐藏）。这些包括我们试图去除的背景模式：色散方向上的1/f噪声变化以及全帧数据中框架顶部和底部的画框效应。然而，清理过程中可能还会因过拟合暗数据而引入低水平伪影。
仔细检查这两幅残余图像，以了解清理过程对您数据的影响。
# 绘制原始数据和清理后的数据，以及残差图。
# 定义清理后的默认掩膜文件路径列表
cleaned_default_masks = [
stage2_nsclean_default_dir + "jw01366003001_04101_00001-seg002_nrs1_nsclean.fits", # NRS1清理后的文件路径
stage2_nsclean_default_dir + "jw01366003001_04101_00001-seg002_nrs2_nsclean.fits", # NRS2清理后的文件路径
]
# 遍历每个关联的rateint数据文件和清理后的文件
for rateint_file, cleaned_file in zip(rateints_names, cleaned_default_masks):
# 调用函数绘制清理后的数据
plot_cleaned_data(
mast_products_dir + rateint_file, # 原始数据文件路径
cleaned_file, # 清理后的数据文件路径
slice_index=100, # 切片索引
layout="rows" # 布局方式
)
默认的NRS2管道掩模会遮盖整个计数率图像，因此在应用NSClean时没有任何差异。
比较从清理后的数据中提取的光谱与从原始计数率文件中提取的光谱。
# 绘制整个光谱。
# 1D 提取的光谱。
x1d_nsclean_skipped = [
# 路径到未清理的光谱数据文件（NRS1）
stage2_nsclean_skipped_dir + "jw01366003001_04101_00001-seg002_nrs1_x1dints.fits",
# 路径到未清理的光谱数据文件（NRS2）
stage2_nsclean_skipped_dir + "jw01366003001_04101_00001-seg002_nrs2_x1dints.fits",
]
x1d_nsclean_default = [
# 路径到默认清理的光谱数据文件（NRS1）
stage2_nsclean_default_dir + "jw01366003001_04101_00001-seg002_nrs1_x1dints.fits",
# 路径到默认清理的光谱数据文件（NRS2）
stage2_nsclean_default_dir + "jw01366003001_04101_00001-seg002_nrs2_x1dints.fits",
]
# 遍历未清理和清理的光谱文件，绘制光谱
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_default):
plot_spectra([original, cleaned], ext_num=100) # 绘制光谱，扩展编号为100
# 感兴趣的波长/通量区域
wavelength_range = {"nrs1": [3.4, 3.45], "nrs2": [4.25, 4.275]} # 定义波长范围
flux_range = {"nrs1": [1100, 1300], "nrs2": [500, 550]} # 定义通量范围
# 再次遍历未清理和清理的光谱文件，绘制指定波长和通量范围的光谱
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_default):
plot_spectra(
[original, cleaned], # 光谱文件列表
ext_num=100, # 扩展编号
scale_percent=0, # 缩放百分比
wavelength_range=wavelength_range, # 波长范围
flux_range=flux_range, # 通量范围
)
注意：
NRS1中的亮谱在清理后几乎没有变化。
在NRS2中，管道动态生成的默认掩模遮盖了整个科学区域，这可能会过度遮掩那些能够增强背景拟合的暗区。这导致NRS2的原始数据与清理后的数据之间没有明显差异。在接下来的部分中，我们将探讨使用管道生成替代掩模的过程，以解决此问题。
6. 使用 NSClean（替代掩膜）清理 1/f 噪声 #
对于某些数据集，掩膜整个科学区域可能会过度掩盖探测器的暗区，这些区域可以用来改善背景拟合。过度掩膜可能会在清理过程中引入一些高频噪声，表现为光谱轨迹上的垂直条纹。
此外，对于某些数据集，探测器上可能存在几个未被 WCS 边界框掩膜的照明区域。在子数组数据中，请注意掩膜并未覆盖光谱轨迹的蓝端。该区域是被照明的，但尚未校准。
在某些情况下，使用替代算法构建掩膜可能会更有利。在这里，我们不使用边界框，而是迭代地掩膜任何超过背景 0.1 sigma 的数据。对于亮源，这样可以在光谱轨迹附近保留更多的暗数据，从而可能改善背景拟合。
然而，请注意，过度清理可能会影响光谱的连续谱水平，如果掩膜中包含的照明数据过多或过少。再次强调，生成的掩膜和输出光谱应仔细检查，以权衡清理的好处与对光谱的影响。
要调整此掩膜中的照明检测，请尝试修改下面的 n_sigma 参数。较高的值将识别较少的照明，较低的值将识别更多的照明。
# 设置用于运行NSClean的目录，使用备用参数。
stage2_nsclean_alternate_dir = "./stage2_nsclean_alternate/"
# 检查目录是否存在
if not os.path.exists(stage2_nsclean_alternate_dir):
# 如果目录不存在，则创建该目录
os.makedirs(stage2_nsclean_alternate_dir)
# 1/f 噪声清理数据（备用 NSClean 管道掩膜）。
# 估计运行时间：66分钟（可能会有所不同）。
start = tt.time() # 记录开始时间
for indx, i in enumerate(rateints_names): # 遍历 rateints_names 列表中的每个元素
print(f"Processing {i}...") # 打印正在处理的文件名
if "nrs1" in i: # 检查文件名中是否包含 "nrs1"
extract1dref = os.path.basename(extract_1d_ref)[:-5] + "_nrs1_bots.json" # 生成 nrs1 的提取文件名
else:
extract1dref = os.path.basename(extract_1d_ref)[:-5] + "_nrs2_bots.json" # 生成 nrs2 的提取文件名
Spec2Pipeline.call( # 调用 Spec2Pipeline 处理数据
mast_products_dir + i, # 输入文件路径
save_results=True, # 保存结果
steps={ # 定义处理步骤
# 为光谱数据分配波长解决方案。
"assign_wcs": {"skip": False, "save_results": True}, # 步骤：分配 WCS
# 从 NIRSpec 图像中去除相关读噪声（1/f 噪声）。
"nsclean": { # 步骤：噪声清理
"skip": False, # 不跳过此步骤
"save_mask": True, # 保存掩膜
"n_sigma": 0.1, # 噪声阈值
"mask_spectral_regions": False, # 不掩盖光谱区域
"save_results": True, # 保存结果
},
"extract_2d": {"skip": False}, # 获取 2D 切片。
"srctype": {"skip": False}, # 分配源类型：'point' 或 'extended'。
# 我们选择跳过此步骤，因为它可能引入虚假像素
# 这会影响光曲线的散射。
"flat_field": {"skip": True}, # 跳过平场校正步骤
# 我们也跳过此步骤，因为 BOTS 系外行星观测是相对的。
"photom": {"skip": True}, # 跳过光度校正步骤
# 估计在 2D 光谱中标记为 DO_NOT_USE 的像素的通量值。
"pixel_replace": { # 步骤：像素替换
"skip": False, # 不跳过此步骤
"n_adjacent_cols": 5, # 邻近列数
"algorithm": "fit_profile", # 使用的算法
},
"extract_1d": { # 步骤：提取 1D 数据
"skip": False, # 不跳过此步骤
"override_extract1d": extract1dref, # 覆盖提取文件名
"apply_apcorr": False, # 不应用光圈校正
"subtract_background": False, # 不减去背景
},
# 'dispaxis': 1, 'xstart': 100, 'xstop': 150, 'extract_width': 10},
},
output_dir=stage2_nsclean_alternate_dir, # 输出目录
)
print(f"Saved {i[:-13]}" + "assign_wcs.fits") # 打印保存的文件信息
print(f"Saved {i[:-13]}" + "mask.fits") # 打印保存的掩膜文件信息
print(f"Saved {i[:-13]}" + "nsclean.fits") # 打印保存的噪声清理文件信息
print(f"Saved {i[:-13]}" + "calints.fits") # 打印保存的校准文件信息
print(f"Saved {i[:-13]}" + "x1dints.fits") # 打印保存的 1D 文件信息
end = tt.time() # 记录结束时间
print("Run time: ", round(end - start, 1) / 60.0, " min") # 打印运行时间
6.1 验证掩模（备用掩模） #
检查掩模与速率数据，以确保它仅保留探测器的暗区。
# 绘制带有屏蔽区域的速率数据。
# 从管道中构建的即时掩膜列表。
nsclean_alternate_masks = [
stage2_nsclean_alternate_dir + "jw01366003001_04101_00001-seg002_nrs1_mask.fits", # NRS1掩膜文件路径
stage2_nsclean_alternate_dir + "jw01366003001_04101_00001-seg002_nrs2_mask.fits", # NRS2掩膜文件路径
]
# 绘制每个相关的速率数据集和掩膜文件。
for rateint_file, mask_file in zip(rateints_names, nsclean_alternate_masks): # 遍历速率数据文件和掩膜文件
plot_dark_data( # 调用绘图函数
mast_products_dir + rateint_file, # 速率数据文件路径
mask_file, # 掩膜文件路径
slice_index=100, # 切片索引
axis=0, # 轴方向
scale=9 # 缩放因子
)
6.2 原始数据与清理数据的比较（备用掩膜） #
# 绘制原始数据和清理后的数据，以及残差图。
# 定义清理后的备用掩膜文件路径列表
cleaned_alternate_masks = [
stage2_nsclean_alternate_dir + "jw01366003001_04101_00001-seg002_nrs1_nsclean.fits", # NRS1清理后的文件路径
stage2_nsclean_alternate_dir + "jw01366003001_04101_00001-seg002_nrs2_nsclean.fits", # NRS2清理后的文件路径
]
# 遍历每个关联的rateint数据和清理后的文件
for rateint_file, cleaned_file in zip(rateints_names, cleaned_alternate_masks):
# 绘制清理后的数据
plot_cleaned_data(
mast_products_dir + rateint_file, # 原始rateint文件路径
cleaned_file, # 清理后的文件路径
slice_index=100, # 切片索引
layout="rows" # 布局方式
)
比较从清理后的数据中提取的光谱与从原始速率文件中提取的光谱。
# 1D 提取的光谱
# 定义未清理的光谱文件路径列表
x1d_nsclean_skipped = [
# 第一段光谱文件路径
stage2_nsclean_skipped_dir + "jw01366003001_04101_00001-seg002_nrs1_x1dints.fits",
# 第二段光谱文件路径
stage2_nsclean_skipped_dir + "jw01366003001_04101_00001-seg002_nrs2_x1dints.fits",
]
# 定义替代的清理光谱文件路径列表
x1d_nsclean_alternate = [
# 第一段替代光谱文件路径
stage2_nsclean_alternate_dir + "jw01366003001_04101_00001-seg002_nrs1_x1dints.fits",
# 第二段替代光谱文件路径
stage2_nsclean_alternate_dir + "jw01366003001_04101_00001-seg002_nrs2_x1dints.fits",
]
# 遍历未清理和清理的光谱文件路径
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_alternate):
# 绘制光谱，传入原始和清理后的光谱文件，扩展编号为100，比例设置为0
plot_spectra([original, cleaned], ext_num=100, scale_percent=0)
# 感兴趣的波长范围
wavelength_range = {"nrs1": [3.4, 3.45], "nrs2": [4.25, 4.275]} # 定义不同通道的波长范围
flux_range = {"nrs1": [1100, 1300], "nrs2": [500, 550]} # 定义不同通道的通量范围
# 遍历原始数据和清理后的数据
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_alternate):
# 绘制光谱
plot_spectra(
[original, cleaned], # 传入原始和清理后的光谱数据
ext_num=100, # 扩展编号
scale_percent=0, # 缩放百分比
wavelength_range=wavelength_range, # 波长范围
flux_range=flux_range, # 通量范围
)
注意： NRS1/NRS2中的亮谱经过清理后几乎没有差异。然而，根据比率，NRS1中的谱线表明清理后的数据可能具有稍高的通量，而在NRS2中，清理后的数据与原始数据相比则显示出稍低的通量。
7. 使用 NSClean 清理 1/f 噪声（手动修改的掩膜） #
在某些情况下，可能需要手动生成掩膜。在这里，我们介绍 一种 从管道输出的默认掩膜开始手动修改掩膜的方法。值得注意的是，使用此方法修改的掩膜不一定会优于之前的两种选项。
# 设置用于运行 NSClean 的目录，并使用用户提供的掩模。
stage2_nsclean_modified_dir = "./stage2_nsclean_modified/" # 定义修改后的 NSClean 目录路径
if not os.path.exists(stage2_nsclean_modified_dir): # 检查目录是否存在
# 如果目录不存在，则创建该目录。
os.makedirs(stage2_nsclean_modified_dir) # 创建目录
# 手动修改某些掩膜区域。
# 定义一个列表来存储修改后的掩膜路径。
nsclean_modified_masks = []
# 遍历原始掩膜列表。
for mask in nsclean_default_masks:
# 新的掩膜文件名。
output_file = os.path.basename(mask)[:-5] + "_modified.fits"
# 打开FITS文件。
with fits.open(mask) as hdul:
# 从科学扩展中提取掩膜数据。
mask_data = hdul["SCI"].data.copy() # 复制数据。
# 步骤1：将默认掩膜区域设置回True。
mask_data[:, :, :] = True
# 步骤2：手动重新定义掩膜区域。
if "nrs1" in mask: # 如果掩膜名称中包含"nrs1"
mask_data[:, 11:29, 400:700] = False # 设置特定区域为False
mask_data[:, 10:27, 700:800] = False # 设置特定区域为False
mask_data[:, 8:26, 800:1000] = False # 设置特定区域为False
mask_data[:, 7:24, 1000:1200] = False # 设置特定区域为False
mask_data[:, 4:23, 1200:1400] = False # 设置特定区域为False
mask_data[:, 4:22, 1400:1600] = False # 设置特定区域为False
mask_data[:, 3:21, 1600:2043] = False # 设置特定区域为False
else: # 如果掩膜名称不包含"nrs1"
mask_data[:, :16, :290] = False # 设置特定区域为False
mask_data[:, 1:18, 290:740] = False # 设置特定区域为False
mask_data[:, 2:22, 740:1290] = False # 设置特定区域为False
mask_data[:, 8:29, 1290:1840] = False # 设置特定区域为False
mask_data[:, 11:, 1840:] = False # 设置特定区域为False
# 更新科学扩展中的数据。
hdul["SCI"].data = mask_data
# 保存修改后的FITS文件。
output_path = os.path.join(stage2_nsclean_modified_dir, output_file)
hdul_modified = hdul.copy() # 复制文件
hdul_modified.writeto(output_path, overwrite=True) # 写入修改后的文件，覆盖同名文件
nsclean_modified_masks.append(output_path) # 将修改后的路径添加到列表中
print(f"Saved modified mask as: {output_path}") # 打印保存的文件路径
7.1 验证掩模（手动修改的掩模） #
检查掩模与速率数据，以确保它仅保留探测器的暗区。
# 绘制带有屏蔽区域的速率数据。
# 修改后的掩膜列表，用于数据处理管道。
nsclean_modified_masks = [
stage2_nsclean_modified_dir # 指定修改后的掩膜目录
+ "jw01366003001_04101_00001-seg002_nrs1_mask_modified.fits", # NRS1的修改掩膜文件
stage2_nsclean_modified_dir # 指定修改后的掩膜目录
+ "jw01366003001_04101_00001-seg002_nrs2_mask_modified.fits", # NRS2的修改掩膜文件
]
# 遍历每个关联的速率数据和掩膜文件进行绘图。
for rateint_file, mask_file in zip(rateints_names, nsclean_modified_masks):
plot_dark_data( # 调用绘图函数
mast_products_dir + rateint_file, # 速率数据文件路径
mask_file, # 掩膜文件路径
slice_index=100, # 切片索引
axis=0, # 绘图轴
scale=9 # 缩放因子
)
# 1/f 噪声清理数据（用户提供的掩码）。
# 估计运行时间：57分钟（可能会有所不同）。
start = tt.time() # 记录开始时间
for indx, i in enumerate(rateints_names): # 遍历每个速率积分名称
print(f"Processing {i}...") # 打印正在处理的文件名
if "nrs1" in i: # 检查文件名中是否包含"nrs1"
extract1dref = os.path.basename(extract_1d_ref)[:-5] + "_nrs1_bots.json" # 设置nrs1的提取文件名
else:
extract1dref = os.path.basename(extract_1d_ref)[:-5] + "_nrs2_bots.json" # 设置nrs2的提取文件名
Spec2Pipeline.call( # 调用Spec2Pipeline进行数据处理
mast_products_dir + i, # 输入文件路径
save_results=True, # 保存结果
steps={ # 定义处理步骤
# 为光谱数据分配波长解决方案。
"assign_wcs": {"skip": False, "save_results": True}, # 赋值WCS
# 从NIRSpec图像中去除相关读噪声（1/f噪声）。
"nsclean": {
"skip": False, # 不跳过此步骤
"save_mask": True, # 保存掩码
"save_results": True, # 保存结果
"user_mask": nsclean_modified_masks[indx], # 使用用户提供的掩码
},
"extract_2d": {"skip": False}, # 获取2D切片
"srctype": {"skip": False}, # 分配源类型：'point'或'extended'
# 我们选择跳过此步骤，因为它可能引入虚假像素
# 这会影响光曲线的散射。
"flat_field": {"skip": True}, # 跳过平场校正步骤
# 我们也跳过此步骤，因为BOTS系外行星观测是相对的。
"photom": {"skip": True}, # 跳过光度校正步骤
# 估计在2D光谱中标记为DO_NOT_USE的像素的通量值。
"pixel_replace": {
"skip": False, # 不跳过此步骤
"n_adjacent_cols": 5, # 邻近列数
"algorithm": "fit_profile", # 使用拟合轮廓算法
},
"extract_1d": { # 进行一维提取
"skip": False, # 不跳过此步骤
"override_extract1d": extract1dref, # 覆盖提取文件名
"apply_apcorr": False, # 不应用光圈校正
"subtract_background": False, # 不减去背景
},
# 'dispaxis': 1, 'xstart': 100, 'xstop': 150, 'extract_width': 10},
},
output_dir=stage2_nsclean_modified_dir, # 输出目录
)
print(f"Saved {i[:-13]}" + "assign_wcs.fits") # 打印保存的assign_wcs文件名
print(f"Saved {i[:-13]}" + "mask.fits") # 打印保存的mask文件名
print(f"Saved {i[:-13]}" + "nsclean.fits") # 打印保存的nsclean文件名
print(f"Saved {i[:-13]}" + "calints.fits") # 打印保存的calints文件名
print(f"Saved {i[:-13]}" + "x1dints.fits") # 打印保存的x1dints文件名
end = tt.time() # 记录结束时间
print("Run time: ", round(end - start, 1) / 60.0, " min") # 打印运行时间
7.2 原始数据与清理数据（手动修改的掩模）比较 #
# 绘制原始数据和清理后的数据，以及残差图。
# 定义清理后的修改掩膜文件路径列表
cleaned_modified_masks = [
stage2_nsclean_modified_dir + "jw01366003001_04101_00001-seg002_nrs1_nsclean.fits", # NRS1清理后的文件
stage2_nsclean_modified_dir + "jw01366003001_04101_00001-seg002_nrs2_nsclean.fits", # NRS2清理后的文件
]
# 遍历每一组速率积分数据和清理后的文件
for rateint_file, cleaned_file in zip(rateints_names, cleaned_modified_masks):
# 绘制清理后的数据
plot_cleaned_data(
mast_products_dir + rateint_file, # 原始速率积分文件路径
cleaned_file, # 清理后的文件路径
slice_index=100, # 切片索引
layout="rows" # 布局方式
)
比较从清理后的数据中提取的光谱与从原始速率文件中提取的光谱。
# 1D 提取的光谱。
# 定义未清理的光谱文件路径列表
x1d_nsclean_skipped = [
# 添加第一段未清理的光谱文件路径
stage2_nsclean_skipped_dir + "jw01366003001_04101_00001-seg002_nrs1_x1dints.fits",
# 添加第二段未清理的光谱文件路径
stage2_nsclean_skipped_dir + "jw01366003001_04101_00001-seg002_nrs2_x1dints.fits",
]
# 定义已清理的光谱文件路径列表
x1d_nsclean_modified = [
# 添加第一段已清理的光谱文件路径
stage2_nsclean_modified_dir + "jw01366003001_04101_00001-seg002_nrs1_x1dints.fits",
# 添加第二段已清理的光谱文件路径
stage2_nsclean_modified_dir + "jw01366003001_04101_00001-seg002_nrs2_x1dints.fits",
]
# 遍历未清理和已清理的光谱文件路径
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_modified):
# 绘制光谱，传入原始和清理后的文件路径，扩展编号为100，比例设置为0
plot_spectra([original, cleaned], ext_num=100, scale_percent=0)
# 感兴趣的波长范围
wavelength_range = {"nrs1": [3.4, 3.45], "nrs2": [4.25, 4.275]} # 定义不同通道的波长范围
flux_range = {"nrs1": [1100, 1300], "nrs2": [500, 550]} # 定义不同通道的通量范围
# 遍历原始数据和清理后的数据
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_modified):
plot_spectra( # 调用绘制光谱的函数
[original, cleaned], # 传入原始和清理后的光谱数据
ext_num=100, # 设置扩展编号
scale_percent=0, # 设置缩放百分比
wavelength_range=wavelength_range, # 传入波长范围
flux_range=flux_range, # 传入通量范围
)
注意： 再次强调，NRS1/NRS2中的亮谱在清理后几乎没有差异。然而，根据比率，NRS1中的光谱表明清理后的数据可能具有稍高的通量，而在NRS2中，清理后的数据与原始数据相比则显示出稍低的通量。
8. 结论 #
下面的最终图表展示了计数率图像和相应的1D提取光谱并排比较不同的清理方法：原始图像（未应用NSClean）、清理后的计数率图像（使用默认的管道掩模）、清理后的计数率图像（使用替代管道掩模），以及最后，清理后的计数率图像（使用手动修改的掩模）。
请注意，本笔记本中呈现的结果可能会因不同数据集（例如，不同亮度、空间范围等的目标）而有所不同。鼓励用户使用不同的掩模方法探索NSClean，以确定最佳结果。
清理算法的输出现在已准备好进行进一步处理。上述Spec2Pipeline运行生成的(_calint.fits)文件可以作为TSO3Pipeline的输入，用于生成最终的组合光谱。
# 未清理数据 vs. 清理数据（默认掩膜） vs. 清理数据（备用掩膜）速率数据
# 读取原始速率数据
original_rate_data = [
fits.open(mast_products_dir + rate_name)[1].data[100, :, :] # 从文件中提取第100层数据
for rate_name in rateints_names # 遍历速率名称列表
]
# 读取使用默认掩膜清理后的速率数据
cleaned_rate_default_data = [
fits.open(cleaned_default_mask)[1].data[100, :, :] # 从文件中提取第100层数据
for cleaned_default_mask in cleaned_default_masks # 遍历默认掩膜列表
]
# 读取使用备用掩膜清理后的速率数据
cleaned_rate_alternate_data = [
fits.open(cleaned_alternate_mask)[1].data[100, :, :] # 从文件中提取第100层数据
for cleaned_alternate_mask in cleaned_alternate_masks # 遍历备用掩膜列表
]
# 读取使用手动修改掩膜清理后的速率数据
cleaned_rate_modified_data = [
fits.open(cleaned_modified_mask)[1].data[100, :, :] # 从文件中提取第100层数据
for cleaned_modified_mask in cleaned_modified_masks # 遍历手动修改掩膜列表
]
# 为可视化绘图处理数据
for data_list in [
original_rate_data, # 原始速率数据
cleaned_rate_default_data, # 清理后的默认掩膜数据
cleaned_rate_alternate_data, # 清理后的备用掩膜数据
cleaned_rate_modified_data, # 清理后的手动修改掩膜数据
]:
for data in data_list:
data[np.isnan(data)] = 0 # 将NaN值替换为0
# 原始数据与清理数据（使用默认掩膜）比较
fig, axs = plt.subplots(2, 4, figsize=(25, 8)) # 创建子图
# 设置y轴标题并绘制数据
titles = [
"Original Rate Data", # 原始速率数据标题
"Cleaned Rate Data (Default Mask)", # 清理后的默认掩膜数据标题
"Cleaned Rate Data (Alternate Mask)", # 清理后的备用掩膜数据标题
"Cleaned Rate Data (Hand-Modified Mask)", # 清理后的手动修改掩膜数据标题
]
# 遍历数据列表和标题
for i, (data_list, title) in enumerate(
zip(
[
original_rate_data, # 原始速率数据
cleaned_rate_default_data, # 清理后的默认掩膜数据
cleaned_rate_alternate_data, # 清理后的备用掩膜数据
cleaned_rate_modified_data, # 清理后的手动修改掩膜数据
],
titles, # 标题列表
)
):
for j, data in enumerate(data_list):
ax = axs[j, i] # 获取当前子图
if j == 0:
ax.set_title(f'{title} \n {"Integration[100,:,:] | NRS1"}', fontsize=12) # 设置标题
else:
ax.set_title(f'{title} \n {"Integration[100,:,:] | NRS2"}', fontsize=12) # 设置标题
im = ax.imshow(data, origin="lower", aspect=30, clim=(-1e-1, 1e-1)) # 显示图像
fig.colorbar(im, ax=ax, pad=0.05, shrink=0.7, label="DN/s") # 添加颜色条
ax.set_xlabel("Pixel Column", fontsize=10) # 设置x轴标签
ax.set_ylabel("Pixel Row", fontsize=10) # 设置y轴标签
plt.tight_layout() # 调整子图布局
plt.show() # 显示图形
# 最终比较
# 绘制第一组光谱数据
plot_spectra(
[
x1d_nsclean_skipped[0], # 跳过的清理数据的第一组
x1d_nsclean_default[0], # 默认清理数据的第一组
x1d_nsclean_alternate[0], # 替代清理数据的第一组
x1d_nsclean_modified[0], # 修改后的清理数据的第一组
],
ext_num=100, # 扩展编号设置为100
wavelength_range=wavelength_range, # 波长范围
flux_range=flux_range, # 通量范围
)
# 绘制第二组光谱数据
plot_spectra(
[
x1d_nsclean_skipped[1], # 跳过的清理数据的第二组
x1d_nsclean_default[1], # 默认清理数据的第二组
x1d_nsclean_alternate[1], # 替代清理数据的第二组
x1d_nsclean_modified[1], # 修改后的清理数据的第二组
],
ext_num=100, # 扩展编号设置为100
wavelength_range=wavelength_range, # 波长范围
flux_range=flux_range, # 通量范围
)
最终说明： 替代的（基于剪辑的）和手动修改的掩膜产生了相似的结果。然而，手动修改的掩膜似乎引入了更少的高频噪声。
关于笔记本 #
作者： Melanie Clarke, Kayli Glidic; NIRSpec仪器团队
更新时间： 2024年2月29日。