运行 spec2 管道#
WFSS 的分散图像将在运行完相应直接图像的 image2 和 image3 管道 后，通过 spec2 管道进行处理。如前一个笔记本所提到的，如果您有一个仅包含您希望查看的源的源目录，这一步将极为有帮助。任何额外的源都会增加校准所需的时间。
用例：在创建自定义源目录后，应在分散的 WFSS 图像上运行 spec2。
数据：JWST/NIRISS 图像和来自程序 2079 观测 004 的光谱。这些数据应存储在一个名为 data 的单一目录中，并可以从笔记本 00_niriss_mast_query_data_setup.ipynb 下载。
工具：astropy, crds, glob, jdaviz, json, jwst, matplotlib, numpy, os, pandas, shutil, urllib, zipfile
跨仪器：NIRISS
内容
作者：Rachel Plesha (rplesha@stsci.edu), Camilla Pacifici (cpacifici@stsci.edu), JWebbinar 笔记本。
首次发布：2024年5月
最后编辑：2024年8月
最后测试：此笔记本最后一次测试是在 JWST 管道版本 1.12.5 和 CRDS 上下文 jwst_1225.pmap。
导入与数据设置#
CRDS 文档#
在这个链接中，您可以找到关于CRDS（Calibration Reference Data System，校准参考数据系统）的详细信息。CRDS是JWST（James Webb Space Telescope，詹姆斯·韦伯太空望远镜）数据处理的重要组成部分，负责管理和分发用于数据校准的参考文件。
主要内容#
CRDS的功能：CRDS的主要功能是确保JWST数据处理所需的所有参考文件都能被有效管理和访问。
参考文件：这些文件包括各种校准数据，如响应曲线、噪声模型等。
数据版本控制：CRDS还提供版本控制，以确保使用的参考文件是最新的，并与科学数据的处理版本相匹配。
请访问上述链接以获取更详细的信息和指导。
# 更新CRDS路径到你的本地目录
%env CRDS_PATH=crds_cache # 设置CRDS缓存路径为'crds_cache'
# 设置CRDS服务器的URL
%env CRDS_SERVER_URL=https://jwst-crds.stsci.edu # 设置CRDS服务器的URL为JWST CRDS服务器
import glob # 导入glob模块，用于文件路径操作
import json # 导入json模块，用于处理JSON数据
import os # 导入os模块，用于操作系统功能
import shutil # 导入shutil模块，用于文件和目录的高级操作
import urllib # 导入urllib模块，用于处理URL
import zipfile # 导入zipfile模块，用于处理ZIP文件
import numpy as np # 导入numpy库，用于数值计算
import pandas as pd # 导入pandas库，用于数据处理和分析
from astropy.table import Table # 从astropy.table导入Table类，用于处理表格数据
from astropy.io import fits # 从astropy.io导入fits模块，用于处理FITS文件
from astropy import constants as const # 从astropy导入常量模块，用于天文常数
from matplotlib import pyplot as plt # 从matplotlib导入pyplot模块，用于绘图
from matplotlib import patches # 从matplotlib导入patches模块，用于绘图中的形状
%matplotlib widget # 设置matplotlib为交互式模式
# %matplotlib inline # 设置matplotlib为内联模式（注释掉的代码）
from jwst.pipeline import Spec2Pipeline # 从jwst.pipeline导入Spec2Pipeline类，用于JWST数据处理
检查您正在使用的JWST（詹姆斯·韦伯太空望远镜）管道版本。要查看可用的最新管道版本或安装以前的版本，请访问 GitHub。同时，请确认您使用的 CRDS（参考数据系统）版本。 CRDS文档 解释了如何在JWST管道中设置特定的上下文。如果这两个值与上面最后测试的说明不同，可能会导致此笔记本中的差异。
import jwst # 导入JWST相关的库
import crds # 导入CRDS（Calibration Reference Data System）库
print('JWST Pipeliene Version:', jwst.__version__) # 打印JWST管道的版本信息
print('CRDS Context:', crds.get_context_name('jwst')) # 获取并打印JWST的CRDS上下文名称
数据设置#
数据目录 data_dir 应该包含所有的关联文件和速率文件，并且放在一个单一的平面目录中。custom_run_image3 应该与笔记本 01_niriss_wfss_image2_image3.ipynb 中的自定义 image3 校准输出目录相匹配。
对于 spec2，我们需要下载的速率文件，以及来自 image3 的分割图和源目录。因此，我们将创建一个新目录（由 custom_run_spec2 定义），切换到该目录，并复制相关数据。
对于常规工作流程，将所有文件保留在一个单一目录中通常比我们在这里创建的多个目录更为简单。
打开我们的数据信息文件，这将使解析速率文件以便复制变得更容易，因为我们只需要分散图像，即那些使用GR150R或GR150C观测的图像。
import os # 导入os模块，用于处理文件和目录
data_dir = 'data' # 定义数据目录名称
# 检查数据目录是否存在
if not os.path.exists(data_dir):
os.mkdir(data_dir) # 如果不存在，则创建数据目录
custom_run_spec2 = 'custom_spec2' # 定义自定义spec2结果保存目录
custom_run_image3 = 'custom_image3_calibrated' # 定义自定义image3校准结果保存目录
# 如果自定义目录不存在，则创建这些目录
for custom_dir in [custom_run_spec2, custom_run_image3]:
if not os.path.exists(os.path.join(data_dir, custom_dir)):
os.mkdir(os.path.join(data_dir, custom_dir)) # 创建自定义目录
# 如果您尚未从笔记本00下载数据或尚未运行笔记本01，请运行此单元格。否则，可以跳过它。
# 从Box下载未校准的数据到数据目录：
boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/niriss_wfss_advanced/niriss_wfss_advanced_02_input.zip' # Box链接
boxfile = os.path.basename(boxlink) # 获取文件名
urllib.request.urlretrieve(boxlink, boxfile) # 下载文件
zf = zipfile.ZipFile(boxfile, 'r') # 打开zip文件
zf.extractall(path=data_dir) # 解压到指定的数据目录
# 将从box文件下载的文件移动到顶层数据目录
box_download_dir = os.path.join(data_dir, boxfile.split('.zip')[0]) # 构建下载目录路径
# 遍历下载目录中的所有文件
for filename in glob.glob(os.path.join(box_download_dir, '*')): # 获取所有文件名
if '.csv' in filename: # 如果是CSV文件
# 移动到当前目录
os.rename(filename, os.path.basename(filename)) # 重命名并移动文件到当前目录
elif '_segm.fits' in filename or '_cat.ecsv' in filename: # 如果是分割图像或分类文件
# 将image2产品移动到适当的目录
os.rename(filename, os.path.join(data_dir, custom_run_spec2, os.path.basename(filename))) # 移动到指定目录
elif '_i2d.fits' in filename: # 如果是图像3产品
# 将image3产品也移动到它们的目录
os.rename(filename, os.path.join(data_dir, custom_run_image3, os.path.basename(filename))) # 移动到指定目录
else: # 其他文件
# 移动到数据目录
os.rename(filename, os.path.join(data_dir, os.path.basename(filename))) # 移动到数据目录
# 现在删除不必要的文件
os.remove(boxfile) # 删除下载的zip文件
os.rmdir(box_download_dir) # 删除下载目录
# 从之前创建的csv文件中，找到我们想要用spec2进行校准的所有光栅观测的列表
listrate_file = 'list_ngdeep_rate.csv' # 定义包含光栅观测数据的CSV文件名
rate_df = pd.read_csv(listrate_file) # 读取CSV文件并将其存储为DataFrame
# 复制所有的光栅率文件
# 获取GR150R滤镜的文件名列表
gr150r_files = list(rate_df[rate_df['FILTER'] == 'GR150R']['FILENAME'])
# 获取GR150C滤镜的文件名列表
gr150c_files = list(rate_df[rate_df['FILTER'] == 'GR150C']['FILENAME'])
# 遍历GR150R和GR150C的文件列表
for grism_rate in gr150r_files + gr150c_files:
# 检查文件是否存在
if os.path.exists(grism_rate):
# 复制文件到指定目录
shutil.copy(grism_rate, os.path.join(data_dir, custom_run_spec2, os.path.basename(grism_rate)))
else:
# 如果文件不存在，打印错误信息
print(f'{grism_rate} does not exist. Not able to copy file.')
# 复制所有的 spec2 asn 文件
for asn in glob.glob(os.path.join(data_dir, '*spec2*asn*.json')): # 遍历匹配的 spec2 asn 文件
if os.path.exists(asn): # 检查文件是否存在
shutil.copy(asn, os.path.join(data_dir, custom_run_spec2, os.path.basename(asn))) # 复制文件到指定目录
else: # 如果文件不存在
print(f'{asn} does not exist. Not able to copy file.') # 输出文件不存在的提示信息
# 复制所有必要的image3输出文件
cats = glob.glob(os.path.join(data_dir, custom_run_image3, '*source*_cat.ecsv')) # 复制source-match和source118目录
segm = glob.glob(os.path.join(data_dir, custom_run_image3, '*_segm.fits')) # 获取所有分割图像文件
# 遍历所有的image3文件，包括目录和分割文件
for image3_file in cats + segm:
if os.path.exists(image3_file): # 检查文件是否存在
# 复制文件到指定的目标目录
shutil.copy(image3_file, os.path.join(data_dir, custom_run_spec2, os.path.basename(image3_file)))
else:
# 如果文件不存在，打印错误信息
print(f'{image3_file} does not exist. Not able to copy file.')
cwd = os.getcwd() # 获取当前工作目录
if cwd != data_dir: # 如果当前目录不是数据目录，则切换到数据目录
try:
os.chdir(data_dir) # 尝试切换到数据目录
except FileNotFoundError: # 如果找不到指定的目录
print(f"Not able to change into: {data_dir}.\nRemaining in: {cwd}") # 打印错误信息，显示无法切换目录
pass # 继续执行后面的代码
if not os.path.exists(custom_run_spec2): # 如果自定义运行目录不存在
os.mkdir(custom_run_spec2) # 创建自定义运行目录
os.chdir(custom_run_spec2) # 切换到自定义运行目录
new_cwd = os.getcwd() # 获取新的当前工作目录
print('Now in:', new_cwd) # 打印当前工作目录
# 为后续校准创建目录
source_outdir = 'new_catalog_calibrated' # 定义校准后源目录名称
if not os.path.exists(source_outdir): # 检查目录是否存在
os.mkdir(source_outdir) # 如果不存在，则创建该目录
param_outdir = 'parameter_input_calibrated' # 定义校准后参数目录名称
if not os.path.exists(param_outdir): # 检查目录是否存在
os.mkdir(param_outdir) # 如果不存在，则创建该目录
自定义 Spec2 管道运行#
由于这是一个自定义的 Spec2 管道运行，第一步是确保使用正确的源目录（source catalog）。这将定义光谱轨迹切割（spectral trace cutouts），从而告知管道在哪里提取光谱。
Spec2 关联文件#
与成像部分的处理流程一样，spec2 也有关联文件。这些文件稍微复杂一些，因为它们需要包含科学数据（WFSS）、直接图像、源目录和分割图作为成员。对于科学数据，速率文件被用作输入，类似于 image2。与 image2 一样，每个观察序列中的每个分散图像抖动位置应该有一个 asn 文件。在这种情况下，这应该与 FILTER=GR150R 或 FILTER=GR150C 的速率文件数量相匹配。对于这个程序和观测，每个光栅都有三个抖动，并且同时使用 GR150R 和 GR150C，总共在每个观察序列中有六个曝光，观察中使用了五个观察序列，使用的阻挡滤光片为 F115W -> F115W -> F150W -> F150W -> F200W。
查看Spec2关联文件#
import numpy as np # 导入NumPy库
import glob # 导入glob库
# 获取当前目录下所有包含'spec2'和'asn'的JSON文件，并按字母顺序排序
spec2_asns = np.sort(glob.glob('*spec2*asn*.json'))
# 打印找到的Spec2 ASN文件的数量
print(len(spec2_asns), 'Spec2 ASN files')
# 检查ASN文件的数量是否与分散图像的数量匹配 -- 对于观测004，应该是30个
print(len(rate_df[(rate_df['FILTER'] == 'GR150R') | (rate_df['FILTER'] == 'GR150C')]), 'Dispersed image rate files')
# 查看其中一个关联文件
# 从指定的关联文件中加载JSON数据
asn_data = json.load(open(spec2_asns[0]))
# 遍历关联数据中的每一个键值对
for key, data in asn_data.items():
# 打印键及其对应的数据
print(f"{key} : {data}")
# 特别是，仔细查看关联文件中的产品文件名：
for product in asn_data['products']: # 遍历关联数据中的每个产品
for key, value in product.items(): # 遍历产品的每个键值对
if key == 'members': # 如果键是'members'
print(f"{key}:") # 打印键名
for member in value: # 遍历成员列表
print(f" {member['expname']} : {member['exptype']}") # 打印每个成员的文件名和类型
else: # 如果键不是'members'
print(f"{key}: {value}") # 打印键名和对应的值
修改关联文件以使用自定义源目录#
当前关联文件中使用的是管道源目录。在之前的笔记本中，我们创建了一个自定义源目录，文件扩展名为 source-match_cat.ecsv，因此我们需要在关联文件中指向该目录。
new_source_ext = 'source-match_cat.ecsv' # 新的源目录扩展名
# loop through all of the spec2 association files in the current directory
for asn in spec2_asns: # 遍历当前目录中的所有spec2关联文件
asn_data = json.load(open(asn)) # 读取关联文件的数据
for product in asn_data['products']: # 对于每个产品，检查其成员
for member in product['members']: # 每个产品有多个成员
if member['exptype'] == 'sourcecat': # 检查成员类型是否为sourcecat
cat_in_asn = member['expname'] # 获取当前成员的文件名
# check that we haven't already updated the source catalog name
if new_source_ext not in cat_in_asn: # 检查是否已经更新过源目录名
new_cat = cat_in_asn.replace('cat.ecsv', new_source_ext) # 替换为新的目录名
# actually update the association file member
if os.path.exists(new_cat): # 检查新的目录文件是否存在
member['expname'] = new_cat # 更新成员的文件名
else:
print(f"{new_cat} does not exist in the current working directory") # 输出文件不存在的提示
# write out the new association file
with open(asn, 'w', encoding='utf-8') as f: # 以写入模式打开关联文件
json.dump(asn_data, f, ensure_ascii=False, indent=4) # 将更新后的数据写入文件
# 双重检查数据是否仍然正常：
asn_check = json.load(open(spec2_asns[0])) # 从第一个ASN文件中加载JSON数据
for product in asn_check['products']: # 遍历产品列表
for key, value in product.items(): # 遍历每个产品的键值对
if key == 'members': # 如果键是'members'
print(f"{key}:") # 打印键名
for member in value: # 遍历成员列表
print(f" {member['expname']} : {member['exptype']}") # 打印每个成员的名称和类型
else: # 如果键不是'members'
print(f"{key}: {value}") # 打印键名和对应的值
或者，在您喜欢的文本编辑器中打开 .json 文件，并手动编辑每个目录名称。
运行 spec2#
与图像处理管道一样，我们可以为管道中的步骤提供自定义参数。
关于在 spec2 阶段运行的所有内容的更多信息可以在 这里 找到。
首先，我们将仅对一个关联文件运行 spec2，因为如果有很多源，运行 spec2 可能会花费较长时间。我们将输出保存到当前目录中，该目录应为上述创建的 custom_spec2 目录。
test_asn = spec2_asns[0] # 从spec2_asns列表中获取第一个归档序列号（ASN）
print(f'Calibrating: {test_asn}') # 打印正在校准的归档序列号
# 检查校准后的文件是否已经存在
asn_data = json.load(open(test_asn)) # 读取ASN文件并解析为JSON格式
x1d_file = f"{asn_data['products'][0]['name']}_x1d.fits" # 构建x1d文件的名称
if os.path.exists(x1d_file): # 检查x1d文件是否存在
print(x1d_file, ': x1d file already exists.') # 如果存在，打印提示信息
else: # 如果文件不存在
spec2 = Spec2Pipeline.call(test_asn, save_results=True) # 调用Spec2Pipeline处理数据并保存结果
检查Spec2的输出#
Spec2的输出文件为cal.fits和x1d.fits。在这里，我们快速查看这些文件的一些重要部分。
# 从指定的ASN文件中加载JSON数据
asn_example = json.load(open(spec2_asns[0]))
# 获取第一个成员的速率文件名
rate_file = asn_example['products'][0]['members'][0]['expname']
# 获取第三个成员的源目录文件名
source_cat = asn_example['products'][0]['members'][2]['expname']
# 将速率文件名中的'rate'替换为'cal'，生成校准文件名
cal_file = rate_file.replace('rate', 'cal')
# 将速率文件名中的'rate'替换为'x1d'，生成一维光谱文件名
x1d_file = rate_file.replace('rate', 'x1d')
# 打开所有文件以进行查看
rate_hdu = fits.open(rate_file) # 打开速率文件
cal_hdu = fits.open(cal_file) # 打开校准文件
x1d_hdu = fits.open(x1d_file) # 打开一维光谱文件
cat = Table.read(source_cat) # 读取源目录表
# 首先查看我们期望从目录中识别出多少个源
print(f'There are {len(cat)} sources identified in the current catalog.\n') # 打印识别到的源数量
# 然后比较校准文件和一维光谱文件的扩展长度
print(f'The x1d file has {len(x1d_hdu)} extensions & the cal file has {len(cal_hdu)} extensions') # 打印扩展数量
请注意，每个文件中的第0个和最后一个扩展不包含科学数据，而其余扩展对应于每个源。x1d 文件包含每个源的提取光谱，而 cal 文件包含七个扩展中的2D切片信息，分别为 SCI、DQ、ERR、WAVELENGTH、VAR_POISSON、VAR_RNOISE 和 VAR_FLAT。
这部分原因就是为什么拥有一个精细化的源目录如此重要。如果有些源对您的研究没有用处，那么就没有必要创建切片并提取它们。
请注意，源的数量超过了文件中的扩展数量。这是因为管道默认只提取100个最亮的源。要更改此行为，请向管道提供参数 wfss_nbright，我们将在 下面 进行说明。
# 打印x1d_hdu对象的信息
print(x1d_hdu.info()) # 输出x1d_hdu的详细信息，包括数据类型、维度等
print(cal_hdu.info()) # 打印cal_hdu对象的信息
x1d 文件是一个 BinTable，因此每个数据扩展中包含了额外的列。x1d 进一步阅读 详细介绍了每一列包含的内容，但我们也可以通过查看其中一个数据列快速了解。
# 快速查看x1d文件中可用的列
print(x1d_hdu[1].data.columns) # 打印x1d文件中第二个HDU的数据列
进一步探索 Spec2 数据#
在Spec2数据中查找源#
每个cal和x1d文件的扩展名在头文件中都有一个源ID。这些值应与源目录中的值相匹配。
def find_source_ext(x1d_hdu, cal_hdu, source_id, info=True):
# 查找给定源ID在x1d和cal扩展中的位置
# x1d扩展首先处理
x1d_source_ids = np.array([x1d_hdu[ext].header['SOURCEID'] for ext in range(len(x1d_hdu))[1:-1]]) # 获取x1d文件中所有源ID，去掉第0个和最后一个数据扩展
wh_x1d = np.where(x1d_source_ids == source_id)[0][0] + 1 # 找到源ID的位置，并加1以考虑主头部
# 也查找cal扩展，但仅在SCI扩展中；
# 对于所有其他扩展填充源ID为-999，以获取正确的扩展值
cal_source_ids = np.array([cal_hdu[ext].header['SOURCEID'] if cal_hdu[ext].header['EXTNAME'] == 'SCI'
else -999 for ext in range(len(cal_hdu))[1:-1]]) # 仅在SCI扩展中获取源ID
wh_cal = np.where(cal_source_ids == source_id)[0][0] + 1 # 找到源ID的位置，并加1以考虑主头部
if info:
# 如果info为真，打印相关信息
print(f"All source IDs in x1d file:\n{x1d_source_ids}\n") # 打印x1d文件中的所有源ID
print(f"Extension {wh_x1d} in {x1d_hdu[0].header['FILENAME']} contains the data for source {source_id} from our catalog") # 打印x1d扩展信息
print(f"Extension {wh_cal} in {cal_hdu[0].header['FILENAME']} contains the data for source {source_id} from our catalog") # 打印cal扩展信息
return wh_x1d, wh_cal # 返回x1d和cal扩展中源ID的位置
# 我们来查找在之前的笔记本中识别的源118。
source_id = 118 # 定义源的ID为118
# 使用find_source_ext函数查找对应的x1d和cal数据的索引
wh_x1d_118, wh_cal_118 = find_source_ext(x1d_hdu, cal_hdu, source_id)
# 从x1d_hdu中提取源118的数据
x1d_data_118 = x1d_hdu[wh_x1d_118].data
# 从cal_hdu中提取源118的数据
cal_data_118 = cal_hdu[wh_cal_118].data
首先，让我们看看提取框在速率图像上的显示。这将帮助我们对管道的提取可靠性有一个整体的了解。
注意：在这个例子中，提取框并没有完全围绕光谱居中。目前正在进行校准工作，以更好地考虑NIRISS探测器上光谱轨迹形状的差异。关于此校准状态的更新可以在NIRISS jdox上找到。
# 用零填充坏像素的nan值，以便更容易查看
rate_data = rate_hdu[1].data # 从rate_hdu中提取数据
rate_data[np.isnan(rate_data)] = 0 # 将nan值替换为0
# 从cal数据的头文件中提取提取框参数：
cal_header = cal_hdu[wh_cal_118].header # 获取指定cal_hdu的头文件
sx_left = cal_header['SLTSTRT1'] # 提取左侧起始位置
swidth = cal_header['SLTSIZE1'] # 提取宽度
sx_right = cal_header['SLTSTRT1'] + swidth # 计算右侧位置
sy_bottom = cal_header['SLTSTRT2'] # 提取底部起始位置
sheight = cal_header['SLTSIZE2'] # 提取高度
sy_top = cal_header['SLTSTRT2'] + sheight # 计算顶部位置
# 绘制速率文件和提取框
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) # 创建1行2列的子图
ax1.imshow(rate_data, origin='lower', vmin=0.2, vmax=1, aspect='auto') # 显示速率数据，可能需要调整缩放
ax2.imshow(rate_data, origin='lower', vmin=0.2, vmax=1, aspect='auto') # 显示速率数据，可能需要调整缩放
rectangle = patches.Rectangle((sx_left, sy_bottom), swidth, sheight, edgecolor='darkorange', facecolor="None", linewidth=1) # 创建提取框
ax1.add_patch(rectangle) # 在ax1中添加提取框
ax1.set_title(rate_file) # 设置ax1的标题
rectangle2 = patches.Rectangle((sx_left, sy_bottom), swidth, sheight, edgecolor='darkorange', facecolor="None", linewidth=2) # 创建第二个提取框
ax2.add_patch(rectangle2) # 在ax2中添加提取框
ax2.set_xlim(sx_left-30, sx_right+30) # 设置ax2的x轴范围
ax2.set_ylim(sy_bottom-30, sy_top+30) # 设置ax2的y轴范围
ax2.set_title(f'Source {source_id} Zoom in') # 设置ax2的标题
plt.suptitle(f"{cal_hdu[0].header['FILTER']} {cal_hdu[0].header['PUPIL']}") # 设置整体标题
我们可以在这个框中查看提取的光谱，包括校准文件（cal file）和一维文件（x1d file）。在下面提取的光谱中，您可以看到来自星系的[OII]和H\(\beta\)发射线。
注意：在一维光谱中看到的上翘边缘效应是由于当前光谱校准时在提取框边缘的插值造成的。这也是一个持续的校准工作的一部分。
附加说明：来自数据处理管道的默认通量单位是Jansky。然而，在这些提取的光谱中，我们显示的单位是erg/s/cm²/Å。要关闭此功能，请在plot_cutout_and_spectrum中设置convert=False
def Fnu_to_Flam(wave_micron, flux_jansky):
# 将Jansky单位的 flux 转换为 erg/s/cm^2/Angstrom，输入波长单位为微米
f_lambda = 1E-21 * flux_jansky * (const.c.value) / (wave_micron**2) # 计算 erg/s/cm^2/Angstrom
return f_lambda # 返回转换后的 flux
def plot_cutout_and_spectrum(cal_data, x1d_data, cal_file, x1d_file, source_id, convert=True):
# 创建一个包含两个子图的图形，大小为11x5英寸
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 5))
# 绘制校准图像
ax1.imshow(cal_data, origin='lower', vmin=0, vmax=np.nanmax(cal_data)*.01, aspect='auto')
ax1.set_title(os.path.basename(cal_file)) # 设置子图标题为校准文件名
# 绘制光谱
wave = x1d_data['WAVELENGTH'].ravel() # 获取波长数据并展平
flux = x1d_data['FLUX'].ravel() # 获取光通量数据并展平
if convert: # 如果需要转换单位
flux = Fnu_to_Flam(wave, flux) # 将光通量从Fnu转换为Flam
fluxunit = 'egs/s/cm^2/Angstrom' # 设置单位为ergs/s/cm^2/Angstrom
else:
fluxunit = 'Jy' # 设置单位为Jy
ax2.plot(wave, flux) # 绘制光谱
ax2.set_title(os.path.basename(x1d_file)) # 设置子图标题为一维文件名
edge_buffer = int(len(flux) * .25) # 计算边缘缓冲区大小
max_flux = np.nanmax(flux[edge_buffer:edge_buffer*-1]) # 计算有效光通量的最大值
ax2.set_ylim(0, max_flux+(max_flux*0.1)) # 设置y轴范围，添加10%的缓冲区
ax2.set_xlabel('Wavelength (Microns)') # 设置x轴标签为波长（微米）
ax2.set_ylabel(f'Flux ({fluxunit})') # 设置y轴标签为光通量及其单位
# 根据滤光片类型设置x轴和y轴标签
if fits.getval(cal_file, 'FILTER') == 'GR150C':
ax1.set_xlabel('X Pixels (<--- dispersion direction)') # 设置x轴标签
ax1.set_ylabel('Y Pixels') # 设置y轴标签
else:
ax1.set_xlabel('X Pixels') # 设置x轴标签
ax1.set_ylabel('Y Pixels (<--- dispersion direction)') # 设置y轴标签
# 设置总标题，包含滤光片和源ID信息
plt.suptitle(f"{fits.getval(cal_file, 'FILTER')} {fits.getval(cal_file, 'PUPIL')} Source {source_id}", x=0.5, y=1)
def plot_cutout_and_spectrum(cal_data, x1d_data, cal_file, x1d_file, source_id):
"""
绘制切片图和光谱图
:param cal_data: 校准后的数据
:param x1d_data: 一维光谱数据
:param cal_file: 校准文件名
:param x1d_file: 一维光谱文件名
:param source_id: 源ID
"""
# 导入必要的库
import matplotlib.pyplot as plt # 导入绘图库
import numpy as np # 导入NumPy库用于数值计算
# 创建一个图形和两个子图
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8)) # 创建2行1列的子图，设置图形大小
# 绘制切片图
ax1.imshow(cal_data, cmap='gray', origin='lower') # 显示校准数据，使用灰度色图
ax1.set_title(f'Cutout for Source ID: {source_id}') # 设置切片图标题
ax1.set_xlabel('X Pixel') # 设置X轴标签
ax1.set_ylabel('Y Pixel') # 设置Y轴标签
# 绘制光谱图
wavelengths = x1d_data['wavelength'] # 获取波长数据
flux = x1d_data['flux'] # 获取光谱通量数据
ax2.plot(wavelengths, flux, color='blue') # 绘制光谱曲线
ax2.set_title(f'Spectrum for Source ID: {source_id}') # 设置光谱图标题
ax2.set_xlabel('Wavelength (microns)') # 设置X轴标签
ax2.set_ylabel('Flux (Jy)') # 设置Y轴标签
# 保存图像
plt.tight_layout() # 自动调整子图参数以填充整个图像区域
plt.savefig(f'{source_id}_cutout_and_spectrum.png') # 保存图像为PNG文件
plt.show() # 显示图像
以上代码在每一行添加了中文注释，以便于理解代码的功能和结构。
限制提取源的数量#
由于提取大量源需要很长时间，让我们看看是否可以减少提取的源的数量。我们将通过参数输入（针对亮点源）以及进一步精炼的源目录来实现这一点。
使用参数输入进行限幅提取#
在此校准中，我们明确指出以下步骤的参数：
extract_2d，在此步骤中我们设置了wfss_nbright，该参数限制了提取的亮源数量，以及wfss_mmag_extract，该参数设置了我们希望提取的最暗星等的限制。进一步阅读
在这种情况下，我们将提取限制为仅提取10个最亮的天体。
# 检查校准文件是否已经存在
asn_data = json.load(open(spec2_asns[0])) # 从指定的ASN文件中加载数据
x1d_file = os.path.join(param_outdir, f"{asn_data['products'][0]['name']}_x1d.fits") # 构建x1d文件的完整路径
if os.path.exists(x1d_file): # 检查x1d文件是否存在
print(x1d_file, ': x1d file already exists.') # 如果存在，打印提示信息
else: # 如果文件不存在
spec2 = Spec2Pipeline.call(spec2_asns[0], # 调用Spec2Pipeline处理数据
steps={'extract_2d': {'wfss_nbright': 10}, }, # 设置提取2D数据的步骤参数
save_results=True, # 保存处理结果
output_dir=param_outdir) # 指定输出目录
# 导入必要的库
import glob # 用于查找符合特定规则的文件路径名
import os # 用于处理文件和目录路径
from astropy.io import fits # 用于处理FITS文件格式
# 打开x1d文件以检查提取了多少个源
x1ds = glob.glob(os.path.join(param_outdir, '*x1d.fits*')) # 查找指定目录下所有符合条件的x1d.fits文件
with fits.open(x1ds[0]) as temp_x1d: # 打开找到的第一个x1d文件
print(f'The x1d file has {len(temp_x1d)-2} extracted sources') # 输出提取的源数量，减去2是因为FITS文件的前两个HDU通常是头信息
使用源目录进行极限提取#
在上一个笔记本中，我们以几种不同的方式限制了目录。在这里，我们将目录限制在一个特定的亮度范围内，然后使用这个新目录提取数据。
# 导入必要的库
source_match_cats = np.sort(glob.glob('*source-match_cat.ecsv')) # 获取所有匹配源的文件并排序
source_match_cat = Table.read(source_cat) # 读取源匹配目录的表格数据
# 查看可能的星等范围
mags = cat['isophotal_vegamag'] # 获取星等数据
min_vegmag = mags.min() # 计算最小星等
max_vegmag = mags.max() # 计算最大星等
print(f"Magnitude range: {min_vegmag} - {max_vegmag}") # 打印星等范围
# 源118的Vega星等应为~21.68
source_id = 118 # 定义要查找的源ID
source_mag = source_match_cat[source_match_cat['label'] == source_id]['isophotal_vegamag'][0] # 获取指定源的星等
print(f"Magnitude for source in previous notebook (source {source_id}) : {source_mag}") # 打印指定源的星等
# 查找在特定星等范围内的源的目录（21.18，以确保包含我们的示例源）
min_mag = 21.1 # 设置最小星等
max_mag = 21.2 # 设置最大星等
# 从源匹配目录中筛选出星等在指定范围内的源
mag_cat = source_match_cat[(source_match_cat['isophotal_vegamag'] >= min_mag) & (source_match_cat['isophotal_vegamag'] <= max_mag)]
# 选择感兴趣的列：标签、x中心、y中心、天空中心、是否扩展、绝对星等、视星等
mag_cat['label', 'xcentroid', 'ycentroid', 'sky_centroid', 'is_extended', 'isophotal_abmag', 'isophotal_vegamag']
mag_source_ext = 'mag-limit_cat.ecsv' # 定义新的文件扩展名
new_cat = Table(mag_cat) # 将行实例转换为数据框
# 保存新的目录，使用唯一的名称
new_cat_name = source_cat.replace('cat.ecsv', mag_source_ext) # 替换文件名中的扩展名
new_cat.write(new_cat_name, overwrite=True) # 写入新的目录，允许覆盖
print('Saved:', new_cat_name) # 打印保存的文件名
一旦我们有了一个已限制为特定源的源目录（source catalog），接下来我们就可以将剩余的目录与这些源进行匹配。
for cat_name in source_match_cats: # 遍历源匹配目录中的每个目录名
if cat_name == source_cat: # 如果当前目录名是源目录名
# 基础目录已经保存过了
continue # 跳过当前循环，继续下一个目录
# 匹配当前目录与基础目录之间的源ID
cat = Table.read(cat_name) # 读取当前目录的数据表
cat.add_index('label') # 为数据表添加索引，索引字段为'label'
new_cat = cat.loc[list(mag_cat['label'])] # 根据基础目录中的'label'筛选当前目录的数据
# 检查确保所有源都存在
print(repr(new_cat['label', 'xcentroid', 'ycentroid', 'sky_centroid', 'is_extended', 'isophotal_abmag', 'isophotal_vegamag'])) # 打印新目录中的相关字段
# 使用唯一名称保存新目录
new_cat_name = cat_name.replace('cat.ecsv', mag_source_ext) # 替换文件扩展名以生成新目录名
new_cat.write(new_cat_name, overwrite=True) # 将新目录写入文件，允许覆盖
print('Saved:', new_cat_name) # 打印保存的文件名
print() # 打印空行以分隔输出
一旦新的目录（catalogs）生成，我们必须更新关联文件（association files）以使用新的目录。注意：如果您希望再次使用之前的源目录进行校准（calibrate），则需要再次更新关联文件。
new_source_ext = mag_source_ext # 将新的源扩展名赋值给变量
# loop through all of the spec2 association files in the current directory
for asn in spec2_asns: # 遍历当前目录中的所有spec2关联文件
asn_data = json.load(open(asn)) # 读取关联文件的内容并解析为JSON格式
for product in asn_data['products']: # 对于每个产品，检查其成员
for member in product['members']: # 每个产品有多个成员
if member['exptype'] == 'sourcecat': # 检查成员类型是否为'sourcecat'
cat_in_asn = member['expname'] # 获取当前成员的扩展名
# check that we haven't already updated the source catalog name
if new_source_ext not in cat_in_asn: # 确保我们还没有更新源目录名称
new_cat = cat_in_asn.replace('cat.ecsv', new_source_ext) # 替换扩展名为新的源扩展名
# actually update the association file member
if os.path.exists(new_cat): # 检查新的目录文件是否存在
member['expname'] = new_cat # 更新成员的扩展名为新的目录文件名
else:
print(f"{new_cat} does not exist in the current working directory") # 如果文件不存在，打印提示信息
# write out the new association file
with open(asn, 'w', encoding='utf-8') as f: # 以写入模式打开关联文件
json.dump(asn_data, f, ensure_ascii=False, indent=4) # 将更新后的数据写回文件，确保中文字符正常显示
# 双重检查源目录是否已更改
for spec2_asn in spec2_asns: # 遍历每个spec2_asn文件
asn_check = json.load(open(spec2_asn)) # 读取并解析JSON文件
for product in asn_check['products']: # 遍历产品列表
for key, value in product.items(): # 遍历每个产品的键值对
if key == 'members': # 如果键是'members'
for member in value: # 遍历成员列表
if member['exptype'] == 'sourcecat': # 如果成员类型是'sourcecat'
print(f" {member['exptype']}: {member['expname']}") # 打印成员类型和名称
else: # 如果键不是'members'
print(f"{key}: {value}") # 打印其他键值对
现在，当我们对所有内容进行校准时，对于单个文件来说，所需时间应该会少很多，因为源的数量有限。然而，我们将在此访问中校准所有文件，因此这个单元格的运行可能需要一些时间。
# 使用新的源目录进行校准
for spec2_asn in spec2_asns: # 遍历所有的spec2关联文件
# 检查校准后的文件是否已经存在
asn_data = json.load(open(spec2_asn)) # 读取关联文件的JSON数据
x1d_file = os.path.join(source_outdir, f"{asn_data['products'][0]['name']}_x1d.fits") # 构建x1d文件的路径
if os.path.exists(x1d_file): # 检查x1d文件是否存在
print(x1d_file, ': x1d file already exists.') # 如果存在，打印提示信息
else: # 如果文件不存在
spec2 = Spec2Pipeline.call(spec2_asn, save_results=True, output_dir=source_outdir) # 调用Spec2Pipeline进行处理并保存结果
最终可视化#
现在所有数据都已校准，查看所有提取的文件是很有用的。来自spec2的cal和x1d文件在每个抖动步骤中被提取，如下所示。然后，Spec3将每个源的单个抖动x1d文件转换为每个光栅和滤光片的单个组合光谱。
请注意，对于GR150R数据，色散方向在-Y方向，而对于GR150C数据，色散方向在-X方向。
查看单个源的所有文件#
# 探索新数据
# 获取所有以"x1d.fits"结尾的文件，并按名称排序
x1ds = np.sort(glob.glob(os.path.join(source_outdir, "*x1d.fits")))
# 从第一个文件中获取所有源ID的列表，以便在此示例中查看
sources = [fits.getval(x1ds[0], 'SOURCEID', ext=ext) for ext in range(len(fits.open(x1ds[0])))[1:-1]]
# 设置要查找的源ID
source_id = 118
# 绘制每个x1d/cal文件
for i, x1d_file in enumerate(x1ds):
# 生成对应的cal文件名
cal_file = x1d_file.replace('x1d.fits', 'cal.fits')
# 打开x1d和cal文件
with fits.open(x1d_file) as x1d_hdu, fits.open(cal_file) as cal_hdu:
try:
# 查找源在x1d和cal文件中的扩展
wh_x1d, wh_cal = find_source_ext(x1d_hdu, cal_hdu, source_id, info=False)
except IndexError:
# 如果源不在此观测中，则跳过
continue
# 获取x1d和cal数据
x1d_data = x1d_hdu[wh_x1d].data
cal_data = cal_hdu[wh_cal].data
# 绘制切片和光谱
plot_cutout_and_spectrum(cal_data, x1d_data, cal_file, x1d_file, source_id)
将这些文件叠加在一起以进行比较。两个光栅（grisms）将使用不同的线条样式，以突出可能由于校准（calibration）、包括污染（contamination）而导致的任何差异，而每个阻挡滤光片（blocking filter）将使用不同的颜色。
# 将不同的光栅图叠加在一起以显示同一源的光谱
x1ds = np.sort(glob.glob(os.path.join(source_outdir, "*x1d.fits"))) # 获取所有x1d文件并排序
# 从第一个文件中获取所有源ID的列表，以便在此示例中查看
sources = [fits.getval(x1ds[0], 'SOURCEID', ext=ext) for ext in range(len(fits.open(x1ds[0])))[1:-1]] # 获取源ID列表
source_id = 118 # 选择要查看的源ID
# 创建一个包含三个面板的图形
src118_fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(11, 9), sharex=True, sharey=True) # 创建子图
ls_dict = {'GR150R': '--', # 定义光栅线型字典
'GR150C': '-',
}
color_dict = {'F115W': '#e1cb00', # 定义滤光片颜色字典
'F150W': '#32b45c',
'F200W': '#099ab4',
}
max_fluxes = [] # 存储最大通量的列表
all_waves = [] # 存储所有波长的列表
# 绘制每个x1d文件
for i, x1d_file in enumerate(x1ds): # 遍历每个x1d文件
with fits.open(x1d_file) as x1d_hdu: # 打开x1d文件
try:
wh_x1d, wh_cal = find_source_ext(x1d_hdu, cal_hdu, source_id, info=False) # 查找源的扩展
except IndexError: # 如果源不在此观测中
continue # 跳过此文件
x1d_data = x1d_hdu[wh_x1d].data # 获取x1d数据
grism = x1d_hdu[0].header['FILTER'] # 获取光栅类型
filter = x1d_hdu[0].header['PUPIL'] # 获取滤光片信息
wave = x1d_data['WAVELENGTH'].ravel() # 获取波长数据并展平
flux = x1d_data['FLUX'].ravel() # 获取通量数据并展平
flux = Fnu_to_Flam(wave, flux) # 将Fnu转换为Flam
fluxunits = 'erg/s/cm^2/Angstrom' # 定义通量单位
ax1.plot(wave, flux, color=color_dict[filter], ls=ls_dict[grism], alpha=0.7) # 在ax1上绘制光谱
if grism == 'GR150C': # 如果光栅是GR150C
ax2.plot(wave, flux, color=color_dict[filter], ls=ls_dict[grism], alpha=0.7) # 在ax2上绘制光谱
else: # 否则
ax3.plot(wave, flux, color=color_dict[filter], ls=ls_dict[grism], alpha=0.7) # 在ax3上绘制光谱
# 保存最大通量以供绘图，去除边缘效应
edge_buffer = int(len(flux) * .25) # 计算边缘缓冲区大小
max_fluxes.append(np.nanmax(flux[edge_buffer:edge_buffer*-1])) # 记录最大通量
all_waves.extend(wave) # 扩展所有波长列表
# 设置绘图限制和标签
ax1.set_ylim(0, np.max(max_fluxes)) # 设置y轴范围
ax1.set_xlim(np.min(all_waves), np.max(all_waves)) # 设置x轴范围
src118_fig.suptitle(f'Source {source_id}') # 设置总标题
ax1.set_title('GR150R & GR150C') # 设置ax1标题
ax2.set_title('GR150C') # 设置ax2标题
ax3.set_title('GR150R') # 设置ax3标题
for ax in [ax1, ax2, ax3]: # 为每个轴设置标签
ax.set_xlabel('Wavelength (Microns)') # 设置x轴标签
ax.set_ylabel(f'Flux ({fluxunits})') # 设置y轴标签
# 为每个滤光片添加标签
for filt, color in color_dict.items(): # 遍历滤光片和颜色字典
ax1.plot(0, 0, color=color, label=filt) # 在ax1上添加标签
ax1.legend(bbox_to_anchor=(1, 1.05)) # 设置图例位置
src118_fig.tight_layout() # 调整布局以避免重叠
查看单个文件的所有源#
请注意，有些源可能实际上并没有提取出任何有趣的内容。如果是这种情况，请返回您的源目录和图像，以确保您已正确识别源，并且目标在切割图像中居中。此笔记本使用简单的示例来创建和限制源目录，因此它可能并不总是显示最具科学价值的源。
import os # 导入os模块，用于处理文件和目录路径
from astropy.io import fits # 导入astropy.io.fits模块，用于处理FITS文件
x1d_file = os.path.join(source_outdir, 'jw02079004002_11101_00002_nis_x1d.fits') # 定义x1d文件的路径
cal_file = x1d_file.replace('x1d.fits', 'cal.fits') # 将x1d文件名替换为cal文件名
with fits.open(x1d_file) as x1d_hdu, fits.open(cal_file) as cal_hdu: # 打开x1d和cal文件
for ext in range(len(x1d_hdu))[1:-1]: # 遍历x1d文件的扩展，跳过第一个和最后一个扩展
source_id = x1d_hdu[ext].header['SOURCEID'] # 获取当前扩展的源ID
try:
wh_x1d, wh_cal = find_source_ext(x1d_hdu, cal_hdu, source_id, info=False) # 查找源在x1d和cal文件中的扩展索引
except IndexError: # 捕获索引错误
# 这意味着该源不在此观测中
continue # 跳过当前循环，继续下一个扩展
x1d_data = x1d_hdu[wh_x1d].data # 获取x1d数据
cal_data = cal_hdu[wh_cal].data # 获取cal数据
plot_cutout_and_spectrum(cal_data, x1d_data, cal_file, x1d_file, source_id) # 绘制cutout和光谱
可视化提取源在分散图像中的位置，以及与直接图像的比较#
# 设置图形
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 5), sharex=True, sharey=True)
# **光栅数据
# 为分散图像绘制
# x1d_file 和 cal_file 使用相同的根名，因为我们正在查看单个源
rate_file = os.path.basename(x1d_file.replace('x1d.fits', 'rate.fits'))
# 用零填充坏像素的nan值，以便更容易查看
with fits.open(rate_file) as rate_hdu:
rate_data = rate_hdu[1].data # 读取速率数据
rate_data[np.isnan(rate_data)] = 0 # 将nan值替换为0
# 绘制速率文件和提取框
ax1.imshow(rate_data, origin='lower', vmin=0, vmax=np.nanmax(rate_data)*0.01, aspect='auto')
with fits.open(x1d_file) as x1d_hdu, fits.open(cal_file) as cal_hdu:
for ext in range(len(x1d_hdu))[1:-1]: # 遍历扩展
source_id = x1d_hdu[ext].header['SOURCEID'] # 获取源ID
try:
wh_x1d, wh_cal = find_source_ext(x1d_hdu, cal_hdu, source_id, info=False) # 查找源的扩展
except IndexError:
# 这意味着源不在此观测中
continue
x1d_data = x1d_hdu[wh_x1d].data # 获取x1d数据
cal_data = cal_hdu[wh_cal].data # 获取校准数据
# 从校准数据的头文件中提取框参数：
cal_header = cal_hdu[wh_cal].header
sx_left = cal_header['SLTSTRT1'] # 提取框左侧
swidth = cal_header['SLTSIZE1'] # 提取框宽度
sx_right = cal_header['SLTSTRT1'] + swidth # 提取框右侧
sy_bottom = cal_header['SLTSTRT2'] # 提取框底部
sheight = cal_header['SLTSIZE2'] # 提取框高度
sy_top = cal_header['SLTSTRT2'] + sheight # 提取框顶部
rectangle = patches.Rectangle((sx_left, sy_bottom), swidth, sheight, edgecolor='darkred', facecolor="None", linewidth=1) # 创建提取框
ax1.add_patch(rectangle) # 将提取框添加到图中
ax1.text(sx_left, sy_top+10, source_id, fontsize=12, color='darkred') # 在提取框上方添加源ID
ax1.set_title(f"{os.path.basename(x1d_file).split('_nis')[0]}: {cal_hdu[0].header['FILTER']} {cal_hdu[0].header['PUPIL']}") # 设置标题
# **成像数据
asn_data = json.load(open(fits.getval(x1d_file, 'ASNTABLE'))) # 读取关联数据
i2d_name = asn_data['products'][0]['members'][1]['expname'] # 获取i2d文件名
cat_name = asn_data['products'][0]['members'][2]['expname'] # 获取目录文件名
with fits.open(os.path.join('../../', data_dir, custom_run_image3, i2d_name)) as i2d:
ax2.imshow(i2d[1].data, origin='lower', aspect='auto', vmin=0, vmax=np.nanmax(i2d[1].data)*0.01) # 绘制i2d数据
ax2.set_title(f"{os.path.basename(i2d_name).split('_i2d')[0]}") # 设置标题
# 也绘制相关目录
cat = Table.read(cat_name) # 读取目录
extended_sources = cat[cat['is_extended'] == 1] # 1为真；即为扩展源
point_sources = cat[cat['is_extended'] == 0] # 0为假；即为点源
for color, sources in zip(['darkred', 'black'], [extended_sources, point_sources]): # 遍历颜色和源
# 绘制源
ax2.scatter(sources['xcentroid'], sources['ycentroid'], s=40, facecolors='None', edgecolors=color, alpha=0.9, lw=2)
# 添加源标签
for i, source_num in enumerate(sources['label']):
ax2.annotate(source_num,
(sources['xcentroid'][i]+0.5, sources['ycentroid'][i]+0.5),
fontsize=12,
color=color) # 在源旁边添加标签
fig.tight_layout() # 调整图形布局
继续深入探索，包括使用spec3阶段的管道！
