使用 NSClean 清理 NIRSpec IFU 产品中的残余 1/f 噪声#
笔记本目标#
本笔记本的目标是通过去除残余的 1/f 噪声来生成清理后的 IFU (_rate.fits) 文件。这些清理后的文件将作为第 3 级 (Spec3Pipeline) 管道的输入。
目录#
1. 介绍 #
JWST NIRSpec 仪器有许多特性和特点，观察者在规划观测和解释数据时应予以注意。NIRSpec 管道产品中一个显著的特征是提取的 1-D 光谱中出现负值和/或过剩的通量，通常伴有不规则的波长依赖性起伏。这种伪影的原因是相关噪声，称为 1/f 噪声，源于低水平的探测器热不稳定性，表现为 2-D 计数率图像中的垂直条带，尤其是在 NRS2 探测器的曝光中。虽然 IRS2 读出模式减少了这种效应，但并未完全消除。
为了解决这个问题，JWST 科学校准管道集成了由 Bernard Rauscher 开发的外部包，称为 NSClean，并在 Spec2Pipeline 中通过 NSCleanStep 使用。该算法利用探测器的暗区在傅里叶空间中拟合背景模型。它需要一个输入掩模来识别探测器的所有暗区。掩模越全面，背景拟合效果越好。
在本笔记本中，我们将使用集成到管道中的 NSClean 算法，利用默认参数动态生成的掩模来去除 1/f 噪声。在某些情况下，该掩模可能不够完整或过于严格，无法实现最佳的噪声去除。为了解决这个问题，我们展示了如何手动修改默认掩模，以及如何通过调整 NSCleanStep 参数 创建替代掩模。如有需要，请参见 NSClean 文档 以获取有关手动创建自定义掩模的一些建议。
本笔记本利用了对类星体 XID2028 的 IFU 观测，使用光栅/滤光片 G140H/F100LP，作为 JWST 早期发布科学计划 ERS-1335 第 4 次观测的示例。
2. 导入库 #
# ---------- 设置CRDS环境变量 ----------
import os # 导入操作系统模块
import 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模块导入自定义函数，用于文件获取和数据绘图
# Hide all log and warning messages. # 隐藏所有日志和警告信息
logging.disable(logging.ERROR) # 禁用错误级别的日志信息
warnings.simplefilter("ignore", RuntimeWarning) # 忽略运行时警告
3. 下载数据 #
输入数据来自于本笔记本，包含了对类星体 XID2028 的 IFU 观测，使用了光栅/滤光片 G140H/F100LP。该数据集是 JWST 早期发布科学计划 ERS-1335 的一部分，具体为观测 4。它由 9 次集成（9 个抖动点）组成，每次集成包含 16 组数据。本笔记本将以第二次抖动曝光（00002）作为示例。然而，需要注意的是，在进入 Spec3Pipeline 之前，所有曝光必须首先通过 Spec2Pipeline 进行处理。
# 定义下载目录
mast_products_dir = "./mast_products/"
# 检查目录是否存在
if not os.path.exists(mast_products_dir):
# 如果目录不存在，则创建该目录
os.makedirs(mast_products_dir)
# 该笔记本专注于第二个抖动曝光。
obs_ids = ["jw01335004001_03101_00002"] # 观测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 文件。
for name in rate_names: # 遍历每个文件名
print(f"Downloading {name}") # 打印正在下载的文件名
get_jwst_file(name, mast_api_token=None, save_directory=mast_products_dir) # 下载指定的FITS文件
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）。
# 预计运行时间：132分钟（可能会有所不同）。
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 # 跳过此步骤，不去除NIRSpec图像中的相关读噪声（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）：
被指定为 IFU 数据的科学区域的像素。
从故障开启的 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 管道掩模）。
# 估计运行时间：105 分钟（可能会有所不同）。
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") # 打印保存的 nsclean 文件名
print(f"Saved {i[:-9]}" + "cal.fits") # 打印保存的校准文件名
print(f"Saved {i[:-9]}" + "x1d.fits") # 打印保存的 x1d 文件名
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 + "jw01335004001_03101_00002_nrs1_mask.fits", # NRS1掩膜文件路径
stage2_nsclean_default_dir + "jw01335004001_03101_00002_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文档以获取有关手动编辑的其他建议。
请注意，在下面的图像中，有一些值与原始速率文件存在较大的相对差异（在下面的相对差异图像中显示）。这些是清理过程的伪影。
此外，还有更广泛的低水平残余背景效应（在右侧的相对差异图像中显示，散布的异常值通过sigma裁剪识别，并被掩模隐藏）。这些包括我们试图去除的背景模式：色散方向上的1/f噪声变化以及全帧数据顶部和底部的画框效应。然而，清理过程中也可能因过拟合暗数据而引入低水平伪影。
仔细检查这两幅残余图像，以了解清理过程对数据的影响。
# 绘制原始数据和清理后的数据，以及残差图。
cleaned_default_masks = [
stage2_nsclean_default_dir + "jw01335004001_03101_00002_nrs1_nsclean.fits", # 第一幅清理后的数据文件路径
stage2_nsclean_default_dir + "jw01335004001_03101_00002_nrs2_nsclean.fits", # 第二幅清理后的数据文件路径
]
# 遍历每一组关联的rateint数据和清理后的文件。
for rate_file, cleaned_file in zip(rate_names, cleaned_default_masks): # 将rate_names和cleaned_default_masks进行配对
plot_cleaned_data( # 调用绘图函数
mast_products_dir + rate_file, cleaned_file, layout="columns", scale=9 # 绘制原始数据和清理后的数据，设置布局和缩放比例
)
比较从清理后的数据中提取的光谱与从原始速率文件中提取的光谱。
# 1D 提取的光谱数据。
x1d_nsclean_skipped = [
# 路径到未清理的光谱数据文件（nrs1）
stage2_nsclean_skipped_dir + "jw01335004001_03101_00002_nrs1_x1d.fits",
# 路径到未清理的光谱数据文件（nrs2）
stage2_nsclean_skipped_dir + "jw01335004001_03101_00002_nrs2_x1d.fits",
]
x1d_nsclean_default = [
# 路径到默认清理的光谱数据文件（nrs1）
stage2_nsclean_default_dir + "jw01335004001_03101_00002_nrs1_x1d.fits",
# 路径到默认清理的光谱数据文件（nrs2）
stage2_nsclean_default_dir + "jw01335004001_03101_00002_nrs2_x1d.fits",
]
# 感兴趣的波长范围。
wavelength_range = {"nrs1": [1.15, 1.25], "nrs2": [1.65, 1.75]}
# 遍历未清理和清理的光谱文件
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_default):
# 绘制光谱
plot_spectra(
[original, cleaned], # 输入的光谱文件列表
scale_percent=4, # 设定缩放百分比
wavelength_range=wavelength_range # 设定波长范围
)
注意：
NRS1在1.2微米附近和NRS2在1.7微米附近的光谱部分，由于1/f噪声，过量的通量得到了减少。
清理后的光谱与原始光谱之间存在几个尖峰。这些尖峰对应于清理过程中引入的散射伪影，如上所述。
6. 使用 NSClean（替代掩膜）清除 1/f 噪声 #
对于某些数据集，掩膜整个科学区域可能会过度掩盖探测器的暗区，而这些区域可以用来改善背景拟合。过度掩膜可能会在清理过程中引入一些高频噪声，这些噪声表现为光谱轨迹上的垂直条纹。此外，对于某些数据集，探测器上可能存在多个未被 IFU 切片边界框掩膜的照明区域。这可能包括瞬态伪影，如宇宙射线或来自饱和源的光晕。
在某些情况下，使用替代算法构建掩膜可能更为有利。在这里，我们不使用边界框，而是迭代地掩膜任何超过背景 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 管道掩膜）。
# 估计运行时间：87分钟（可能会有所不同）。
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, # 保存结果设置为 True
steps={ # 定义处理步骤
"nsclean": { # NSClean 步骤
"skip": False, # 不跳过此步骤
"save_mask": True, # 保存掩膜设置为 True
"n_sigma": 1, # 设置 sigma 值为 1
"mask_spectral_regions": False, # 不掩膜光谱区域
"save_results": True, # 保存结果设置为 True
}, # 从 NIRSpec 图像中去除相关的读噪声（1/f 噪声）。
},
output_dir=stage2_nsclean_alternate_dir, # 输出目录
)
print(f"Saved {i[:-9]}" + "mask.fits") # 打印保存的掩膜文件名
print(f"Saved {i[:-9]}" + "nsclean.fits") # 打印保存的 nsclean 文件名
print(f"Saved {i[:-9]}" + "cal.fits") # 打印保存的校准文件名
print(f"Saved {i[:-9]}" + "x1d.fits") # 打印保存的 x1d 文件名
end = tt.time() # 记录结束时间
print("Run time: ", round(end - start, 1) / 60.0, " min") # 打印运行时间
6.1 验证掩模（备用掩模） #
检查掩模与速率数据，以确保它仅保留探测器的暗区。
# 绘制被遮挡区域的速率数据。
# 从管道中构建的按需掩膜列表。
nsclean_alternate_masks = [
stage2_nsclean_alternate_dir + "jw01335004001_03101_00002_nrs1_mask.fits", # NRS1掩膜文件路径
stage2_nsclean_alternate_dir + "jw01335004001_03101_00002_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 + "jw01335004001_03101_00002_nrs1_nsclean.fits", # NRS1清理后的文件路径
stage2_nsclean_alternate_dir + "jw01335004001_03101_00002_nrs2_nsclean.fits", # NRS2清理后的文件路径
]
# 遍历每一组相关的rateint数据和清理后的文件
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 + "jw01335004001_03101_00002_nrs1_x1d.fits",
# 添加NRS2的X1D文件路径
stage2_nsclean_alternate_dir + "jw01335004001_03101_00002_nrs2_x1d.fits",
]
# 定义感兴趣的波长范围
wavelength_range = {"nrs1": [1.15, 1.25], "nrs2": [1.65, 1.75]}
# 遍历原始和清理后的光谱数据
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_alternate):
# 绘制光谱图
plot_spectra(
# 传入原始和清理后的光谱数据
[original, cleaned], scale_percent=4, wavelength_range=wavelength_range
)
注意：
NRS1光谱的整体连续体水平在清理后发生了变化。这是因为sigma-clipping方法与默认方法不同，仅掩蔽明亮的异常值，而不是整个科学区域。因此，一些天空区域被包含在从数据中减去的背景模型中，这可能导致连续体水平的变化。
NRS2光谱中接近1.7微米的部分，由于1/f噪声导致的过量通量已被减少。
清理后的光谱与原始光谱之间的差异中存在几个尖峰。这些尖峰对应于清理过程中引入的散射伪影，如上所述。
7. 使用 NSClean 清理 1/f 噪声（手动修改的掩模） #
在某些情况下，可能需要手动生成掩模。在这里，我们介绍 一种 手动修改掩模的方法（排除 NRS1/NRS2 中的一些大雪球），该方法以管道输出的默认掩模为起点。值得注意的是，使用此方法修改的掩模不一定会优于之前的两种选项。
# 设置用于运行 NSClean 的目录，并使用用户提供的掩模。
stage2_nsclean_modified_dir = "./stage2_nsclean_modified/"
if not os.path.exists(stage2_nsclean_modified_dir): # 检查目录是否存在
# 如果目录不存在，则创建该目录。
os.makedirs(stage2_nsclean_modified_dir) # 创建目录
# 手动修改某些掩膜区域
# NRS1/2中的雪球区域。
# 定义一个列表来存储修改后的掩膜路径。
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 "nrs1" in mask:
# 针对NRS1掩膜，遮罩雪球区域。
mask_data[10:50, 180:215] = False
mask_data[560:580, 450:465] = False
mask_data[1720:1810, 1150:1240] = False
mask_data[1540:1580, 535:580] = False
mask_data[1865:1910, 725:770] = False
else:
# 针对NRS2掩膜，遮罩雪球区域。
mask_data[200:250, 130:190] = False
mask_data[1400:1420, 40:80] = False
mask_data[1700:1730, 150:180] = False
mask_data[630:710, 1100:1170] = False
mask_data[520:570, 1770:1820] = 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 + "jw01335004001_03101_00002_nrs1_mask_modified.fits", # NRS1的修改掩膜文件路径
stage2_nsclean_modified_dir + "jw01335004001_03101_00002_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) # 绘制速率数据和掩膜
# 1/f 噪声清理数据（用户提供的掩膜）。
# 估计运行时间：87分钟（可能会有所不同）。
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") # 打印保存的nsclean文件名
print(f"Saved {i[:-9]}" + "cal.fits") # 打印保存的cal文件名
print(f"Saved {i[:-9]}" + "x1d.fits") # 打印保存的x1d文件名
end = tt.time() # 记录结束时间
print("Run time: ", round(end - start, 1) / 60.0, " min") # 打印运行时间
7.2 原始数据与清理数据（手动修改的掩膜）比较 #
# 绘制原始数据和清理后的数据，以及残差图。
cleaned_modified_masks = [
stage2_nsclean_modified_dir + "jw01335004001_03101_00002_nrs1_nsclean.fits", # 第一组清理后的数据文件路径
stage2_nsclean_modified_dir + "jw01335004001_03101_00002_nrs2_nsclean.fits", # 第二组清理后的数据文件路径
]
# 遍历每一对关联的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 数据文件列表
x1d_nsclean_modified = [
# 添加第一幅图像的文件路径
stage2_nsclean_modified_dir + "jw01335004001_03101_00002_nrs1_x1d.fits",
# 添加第二幅图像的文件路径
stage2_nsclean_modified_dir + "jw01335004001_03101_00002_nrs2_x1d.fits",
]
# 定义感兴趣的波长范围
wavelength_range = {"nrs1": [1.15, 1.25], "nrs2": [1.65, 1.75]}
# 遍历原始数据和清理后的数据
for original, cleaned in zip(x1d_nsclean_skipped, x1d_nsclean_modified):
# 绘制光谱图
plot_spectra(
# 传入原始和清理后的数据
[original, cleaned], scale_percent=4, wavelength_range=wavelength_range
)
注意：
即使在遮蔽我们遇到的少量雪球时，光谱仍然与使用默认遮罩清理的光谱相似。
对于NRS1，接近1.2微米的光谱部分，以及对于NRS2，接近1.7微米的部分，由于1/f噪声导致的过量通量有所减少。
清理后的光谱与原始光谱之间存在几个尖峰。这些尖峰对应于清理过程中引入的散射伪影，如上所述。
8. 结论 #
下面的最终图表展示了计数率图像和相应的1D提取光谱并排比较不同的清理方法：原始图像（未应用NSClean）、清理后的计数率图像（使用默认的管道掩模）、清理后的计数率图像（使用替代的管道掩模），以及最后，清理后的计数率图像（使用手动修改的掩模）。
请注意，本笔记本中呈现的结果可能因不同数据集而异（例如，不同亮度、空间范围等的目标）。鼓励用户使用不同的掩模方法探索NSClean，以确定最佳结果。
清理算法的输出现在已准备好进行进一步处理。上述Spec2Pipeline运行生成的(_cal.fits)文件可作为Spec3Pipeline的输入，用于生成最终的组合光谱立方体和提取光谱。
# 未清理与清理后的（默认掩膜）与清理后的（备用掩膜）速率数据
# 读取原始速率数据
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-3, 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() # 显示图形
# 最终比较
# 感兴趣的波长范围
wavelength_range = {"nrs1": [1.15, 1.25], "nrs2": [1.65, 1.75]} # 定义波长范围字典
# 绘制第一个光谱
plot_spectra(
[
x1d_nsclean_skipped[0], # 跳过的清理数据的第一个光谱
x1d_nsclean_default[0], # 默认清理数据的第一个光谱
x1d_nsclean_alternate[0], # 替代清理数据的第一个光谱
x1d_nsclean_modified[0], # 修改后的清理数据的第一个光谱
],
wavelength_range=wavelength_range, # 使用定义的波长范围
scale_percent=4, # 设置缩放百分比
)
# 绘制第二个光谱
plot_spectra(
[
x1d_nsclean_skipped[1], # 跳过的清理数据的第二个光谱
x1d_nsclean_default[1], # 默认清理数据的第二个光谱
x1d_nsclean_alternate[1], # 替代清理数据的第二个光谱
x1d_nsclean_modified[1], # 修改后的清理数据的第二个光谱
],
wavelength_range=wavelength_range, # 使用定义的波长范围
scale_percent=4, # 设置缩放百分比
)
注意： 使用替代掩模进行清理仍然会去除一些波长依赖的变化，但在NRS2中留下了残余的背景变化。还要注意，光谱的整体连续体水平发生了变化，尤其是对于NRS1。这是因为sigma-clipping（σ裁剪）与默认方法不同，仅对亮的异常值进行掩蔽，而不是整个科学区域。因此，一些天空区域被包含在从数据中减去的背景模型中，这可能导致连续体水平的变化。在这种情况下，原始算法（对每个IFU切片屏蔽整个科学区域）或手动修改的方法（排除雪球）比通过sigma-clipping创建掩模更为可取。
关于笔记本 #
作者： Melanie Clarke, Kayli Glidic; NIRSpec仪器团队
更新时间： 2024年2月29日。