运行图像处理管道并创建源目录#
在本次运行中，我们从已经通过探测器1管道校准的速率文件开始。这将节省时间，因为我们不需要编辑探测器1过程中执行的任何步骤。因此，WFSS（宽场分光光谱）运行的第一步校准是将速率文件的直接图像通过JWST管道的Image2和Image3步骤。这包括创建源目录，通常需要根据管道的默认值进行调整。没有良好的源目录将导致在分散的WFSS图像中源的提取不够优化。
用例：管道的默认参数未能提取预期的源，因此需要设置自定义参数以获得新的合成图像和源目录。
数据：JWST/NIRISS图像和来自程序2079观察004的光谱。这些数据应存储在一个名为data的单一目录中，可以从之前的笔记本00_niriss_mast_query_data_setup.ipynb下载。
工具：astropy, crds, glob, jdaviz, json, jwst, matplotlib, numpy, os, pandas, urllib, warnings, zipfile
跨仪器：NIRISS
内容
作者：Rachel Plesha (rplesha@stsci.edu), Camilla Pacifici (cpacifici@stsci.edu), JWebbinar笔记本。
首次发布：2024年5月
最后测试：此笔记本最后一次测试是在JWST管道版本1.12.5和CRDS上下文jwst_1225.pmap下进行的。
导入与数据设置#
CRDS 文档#
此文档提供了关于 CRDS（Calibration Reference Data System）的详细信息，特别是在 JWST（James Webb Space Telescope）数据处理中的应用。CRDS 是一个用于管理和分发天文数据处理所需的参考文件的系统。
主要内容#
CRDS 的功能：CRDS 负责确保在数据处理过程中使用的所有参考文件都是最新的，并且与特定的科学目标和观测条件相匹配。
参考文件的类型：CRDS 管理多种类型的参考文件，包括但不限于：
校准文件（Calibration Files）
掩蔽文件（Mask Files）
增益文件（Gain Files）
如何使用 CRDS：用户可以通过命令行工具或编程接口访问 CRDS，以获取所需的参考文件。
更新和维护：CRDS 定期更新，以确保用户能够访问最新的校准数据和参考文件。
有关更详细的信息和使用指南，请访问 CRDS Documentation。
# 更新CRDS路径到你的本地目录
%env CRDS_PATH=crds_cache # 设置CRDS缓存路径为'crds_cache'
# 设置CRDS服务器的URL
%env CRDS_SERVER_URL=https://jwst-crds.stsci.edu # 设置CRDS服务器的URL
import os # 导入操作系统模块
import glob # 导入用于文件路径操作的模块
import json # 导入JSON处理模块
import warnings # 导入警告模块
import urllib # 导入用于URL处理的模块
import zipfile # 导入用于处理ZIP文件的模块
import numpy as np # 导入NumPy库，用于数值计算
import pandas as pd # 导入Pandas库，用于数据处理和分析
import astropy.units as u # 导入Astropy单位模块
from astropy.io import fits # 从Astropy导入FITS文件处理模块
from astropy.coordinates import SkyCoord # 从Astropy导入天文坐标模块
from astropy.table import Table # 从Astropy导入表格模块
from jdaviz import Imviz # 从jdaviz导入Imviz模块，用于图像可视化
from matplotlib import pyplot as plt # 从Matplotlib导入绘图模块
%matplotlib widget # 使用交互式绘图
# %matplotlib inline # 注释掉的代码，用于在Jupyter Notebook中内嵌绘图
from jwst.pipeline import Image2Pipeline # 从JWST管道导入Image2Pipeline模块
from jwst.pipeline import Image3Pipeline # 从JWST管道导入Image3Pipeline模块
检查您正在使用的JWST（詹姆斯·韦伯太空望远镜）管道版本。要查看可用的最新管道版本或安装以前的版本，请访问 GitHub。同时，请确认您使用的 CRDS（参考数据服务）版本。 CRDS文档 解释了如何在JWST管道中设置特定的上下文。如果这两个值与上述最后测试的说明不同，则此笔记本可能会存在差异。
import jwst # 导入JWST相关库
import crds # 导入CRDS（天文数据记录系统）库
print('JWST Pipeliene Version:', jwst.__version__) # 打印JWST管道的版本信息
print('CRDS Context:', crds.get_context_name('jwst')) # 打印当前使用的CRDS上下文
数据目录 data_dir 应该包含所有的关联文件和速率文件，并且放在一个单一的平面目录中。default_run_image3 和 custom_run_image3 是我们稍后将用于校准数据的目录。它们被分开是为了便于比较这两种输出。
data_dir = 'data' # 数据目录，存放处理后的数据
default_run_image3 = 'default_image3_calibrated' # 默认图像3运行结果保存路径（在data_dir内部）
custom_run_image3 = 'custom_image3_calibrated' # 自定义图像3运行结果保存路径（在data_dir内部）
关联文件期望满足以下条件：1) 所有数据都在同一目录中；2) 您在该目录中执行管道调用。由于这个原因，我们需要切换到数据目录以运行成像管道。
# 如果您尚未从笔记本00下载数据，请运行此单元格。否则，可以跳过它。
# 从Box下载未校准的数据到数据目录：
boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/niriss_wfss_advanced/niriss_wfss_advanced_01_input.zip' # 数据链接
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)) # 重命名并移动到当前目录
else: # 如果文件不是CSV格式
# 移动到数据目录
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
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 # 继续执行后面的代码
for temp_dir in [default_run_image3, custom_run_image3]: # 遍历默认和自定义的运行图像目录
if not os.path.exists(temp_dir): # 检查目录是否存在
os.mkdir(temp_dir) # 如果不存在，则创建该目录
默认成像处理管道运行#
首先，对所有使用WFSS数据观测的直接图像运行默认的image2和image3步骤。
运行默认图像2#
图像2是在直接图像速率文件上运行的。虽然您的程序应该有有效的关联文件可以从MAST下载，但如果出于任何原因您需要自己制作关联文件，请参见创建自定义ASN文件。
查看 Level 2 图像关联文件#
首先，查看关联（ASN）文件，以更好地理解其中包含的所有内容。
对于 image2 关联文件，每个观察序列中的每个抖动位置应该有一个 ASN 文件，这由 曝光策略 设置。在这种情况下，ASN 文件的数量应该与 rate_df 中的直接图像数量相匹配（FILTER=CLEAR），因为每个直接图像在观察序列中处于独特的抖动位置（XOFFSET, YOFFSET）。对于该程序和观测，在更换光栅之前有一张只有一个抖动的直接图像，在更换光栅之间有一张带有四个抖动的直接图像，以及在序列末尾有一张带有三个抖动的直接图像。这导致每个观察序列总共有八张图像，而在使用阻挡滤光片 F115W -> F115W -> F150W -> F150W -> F200W 的观测中，有五个观察序列。
import glob # 导入glob模块，用于文件路径匹配
# 查找当前目录下所有包含'image2'和'asn'的json文件
image2_asns = glob.glob('*image2*asn*.json')
# 打印找到的'image2' ASN文件的数量
print(len(image2_asns), 'Image2 ASN files') # 应该有40个image2的asn文件
# 打印与'CLEAR'过滤器匹配的速率文件数量
print(len(rate_df[rate_df['FILTER'] == 'CLEAR']), 'Direct Image rate files') # 速率文件的数量应与关联文件的数量匹配
# 查看其中一个关联文件
# 加载第一个关联文件的数据
asn_data = json.load(open(image2_asns[0]))
# 遍历关联数据中的每个键值对
for key, data in asn_data.items():
# 打印键及其对应的数据
print(f"{key} : {data}")
从这个关联中，我们可以得出关于观测的许多信息：
从
asn_type和asn_rule，我们可以看到这是一个 image2 关联。从
degraded_status，我们可以看到没有需要排除在校准之外的曝光。从
constraints，我们可以看到这不是一个时间序列观测（TSO），该观测是程序 2079 的一部分，使用 NIRISS 进行观测，采用了 CLEAR（即 WFSS 的成像）和 F150W 滤光片。从
products，我们可以看到只有一个关联的曝光。这在 image2 中是典型的，因为通常每个观测序列的每个微调只有一个曝光。
# 特别是，查看关联文件中的产品文件名：
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}") # 打印键名和对应的值
运行 image2#
rate.fits 产品将被校准为 cal.fits 文件。关于在图像处理（Image2）部分执行的步骤的更多信息，可以在 Image2 管道文档 中找到。
在这种情况下，我们将输出保存到运行管道的同一目录中，以便随后可以使用输出的 cal 文件来运行图像处理（Image3）管道。
for img2_asn in image2_asns: # 遍历所有的image2关联文件
# check if the calibrated file already exists # 检查校准文件是否已经存在
asn_data = json.load(open(img2_asn)) # 读取关联文件的内容
cal_file = f"{asn_data['products'][0]['name']}_cal.fits" # 构建校准文件的名称
if os.path.exists(cal_file): # 如果校准文件已存在
print(cal_file, 'cal file already exists.') # 输出文件已存在的提示
continue # 跳过当前循环，继续下一个关联文件
# if not, calibrated with image2 # 如果不存在，则进行校准
img2 = Image2Pipeline.call(img2_asn, save_results=True) # 调用Image2Pipeline进行校准并保存结果
运行默认图像3#
查看 Level 3 关联文件#
内容与 image2 非常相似，但请注意，现在有更多的成员被关联在一起，并且它们使用的是来自 image2 的单个指向校准文件。image3 重新采样并组合来自所有抖动和观测序列的相同阻挡滤光片（NIRISS WFSS 的 PUPIL）图像，以形成单个图像，这导致了更少的 image3 关联文件。
import glob # 导入glob模块，用于文件路径操作
import numpy as np # 导入numpy库，用于数组和数值计算
# 查找当前目录下所有包含'image3'和'asn'的json文件
image3_asns = glob.glob('*image3*asn*.json')
# 打印找到的image3 ASN文件的数量，应该是3个
print(len(image3_asns), 'Image3 ASN files') # 应该有3个image3关联文件
# 计算使用的唯一阻挡滤光片的数量，应该与image3关联文件数量匹配
uniq_filters = np.unique(rate_df[rate_df['FILTER'] == 'CLEAR']['PUPIL'])
# 打印使用的唯一滤光片数量及其名称
print(f"{len(uniq_filters)} unique filters used: {uniq_filters}")
# 查看其中一个关联文件
# 从第一个关联文件中加载数据
image3_asn_data = json.load(open(image3_asns[0]))
# 遍历关联文件中的每个键值对
for key, data in image3_asn_data.items():
# 打印键和值
print(f"{key} : {data}")
# 特别地，查看关联文件中的产品文件名：
for product in image3_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}") # 打印键名和对应的值
运行 image3#
在Image3中，cal.fits个别指向文件将被校准为一个合并的i2d.fits图像。Image3步骤也是最终源目录创建的地方，因此我们可以更改一些输入参数以获得更精细的输出源目录。这将在下面的自定义成像管道运行部分中完成。然而，我们将首先使用默认参数对数据进行校准。有关在管道的Image3部分执行的步骤的更多信息，请参阅Image3管道文档。
注意：Image3的运行可能需要一些时间
for img3_asn in image3_asns: # 遍历所有的image3关联文件
# check if the calibrated file already exists # 检查校准文件是否已存在
asn_data = json.load(open(img3_asn)) # 读取关联文件的JSON数据
cal_file = os.path.join(default_run_image3, f"{asn_data['products'][0]['name']}_i2d.fits") # 构建校准文件的路径
if os.path.exists(cal_file): # 如果校准文件已经存在
print(cal_file, 'cal file already exists.') # 打印文件已存在的信息
continue # 跳过当前循环，继续下一个
# if not, calibrated with image3 # 如果不存在，则进行image3的校准
img3 = Image3Pipeline.call(img3_asn, save_results=True, output_dir=default_run_image3) # 调用Image3Pipeline进行校准并保存结果
# 移除不必要的文件以节省磁盘空间
for crf in glob.glob(os.path.join(default_run_image3, '*crf.fits')): # 遍历匹配路径下所有以'crf.fits'结尾的文件
os.remove(crf) # 删除找到的文件
检查默认结果#
# 这些都是来自Image3管道的结果
image3_i2d = np.sort(glob.glob(os.path.join(default_run_image3, '*i2d.fits'))) # 获取多个dither/mosaic的组合图像并排序
image3_segm = np.sort(glob.glob(os.path.join(default_run_image3, '*segm.fits'))) # 获取定义源范围的分割图并排序
image3_cat = np.sort(glob.glob(os.path.join(default_run_image3, '*cat.ecsv'))) # 获取定义特定像素下源的RA/Dec的源目录并排序
Matplotlib#
Matplotlib 在某些方面存在局限性，而 ImViz 可能更适合您的需求——特别是如果您喜欢以 WCS 坐标查看事物的话。为了笔记本的目的，我们将重点介绍一些使用 matplotlib 包的关键领域。
使用 i2d 组合图像和由 Image3 生成的源目录，我们可以直观地检查管道找到的源位置是否令人满意。在以下图中，管道定义为扩展源的部分以橙红色显示，而管道定义为点源的部分则以灰色显示。这一定义会影响 WFSS 图像中的提取框。
fig = plt.figure(figsize=(10, 10)) # 创建一个10x10英寸的图形
cols = 2 # 设置每行的子图数量为2
rows = int(np.ceil(len(image3_i2d) / cols)) # 计算所需的行数，向上取整
for plt_num, img in enumerate(image3_i2d): # 遍历每个图像及其索引
# 确定子图的位置
xpos = (plt_num % 40) % cols # 计算当前子图在列中的位置
ypos = ((plt_num % 40) // cols) # 计算当前子图在行中的位置，使用整除以获得整数
# 创建子图
ax = plt.subplot2grid((rows, cols), (ypos, xpos)) # 在指定的网格位置创建子图
# 绘制图像
with fits.open(img) as hdu: # 打开FITS文件
ax.imshow(hdu[1].data, vmin=0, vmax=0.3, origin='lower') # 显示图像数据，设置颜色范围和原点位置
ax.set_title(f"obs{hdu[0].header['OBSERVTN']} {hdu[0].header['PUPIL']}") # 设置子图标题
# 同时绘制相关的目录
cat = Table.read(img.replace('i2d.fits', 'cat.ecsv')) # 读取对应的目录文件
extended_sources = cat[cat['is_extended'] == 1] # 筛选出扩展源，1表示为真
point_sources = cat[cat['is_extended'] == 0] # 筛选出点源，0表示为假
# 绘制扩展源的散点图
ax.scatter(extended_sources['xcentroid'], extended_sources['ycentroid'], s=20, facecolors='None', edgecolors='orangered', alpha=0.9)
# 绘制点源的散点图
ax.scatter(point_sources['xcentroid'], point_sources['ycentroid'], s=20, facecolors='None', edgecolors='dimgrey', alpha=0.9)
# 帮助调整坐标轴以避免重叠；如果这不起作用，可以手动设置
plt.tight_layout() # 自动调整子图参数以使其适合图形区域
分割图（segmentation maps）也是Image3流程的产物，它们用于帮助确定源目录（source catalog）。让我们来看看这些图，以确保我们对它所定义的源感到满意。
在分割图中，每个蓝色斑块（blob）应该对应一个物理目标（physical target）。有些情况下，源可能会发生混叠（blended），在这种情况下，制作分割图和源目录的参数应该进行调整。下面的观测004 F200W滤光片图像中可以看到一个例子，其中两个星系在坐标大约为(1600, 1300)的位置混叠成一个源。关于这一点将在自定义成像流程运行中详细讨论。
请注意，由于滤光片偏移差异，下面切割图所显示的视场（field of views）在每个滤光片下是不同的。
# 我们将多次查看这个，因此将其定义为一个函数
def plot_image_and_segmentation_map(i2d_images, segm_images, xmin=1250, xmax=1750, ymin=1250, ymax=1750, cat_suffix='cat.ecsv'):
cols = 2 # 设置每行的列数为2
rows = len(i2d_images) # 获取i2d_images的数量作为行数
fig = plt.figure(figsize=(10, 10*(rows/2))) # 创建一个图形对象，设置大小
for plt_num, img in enumerate(np.sort(np.concatenate([segm_images, i2d_images]))):
# 确定子图的位置
xpos = (plt_num % 40) % cols # 计算当前子图的x坐标
ypos = ((plt_num % 40) // cols) # 计算当前子图的y坐标，使用整除以获得整数
# 创建子图
ax = plt.subplot2grid((rows, cols), (ypos, xpos))
if 'i2d' in img: # 如果图像是i2d类型
cat = Table.read(img.replace('i2d.fits', cat_suffix)) # 读取对应的目录
cmap = 'gist_gray' # 设置颜色映射为灰度
else:
cmap = 'tab20c_r' # 设置颜色映射为分类颜色
# 绘制图像
with fits.open(img) as hdu: # 打开FITS文件
ax.imshow(hdu[1].data, vmin=0, vmax=0.3, origin='lower', cmap=cmap) # 显示图像数据
title = f"obs{hdu[0].header['OBSERVTN']} {hdu[0].header['PUPIL']}" # 获取标题信息
# 也绘制相关的目录
extended_sources = cat[cat['is_extended'] == 1] # 选择扩展源
point_sources = cat[cat['is_extended'] == 0] # 选择点源
for color, sources in zip(['darkred', 'black'], [extended_sources, point_sources]):
# 绘制源
ax.scatter(sources['xcentroid'], sources['ycentroid'], s=20, facecolors='None', edgecolors=color, alpha=0.9)
# 添加源标签
for i, source_num in enumerate(sources['label']):
ax.annotate(source_num,
(sources['xcentroid'][i]+0.5, sources['ycentroid'][i]+0.5),
fontsize=8,
color=color) # 在源位置添加标签
if 'i2d' in img: # 如果是i2d图像
ax.set_title(f"{title} combined image") # 设置标题为组合图像
else:
ax.set_title(f"{title} segmentation map") # 设置标题为分割图
# 放大到较小的区域
ax.set_xlim(xmin, xmax) # 设置x轴范围
ax.set_ylim(ymin, ymax) # 设置y轴范围
# 帮助使坐标轴不重叠；如果这不起作用，您也可以手动设置
plt.tight_layout() # 调整布局以避免重叠
return fig # 返回图形对象
# 调用函数绘制图像及其分割图
default_fig = plot_image_and_segmentation_map(image3_i2d, image3_segm)
ImViz#
与DS9类似，ImViz允许您交互式地查看这些图像及其对应的源目录。
imviz = Imviz() # 创建 Imviz 实例
viewer = imviz.default_viewer # 获取默认的查看器
# 绘制每个 i2d 图像
catalogs = [] # 用于绘制目录的列表
labels = [] # 用于绘制目录的标签列表
for img in image3_i2d: # 遍历每个 i2d 图像
print(f'Plotting: {img}') # 打印当前绘制的图像名称
# 获取观测编号和光学元件信息作为标签
label = f"obs{fits.getval(img, 'OBSERVTN')} {fits.getval(img, 'PUPIL')}"
with warnings.catch_warnings(): # 捕获警告
warnings.simplefilter('ignore') # 忽略警告
imviz.load_data(img, data_label=label) # 加载图像数据
# 保存信息以便后续绘制目录
catalogs.append(img.replace('i2d.fits', 'cat.ecsv')) # 生成目录文件名并添加到列表
labels.append(label) # 添加标签到列表
# 这将图像对齐以使用 WCS 坐标；
# 图像需要先加载，但在添加标记之前
linking = imviz.plugins['Orientation'] # 获取方向插件
linking.link_type = 'WCS' # 设置链接类型为 WCS
# 还要绘制相关的目录
# 由于在使用天空坐标时需要在 imviz 中链接，因此需要一个单独的循环
for catname, label in zip(catalogs, labels): # 遍历目录名称和标签
cat = Table.read(catname) # 读取目录文件
# 将表格格式化为 imviz 期望的格式
sky_coords = Table({'coord': [SkyCoord(ra=cat['sky_centroid'].ra.degree, # 创建天空坐标
dec=cat['sky_centroid'].dec.degree,
unit="deg")]})
viewer.marker = {'color': 'orange', 'alpha': 1, 'markersize': 20, 'fill': False} # 设置标记样式
viewer.add_markers(sky_coords, use_skycoord=True, marker_name=f"{label} catalog") # 添加标记到查看器
# 这将改变所有图像的拉伸
plotopt = imviz.plugins['Plot Options'] # 获取绘图选项插件
plotopt.select_all(viewers=True, layers=True) # 选择所有查看器和图层
plotopt.stretch_preset = '99.5%' # 设置拉伸预设为 99.5%
imviz.show() # 显示 Imviz 窗口
自定义成像处理流程运行#
图像3#
尝试编辑几个参数，并将结果与上述默认运行进行比较，最初针对单个文件。
当我们调用 image3 管道时，可以对管道的特定步骤进行修改。在这种情况下，我们将编辑管道的 source_catalog 和 tweakreg 步骤。关于可以调整的不同参数的解释，请参见下面的进一步信息，而默认值是该处列出的默认管道值和提供的参数参考文件的组合。
# 获取所有包含'image3'和'asn'的JSON文件，并按名称排序
image3_asns = np.sort(glob.glob('*image3*asn*.json'))
# 选择第二个asn文件进行处理
test_asn = image3_asns[1]
# 检查校准后的文件是否已经存在
asn_data = json.load(open(test_asn)) # 读取asn文件内容
# 构建i2d文件的路径
i2d_file = os.path.join(custom_run_image3, f"{asn_data['products'][0]['name']}_i2d.fits")
# 如果i2d文件已经存在，打印提示信息
if os.path.exists(i2d_file):
print(i2d_file, 'i2d file already exists.') # 输出文件已存在的信息
else:
# 调用image3管道，进行处理，并添加一些新的修改
cust_img3 = Image3Pipeline.call(test_asn,
steps={
'source_catalog': {'kernel_fwhm': 5.0, # 设置源目录的FWHM
'snr_threshold': 10.0, # 设置信噪比阈值
'npixels': 50, # 设置像素数量
'deblend': True, # 启用去混叠
},
'tweakreg': {'snr_threshold': 20, # 设置信噪比阈值
'abs_refcat': 'GAIADR3', # 设置绝对参考目录
'save_catalogs': True, # 保存目录
'searchrad': 3.0, # 设置搜索半径
'kernel_fwhm': 2.302, # 设置核的FWHM
'fitgeometry': 'shift', # 设置拟合几何
},
},
save_results=True, # 保存结果
output_dir=custom_run_image3) # 指定输出目录
# 导入必要的库
import os # 用于处理文件和目录
import glob # 用于查找符合特定规则的文件路径名
# 删除不必要的文件以节省磁盘空间
for crf in glob.glob(os.path.join(custom_run_image3, '*crf.fits')): # 查找所有以'crf.fits'结尾的文件
os.remove(crf) # 删除找到的文件
检查自定义结果#
# 生成默认的i2d文件路径，将默认运行图像路径与i2d文件名结合
default_i2d = os.path.join(default_run_image3, os.path.basename(i2d_file))
# 创建一个列表，包含待比较的i2d文件
compare_i2ds = [i2d_file, default_i2d]
# 创建一个列表，包含待比较的分割文件，替换文件名中的'i2d.fits'为'segm.fits'
compare_segm = [i2d_file.replace('i2d.fits', 'segm.fits'), default_i2d.replace('i2d.fits', 'segm.fits')]
# 调用函数绘制图像和分割图，并将比较结果存储在compare_fig中
compare_fig = plot_image_and_segmentation_map(compare_i2ds, compare_segm)
下面的单元格使用Imviz来可视化类似的信息。
imviz = Imviz() # 创建Imviz实例
viewer = imviz.default_viewer # 获取默认查看器
# 遍历图像文件和标签
for img, label in zip([i2d_file, os.path.join(default_run_image3, os.path.basename(i2d_file))], ['custom', 'default']):
print(f'Plotting: {img}') # 打印当前绘制的图像文件名
# 创建标题，包含观测编号和光学元件信息
title = f"{label} obs{fits.getval(img, 'OBSERVTN')} {fits.getval(img, 'PUPIL')}"
with warnings.catch_warnings(): # 捕获警告
warnings.simplefilter('ignore') # 忽略警告
imviz.load_data(img, data_label=title) # 加载数据到Imviz中
# 这将图像对齐以使用WCS坐标
linking = imviz.plugins['Orientation'] # 获取方向插件
linking.link_type = 'WCS' # 设置链接类型为WCS
# 还绘制相关的目录
cat = Table.read(img.replace('i2d.fits', 'cat.ecsv')) # 读取对应的目录文件
# 将表格格式化为imviz所需的格式
t_xy = Table({'x': cat['xcentroid'], # 提取x中心坐标
'y': cat['ycentroid']}) # 提取y中心坐标
viewer.marker = {'color': 'orange', 'alpha': 1, 'markersize': 20, 'fill': False} # 设置标记样式
viewer.add_markers(t_xy, marker_name=f"{label} catalog") # 添加标记到查看器中
# 这将改变所有图像的拉伸
plotopt = imviz.plugins['Plot Options'] # 获取绘图选项插件
plotopt.select_all(viewers=True, layers=True) # 选择所有查看器和图层
plotopt.stretch_preset = '99.5%' # 设置拉伸预设为99.5%
imviz.show() # 显示Imviz界面
如果您对上述结果满意，请对剩余图像进行校准。
# 获取所有包含'image3'和'asn'的JSON文件，并按名称排序
image3_asns = np.sort(glob.glob('*image3*asn*.json'))
# 遍历每个找到的'image3'的asn文件
for img3_asn in image3_asns:
# 检查校准后的文件是否已经存在
asn_data = json.load(open(img3_asn)) # 读取asn文件内容为JSON格式
# 构建输出文件的路径
i2d_file = os.path.join(custom_run_image3, f"{asn_data['products'][0]['name']}_i2d.fits")
# 如果输出文件已存在，打印提示并跳过该文件
if os.path.exists(i2d_file):
print(i2d_file, 'cal file already exists.') # 输出已存在的文件信息
continue # 跳过后续处理，继续下一个文件
# 调用image3处理管道，进行图像处理，并添加一些新的修改
cust_img3 = Image3Pipeline.call(img3_asn,
steps={ # 定义处理步骤
'source_catalog': {'kernel_fwhm': 5.0, # 源目录步骤的参数设置
'snr_threshold': 10.0, # 信噪比阈值
'npixels': 50, # 最小像素数
'deblend': True, # 是否去混叠
},
'tweakreg': {'snr_threshold': 20, # 调整注册步骤的信噪比阈值
'abs_refcat': 'GAIADR3', # 绝对参考目录
'save_catalogs': True, # 是否保存目录
'searchrad': 3.0, # 搜索半径
'kernel_fwhm': 2.302, # 核函数的全宽半高
'fitgeometry': 'shift', # 拟合几何形状
},
},
save_results=True, # 保存结果
output_dir=custom_run_image3) # 指定输出目录
# 导入所需的库
import os # 用于文件和目录操作
import glob # 用于文件路径匹配
# 移除不必要的文件以节省磁盘空间
for crf in glob.glob(os.path.join(custom_run_image3, '*crf.fits')): # 查找匹配的文件
os.remove(crf) # 删除找到的文件
# 这些都是来自Image3管道的结果
cust_image3_i2d = np.sort(glob.glob(os.path.join(custom_run_image3, '*i2d.fits'))) # 获取多个dither/mosaic合成图像的文件路径并排序
cust_image3_segm = np.sort(glob.glob(os.path.join(custom_run_image3, '*segm.fits'))) # 获取分割图的文件路径并排序，该图定义了源的范围
custom_fig = plot_image_and_segmentation_map(cust_image3_i2d, cust_image3_segm) # 绘制图像及其分割图
进一步完善源目录#
在上述情况下，我们直接使用管道（pipeline）修改了结果。也许您想使用其他人的自定义源目录（source catalog），或者进一步修改管道输出的源目录。在这些情况下，我们需要修改spec2 ASN文件，以指向新的源目录，这将在spec2笔记本中讨论。此外，匹配不同目录中的所有源ID（source ID）也是很有用的。在这种情况下，管道创建了三个不同的目录，它们识别相同的RA/Dec或X/Y像素位置，但源ID不同，因此我们将编辑这些目录，使它们在源ID值上保持一致。这些额外的步骤并不总是必要的，但在NIRISS WFSS数据的分析中可能会有所帮助。
在目录中匹配源ID#
在上面的图中，您可以看到同一个源有多个源ID。在这里，我们希望匹配所有观测中的源ID，以确保无论我们查看哪个滤光片或观测，都在讨论同一个源。为此，我们使用 astropy 的 match_to_catalog_3d 函数并重新基准化标签。
第一步是决定一个基础目录，以便与所有其他目录进行匹配。在这里，我们将使用观测004 F115W滤光片。
custom_cats = np.sort(glob.glob(os.path.join(custom_run_image3, '*niriss_clear-f????_cat.ecsv'))) # 获取所有符合格式的cat文件路径并排序
print("All image3 catalogs:\n", custom_cats) # 打印所有找到的image3目录
base_cat = Table.read(custom_cats[0]) # 读取第一个catalog文件作为基础目录
# save out the base catalog with a new name to be consistent
base_cat_name = custom_cats[0].replace('cat.ecsv', 'source-match_cat.ecsv') # 替换文件名以生成新的基础目录名称
base_cat.write(base_cat_name, overwrite=True) # 将基础目录写入新文件，允许覆盖
print("\nBase catalog:", base_cat_name) # 打印新基础目录的名称
遍历剩余的目录，以基于天球坐标匹配 ID，匹配精度为 1 角秒以内。
max_sep = 1 * u.arcsec # 设置最大分离角度为1弧秒，如果需要可以调整
base_sky = base_cat['sky_centroid'] # 获取基础目录中的天空中心坐标
for to_match_cat in custom_cats[1:]: # 遍历自定义目录中的每个目录，从第二个开始
# 读取目录
other_cat = Table.read(to_match_cat) # 读取当前匹配的目录
other_sky = other_cat['sky_centroid'] # 获取当前目录中的天空中心坐标
# 根据天空坐标找到两个目录之间的匹配源
idx, d2d, d3d = base_sky.match_to_catalog_3d(other_sky) # 进行3D匹配，返回索引和距离
sep_constraint = d2d < max_sep # 根据最大分离角度筛选匹配源
base_matches = base_cat[sep_constraint] # 获取基础目录中匹配的源
other_matches = other_cat[idx[sep_constraint]] # 获取其他目录中匹配的源
# 重新设置ID值，使其相同
other_matches['label'] = base_matches['label'] # 将匹配的标签从基础目录复制到其他目录
# 保存新的目录
match_cat_name = to_match_cat.replace('cat.ecsv', 'source-match_cat.ecsv') # 修改文件名以保存匹配目录
other_matches.write(match_cat_name, overwrite=True) # 写入匹配目录，覆盖同名文件
print('Saved:', match_cat_name) # 打印保存的文件名
查看新的源标签编号。它们应该匹配！
# 调用函数绘制图像和分割图，并保存结果
new_cat_fig = plot_image_and_segmentation_map(
cust_image3_i2d, # 输入的图像数据
cust_image3_segm, # 输入的分割数据
cat_suffix='source-match_cat.ecsv', # 输出文件的后缀名
xmin=1500, # x轴最小值
xmax=2000, # x轴最大值
ymin=800, # y轴最小值
ymax=1300 # y轴最大值
)
手动编辑源目录#
展望WFSS（宽场光谱仪）提取，我们可能在某一特定时刻只关心一到两个源。在这种情况下，缩减源目录至仅包含这些源会更有助于查看WFSS数据中提取的1-D光谱。
要开始，请查看当前某个滤光器的自定义源目录，以了解目录中包含的内容。
# 首先，查看F200W滤光器的当前自定义源目录
# 获取自定义运行目录中所有匹配的源目录文件，并按字母顺序排序，选择第三个文件
cat_name = np.sort(glob.glob(os.path.join(custom_run_image3, '*source-match_cat.ecsv')))[2]
# 读取源目录文件，存储为Table对象
cat = Table.read(cat_name)
# 打印源目录文件的名称
print(cat_name)
# 显示源目录的内容
cat
寻找源（source）可能有多种方法，以下是三种选项：
使用已知的赤经/赤纬（RA/Dec）来定位一个对象
使用已知的x/y位置来定位一个对象
使用对象的源ID（source ID）。请注意，我们在这里使用的是重新基准化的源目录（rebased source catalogs）中的ID。
这个概念可以扩展到通过源目录中包含的任何列进行过滤。
# 已知的赤经/赤纬
desired_ra = 53.15437048946369 # 目标赤经
desired_dec = -27.771689847051736 # 目标赤纬
# 创建一个SkyCoord对象，使用目标的赤经和赤纬
c = SkyCoord(ra=desired_ra*u.degree, dec=desired_dec*u.degree)
# 将目标坐标与目录中的天体坐标进行匹配，返回最近的天体ID和距离
nearest_id, distance_2d, distance_3d = c.match_to_catalog_sky(cat['sky_centroid'])
# 从目录中提取与最近天体ID对应的相关信息
cat[['label', 'xcentroid', 'ycentroid', 'sky_centroid', 'is_extended', 'isophotal_abmag', 'isophotal_vegamag']][nearest_id]
# 在F200W图像中使用已知的X/Y像素位置（基于你定义的cat）
known_x = 1880 # 定义已知的X坐标
known_y = 1100 # 定义已知的Y坐标
# 计算每个天体到已知位置的距离
nearest_pos = [np.sqrt((x-known_x)**2 + (y-known_y)**2) for x, y in zip(cat['xcentroid'], cat['ycentroid'])]
# 找到距离最近的天体的索引
wh_nearest = np.where(np.array(nearest_pos) == min(nearest_pos))[0][0]
# 输出距离最近的天体的相关信息
cat[['label', 'xcentroid', 'ycentroid', 'sky_centroid', 'is_extended', 'isophotal_abmag', 'isophotal_vegamag']][wh_nearest]
# alternatively with a known source number
# 另外一种方法是使用已知的源编号
# note that this number might be different for you depending on pipeline versions and parameters changed.
# 请注意，这个编号可能会因管道版本和参数更改而有所不同。
# source = 109 # if you want to specify the number
# source = 109 # 如果你想指定一个编号
source = cat['label'][nearest_id] # to match the one chosen in this notebook
# source = cat['label'][nearest_id] # 从分类目录中获取最近的源编号，以匹配本笔记本中选择的编号
wh_source = np.where(np.array(cat['label'] == source))[0][0]
# wh_source = np.where(np.array(cat['label'] == source))[0][0] # 找到与源编号匹配的索引位置
cat[['label', 'xcentroid', 'ycentroid', 'sky_centroid', 'is_extended', 'isophotal_abmag', 'isophotal_vegamag']][wh_source]
# 从分类目录中提取相关信息，包括标签、中心坐标、天空坐标、是否扩展、等光度绝对星等和等光度维根星等
使用上述三种方法中的任意一种，编写新的目录。
使用 RA/Dec:
new_cat = Table(cat[nearest_id])使用 x/y:
new_cat = Table(cat[wh_nearest])使用源 ID:
new_cat = Table(cat[wh_source])
new_cat = Table(cat[wh_source]) # 将行实例转换回数据框
# 保存新目录，使用唯一名称
new_cat_name = cat_name.replace('cat.ecsv', f'source{source}_cat.ecsv') # 替换文件名中的部分以生成新名称
new_cat.write(new_cat_name, overwrite=True) # 将新目录写入文件，允许覆盖
print('Saved:', new_cat_name) # 打印保存的文件名
一旦我们有了一个满意的更新源目录，我们就可以进入管道的 spec2 步骤。在运行 spec2 后，可能需要返回到这一步。让我们快速查看一下在接下来的笔记本中，我们将使用 spec2 提取的源。
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 6)) # 创建一个1行2列的子图，设置图形大小为9x6
img = cust_image3_i2d[-1] # 获取最后一幅图像
cat = Table.read(new_cat_name) # 读取新的目录数据
for ax in [ax1, ax2]: # 遍历两个子图
# plot the image # 绘制图像
with fits.open(img) as hdu: # 打开图像文件
ax.imshow(hdu[1].data, vmin=0, vmax=0.3, origin='lower') # 显示图像数据，设置颜色范围和原点位置
title = f"obs{hdu[0].header['OBSERVTN']} {hdu[0].header['PUPIL']}" # 获取观测和光学元件信息
# also plot the associated catalog # 也绘制相关的目录
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]): # 遍历颜色和源
# plotting the sources # 绘制源
ax.scatter(sources['xcentroid'], sources['ycentroid'], s=20, facecolors='None', edgecolors=color, alpha=0.9) # 绘制源的散点图
# adding source labels # 添加源标签
for i, source_num in enumerate(sources['label']): # 遍历源标签
ax.annotate(source_num, # 注释源标签
(sources['xcentroid'][i]+0.5, sources['ycentroid'][i]+0.5), # 设置标签位置
fontsize=8, # 设置字体大小
color=color) # 设置标签颜色
fig.suptitle("Speicifc Source to Extract with Spec2") # 设置整个图形的标题
# zooming in on a smaller region # 放大到一个较小的区域
ax2.set_xlim(known_x-50, known_x+50) # 设置x轴范围
ax2.set_ylim(known_y-50, known_y+50) # 设置y轴范围
fig.tight_layout() # 调整子图间距
