WFSS光谱 第1部分：合并和归一化1D光谱#
用例： 最优提取光栅光谱；红移测量；发射线图。简化版的 JDox 科学用例 # 33。
数据： JWST 模拟的 NIRISS 图像来自 MIRAGE，经过 JWST 校准管道 处理；星系团。
工具： specutils, astropy, pandas, emcee, lmfit, corner, h5py。
跨仪器： NIRSpec
文档： 本笔记本是 STScI 更大 后处理管道数据分析工具生态系统 的一部分。
引言#
本笔记本是关于 NIRISS WFSS 数据的一组 4 个笔记本中的第 2 个：
1D 最优提取，因为 JWST 管道仅提供盒子提取。最优提取提高了微弱源的光谱信噪比（S/N）。
合并并归一化 1D 光谱。
与模板进行交叉相关以获取红移。
空间分辨的发射线图。
本笔记本将首先在不同的抖动位置合并最优提取的 1D 光谱（来自第 1 个笔记本）。然后，合并的光谱将通过从直接图像（即宽带光度）估算的通量进行归一化。
本笔记本将从前一个笔记本中保存的最优提取的 1D 光谱开始，文件名类似于 “l3_nis_f115w_G150C_s00002_ndither0_1d_opt.fits”。来自每个抖动的各种一维光谱将被合并并保存到一个单一文件中，命名为 “l3_nis_f115w_G150C_s00004_combine_1d_opt.fits”。
注意： 该过程并不打算合并来自两种不同方向（C&R）的光谱，因为每种方向的光谱分辨率因源形态而异。
注意： 将光栅光谱归一化到每个宽带滤光器的光度将改善零点校准，这在多个滤光器的连续拟合中至关重要。
%matplotlib inline # 在Jupyter Notebook中内联显示Matplotlib图形
import os # 导入操作系统模块
import numpy as np # 导入NumPy库，用于数组和数学运算
from scipy.integrate import simpson # 从SciPy库导入辛普森积分方法
from urllib.request import urlretrieve # 从urllib库导入urlretrieve，用于下载文件
import tarfile # 导入tarfile模块，用于处理tar文件
from astropy.io import fits # 从Astropy库导入fits模块，用于处理FITS文件
import astropy.units as u # 导入Astropy单位模块
from astropy.io import ascii # 从Astropy库导入ascii模块，用于处理ASCII文件
import astropy # 导入Astropy库
print('astropy', astropy.__version__) # 打印Astropy库的版本
import matplotlib.pyplot as plt # 导入绘图库 matplotlib 的 pyplot 模块
import matplotlib as mpl # 导入 matplotlib 库
# 这些操作是为了确保图形在内联和笔记本版本中大小一致
%config InlineBackend.print_figure_kwargs = {'bbox_inches': None} # 配置内联后端的图形输出参数
mpl.rcParams['savefig.dpi'] = 80 # 设置保存图形时的分辨率为 80 dpi
mpl.rcParams['figure.dpi'] = 80 # 设置图形显示时的分辨率为 80 dpi
import specutils # 导入specutils库，用于光谱数据处理
from specutils.spectra.spectrum1d import Spectrum1D # 从specutils中导入一维光谱类
from astropy.nddata import StdDevUncertainty # 从astropy导入标准差不确定性类
print("Specutils: ", specutils.__version__) # 打印specutils库的版本
0. 下载笔记本 01 产品#
这些也可以通过运行笔记本获得。
if not os.path.exists('./output'): # 检查输出目录是否存在
import zipfile # 导入zipfile模块以处理zip文件
import urllib.request # 导入urllib.request模块以处理URL请求
# 定义要下载的zip文件的链接
boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/NIRISS_lensing_cluster/output.zip'
boxfile = './output.zip' # 定义下载后保存的zip文件路径
urllib.request.urlretrieve(boxlink, boxfile) # 从指定链接下载zip文件并保存到本地
zf = zipfile.ZipFile(boxfile, 'r') # 打开下载的zip文件以进行读取
list_names = zf.namelist() # 获取zip文件中所有文件的名称列表
for item in list_names: # 遍历zip文件中的每个文件
zf.extract(member=item) # 提取当前文件，使用extract比extractall更安全，避免绝对路径或相对路径问题
else:
print('Already exists') # 如果输出目录已存在，则打印提示信息
# 哪个数据集？
DIR_OUT = './output/' # 输出目录
filt = 'f200w' # 选择的滤光片
grism = 'G150C' # 选择的光栅
id = '00004' # 数据集的ID
ndither = 2 # 像素移动的数量
1. 合并来自不同偏移（dithers）的光谱；#
dithers = np.arange(0, ndither, 1) # 创建一个从0到ndither的数组，步长为1
for dither in dithers: # 遍历每个dither
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{dither}_1d_opt.fits' # 构建文件名
fd = fits.open(file_1d)[1].data # 打开FITS文件并读取数据
if dither == 0: # 如果是第一个dither
wave = fd['wavelength'] # 获取波长数据
flux = np.zeros((ndither, len(wave)), 'float') # 初始化flux数组，形状为(ndither, 波长长度)
flux_err = np.zeros((ndither, len(wave)), 'float') # 初始化flux_err数组，形状为(ndither, 波长长度)
flux[dither, :] = fd['flux'] # 将第一个dither的flux数据存入数组
flux_err[dither, :] = fd['uncertainty'] # 将第一个dither的flux误差数据存入数组
else: # 如果不是第一个dither
wave_tmp = fd['wavelength'] # 获取当前dither的波长数据
flux_tmp = fd['flux'] # 获取当前dither的flux数据
flux_err_tmp = fd['uncertainty'] # 获取当前dither的flux误差数据
flux[dither, :] = np.interp(wave, wave_tmp, flux_tmp) # 对flux进行插值
flux_err[dither, :] = np.sqrt(np.interp(wave, wave_tmp, flux_err_tmp[:]**2)) # 对flux误差进行插值并开方
plt.plot(wave, flux[0, :], color='r', label='dither1') # 绘制第一组数据，使用红色表示
plt.plot(wave, flux[1, :], color='b', label='dither2') # 绘制第二组数据，使用蓝色表示
plt.legend(loc=2) # 添加图例，位置设置为左上角
对来自不同抖动位置的两个光谱进行加权平均；#
flux_combine = np.zeros(len(wave), 'float') # 初始化一个与波长相同长度的数组，用于存储组合后的flux值
flux_err_combine = np.zeros(len(wave), 'float') # 初始化一个与波长相同长度的数组，用于存储组合后的flux误差值
for ii in range(len(wave)): # 遍历每个波长点
# 计算组合后的flux值，使用加权平均法，权重为flux误差的倒数平方
flux_combine[ii] = np.sum(flux[:, ii] * 1 / flux_err[:, ii]**2) / np.sum(1 / flux_err[:, ii]**2)
# 计算组合后的flux误差，取所有flux误差的平方和的平方根，然后除以样本数量
flux_err_combine[ii] = np.sqrt(np.sum(flux_err[:, ii]**2)) / len(flux_err[:, ii])
# 绘制第一个dither的flux曲线，颜色为灰色
plt.plot(wave, flux[0, :], color='gray', label='dither1')
# 绘制第二个dither的flux曲线，颜色为灰色
plt.plot(wave, flux[1, :], color='gray', label='dither2')
# 绘制组合后的flux曲线，颜色为红色
plt.plot(wave, flux_combine[:], color='r', label='Combined')
# 添加图例，位置为左上角
plt.legend(loc=2)
2. 连续体归一化 (Continuum normalization)#
由于 NIRISS 光谱可能受到背景减法 (background subtraction) 和污染 (contamination) 的影响，这可能导致滤光片 (filters) 之间的不匹配，因此我们旨在将光谱 (spectra) 从滤光片归一化到其宽带光度 (broadband magnitude)。
# 从Notebook 01a中打开宽带通量目录；
# 通量已经以Fnu形式给出，magzp = 25.0;
file = DIR_OUT + 'l3_nis_flux.cat' # 定义文件路径
fd_cat = ascii.read(file) # 读取目录文件
id_cat = fd_cat['id'] # 提取id列
magzp = 25.0 # 设置零点星等
fd_cat # 显示目录内容
# 获取NIRISS图像的滤波器曲线;
url = 'https://jwst-docs.stsci.edu/files/97978948/97978953/1/1596073225227/NIRISS_Filters_2017May.tar.gz' # 定义要下载的文件URL
filename = 'tmp.tar.gz' # 定义下载后保存的文件名
urlretrieve(url, filename) # 从指定URL下载文件并保存为filename
my_tar = tarfile.open(filename) # 打开下载的tar.gz文件
my_tar.extractall('./', filter='data') # 解压缩文件到当前目录，使用filter='data'以跳过绝对路径或相对路径
# 加载滤光器响应数据：
DIR_FIL = './NIRISS_Filters/DATA/' # 定义滤光器数据文件夹路径
# 数字对应EAZY滤光器响应中的F200W；
eazy_filt = 311 # 设置EAZY滤光器编号
# 读取传输数据；
filt_data = ascii.read(f'{DIR_FIL}/NIRISS_{filt.upper()}.txt') # 从指定路径读取滤光器数据
print(filt_data) # 打印读取的数据
wave_filt = filt_data['Wavelength'] # 获取波长数据
flux_filt = filt_data['FilterTrans'] # 获取滤光器传输数据
plt.plot(wave_filt, flux_filt, ls='-', label=f'NIRISS {filt.upper()}') # 绘制波长与响应的曲线
plt.xlabel('wavelength') # 设置x轴标签为'波长'
plt.ylabel('response') # 设置y轴标签为'响应'
plt.legend(loc=1) # 显示图例，位置在右上角
# 定义一个用于滤波器的通量卷积的小函数
def filconv(lfil, ffil, l0, f0, DIR='', c=3e18):
'''
lfil : 滤波器响应曲线的波长数组。
ffil : 滤波器响应曲线的通量数组。
l0 : 光谱的波长数组，单位为 f_nu，而不是 f_lam
f0 : 光谱的通量数组，单位为 f_nu，而不是 f_lam
'''
fhalf = np.max(ffil)/2.0 # 计算滤波器响应的半最大值
con = (ffil > fhalf) # 找到响应大于半最大值的区域
lfwhml = np.min(lfil[con]) # 计算左半宽的波长
lfwhmr = np.max(lfil[con]) # 计算右半宽的波长
lcen = (lfwhmr + lfwhml)/2. # 计算中心波长
lamS, spec = l0, f0 # 两列数据：波长和通量密度
lamF, filt = lfil, ffil # 两列数据：波长和响应范围 [0,1]
filt_int = np.interp(lamS, lamF, filt) # 将滤波器插值到共同的（光谱）波长轴上
if len(lamS) > 0: # 检查光谱波长数组是否为空
I1 = simpson(spec / lamS**2 * c * filt_int * lamS, x=lamS) # 计算 Fnu 的分母
I2 = simpson(filt_int/lamS, x=lamS) # 计算分子
fnu = I1/I2/c # 计算平均通量密度
else:
I1 = 0 # 如果光谱波长数组为空，分母为0
I2 = 0 # 如果光谱波长数组为空，分子为0
fnu = 0 # 如果光谱波长数组为空，平均通量密度为0
return lcen, fnu # 返回中心波长和平均通量密度
# 对观测光谱进行卷积处理，使用滤光器响应。
# 输入的通量需要以Fnu为单位，且magzp=25；
# 截取观测通量的边缘部分；
con = (wave > 1.75) & (wave < 2.25) # 选择波长在1.75到2.25之间的部分
lcen, fnu = filconv(wave_filt, flux_filt, wave[con], flux_combine[con], DIR='./') # 调用filconv函数进行卷积，计算中心波长和总通量
print('Central wavelength and total flux are ;', lcen, fnu) # 输出中心波长和总通量
# 获取归一化因子;
iix = np.where(id_cat[:] == int(id)) # 找到与给定id匹配的索引
Cnorm = fd_cat[f'F{eazy_filt}'][iix] / fnu # 计算归一化因子
Cnorm # 返回归一化因子
# 写入归一化的光谱：
# 将其转换为Spectrum1D实例。
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{dither}_1d_opt.fits' # 定义输出文件名
# 创建Spectrum1D实例，包含光谱轴和归一化的光谱数据
obs = Spectrum1D(spectral_axis=wave*u.um, # 光谱轴，单位为微米
flux=flux_combine * Cnorm * u.dimensionless_unscaled, # 归一化的光谱数据
uncertainty=StdDevUncertainty(flux_err_combine * Cnorm), # 归一化的光谱误差
unit='') # 单位为空
# 将光谱数据写入文件，格式为tabular-fits，若文件已存在则覆盖
obs.write(file_1d, format='tabular-fits', overwrite=True)
对其他滤光片重复步骤1和2；#
DIR_FIL = './NIRISS_Filters/DATA/' # 数据目录
filts = ['f115w', 'f150w', 'f200w'] # NIRISS滤光片列表
# 对应EAZY滤光响应的NIRISS滤光片编号
eazy_filts = [309, 310, 311]
# 边缘问题的通量掩码
mask_lw = [1.05, 1.35, 1.75] # 低波长掩码
mask_uw = [1.25, 1.65, 2.23] # 高波长掩码
for ff, filt in enumerate(filts): # 遍历每个滤光片
for dither in dithers: # 遍历每个抖动
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{dither}_1d_opt.fits' # 生成文件名
fd = fits.open(file_1d)[1].data # 打开FITS文件并读取数据
if dither == 0: # 如果是第一个抖动
wave = fd['wavelength'] # 获取波长数据
flux = np.zeros((ndither, len(wave)), 'float') # 初始化通量数组
flux_err = np.zeros((ndither, len(wave)), 'float') # 初始化通量误差数组
flux[dither, :] = fd['flux'] # 存储通量数据
flux_err[dither, :] = fd['uncertainty'] # 存储通量误差数据
else: # 如果不是第一个抖动
wave_tmp = fd['wavelength'] # 获取临时波长数据
flux_tmp = fd['flux'] # 获取临时通量数据
flux_err_tmp = fd['uncertainty'] # 获取临时通量误差数据
# 使用插值方法填充通量和通量误差
flux[dither, :] = np.interp(wave, wave_tmp, flux_tmp)
flux_err[dither, :] = np.sqrt(np.interp(wave, wave_tmp, flux_err_tmp[:]**2))
# 合并通量和误差
flux_combine = np.zeros(len(wave), 'float') # 初始化合并通量数组
flux_err_combine = np.zeros(len(wave), 'float') # 初始化合并通量误差数组
for ii in range(len(wave)): # 遍历每个波长点
# 计算加权平均通量
flux_combine[ii] = np.sum(flux[:, ii] * 1 / flux_err[:, ii]**2) / np.sum(1 / flux_err[:, ii]**2)
# 计算合并通量误差
flux_err_combine[ii] = np.sqrt(np.sum(flux_err[:, ii]**2)) / len(flux_err[:, ii])
# 归一化
filt_data = ascii.read(f'{DIR_FIL}/NIRISS_{filt.upper()}.txt') # 读取滤光片数据
wave_filt = filt_data['Wavelength'] # 获取滤光片波长
flux_filt = filt_data['FilterTrans'] # 获取滤光片透过率
con = (wave > mask_lw[ff]) & (wave < mask_uw[ff]) # 选择有效波长范围
lcen, fnu = filconv(wave_filt, flux_filt, wave[con], flux_combine[con]) # 进行卷积计算
iix = np.where(id_cat[:] == int(id)) # 查找对应的ID索引
Cnorm = fd_cat[f'F{eazy_filts[ff]}'][iix] / fnu # 计算归一化因子
# 写入结果
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits' # 生成输出文件名
# 创建光谱对象并写入文件
obs = Spectrum1D(spectral_axis=wave*u.um,
flux=flux_combine * Cnorm * u.dimensionless_unscaled,
uncertainty=StdDevUncertainty(flux_err_combine * Cnorm), unit='')
obs.write(file_1d, format='tabular-fits', overwrite=True) # 写入FITS文件，覆盖已有文件
结果;#
光谱彼此对齐，因为它们与黑体（BB）光度匹配。
filts = ['f115w', 'f150w', 'f200w'] # 定义滤光片列表
cols = ['lightblue', 'orange', 'purple'] # 定义颜色列表
for ff, filt in enumerate(filts): # 遍历滤光片列表及其索引
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits' # 构建文件名
fd = fits.open(file_1d)[1].data # 打开FITS文件并读取数据
wave = fd['wavelength'] # 提取波长数据
flux = fd['flux'] # 提取光通量数据
flux_err = fd['uncertainty'] # 提取不确定性数据
con = (wave > mask_lw[ff]) & (wave < mask_uw[ff]) # 创建波长范围的掩膜
plt.plot(wave[con], flux[con], ls='-', label=f'{grism} {filts[ff]}', color=cols[ff]) # 绘制光通量曲线
filt_data = ascii.read(f'{DIR_FIL}/NIRISS_{filt.upper()}.txt') # 读取滤光片传输数据
wave_filt = filt_data['Wavelength'] # 提取滤光片波长
flux_filt = filt_data['FilterTrans'] # 提取滤光片透过率
con = (wave > mask_lw[ff]) & (wave < mask_uw[ff]) # 再次创建波长范围的掩膜
lcen, fnu = filconv(wave_filt, flux_filt, wave[con], flux[con]) # 进行滤波卷积计算
iix = np.where(id_cat[:] == int(id)) # 找到对应的ID索引
plt.scatter(lcen, fd_cat[f'F{eazy_filts[ff]}'][iix], marker='s', s=50, edgecolor='k', color=cols[ff], label=f'BB {filts[ff]}') # 绘制散点图
plt.xlabel('wavelength') # 设置x轴标签为波长
plt.ylabel('flux') # 设置y轴标签为光通量
plt.ylim(0, 60) # 设置y轴范围
plt.legend(loc=0) # 显示图例
对另一个目标重复操作；#
id = '00003' # 定义目标ID
filts = ['f115w', 'f150w', 'f200w'] # 定义滤光片列表
eazy_filts = [309, 310, 311] # 定义EAZY滤光片索引
DIR_FIL = './NIRISS_Filters/DATA/' # 定义数据目录
# 边缘问题的通量掩膜;
mask_lw = [1.05, 1.35, 1.75] # 定义低波长掩膜
mask_uw = [1.25, 1.65, 2.23] # 定义高波长掩膜
for ff, filt in enumerate(filts): # 遍历每个滤光片
for dither in dithers: # 遍历每个抖动
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{dither}_1d_opt.fits' # 定义文件名
fd = fits.open(file_1d)[1].data # 打开FITS文件并读取数据
if dither == 0: # 如果是第一个抖动
wave = fd['wavelength'] # 获取波长数据
flux = np.zeros((ndither, len(wave)), 'float') # 初始化通量数组
flux_err = np.zeros((ndither, len(wave)), 'float') # 初始化通量误差数组
flux[dither, :] = fd['flux'] # 存储通量数据
flux_err[dither, :] = fd['uncertainty'] # 存储通量误差数据
else: # 如果不是第一个抖动
wave_tmp = fd['wavelength'] # 获取临时波长数据
flux_tmp = fd['flux'] # 获取临时通量数据
flux_err_tmp = fd['uncertainty'] # 获取临时通量误差数据
flux[dither, :] = np.interp(wave, wave_tmp, flux_tmp) # 插值计算通量
flux_err[dither, :] = np.sqrt(np.interp(wave, wave_tmp, flux_err_tmp[:]**2)) # 插值计算通量误差
# 合并通量;
flux_combine = np.zeros(len(wave), 'float') # 初始化合并通量数组
flux_err_combine = np.zeros(len(wave), 'float') # 初始化合并通量误差数组
for ii in range(len(wave)): # 遍历每个波长
flux_combine[ii] = np.sum(flux[:, ii] * 1 / flux_err[:, ii]**2) / np.sum(1 / flux_err[:, ii]**2) # 计算加权合并通量
flux_err_combine[ii] = np.sqrt(np.sum(flux_err[:, ii]**2)) / len(flux_err[:, ii]) # 计算合并通量误差
# 归一化;
filt_data = ascii.read(f'{DIR_FIL}/NIRISS_{filt.upper()}.txt') # 读取滤光片数据
wave_filt = filt_data['Wavelength'] # 获取滤光片波长
flux_filt = filt_data['FilterTrans'] # 获取滤光片透过率
con = (wave > mask_lw[ff]) & (wave < mask_uw[ff]) # 选择有效波长范围
lcen, fnu = filconv(wave_filt, flux_filt, wave[con], flux_combine[con]) # 进行滤波卷积
iix = np.where(id_cat[:] == int(id)) # 查找目标ID的索引
Cnorm = fd_cat[f'F{eazy_filts[ff]}'][iix] / fnu # 计算归一化因子
# 写入文件:
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits' # 定义输出文件名
obs = Spectrum1D(spectral_axis=wave*u.um, # 创建光谱数据对象
flux=flux_combine * Cnorm * u.dimensionless_unscaled, # 设置通量
uncertainty=StdDevUncertainty(flux_err_combine * Cnorm), unit='') # 设置不确定性
obs.write(file_1d, format='tabular-fits', overwrite=True) # 写入FITS文件，允许覆盖
filts = ['f115w', 'f150w', 'f200w'] # 定义滤光片列表
cols = ['lightblue', 'orange', 'purple'] # 定义颜色列表
for ff, filt in enumerate(filts): # 遍历滤光片列表及其索引
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits' # 构建文件名
fd = fits.open(file_1d)[1].data # 打开FITS文件并读取数据
wave = fd['wavelength'] # 提取波长数据
flux = fd['flux'] # 提取光通量数据
flux_err = fd['uncertainty'] # 提取不确定性数据
con = (wave > mask_lw[ff]) & (wave < mask_uw[ff]) # 创建波长范围的掩膜
plt.plot(wave[con], flux[con], ls='-', label=f'{grism} {filts[ff]}', color=cols[ff]) # 绘制光通量曲线
filt_data = ascii.read(f'{DIR_FIL}/NIRISS_{filt.upper()}.txt') # 读取滤光片传输数据
wave_filt = filt_data['Wavelength'] # 提取滤光片波长
flux_filt = filt_data['FilterTrans'] # 提取滤光片传输率
con = (wave > mask_lw[ff]) & (wave < mask_uw[ff]) # 创建波长范围的掩膜
lcen, fnu = filconv(wave_filt, flux_filt, wave[con], flux[con]) # 进行滤波卷积计算
iix = np.where(id_cat[:] == int(id)) # 找到与当前ID匹配的索引
plt.scatter(lcen, fd_cat[f'F{eazy_filts[ff]}'][iix], marker='s', s=50, edgecolor='k', color=cols[ff], label=f'BB {filts[ff]}') # 绘制散点图
plt.xlabel('wavelength') # 设置x轴标签为波长
plt.ylabel('flux') # 设置y轴标签为光通量
plt.ylim(0, 60) # 设置y轴范围
plt.legend(loc=0) # 显示图例
看起来您没有提供任何代码。如果您能提供需要添加注释的Python代码，我将很乐意帮助您添加行级中文注释。请将代码粘贴到这里。