使用 NSClean 清理 NIRSpec MOS 产品中的残余 1/f 噪声#
Notebook 目标#
本 Notebook 的目标是通过去除残余的 1/f 噪声来生成清理后的 MOS (_rate.fits) 文件。这些清理后的文件将作为第 3 级 (Spec3Pipeline) 管道的输入。
目录#
1. 介绍 #
JWST NIRSpec 仪器具有一些观察者在规划观测和解释数据时应注意的特性和特点。其中一个显著特征是在提取的 1-D 光谱中出现负值和/或多余的通量，通常伴随不规则的波长依赖性起伏。这种伪影的原因是相关噪声，称为 1/f 噪声，源于低水平的探测器热不稳定性，在 2-D 计数率图像中表现为垂直条带，特别是在 NRS2 探测器的曝光中。虽然 IRS2 读出模式减少了这种效应，但并未完全消除。
为了应对这一问题，JWST 科学校准管道在 Spec2Pipeline 中集成了由 Bernard Rauscher 开发的外部包，称为 NSClean，并在 NSCleanStep 下使用。该算法利用探测器的暗区在傅里叶空间中拟合背景模型。它需要一个输入掩模来识别探测器的所有暗区。这个掩模越全面，背景拟合效果越好。
在本 Notebook 中，我们将使用集成到管道中的 NSClean 算法，利用默认参数动态生成的掩模来去除 1/f 噪声。在某些情况下，这个掩模可能不够完整或过于严格，无法实现最佳的噪声去除。为了解决这个问题，我们演示了如何手动修改默认掩模，以及如何通过调整 NSCleanStep 参数 创建替代掩模。如有需要，请参阅 NSClean 文档，获取有关手动创建自定义掩模的一些建议。
本 Notebook 使用来自 CEERS 项目的 MOS 观测数据，光栅/滤光片为 G395M/F290LP，作为 JWST 早期发布科学计划 ERS-1345 的第 61 次观测示例。
2. 导入库 #
# ---------- 设置CRDS环境变量 ----------
import os # 导入操作系统模块
import jwst # 导入JWST模块
# 设置CRDS上下文为指定的.pmap文件
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 # 导入时间库，用于时间相关操作
import logging # 导入日志库，用于记录日志信息
import 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模块导入自定义函数
# Hide all log and warning messages. # 隐藏所有日志和警告信息
logging.disable(logging.ERROR) # 禁用错误级别的日志信息
warnings.simplefilter("ignore", RuntimeWarning) # 忽略运行时警告
3. 下载数据 #
输入数据来自于CEERS项目的MOS观测，使用了光栅/滤光片G395H/F290LP。该数据集是JWST早期发布科学计划 ERS-1345的一部分，具体为观测61。它包含3次积分（3个抖动点；3-SHUTTER-SLITLET），每次积分有9个组。此笔记本将以第三次抖动曝光（00003）作为示例。然而，需要注意的是，在进行Spec3Pipeline之前，所有曝光必须首先通过Spec2Pipeline进行处理。
# 定义下载目录
mast_products_dir = "./mast_products/"
# 检查目录是否存在
if not os.path.exists(mast_products_dir):
# 如果目录不存在，则创建该目录
os.makedirs(mast_products_dir)
# 该笔记本专注于第三次抖动曝光。
obs_ids = ["jw01345061001_07101_00003"] # 观测ID列表
detectors = [1, 2] # 两个探测器 NRS1 和 NRS2。
# 指定计数率产品
rate_names = [] # 初始化计数率文件名列表
for obs_id in obs_ids: # 遍历每个观测ID
for detector in detectors: # 遍历每个探测器
rate_names.append(f"{obs_id}_nrs{detector}_rate.fits") # 生成并添加计数率文件名
# 下载所有的FITS文件和相关的MSA文件。
msa_names = [] # 初始化MSA文件名列表
for name in rate_names: # 遍历每个计数率文件名
print(f"Downloading {name}") # 打印正在下载的文件名
get_jwst_file( # 调用函数下载FITS文件
name, # 文件名
mast_api_token="75f915f251544013bb127b7f6f6edc80", # MAST API令牌
save_directory=mast_products_dir, # 保存目录
)
# 从下载的文件头中检索MSA文件名。
msa_name = fits.getval(mast_products_dir + name, "MSAMETFL") # 获取MSA文件名
if not os.path.isfile(msa_name): # 如果MSA文件不存在
print(f"Downloading {msa_name}") # 打印正在下载的MSA文件名
get_jwst_file( # 调用函数下载MSA文件
msa_name, # MSA文件名
mast_api_token=None, # 不使用API令牌
save_directory=mast_products_dir, # 保存目录
)
msa_names.append(msa_name) # 将MSA文件名添加到列表中
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}...") # 打印正在处理的文件名
if "nrs1" in i: # 检查文件名中是否包含"nrs1"
slit_name = "80" # 如果包含，设置slit_name为"80"
else:
slit_name = "11" # 否则，设置slit_name为"11"
Spec2Pipeline.call( # 调用Spec2Pipeline进行数据处理
mast_products_dir + i, # 输入文件路径
save_results=True, # 设置保存结果为True
steps={"nsclean": {"skip": True}, "extract_2d": {"slit_name": slit_name}}, # 指定处理步骤，跳过nsclean
output_dir=stage2_nsclean_skipped_dir, # 输出目录
)
print(f"Saved {i[:-9]}" + "cal.fits") # 打印保存的cal.fits文件名
print(f"Saved {i[:-9]}" + "x1d.fits") # 打印保存的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）：
被指定为此 MOS 观测的开放狭缝的像素。
从失败开启的 MSA 快门中产生的轨迹。
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 流水线掩膜）。
# 估计运行时间：6分钟。
start = tt.time() # 记录开始时间
for i in rate_names: # 遍历速率名称列表
print(f"Processing {i}...") # 打印正在处理的文件名
if "nrs1" in i: # 检查文件名中是否包含 "nrs1"
slit_name = "80" # 如果包含，设置缝隙名称为 "80"
else:
slit_name = "11" # 否则，设置缝隙名称为 "11"
Spec2Pipeline.call( # 调用 Spec2Pipeline 处理数据
mast_products_dir + i, # 输入文件路径
save_results=True, # 保存结果设置为 True
steps={ # 处理步骤配置
"nsclean": {"skip": False, "save_mask": True, "save_results": True}, # NSClean 步骤配置
"extract_2d": {"slit_name": slit_name}, # 提取 2D 数据时使用的缝隙名称
},
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 验证掩膜（默认管道掩膜） #
检查掩膜与速率数据，以确保它仅保留探测器的暗区。
请注意，仍然存在一些被照亮的区域，主要是由于瞬态伪影，例如宇宙射线和雪球。
# 绘制带有被遮挡区域的速率数据。
# 从管道中构建的即时掩膜列表。
nsclean_default_masks = [
stage2_nsclean_default_dir + "jw01345061001_07101_00003_nrs1_mask.fits", # NRS1掩膜文件路径
stage2_nsclean_default_dir + "jw01345061001_07101_00003_nrs2_mask.fits", # NRS2掩膜文件路径
]
# 绘制每个相关的速率数据集和掩膜文件。
for rate_file, mask_file in zip(rate_names, nsclean_default_masks): # 遍历速率文件和掩膜文件的组合
plot_dark_data(mast_products_dir + rate_file, mask_file, layout="columns", scale=9) # 调用绘图函数，传入速率文件、掩膜文件、布局和缩放参数
5.2 比较原始数据与清理后的数据（默认管道掩模） #
现在我们可以将清理后的数据（使用默认管道掩模）与原始速率文件进行比较，并验证1/f噪声是否已被降低。
在许多情况下，清理过程会向速率文件引入新的伪影。这些伪影应仔细检查，并与噪声降低的好处进行权衡。如果瞬态伪影，如雪球，干扰了清理过程，手动编辑掩模以将这些区域从背景拟合中排除可能是有益的。为此，可以尝试改变异常值检测阈值或直接编辑掩模数组中的特定像素（将在接下来的几个部分中探讨）。否则，请参考NSClean文档以获取有关手动编辑的其他建议。
请注意，在下面的图像中，有一些散布的值与原始速率文件存在较大的相对差异（在下面的相对差异图像中显示）。这些是清理过程的伪影。
还有更广泛的低水平残余背景效应（在右侧的相对差异图像中显示，散布的异常值通过σ裁剪识别，并被掩模隐藏）。这些包括我们试图去除的背景模式：在色散方向上的1/f噪声变化以及全帧数据顶部和底部的画框效应。然而，清理过程中可能还会因过拟合暗数据而引入低水平伪影。
仔细检查这两幅残余图像，以了解清理过程对数据的影响。
# 绘制原始数据和清理后的数据，以及残差图。
# 定义清理后的默认掩膜文件路径
cleaned_default_masks = [
stage2_nsclean_default_dir + "jw01345061001_07101_00003_nrs1_nsclean.fits", # NRS1清理后的文件路径
stage2_nsclean_default_dir + "jw01345061001_07101_00003_nrs2_nsclean.fits", # NRS2清理后的文件路径
]
# 遍历每一组关联的rateint数据和清理后的文件
for rate_file, cleaned_file in zip(rate_names, cleaned_default_masks):
# 绘制清理后的数据
plot_cleaned_data(
mast_products_dir + rate_file, # 原始数据文件路径
cleaned_file, # 清理后的数据文件路径
layout="columns", # 布局方式为列
scale=9 # 缩放比例设置为9
)
# 绘制 NRS2 中 x=1000, y=600 附近的重度掩蔽区域
# 掩蔽过程可能引入一些高频噪声，这些噪声在光谱轨迹上表现为垂直条纹
rate_file = mast_products_dir + rate_names[1] # NRS2 文件路径
original_rate_data = fits.open(rate_file)[1].data[550:651, 800:1401] # 原始数据，裁剪区域
cleaned_rate_data = fits.open(cleaned_default_masks[1])[1].data[550:651, 800:1401] # 清理后的数据，裁剪区域
# 为了绘图可视化处理 NaN 值
original_rate_data[np.isnan(original_rate_data)] = 0 # 将原始数据中的 NaN 替换为 0
cleaned_rate_data[np.isnan(cleaned_rate_data)] = 0 # 将清理后的数据中的 NaN 替换为 0
# 设置颜色映射的最小值和最大值
vmin = np.nanpercentile(cleaned_rate_data, 5) # 清理后数据的 5 百分位
vmax = np.nanpercentile(cleaned_rate_data, 100-9) # 清理后数据的 91 百分位
# 原始数据与清理数据的对比图
fig, axs = plt.subplots(1, 2, figsize=(15, 4)) # 创建 1 行 2 列的子图
# 绘制原始数据
fig.colorbar(
axs[0].imshow(
original_rate_data, # 原始数据
cmap="viridis", # 颜色映射
aspect=4, # 纵横比
origin="lower", # 原点位置
clim=(vmin, vmax), # 颜色范围
),
ax=axs[0], # 关联的子图
pad=0.05, # 颜色条与子图的间距
shrink=0.7, # 颜色条缩放
label="DN/s", # 颜色条标签
)
# 绘制清理后的数据
fig.colorbar(
axs[1].imshow(
cleaned_rate_data, # 清理后的数据
cmap="viridis", # 颜色映射
aspect=4, # 纵横比
origin="lower", # 原点位置
clim=(vmin, vmax), # 颜色范围
),
ax=axs[1], # 关联的子图
pad=0.05, # 颜色条与子图的间距
shrink=0.7, # 颜色条缩放
label="DN/s", # 颜色条标签
)
# 设置子图的标题、刻度值、x轴标签和y轴标签
for ax, title in zip(axs, ["Original Rate Data (NRS2)", "Cleaned Rate Data (NRS2)"]):
ax.set(title=title, xlabel="Pixel Column", ylabel="Pixel Row") # 设置标题和轴标签
ax.set_xticklabels([700, 800, 900, 1000, 1100, 1200, 1300, 1400]) # 设置 x 轴刻度标签
ax.set_yticklabels([475, 500, 525, 575, 600, 625, 650]) # 设置 y 轴刻度标签
比较从清理后的数据中提取的光谱与从原始速率文件中提取的光谱。
# 1D 提取的光谱。
# 定义跳过清理的光谱文件列表
x1d_nsclean_skipped = [
stage2_nsclean_skipped_dir + "jw01345061001_07101_00003_nrs1_x1d.fits", # NRS1光谱文件路径
stage2_nsclean_skipped_dir + "jw01345061001_07101_00003_nrs2_x1d.fits", # NRS2光谱文件路径
]
# 定义默认清理的光谱文件列表
x1d_nsclean_default = [
stage2_nsclean_default_dir + "jw01345061001_07101_00003_nrs1_x1d.fits", # NRS1光谱文件路径
stage2_nsclean_default_dir + "jw01345061001_07101_00003_nrs2_x1d.fits", # NRS2光谱文件路径
]
# 感兴趣的波长区域。
# 遍历跳过清理和默认清理的光谱文件
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_default):
plot_spectra([original, cleaned], scale_percent=9) # 绘制光谱，缩放百分比为9
注意：
在NRS1的80号狭缝和NRS2的11号狭缝中，整体连续体水平已因清理过程而略有改变。
在NRS2的11号狭缝中，原始光谱中的一些负通量已被修正。然而，默认掩膜中的过度掩膜引入了高频噪声，特别是在4.5 - 5.0微米之间影响了11号狭缝。这导致从清理后的数据提取的光谱与原始原始数据（如上所示）相比，出现了明显的下降。
6. 使用 NSClean（替代掩膜）清理 1/f 噪声 #
对于这个数据集，请注意，由于重叠的光栅区域，探测器的某些区域被严重掩膜。例如，查看上面 NRS2 的清理速率数据，在 x=1000，y=600 附近。
在这个区域，清理过程引入了一些高频噪声，表现为光谱轨迹上的垂直条纹。光栅 11 是从这个区域提取的，与原始速率数据（上面）相比，从清理后的数据中提取的光谱显示出明显的下降。
还要注意，对于 MOS 数据，探测器上可能存在几个未被光栅边界框掩膜的照明区域。对于 M 光栅，零级光谱可能出现在探测器上，并且不易定位。对于长通滤光片，仍然有一些光线穿过光栅边界框的红色截止。
在这种情况下，使用替代算法构建掩膜可能是有益的。在这里，我们不使用光栅边界框，而是迭代地掩膜任何超过背景 1 个标准差的数据。这在光谱轨迹之间留下了更多的暗数据，并改善了问题区域的背景拟合。
然而，请注意，过度清理可能会影响光谱的连续性水平，如果掩膜中包含的照明数据过多或过少。再次强调，生成的掩膜和输出光谱应仔细检查，以权衡清理的好处与对光谱的影响。
要调整此掩膜中的照明检测，请尝试修改下面的 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 管道掩膜）。
# 估计运行时间：7分钟。
start = tt.time() # 记录开始时间
for indx, i in enumerate(rate_names): # 遍历 rate_names 列表，获取索引和文件名
print(f"Processing {i}...") # 打印正在处理的文件名
if "nrs1" in i: # 判断文件名中是否包含 "nrs1"
slit_name = "80" # 如果包含，设置光谱条带名称为 "80"
else:
slit_name = "11" # 否则，设置光谱条带名称为 "11"
Spec2Pipeline.call( # 调用 Spec2Pipeline 的 call 方法
mast_products_dir + i, # 输入文件路径
save_results=True, # 保存结果
steps={ # 处理步骤
"nsclean": { # NSClean 步骤
"skip": False, # 不跳过此步骤
"save_mask": True, # 保存掩膜
"save_results": True, # 保存结果
"n_sigma": 1, # 设置 sigma 值
"mask_spectral_regions": False, # 不掩膜光谱区域
},
"extract_2d": {"slit_name": slit_name}, # 提取 2D 数据，使用指定的光谱条带名称
},
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") # 保存的 1D 文件
end = tt.time() # 记录结束时间
print("Run time: ", round(end - start, 1) / 60.0, " min") # 打印运行时间
6.1 验证掩模（备用掩模） #
检查掩模与速率数据，以确保它仅保留探测器的暗区。
# 绘制带有屏蔽区域的速率数据。
# 从管道中构建的即时掩膜列表。
nsclean_alternate_masks = [
stage2_nsclean_alternate_dir + "jw01345061001_07101_00003_nrs1_mask.fits", # NRS1掩膜文件路径
stage2_nsclean_alternate_dir + "jw01345061001_07101_00003_nrs2_mask.fits", # NRS2掩膜文件路径
]
# 绘制每个相关的速率数据集和掩膜文件。
for rate_file, mask_file in zip(rate_names, nsclean_alternate_masks): # 遍历速率文件和掩膜文件
plot_dark_data(mast_products_dir + rate_file, mask_file, layout="columns", scale=9) # 调用绘图函数
6.2 原始数据与清理数据的比较（备用掩模） #
# 绘制原始数据和清理后的数据，以及残差图。
# 定义清理后的备用掩膜文件路径列表
cleaned_alternate_masks = [
stage2_nsclean_alternate_dir + "jw01345061001_07101_00003_nrs1_nsclean.fits", # NRS1清理后的文件路径
stage2_nsclean_alternate_dir + "jw01345061001_07101_00003_nrs2_nsclean.fits", # NRS2清理后的文件路径
]
# 遍历每一组速率数据文件和清理后的文件
for rate_file, cleaned_file in zip(rate_names, cleaned_alternate_masks):
# 绘制清理后的数据
plot_cleaned_data(
mast_products_dir + rate_file, # 原始速率数据文件路径
cleaned_file, # 清理后的数据文件路径
layout="columns", # 布局设置为列
scale=9 # 缩放因子设置为9
)
比较从清理后的数据中提取的光谱与从原始速率文件中提取的光谱。
# 定义包含清理后的X1D数据文件路径的列表
x1d_nsclean_alternate = [
# 添加NRS1的清理后X1D数据文件路径
stage2_nsclean_alternate_dir + "jw01345061001_07101_00003_nrs1_x1d.fits",
# 添加NRS2的清理后X1D数据文件路径
stage2_nsclean_alternate_dir + "jw01345061001_07101_00003_nrs2_x1d.fits",
]
# 遍历原始数据和清理后的数据，进行光谱绘制
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_alternate):
# 绘制原始和清理后的光谱，设置缩放百分比为9%
plot_spectra([original, cleaned], scale_percent=9)
注意：
在NRS1的80号光缝和NRS2的11号光缝中，整体连续体水平因清理过程（类似于默认掩蔽）而略有改变。
在NRS2的11号光缝中，原始光谱中的一些负通量已被修正。默认掩蔽引入的高频噪声影响了11号光缝在4.5 - 5.0微米之间的表现（导致光谱中出现负峰），而替代掩蔽改善了这一情况。
7. 使用 NSClean 清理 1/f 噪声（手动修改的掩模） #
在某些情况下，可能需要手动生成掩模。在这里，我们介绍 一种 手动修改掩模的方法（在 NRS2 中编辑 x=1000, y=600 附近的过度掩模区域，以包含更多背景；排除 NRS1 中的一些大雪球），从管道的默认掩模输出开始。值得注意的是，使用这种方法修改的掩模不一定会优于之前的两种选项。
# 设置用于运行 NSClean 的目录，并使用用户提供的掩膜。
stage2_nsclean_modified_dir = "./stage2_nsclean_modified/" # 定义修改后的 NSClean 目录路径
if not os.path.exists(stage2_nsclean_modified_dir): # 检查目录是否存在
# 如果目录不存在，则创建该目录。
os.makedirs(stage2_nsclean_modified_dir) # 创建目录
# 手动修改某些掩膜区域。
# 特别是修改NRS2中x=1000, y=600附近的区域。
# 定义一个列表来存储修改后的掩膜路径。
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 "nrs2" in mask:
# 步骤1：将默认掩膜区域设置回True。
mask_data[550:651, :1300] = True
# 步骤2：手动重新定义掩膜区域。
mask_data[550:575, :1300] = False # NRS2中的拥挤区域。
mask_data[590:615, 150:1720] = False
mask_data[622:647, 100:1620] = False
else:
mask_data[50:130, 780:850] = False # 雪球区域
mask_data[110:140, 920:945] = False
mask_data[820:940, 2000:2040] = False
mask_data[1650:1700, 1900:1960] = False
mask_data[1800:1900, 330:410] = False
# 更新科学扩展中的数据。
hdul["SCI"].data = mask_data
# 保存修改后的FITS文件
output_path = os.path.join(stage2_nsclean_modified_dir, output_file)
hdul_modified = hdul.copy() # 复制hdul。
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 + "jw01345061001_07101_00003_nrs1_mask_modified.fits", # NRS1的修改掩膜文件路径
stage2_nsclean_modified_dir + "jw01345061001_07101_00003_nrs2_mask_modified.fits", # NRS2的修改掩膜文件路径
]
# 遍历每个相关的速率数据文件和掩膜文件进行绘图。
for rate_file, mask_file in zip(rate_names, nsclean_modified_masks):
plot_dark_data(mast_products_dir + rate_file, mask_file, layout="columns", scale=9) # 绘制速率数据，使用指定的掩膜文件和布局
注意： 在修改NRS2的默认掩膜时，我们选择了一个区域进行解掩膜，然后以不同的方式重新掩膜。然而，这个过程不小心导致了一些之前被掩膜的NaN值（在NRS2的暗数据图中看到的白色像素）被解掩膜。因此，在修改掩膜时应谨慎。此外，尽管我们掩膜了各种雪球，但根据我们提取的狭缝，1D提取光谱可能不会观察到差异。
# 1/f 噪声清理数据（用户提供的掩膜）。
# 估计运行时间：5分钟。
start = tt.time() # 记录开始时间
for indx, i in enumerate(rate_names): # 遍历速率名称列表
print(f"Processing {i}...") # 打印正在处理的文件名
if "nrs1" in i: # 判断文件名中是否包含"nrs1"
slit_name = "80" # 如果包含，设置光谱缝名称为"80"
else:
slit_name = "11" # 否则，设置光谱缝名称为"11"
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], # 使用用户提供的掩膜
},
"extract_2d": {"slit_name": slit_name}, # 提取2D数据，使用指定的光谱缝名称
},
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") # 保存的1D提取文件
end = tt.time() # 记录结束时间
print("Run time: ", round(end - start, 1) / 60.0, " min") # 打印运行时间
7.2 原始数据与清理数据（手动修改的掩模）比较 #
# 绘制原始数据和清理后的数据，以及残差图。
# 定义清理后的文件路径列表
cleaned_modified_masks = [
stage2_nsclean_modified_dir + "jw01345061001_07101_00003_nrs1_nsclean.fits", # NRS1清理后的文件路径
stage2_nsclean_modified_dir + "jw01345061001_07101_00003_nrs2_nsclean.fits", # NRS2清理后的文件路径
]
# 遍历每个关联的rateint数据和清理后的文件
for rate_file, cleaned_file in zip(rate_names, cleaned_modified_masks):
# 绘制清理后的数据
plot_cleaned_data(
mast_products_dir + rate_file, # 原始数据文件路径
cleaned_file, # 清理后的数据文件路径
layout="columns", # 布局方式为列
scale=9 # 缩放比例设置为9
)
比较从清理后的数据中提取的光谱与从原始速率文件中提取的光谱。
# 定义修改后的 x1d_nsclean 文件路径列表
x1d_nsclean_modified = [
# 添加第一个文件的路径
stage2_nsclean_modified_dir + "jw01345061001_07101_00003_nrs1_x1d.fits",
# 添加第二个文件的路径
stage2_nsclean_modified_dir + "jw01345061001_07101_00003_nrs2_x1d.fits",
]
# 感兴趣的波长区域
# 遍历原始和清理后的光谱文件
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_modified):
# 绘制原始和清理后的光谱，缩放百分比设置为9%
plot_spectra([original, cleaned], scale_percent=9)
注意事项:
在NRS1的狭缝80上，整体连续谱水平因清理过程而略有改变（类似于默认掩模和替代掩模）。
在NRS2的狭缝11上，由于使用手动修改的掩模进行清理过程，4-4.55微米之间的通量有所下降。原始光谱中（在更短和更长波长处）的一些负通量已被修正。
8. 结论 #
下面的最终图表展示了计数率图像和相应的1D提取光谱并排比较不同的清理方法：原始图像（未应用NSClean）、清理后的计数率图像（使用默认的管道掩模）、清理后的计数率图像（使用替代的管道掩模），以及最后，清理后的计数率图像（使用手动修改的掩模）。
请注意，本笔记本中呈现的结果可能会因不同的数据集而有所不同（例如，不同亮度、空间范围等的目标）。鼓励用户使用不同的掩模方法探索NSClean，以确定最佳结果。
清理算法的输出现在已准备好进行进一步处理。由上述Spec2Pipeline运行生成的(_cal.fits)文件可用作Spec3Pipeline的输入，以生成最终的组合光谱。
# 未清理数据 vs. 清理数据（默认掩膜） vs. 清理数据（备用掩膜）速率数据
# 读取原始速率数据
original_rate_data = [
fits.open(mast_products_dir + rate_name)[1].data for rate_name in rate_names # 从文件中打开并提取数据
]
# 读取清理后的默认掩膜速率数据
cleaned_rate_default_data = [
fits.open(cleaned_default_mask)[1].data # 从清理后的默认掩膜文件中打开并提取数据
for cleaned_default_mask in cleaned_default_masks
]
# 读取清理后的备用掩膜速率数据
cleaned_rate_alternate_data = [
fits.open(cleaned_alternate_mask)[1].data # 从清理后的备用掩膜文件中打开并提取数据
for cleaned_alternate_mask in cleaned_alternate_masks
]
# 读取清理后的手动修改掩膜速率数据
cleaned_rate_modified_data = [
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:
data[np.isnan(data)] = 0 # 将NaN值替换为0
# 原始数据与清理数据（使用默认掩膜）进行比较
fig, axs = plt.subplots(2, 4, figsize=(25, 12)) # 创建子图
# 设置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] # 获取当前子图
ax.set_title(f'{title} \n {"NRS1" if j == 0 else "NRS2"}', fontsize=12) # 设置标题
im = ax.imshow(data, origin="lower", clim=(-1e-2, 1e-2)) # 绘制图像
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], # 修改后的清理数据的第一组光谱
],
scale_percent=4, # 设置缩放百分比为4%
)
# 绘制第二组光谱数据
plot_spectra(
[
x1d_nsclean_skipped[1], # 跳过的清理数据的第二组光谱
x1d_nsclean_default[1], # 默认清理数据的第二组光谱
x1d_nsclean_alternate[1], # 替代清理数据的第二组光谱
x1d_nsclean_modified[1], # 修改后的清理数据的第二组光谱
],
scale_percent=4, # 设置缩放百分比为4%
)
最终说明：
默认掩膜引入的NRS2高频噪声影响了在4.55微米附近的11号狭缝的1D提取光谱，但使用替代掩膜后该噪声不再出现（手动修改的掩膜略有改善），显著改善了11号狭缝的光谱。
11号狭缝的负通量（在短波长和长波长处）已通过任一掩膜进行了修正。然而，手动修改的掩膜似乎降低了11号狭缝在4-4.55微米之间的光谱连续性水平。在这种情况下，替代掩膜（基于剪辑的）算法优于对每个MSA快门阻塞整个科学区域。
关于笔记本 #
作者： Melanie Clarke, Kayli Glidic; NIRSpec仪器团队
更新时间： 2024年2月29日。