使用 NSClean 清理 NIRSpec FS 产品中的残余 1/f 噪声#
笔记本目标#
本笔记本的目标是通过去除残余的 1/f 噪声生成清理后的 FS (_rate.fits) 文件。这些清理后的文件将作为第 3 级 (Spec3Pipeline) 管道的输入。
目录#
1. 介绍 #
JWST NIRSpec 仪器有许多特性和特点，观察者在规划观测和解释数据时应当了解。其中一个显著特征是在 NIRSpec 管道产品中提取的 1-D 光谱中出现负值和/或多余的通量，通常伴有不规则的波长依赖性起伏。这种伪影的原因是相关噪声，称为 1/f 噪声，源于低水平的探测器热不稳定性，在 2-D 计数率图像中表现为垂直条纹，尤其是在 NRS2 探测器的曝光中。虽然 IRS2 读出模式减少了这种效应，但并未完全消除。
为了解决这个问题，JWST 科学校准管道在 Spec2Pipeline 中集成了由 Bernard Rauscher 开发的外部包，称为 NSClean，在 NSCleanStep 下使用。该算法利用探测器的暗区在傅里叶空间中拟合背景模型。它需要一个输入掩膜来识别探测器的所有暗区。这个掩膜越全面，背景拟合效果越好。
在本笔记本中，我们将使用集成到管道中的 NSClean 算法，利用默认参数动态生成的掩膜来去除 1/f 噪声。在某些情况下，这个掩膜可能不够完整或过于严格，无法实现最佳的噪声去除。为了解决这个问题，我们演示了如何手动修改默认掩膜，以及如何通过调整 NSCleanStep 参数 创建替代掩膜。如有需要，请参阅 NSClean 文档 获取有关手动创建自定义掩膜的一些建议。
本笔记本提供了来自两个 GTO 观测的示例：来自程序 2757 的 F 矮星 GSC 03162-00665 的固定狭缝观测，以及来自程序 1189 的棕矮星 J03480772-6022270 的观测。
2. 导入库 #
# ---------- 设置CRDS环境变量 ----------
import os # 导入os模块，用于操作系统相关功能
import jwst # 导入jwst模块，用于JWST相关功能
# 设置CRDS上下文为指定的版本
os.environ['CRDS_CONTEXT'] = 'jwst_1210.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库，用于处理警告
# ------ JWST Calibration Pipeline Imports ------
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 # 从utils模块导入自定义函数，用于获取JWST文件和绘图
# Hide all log and warning messages. # 隐藏所有日志和警告信息
logging.disable(logging.ERROR) # 禁用错误级别的日志信息
warnings.simplefilter("ignore", RuntimeWarning) # 忽略运行时警告
3. 下载数据 #
本笔记本的输入数据包括使用 S200A1 子阵列 (SUBS200A1) 的固定缝隙 (FS) 观测和以 S200A1 作为主要缝隙的全帧观测。FS 观测的 F 型矮星 GSC 03162-00665 使用子阵列 SUBS200A1 和光栅/滤光片 G395M/F290LP，属于 GTO 计划 2757，具体为观测 2。该观测包含 15 次积分（每次 5 次积分/曝光；3-POINT-NOD），每次有 4 个组。棕矮星 J03480772-6022270 的 FS 观测使用全帧读出和光栅/滤光片 G140M/F100LP，属于 GTO 计划 1189，具体为观测 1。该观测包含 3 次积分（3-POINT-NOD），每次有 11 个组。
本笔记本专注于每个数据集的第一个抖动曝光 (00001)。然而，重要的是要注意，在继续进行 Spec3Pipeline 之前，所有曝光必须首先通过 Spec2Pipeline 处理。
# 定义下载目录
mast_products_dir = "./mast_products/"
# 检查目录是否存在
if not os.path.exists(mast_products_dir):
# 如果目录不存在，则创建该目录
os.makedirs(mast_products_dir)
# 子数组，S200A1；该笔记本专注于第一个抖动/节点。
subarray_obs_ids = ["jw02757002001_03104_00001"] # 子数组观测ID
subarray_detectors = [1] # 仅使用探测器 NRS1。
# 全帧，S200A1 主体；该笔记本专注于第一个抖动/节点。
full_obs_ids = ["jw01189001001_04101_00001"] # 全帧观测ID
full_detectors = [1, 2] # 使用探测器 NRS1 和 NRS2。
# 指定计数率数据产品。
rate_names = [] # 存储计数率文件名的列表
rate_types = [] # 存储计数率类型的列表
for obs_id in subarray_obs_ids: # 遍历子数组观测ID
for detector in subarray_detectors: # 遍历子数组探测器
rate_names.append(f"{obs_id}_nrs{detector}_rate.fits") # 添加计数率文件名
rate_types.append("subarray") # 添加计数率类型为子数组
for obs_id in full_obs_ids: # 遍历全帧观测ID
for detector in full_detectors: # 遍历全帧探测器
rate_names.append(f"{obs_id}_nrs{detector}_rate.fits") # 添加计数率文件名
rate_types.append("full") # 添加计数率类型为全帧
# 下载所有的 FITS 文件。
for name in rate_names: # 遍历所有计数率文件名
print(f"Downloading {name}") # 打印正在下载的文件名
get_jwst_file(name, mast_api_token=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) # 如果目录不存在，则创建该目录
# 原始数据（未应用NSClean）。
# 预计运行时间：1-2分钟。
start = tt.time() # 记录开始时间
for i in rate_names: # 遍历每个速率名称
print(f"Processing {i}...") # 打印正在处理的文件名
Spec2Pipeline.call( # 调用Spec2Pipeline处理数据
mast_products_dir + i, # 输入文件路径
save_results=True, # 保存结果
steps={ # 处理步骤
"nsclean": { # NSClean步骤
"skip": True # 跳过此步骤，避免去除相关读噪声（1/f噪声）
},
},
output_dir=stage2_nsclean_skipped_dir, # 输出目录
)
print(f"Saved {i[:-9]}" + "cal.fits") # 打印保存的校正文件名
print(f"Saved {i[:-9]}" + "x1d.fits") # 打印保存的一维光谱文件名
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 管道掩模）。
# 预计运行时间：7分钟。
start = tt.time() # 记录开始时间
for i in rate_names: # 遍历每个速率名称
print(f"Processing {i}...") # 打印正在处理的文件名
Spec2Pipeline.call( # 调用 Spec2Pipeline 进行数据处理
mast_products_dir + i, # 输入文件路径
save_results=True, # 保存结果
steps={ # 处理步骤
"nsclean": { # NSClean 步骤
"skip": False, # 不跳过此步骤
"save_mask": True, # 保存掩模
"save_results": True, # 保存处理结果
}, # 从 NIRSpec 图像中去除相关的读噪声（1/f 噪声）。
},
output_dir=stage2_nsclean_default_dir, # 输出目录
)
# 打印保存的文件名
print(f"Saved {i[:-9]}" + "mask.fits") # 保存的掩模文件
print(f"Saved {i[:-9]}" + "nsclean.fits") # 保存的清理后的文件
print(f"Saved {i[:-9]}" + "cal.fits") # 保存的校准文件
print(f"Saved {i[:-9]}" + "x1d.fits") # 保存的一维光谱文件
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 + "jw02757002001_03104_00001_nrs1_mask.fits", # NRS1掩膜文件
stage2_nsclean_default_dir + "jw01189001001_04101_00001_nrs1_mask.fits", # NRS1掩膜文件
stage2_nsclean_default_dir + "jw01189001001_04101_00001_nrs2_mask.fits", # NRS2掩膜文件
]
# 绘制每个相关的速率数据集和掩膜文件。
for rate_file, rate_type, mask_file in zip(
rate_names, rate_types, nsclean_default_masks # 将速率文件、类型和掩膜文件组合在一起
):
layout = "rows" if rate_type == "subarray" else "cols" # 根据速率类型决定布局
plot_dark_data(mast_products_dir + rate_file, mask_file, layout=layout, scale=8) # 绘制数据
5.2 原始数据与清理数据的比较（默认管道掩模） #
现在我们可以将清理后的数据（使用默认管道掩模）与原始速率文件进行比较，并验证1/f噪声是否已被减少。
在许多情况下，清理过程会在速率文件中引入新的伪影。这些伪影应仔细检查，并与噪声减少的好处进行权衡。如果瞬态伪影，如雪球，干扰了清理过程，手动编辑掩模以排除这些区域在背景拟合中的考虑可能是有益的。为此，可以尝试调整异常值检测阈值或直接编辑掩模数组中的特定像素（将在接下来的几节中探讨）。否则，请参考NSClean文档以获取有关手动编辑的其他建议。
请注意，在下面的图像中，有一些值与原始速率文件存在较大的相对差异（在下面的相对差异图像中显示）。这些是清理过程的伪影。
此外，还有更广泛的低水平残余背景效应（在右侧的相对差异图像中显示，散布的异常值通过sigma裁剪识别，并被掩模隐藏）。这些包括我们试图去除的背景模式：在色散方向上的1/f噪声变化以及全帧数据中图像顶部和底部的画框效应。然而，清理过程中可能还会因过拟合暗数据而引入低水平伪影。
仔细检查这两幅残余图像，以了解清理过程对数据的影响。
# 绘制原始数据和清理后的数据，以及残差图。
# 定义清理后的默认掩膜文件路径列表
cleaned_default_masks = [
stage2_nsclean_default_dir + "jw02757002001_03104_00001_nrs1_nsclean.fits", # 第一个清理后的文件路径
stage2_nsclean_default_dir + "jw01189001001_04101_00001_nrs1_nsclean.fits", # 第二个清理后的文件路径
stage2_nsclean_default_dir + "jw01189001001_04101_00001_nrs2_nsclean.fits", # 第三个清理后的文件路径
]
# 遍历每个关联的rateint数据和清理后的文件
for rate_file, cleaned_file in zip(rate_names, cleaned_default_masks):
layout = "rows" if rate_type == "subarray" else "cols" # 根据rate_type决定布局方式
scale = 8 if "2757" in rate_file else 1 # 根据文件名决定缩放比例
# 调用函数绘制清理后的数据
plot_cleaned_data(
mast_products_dir + rate_file, cleaned_file, layout=layout, scale=scale # 绘制原始数据和清理后的数据
)
比较从清理后的数据中提取的光谱与从原始速率文件中提取的光谱。请注意，以下针对全帧数据的图表显示了所有的FS（全谱），但主要来源于FS S200A1。
# 全帧1D提取光谱（绘制所有FS）
# 定义未清理的光谱文件路径列表
x1d_nsclean_skipped = [
stage2_nsclean_skipped_dir + "jw02757002001_03104_00001_nrs1_x1d.fits", # 文件1路径
stage2_nsclean_skipped_dir + "jw01189001001_04101_00001_nrs1_x1d.fits", # 文件2路径
stage2_nsclean_skipped_dir + "jw01189001001_04101_00001_nrs2_x1d.fits", # 文件3路径
]
# 定义默认清理的光谱文件路径列表
x1d_nsclean_default = [
stage2_nsclean_default_dir + "jw02757002001_03104_00001_nrs1_x1d.fits", # 文件1路径
stage2_nsclean_default_dir + "jw01189001001_04101_00001_nrs1_x1d.fits", # 文件2路径
stage2_nsclean_default_dir + "jw01189001001_04101_00001_nrs2_x1d.fits", # 文件3路径
]
# 绘制每个全帧中的FS
for original, cleaned, rate_type in zip(
x1d_nsclean_skipped, x1d_nsclean_default, rate_types # 遍历未清理和清理的光谱文件及其类型
):
if rate_type == "subarray": # 如果光谱类型为子阵列，则跳过
continue
for slit in range(5): # 在全帧中有5个FS
plot_spectra([original, cleaned], ext_num=slit, scale_percent=1) # 绘制光谱
注意： 主光谱缝隙 S200A1 中的亮谱在清理后几乎没有变化（在某些波长下，清理后的光谱根据比率显示出微弱的通量增加）。次光谱缝隙（S200A2、S400A1、S1600A1、S200B1）中的微弱背景光谱已针对负通量进行了修正（通过接近零的比率指示），并考虑了波长依赖的变化。
# Subarray SUBS200A1 1D extracted spectra.
for original, cleaned, rate_type in zip(
x1d_nsclean_skipped, x1d_nsclean_default, rate_types # 将原始、清理后的数据和速率类型打包在一起
):
if rate_type != "subarray": # 如果速率类型不是子阵列，则跳过
continue
for slit in range(5): # 处理5个狭缝，针对S200A1
plot_spectra([original, cleaned], ext_num=slit, scale_percent=4) # 绘制光谱，设置扩展编号和缩放百分比
# Subarray SUBS200A1 1D extracted Spectra -- smaller wavelength region of interest.
wavelength_range = {"nrs1": [4.3, 4.7]} # 定义感兴趣的波长范围
flux_range = {"nrs1": [0.004, 0.007]} # 定义感兴趣的通量范围
for original, cleaned, rate_type in zip(
x1d_nsclean_skipped, x1d_nsclean_default, rate_types # 将原始、清理后的数据和速率类型打包在一起
):
if rate_type != "subarray": # 如果速率类型不是子阵列，则跳过
continue
for slit in range(5): # 处理5个狭缝，针对S200A1
plot_spectra(
[original, cleaned], # 绘制原始和清理后的光谱
ext_num=slit, # 设置扩展编号
scale_percent=5, # 设置缩放百分比
wavelength_range=wavelength_range, # 设置波长范围
flux_range=flux_range, # 设置通量范围
)
注意：
对于子阵列数据中较弱的光谱，我们可以看到清理过程中连续谱水平的微小变化。
原始数据中约4.55微米的过量通量在清理后得到了修正。
6. 使用 NSClean（替代掩膜）清除 1/f 噪声 #
对于某些数据集，掩膜整个科学区域可能会过度掩盖探测器的暗区，这些区域可以用来改善背景拟合。过度掩膜可能会在清理过程中引入一些高频噪声，表现为光谱轨迹上的垂直条纹。例如，请参见第 6.2 节中子阵列（PID 2757）的清理速率数据。由于该区域完全被掩膜，一些残余伪影被引入。
此外，对于某些数据集，探测器上可能存在多个未被 WCS 边界框掩膜的照明区域。在子阵列数据中，请注意掩膜未覆盖光谱轨迹的蓝端。这一区域是被照明的，但尚未校准。
在某些情况下，使用替代算法构建掩膜可能会更有利。在这里，我们不使用边界框，而是迭代地掩膜任何超过背景 2 个标准差的数据。对于亮源，这样可以在光谱轨迹附近保留更多的暗数据，从而可能改善背景拟合。
然而，请注意，过度清理可能会影响光谱的连续性水平，如果掩膜中包含了过多或过少的照明数据。再次强调，生成的掩膜和输出光谱应仔细检查，以权衡清理的好处与对光谱的影响。
要调整此掩膜中的照明检测，请尝试修改下面的 n_sigma 参数。较高的值将识别较少的照明，较低的值将识别更多的照明。
# 设置用于运行NSClean的目录，并使用用户提供的掩模。
stage2_nsclean_alternate_dir = "./stage2_nsclean_alternate/" # 定义NSClean的备用目录路径
if not os.path.exists(stage2_nsclean_alternate_dir): # 检查目录是否存在
os.makedirs( # 如果目录不存在，则创建该目录
stage2_nsclean_alternate_dir # 使用定义的路径创建目录
) # 创建目录，如果它不存在。
# 1/f 噪声清理数据（备用 NSClean 管道掩膜）。
# 估计运行时间：9 分钟。
start = tt.time() # 记录开始时间
for indx, i in enumerate(rate_names): # 遍历 rate_names 列表中的每个元素及其索引
print(f"Processing {i}... ") # 打印正在处理的文件名
Spec2Pipeline.call( # 调用 Spec2Pipeline 的 call 方法
mast_products_dir + i, # 输入文件路径
save_results=True, # 保存结果
steps={ # 定义处理步骤
"nsclean": { # NSClean 步骤
"skip": False, # 不跳过此步骤
"save_mask": True, # 保存掩膜
"n_sigma": 2, # 设置 sigma 值
"mask_spectral_regions": False, # 不掩盖光谱区域
"save_results": True, # 保存结果
}, # 从 NIRSpec 图像中去除相关的读噪声（1/f 噪声）
},
output_dir=stage2_nsclean_alternate_dir, # 输出目录
)
# 打印保存的文件名
print(f"Saved {i[:-9]}" + "mask.fits") # 保存的掩膜文件
print(f"Saved {i[:-9]}" + "nsclean.fits") # 保存的清理后文件
print(f"Saved {i[:-9]}" + "cal.fits") # 保存的校准文件
print(f"Saved {i[:-9]}" + "x1d.fits") # 保存的一维光谱文件
end = tt.time() # 记录结束时间
print("Run time: ", round(end - start, 1) / 60.0, " min") # 打印运行时间
6.1 验证掩模（备用掩模） #
检查掩模与速率数据，以确保它仅保留探测器的暗区。
# 绘制带有被屏蔽区域的速率数据。
# 从管道中构建的即时掩膜列表。
nsclean_alternate_masks = [
stage2_nsclean_alternate_dir + "jw02757002001_03104_00001_nrs1_mask.fits", # 第一幅掩膜文件路径
stage2_nsclean_alternate_dir + "jw01189001001_04101_00001_nrs1_mask.fits", # 第二幅掩膜文件路径
stage2_nsclean_alternate_dir + "jw01189001001_04101_00001_nrs2_mask.fits", # 第三幅掩膜文件路径
]
# 绘制每组相关的速率数据和掩膜文件。
for rate_file, rate_type, mask_file in zip(
rate_names, rate_types, nsclean_alternate_masks # 将速率文件名、类型和掩膜文件组合在一起
):
layout = "rows" if rate_type == "subarray" else "cols" # 根据速率类型确定布局方式
plot_dark_data(mast_products_dir + rate_file, mask_file, layout=layout, scale=8) # 绘制数据
6.2 原始数据与清理数据的比较（替代掩膜） #
# 绘制原始数据和清理后的数据，以及残差图。
# 定义清理后的备用掩膜文件路径列表
cleaned_alternate_masks = [
stage2_nsclean_alternate_dir + "jw02757002001_03104_00001_nrs1_nsclean.fits", # 第一个清理后的文件路径
stage2_nsclean_alternate_dir + "jw01189001001_04101_00001_nrs1_nsclean.fits", # 第二个清理后的文件路径
stage2_nsclean_alternate_dir + "jw01189001001_04101_00001_nrs2_nsclean.fits", # 第三个清理后的文件路径
]
# 遍历每个关联的速率数据文件和清理后的文件
for rate_file, cleaned_file in zip(rate_names, cleaned_alternate_masks):
layout = "rows" if rate_type == "subarray" else "cols" # 根据速率类型设置布局
scale = 8 if "2757" in rate_file else 1 # 根据文件名设置缩放比例
# 调用函数绘制清理后的数据
plot_cleaned_data(
mast_products_dir + rate_file, cleaned_file, layout=layout, scale=scale # 绘制原始数据和清理后的数据
)
比较从清理后的数据中提取的光谱与从原始速率文件中提取的光谱。请注意，以下针对全帧数据的图表显示了所有的FS（全帧），但主要来源于FS S200A1。
# 全帧 1D 提取的光谱（绘制所有 FS）
# 定义清理后的光谱文件路径列表
x1d_nsclean_alternate = [
stage2_nsclean_alternate_dir + "jw02757002001_03104_00001_nrs1_x1d.fits", # 文件1路径
stage2_nsclean_alternate_dir + "jw01189001001_04101_00001_nrs1_x1d.fits", # 文件2路径
stage2_nsclean_alternate_dir + "jw01189001001_04101_00001_nrs2_x1d.fits", # 文件3路径
]
# 绘制每个全帧中的 FS
for original, cleaned, rate_type in zip(
x1d_nsclean_skipped, x1d_nsclean_alternate, rate_types # 将原始、清理后的光谱和速率类型打包
):
if rate_type == "subarray": # 如果速率类型是子阵列，则跳过
continue
for slit in range(5): # 在全帧中有 5 个 FS
plot_spectra([original, cleaned], ext_num=slit, scale_percent=1) # 绘制光谱，指定扩展号和缩放百分比
注意： 使用此掩模进行清理时，主光谱 S200A1 中的亮谱在清理后仍然显示出较小的差异（在某些波长下，清理后的光谱根据比率显示出微妙的通量增加）。次光谱中的微弱背景光谱（S200A2, S400A1, S1600A1, S200B1）已针对负通量进行了修正（通过接近零的比率指示）以及波长依赖的变化。
# Subarray SUBS200A1 1D extracted spectra.
# 在全帧中绘制每个FS（光谱）。
for original, cleaned, rate_type in zip(
x1d_nsclean_skipped, x1d_nsclean_alternate, rate_types
):
if rate_type != "subarray": # 如果不是子阵列，则跳过
continue
for slit in range(5): # 在全帧中有5个FS
plot_spectra([original, cleaned], ext_num=slit, scale_percent=4) # 绘制光谱，设置扩展编号和缩放百分比
# Subarray SUBS200A1 1D extracted spectra -- smaller wavelength region of interest.
# 定义感兴趣的较小波长范围
wavelength_range = {"nrs1": [4.3, 4.7]} # 波长范围
flux_range = {"nrs1": [0.004, 0.007]} # 流量范围
for original, cleaned, rate_type in zip(
x1d_nsclean_skipped, x1d_nsclean_alternate, rate_types
):
if rate_type != "subarray": # 如果不是子阵列，则跳过
continue
for slit in range(5): # 处理1个狭缝 S200A1
plot_spectra(
[original, cleaned], # 绘制原始和清理后的光谱
ext_num=slit, # 设置扩展编号
scale_percent=5, # 设置缩放百分比
wavelength_range=wavelength_range, # 设置波长范围
flux_range=flux_range, # 设置流量范围
)
注意： 使用此掩模进行清理仍然会去除接近4.55微米的波长依赖性变化。与原始光谱相比，基于剪切的掩模可能会引入稍微少一些的高频噪声。
7. 使用 NSClean 清理 1/f 噪声（手动修改的掩膜） #
在某些情况下，可能需要手动生成掩膜。在这里，我们介绍 一种 方法来手动修改两个数据集的掩膜（缩小一些被掩膜的 FS 区域以包含更多背景），从管道的默认掩膜输出开始。值得注意的是，使用此方法修改的掩膜不一定会优于之前的两个选项。
# 设置用于运行 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() # 复制数据。
if "2757" in mask:
# 步骤1：将默认掩膜区域设置回True。
mask_data[:, :] = True
# 步骤2：手动重新定义掩膜区域。
mask_data[26:45, 680:] = False
elif "1189" and "nrs1" in mask:
# 将特定区域设置为True。
mask_data[1054:1150, 650:] = True
# 将特定区域设置为False。
mask_data[1077:1098, 650:] = 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 + "jw02757002001_03104_00001_nrs1_mask_modified.fits", # 第一个掩膜文件路径
stage2_nsclean_modified_dir + "jw01189001001_04101_00001_nrs1_mask_modified.fits", # 第二个掩膜文件路径
stage2_nsclean_modified_dir + "jw01189001001_04101_00001_nrs2_mask_modified.fits", # 第三个掩膜文件路径
]
# 遍历每个相关的速率数据和掩膜文件进行绘图。
for rate_file, rate_type, mask_file in zip(
rate_names, rate_types, nsclean_modified_masks # 将速率文件名、类型和掩膜文件组合在一起
):
layout = "rows" if rate_type == "subarray" else "cols" # 根据速率类型决定布局方式
plot_dark_data(mast_products_dir + rate_file, mask_file, layout=layout, scale=8) # 绘制速率数据，使用指定的掩膜和布局
注意： 在修改PID 1189（NRS1）的掩膜时，我们选择了一个区域进行解掩膜，然后以不同的方式重新掩膜。然而，这一过程不小心导致了一些之前被掩膜的NaN值（在NRS1的暗数据图中看到的白色像素）被解掩膜。因此，在修改掩膜时应谨慎行事。
# 1/f 噪声清理数据（用户提供的掩码）。
# 预计运行时间：8分钟。
start = tt.time() # 记录开始时间
for indx, i in enumerate(rate_names): # 遍历每个速率名称及其索引
print(f"Processing {i}... ") # 打印正在处理的文件名
Spec2Pipeline.call( # 调用Spec2Pipeline进行处理
mast_products_dir + i, # 输入文件路径
save_results=True, # 保存结果
steps={ # 处理步骤
"nsclean": { # 噪声清理步骤
"skip": False, # 不跳过此步骤
"save_mask": True, # 保存掩码
"save_results": True, # 保存处理结果
"user_mask": nsclean_modified_masks[indx], # 使用用户提供的掩码
}, # 从NIRSpec图像中去除相关的读取噪声（1/f噪声）。
},
output_dir=stage2_nsclean_modified_dir, # 输出目录
)
# 打印保存的文件名
print(f"Saved {i[:-9]}" + "mask.fits") # 保存的掩码文件
print(f"Saved {i[:-9]}" + "nsclean.fits") # 保存的清理后的文件
print(f"Saved {i[:-9]}" + "cal.fits") # 保存的校准文件
print(f"Saved {i[:-9]}" + "x1d.fits") # 保存的一维光谱文件
end = tt.time() # 记录结束时间
print("Run time: ", round(end - start, 1) / 60.0, " min") # 打印运行时间
7.2 原始数据与清理数据（手动修改的掩膜）比较 #
# 绘制原始数据、清理后的数据以及残差图。
cleaned_modified_masks = [
stage2_nsclean_modified_dir + "jw02757002001_03104_00001_nrs1_nsclean.fits", # 添加第一个清理后的文件路径
stage2_nsclean_modified_dir + "jw01189001001_04101_00001_nrs1_nsclean.fits", # 添加第二个清理后的文件路径
stage2_nsclean_modified_dir + "jw01189001001_04101_00001_nrs2_nsclean.fits", # 添加第三个清理后的文件路径
]
# 绘制每组相关的rateint数据和清理后的文件。
for rate_file, cleaned_file in zip(rate_names, cleaned_modified_masks): # 遍历rate_names和cleaned_modified_masks的配对
layout = "rows" if rate_type == "subarray" else "cols" # 根据rate_type决定布局方式
scale = 8 if "2757" in rate_file else 1 # 根据文件名决定缩放比例
plot_cleaned_data( # 调用绘图函数
mast_products_dir + rate_file, cleaned_file, layout=layout, scale=scale # 传入原始文件路径、清理后的文件路径、布局和缩放比例
)
比较从清理后的数据中提取的光谱，并将其与从原始速率文件中提取的光谱进行比较。
# 全帧 1D 提取光谱（绘制所有 FS）
# 定义修改后的无噪声清理的 1D 光谱文件路径
x1d_nsclean_modified = [
stage2_nsclean_modified_dir + "jw02757002001_03104_00001_nrs1_x1d.fits", # 文件路径 1
stage2_nsclean_modified_dir + "jw01189001001_04101_00001_nrs1_x1d.fits", # 文件路径 2
stage2_nsclean_modified_dir + "jw01189001001_04101_00001_nrs2_x1d.fits", # 文件路径 3
]
# 在全帧中绘制每个 FS
for original, cleaned, rate_type in zip(
x1d_nsclean_skipped, x1d_nsclean_modified, rate_types # 逐个获取原始、清理后的光谱和速率类型
):
if rate_type == "subarray": # 如果速率类型是子阵列，则跳过
continue
for slit in range(5): # 在全帧中有 5 个 FS
plot_spectra([original, cleaned], ext_num=slit, scale_percent=1) # 绘制光谱，设置扩展号和缩放百分比
注意： 再次使用这个修改过的掩模进行清理，主缝隙 S200A1 中的亮谱在清理后仍然显示出很小的差异（在某些波长下，清理后的光谱根据比率显示出微妙的通量增加）。应用于次缝隙（S200A2, S400A1, S1600A1, S200B1）中微弱背景光谱的掩模与管道提供的默认掩模保持一致。因此，使用该掩模获得的这些缝隙的结果与默认掩模相比保持不变。
# Subarray SUBS200A1 1D Extracted Spectra
# 遍历每个原始和清理后的光谱数据以及对应的速率类型
for original, cleaned, rate_type in zip(
x1d_nsclean_skipped, x1d_nsclean_modified, rate_types
):
# 仅处理速率类型为"subarray"的情况
if rate_type != "subarray":
continue
# 遍历全帧中的5个光谱切口
for slit in range(5): # 5 FS in full-frame
# 绘制原始和清理后的光谱
plot_spectra([original, cleaned], ext_num=slit, scale_percent=4)
# Subarray SUBS200A1 1D Extracted Spectra -- smaller wavelength region of interest
# 定义感兴趣的波长范围
wavelength_range = {"nrs1": [4.3, 4.7]}
# 定义感兴趣的通量范围
flux_range = {"nrs1": [0.004, 0.007]}
# 遍历每个原始和清理后的光谱数据以及对应的速率类型
for original, cleaned, rate_type in zip(
x1d_nsclean_skipped, x1d_nsclean_modified, rate_types
):
# 仅处理速率类型为"subarray"的情况
if rate_type != "subarray":
continue
# 遍历1个切口 S200A1
for slit in range(5): # 1 slit S200A1
# 绘制原始和清理后的光谱，指定波长和通量范围
plot_spectra(
[original, cleaned],
ext_num=slit,
scale_percent=5,
wavelength_range=wavelength_range,
flux_range=flux_range,
)
注意：对于子阵列数据中的较弱光谱，我们可以看到清理过程对连续谱水平的微小变化。使用此掩模进行清理仍然会去除接近4.55微米的波长依赖性变化。
8. 结论 #
下面的最终图表展示了计数率图像和相应的1D提取光谱并排比较不同的清理方法：原始图像（未应用NSClean）、清理后的计数率图像（使用默认的管道掩模）、清理后的计数率图像（使用替代的管道掩模），以及最后，清理后的计数率图像（使用手动修改的掩模）。
请注意，本笔记本中呈现的结果可能会因不同的数据集而有所不同（例如，不同亮度、空间范围等的目标）。鼓励用户使用不同的掩模方法探索NSClean，以确定最佳结果。
清理算法的输出现在可以进行进一步处理。上述Spec2Pipeline运行生成的(_cal.fits)文件可以作为输入用于Spec3Pipeline，以生成最终的组合光谱。
# 加载速率数据和清理后的掩膜
original_rate_data = [
# 打开每个速率名称对应的FITS文件，并提取数据部分
fits.open(mast_products_dir + rate_name)[1].data for rate_name in rate_names
]
cleaned_rate_default_data = [
# 打开每个清理后的默认掩膜对应的FITS文件，并提取数据部分
fits.open(cleaned_default_mask)[1].data
for cleaned_default_mask in cleaned_default_masks
]
cleaned_rate_alternate_data = [
# 打开每个清理后的替代掩膜对应的FITS文件，并提取数据部分
fits.open(cleaned_alternate_mask)[1].data
for cleaned_alternate_mask in cleaned_alternate_masks
]
cleaned_rate_modified_data = [
# 打开每个清理后的手动修改掩膜对应的FITS文件，并提取数据部分
fits.open(cleaned_modified_mask)[1].data
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:
# 将数据中的NaN值替换为0
data[np.isnan(data)] = 0
# 设置颜色范围(clim)值
clim_first_row = (-0.5, 2) # 第一行的颜色范围
clim_subsequent_rows = (-1e-1, 1e-1) # 后续行的颜色范围
# 原始数据与清理后的数据（使用默认掩膜）
fig, axs = plt.subplots(3, 4, figsize=(25, 15), gridspec_kw={"hspace": 0.9})
# 设置标题并绘制数据
titles = [
"Original Rate Data", # 原始速率数据
"Cleaned Rate Data (Default Mask)", # 清理后的速率数据（默认掩膜）
"Cleaned Rate Data (Alternate Mask)", # 清理后的速率数据（替代掩膜）
"Cleaned Rate Data (Hand-Modified Mask)", # 清理后的速率数据（手动修改掩膜）
]
for i, title in enumerate(titles):
for j, data in enumerate(
[
original_rate_data,
cleaned_rate_default_data,
cleaned_rate_alternate_data,
cleaned_rate_modified_data,
][i]
):
ax = axs[j, i] # 获取当前子图
aspect = 10 if j == 0 else "auto" # 设置纵横比
clim = (
clim_first_row if j == 0 else clim_subsequent_rows
) # 根据行调整clim值
# 获取速率名称，处理标题
rate_name = rate_names[0][:-10] if i == 1 and j == 0 else rate_names[1][:-10]
ax.set_title(f"{rate_name} \n {title}", fontsize=12) # 设置子图标题
im = ax.imshow(data, origin="lower", aspect=aspect, clim=clim) # 绘制图像
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() # 显示图形
def plot_spectra(
spec_list,
ext_num=1,
scale_percent=2.0,
wavelength_range=None,
flux_range=None,
xlim_low=None,
):
"""
绘制清理后的提取1D光谱，并与原始（未应用NSClean）提取的1D光谱进行比较。
参数：
----------
spec_list: list
包含原始/清理后的1D提取数据的FITS文件路径列表。
scale_percent: float
用于确定可视化强度范围的缩放因子。
ext_num: int
要绘制的切片的索引/扩展。
EXTVER头部值。默认值为2D数据的1。
wavelength_range: dict
波长范围（x轴） {'nrs1': [3.6, 3.65], 'nrs2': [1.65, 1.75]}。
flux_range: dict
流量范围（y轴） {'nrs1': [1, 2], 'nrs2': [1, 2]}。
xlim_low: int
定义1D光谱的下波长端。对BOTS数据有帮助。
"""
if wavelength_range is None:
wavelength_range = {}
# 打开FITS文件。
original_hdul = fits.open(spec_list[0]) # 打开原始光谱文件
cleaned_hdul = fits.open(spec_list[1]) # 打开清理后的光谱文件
if len(spec_list) == 3:
alternate_hdul = fits.open(spec_list[2]) # 打开备用光谱文件
elif len(spec_list) == 4:
alternate_hdul = fits.open(spec_list[2]) # 打开备用光谱文件
handmod_hdul = fits.open(spec_list[3]) # 打开手动修改的光谱文件
else:
alternate_hdul = None # 如果没有备用文件，则设置为None
handmod_hdul = None # 如果没有手动修改文件，则设置为None
# 查找光谱扩展（EXTRACT1D）。
for extnum in range(len(original_hdul)):
hdu = original_hdul[extnum] # 获取当前扩展
if hdu.name != "EXTRACT1D": # 如果扩展名称不是EXTRACT1D，则跳过
continue
slit_name = hdu.header["SLTNAME"] # 获取光谱缝隙名称
if hdu.header["EXTVER"] == ext_num: # 如果扩展版本与指定版本匹配
# 将原始光谱和清理后的光谱一起绘制。
fig, ax = plt.subplots(1, 2, figsize=(10, 5)) # 创建子图
plt.suptitle(
f"Compare 1D Spectra for {os.path.basename(spec_list[0])};"
f"EXP_TYPE/Slit = {slit_name}",
fontsize=15,
)
if "nrs1" in spec_list[0] and xlim_low is not None:
xlim_low = xlim_low # 如果存在nrs1并且xlim_low不为None，则保持不变
else:
xlim_low = 0 # 否则设置为0
ax[0].step(
hdu.data["WAVELENGTH"][xlim_low:], # 绘制原始光谱的波长
hdu.data["FLUX"][xlim_low:], # 绘制原始光谱的流量
linewidth=1,
label="Original", # 标签为“原始”
)
ax[0].step(
cleaned_hdul[extnum].data["WAVELENGTH"][xlim_low:], # 绘制清理后的光谱的波长
cleaned_hdul[extnum].data["FLUX"][xlim_low:], # 绘制清理后的光谱的流量
linewidth=1,
label="Cleaned", # 标签为“清理后的”
)
ax[0].set_xlabel(f"Wavelength ({hdu.header['TUNIT1']})", fontsize=12) # 设置x轴标签
ax[0].set_ylabel(f"Flux ({hdu.header['TUNIT2']})", fontsize=12) # 设置y轴标签
ax[0].set_title("1D Extracted Spectra", fontsize=12) # 设置子图标题
# 绘制光谱之间的差异作为比率。
diff = (
cleaned_hdul[extnum].data["FLUX"][xlim_low:] # 清理后的流量
/ hdu.data["FLUX"][xlim_low:] # 原始流量
)
ax[1].step(
hdu.data["WAVELENGTH"][xlim_low:], # 绘制波长
diff, # 绘制差异
linewidth=1,
label="Cleaned/Original", # 标签为“清理后的/原始”
)
ax[1].set_xlabel(f"Wavelength ({hdu.header['TUNIT1']})", fontsize=12) # 设置x轴标签
ax[1].set_ylabel("Cleaned/Original", fontsize=12) # 设置y轴标签
ax[1].set_title("Difference After Cleaning", fontsize=12) # 设置子图标题
# print("Average Difference = {}".format(np.nanmean(diff))) # 可选：打印平均差异
# 如果可用，还绘制备用光谱。
if alternate_hdul is not None:
ax[0].step(
alternate_hdul[extnum].data["WAVELENGTH"][xlim_low:], # 绘制备用光谱的波长
alternate_hdul[extnum].data["FLUX"][xlim_low:], # 绘制备用光谱的流量
linewidth=1,
label="Cleaned (alternate mask)", # 标签为“清理（备用掩模）”
)
diff2 = (
alternate_hdul[extnum].data["FLUX"][xlim_low:] # 备用流量
/ hdu.data["FLUX"][xlim_low:] # 原始流量
)
ax[1].step(
hdu.data["WAVELENGTH"][xlim_low:], # 绘制波长
diff2, # 绘制备用光谱的差异
linewidth=1,
label="Cleaned (alternate mask)/Original", # 标签为“清理（备用掩模）/原始”
)
if handmod_hdul is not None:
ax[0].step(
handmod_hdul[extnum].data["WAVELENGTH"][xlim_low:], # 绘制手动修改光谱的波长
handmod_hdul[extnum].data["FLUX"][xlim_low:], # 绘制手动修改光谱的流量
linewidth=1,
label="Cleaned (hand-modified mask)", # 标签为“清理（手动修改掩模）”
color="Purple" # 设置颜色为紫色
)
diff3 = (
handmod_hdul[extnum].data["FLUX"][xlim_low:] # 手动修改流量
/ hdu.data["FLUX"][xlim_low:] # 原始流量
)
ax[1].step(
hdu.data["WAVELENGTH"][xlim_low:], # 绘制波长
diff3, # 绘制手动修改光谱的差异
linewidth=1,
label="Cleaned (hand-modified mask)/Original", # 标签为“清理（手动修改掩模）/原始”
)
# 如果需要，设置绘图的y范围。
if flux_range is not None:
for key, y_range in flux_range.items():
if key.lower() in spec_list[0].lower() and y_range is not None:
if len(y_range) == 2:
ax[0].set_ylim(y_range) # 设置y轴范围
ax[1].set_ylim(
[
np.nanpercentile(diff, scale_percent), # 设置y轴范围
np.nanpercentile(diff, 100 - scale_percent),
]
)
# ax[1].set_ylim(0.5, 1.5) # 可选：设置y轴范围
else:
all_flux = [
hdu.data["FLUX"], # 原始流量
cleaned_hdul[extnum].data["FLUX"], # 清理后的流量
]
y_range_ax0 = [
np.nanpercentile(all_flux[0], scale_percent), # 计算y轴范围
np.nanpercentile(all
# 最终比较 - 全帧 (所有FS)
# 在全帧中绘制每个FS
for original, cleaned_og, cleaned_alt, cleaned_mod, rate_type in zip(
x1d_nsclean_skipped, # 原始数据，跳过的清理数据
x1d_nsclean_default, # 默认清理数据
x1d_nsclean_alternate, # 替代清理数据
x1d_nsclean_modified, # 修改后的清理数据
rate_types, # 数据类型
):
if rate_type == "subarray": # 如果数据类型是子数组，则跳过
continue
for slit in range(5): # 在全帧中有5个FS
plot_spectra( # 绘制光谱
[original, cleaned_og, cleaned_alt, cleaned_mod], # 要绘制的光谱数据
ext_num=slit, # 当前光谱的扩展编号
scale_percent=4, # 缩放百分比
)
# 比较 -- 子阵列 SUBS200A1 1D 提取的光谱
# 感兴趣的波长范围
wavelength_range = {"nrs1": [4.3, 4.7]} # 定义波长范围
flux_range = {"nrs1": [0.004, 0.007]} # 定义光通量范围
# 遍历原始数据和清理后的数据
for original, cleaned_og, cleaned_alt, cleaned_mod, rate_type in zip(
x1d_nsclean_skipped, # 原始未清理数据
x1d_nsclean_default, # 默认清理后的数据
x1d_nsclean_alternate, # 替代清理后的数据
x1d_nsclean_modified, # 修改后的清理数据
rate_types, # 数据类型
):
if rate_type != "subarray": # 检查数据类型是否为子阵列
continue # 如果不是，跳过当前循环
for slit in range(5): # 遍历5个光纤（FS）在全帧中
plot_spectra( # 绘制光谱
[original, cleaned_og, cleaned_alt, cleaned_mod], # 传入不同的数据集
ext_num=slit, # 当前光纤编号
wavelength_range=wavelength_range, # 波长范围
flux_range=flux_range, # 光通量范围
)
最终说明：
总体而言，在使用每种掩模（默认、替代、手动修改）进行清理后，S200A1的全帧数据中的亮谱显示出最小的差异。
全帧数据的次级狭缝中的背景光谱已针对负通量进行了修正，无论使用何种掩模。
关于S200A1子阵列数据中的微弱源，我们观察到清理过程中，使用任何掩模时连续谱水平有轻微变化。每种掩模有效地去除了接近4.55微米的波长依赖性变化。然而，默认掩模在清理过程中引入了一些高频噪声，导致光谱轨迹上出现垂直条纹。相比之下，手动修改的掩模引入的高频噪声较少，而替代（基于剪辑的）掩模引入的高频噪声最少。
关于笔记本 #
作者： Melanie Clarke, Kayli Glidic; NIRSpec仪器团队
更新时间： 2024年2月29日。