WFSS光谱 第0部分：最佳提取#
用例： 最优提取光谱；红移测量；发射线图。简化版的 JDox 科学用例 # 33。
数据： JWST 模拟的 NIRISS 图像来自 MIRAGE，经过 JWST 校准管道 处理；星系团。
工具： specutils, astropy, pandas, emcee, lmfit, corner, h5py。
跨仪器： NIRSpec
文档： 本笔记本是 STScI 更大 后处理管道数据分析工具生态系统 的一部分。
引言#
本笔记本是关于 NIRISS WFSS 数据的一组 4 个笔记本中的第 1 个：
1D 最优提取，因为 JWST 管道仅提供盒式提取。最优提取提高了微弱源的光谱信噪比（S/N）。
合并并归一化 1D 光谱。
将星系与模板进行交叉相关以获取红移。
空间分辨的发射线图。
本笔记本将从 NIRISS WFSS 的 后处理管道 产品开始，使用来自 spec level3 的 2D 矫正光谱。
最优提取需要在交叉色散方向上的源形态，这将通过与 WFSS 观测同时拍摄的直接图像获取。色散方向上的形态对于推断光谱分辨率也至关重要，这将通过模板拟合来获取红移和星族信息，在本组的第 3 个笔记本中进行。
注意： 我们假设 2D 矫正光谱的处理已经达到合理水平，即目标 2D 光谱上没有来自其他源的污染光通量，并且背景已经被减去。
%matplotlib inline # 在Jupyter Notebook中内联显示Matplotlib绘图
import os # 导入操作系统模块
import numpy as np # 导入NumPy库，通常用于数组和数值计算
from scipy.ndimage import rotate # 从SciPy库导入旋转函数，用于图像处理
from scipy.optimize import curve_fit # 从SciPy库导入曲线拟合函数
from astropy.io import fits # 从Astropy库导入FITS文件处理模块
import astropy.units as u # 导入Astropy单位模块，用于物理量的单位处理
import astropy.wcs as wcs # 导入Astropy的WCS模块，用于天文坐标转换
from astropy.io import ascii # 从Astropy库导入ASCII文件处理模块
from specutils import Spectrum1D # 从specutils库导入一维光谱类
from astropy.nddata import StdDevUncertainty # 从Astropy库导入标准偏差不确定性类
import specutils # 导入specutils库
print('specutils', specutils.__version__) # 打印specutils库的版本
import astropy # 导入Astropy库
print('astropy', astropy.__version__) # 打印Astropy库的版本
版本应为#
specutils 1.11.0
astropy 5.3.4
import matplotlib.pyplot as plt # 导入绘图库 matplotlib 的 pyplot 模块
import matplotlib as mpl # 导入 matplotlib 库
# 设置保存图像时的分辨率为 80 dpi
mpl.rcParams['savefig.dpi'] = 80
# 设置图像显示时的分辨率为 80 dpi
mpl.rcParams['figure.dpi'] = 80
0. 下载和加载数据：#
这些包括NIRISS的管道处理数据，以及来自图像步骤3的光度目录。
if not os.path.exists('./pipeline_products'): # 检查'./pipeline_products'目录是否存在
import zipfile # 导入zipfile模块，用于处理zip文件
import urllib.request # 导入urllib.request模块，用于下载文件
# 定义要下载的zip文件的链接
boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/NIRISS_lensing_cluster/pipeline_products.zip'
boxfile = './pipeline_products.zip' # 定义下载后保存的zip文件路径
# 下载zip文件
urllib.request.urlretrieve(boxlink, boxfile)
zf = zipfile.ZipFile(boxfile, 'r') # 打开下载的zip文件
list_names = zf.namelist() # 获取zip文件内所有文件的名称列表
for item in list_names: # 遍历zip文件内的每个文件
zf.extract(member=item) # 使用extract提取文件，避免使用extractall以防止路径问题（绝对路径或相对路径）
import os # 导入操作系统模块
from astropy.io import fits # 导入Astropy库中的fits模块，用于处理FITS文件
from astropy import wcs # 导入Astropy库中的wcs模块，用于处理天文坐标系统
import numpy as np # 导入NumPy库，用于数值计算
DIR_DATA = './pipeline_products/' # 数据目录路径
# 输出目录;
DIR_OUT = './output/' # 输出目录路径
if not os.path.exists(DIR_OUT): # 检查输出目录是否存在
os.mkdir(DIR_OUT) # 如果不存在，则创建输出目录
# 检测和科学图像的滤波器;
filt_det = 'f200w' # 设置用于检测的滤波器
# 从直接图像中获取图像数组。用于最佳提取和掩膜。
# 此图像应已进行天空减法；否则，您将遇到最佳提取的错误结果。
infile = f'{DIR_DATA}l3_nis_{filt_det}_i2d_skysub.fits' # 输入文件路径，包含天空减法后的数据
hdu = fits.open(infile) # 打开FITS文件
# 这只是用于误差数组;
infile = f'{DIR_DATA}l3_nis_{filt_det}_i2d.fits' # 输入文件路径，包含原始数据
hdu_err = fits.open(infile) # 打开误差FITS文件
data = hdu[0].data # 获取图像数据
imwcs = wcs.WCS(hdu[0].header, hdu) # 创建WCS对象，用于处理图像的天文坐标
err = hdu_err[2].data # 获取误差数据
weight = 1/np.square(err) # 计算权重，权重为误差的平方的倒数
# 分割图
# 如果管道没有生成分割图，可以通过运行Photutils来准备。
segfile = f'{DIR_DATA}l3_nis_{filt_det}_i2d_seg.fits' # 输入分割图文件路径
seghdu = fits.open(segfile) # 打开分割图FITS文件
segdata = seghdu[0].data # 获取分割图数据
# 从图像的第三级加载目录;
# 以获取源在像素坐标中的位置。
catfile = f'{DIR_DATA}l3_nis_{filt_det}_cat.ecsv' # 构建目录文件的路径
fd = ascii.read(catfile) # 读取目录文件
看起来您输入的内容不完整。如果您有特定的Python代码需要我添加中文注释，请提供代码段，我将很乐意帮助您。
制作宽带通量目录。#
出于方便的考虑，我们将目录汇编成一个通量目录，该目录将在后续笔记本（01b）中使用。
要运行此单元，您需要一个源的光度目录，该目录列出了每个滤光片的源位置和通量。目前，我使用的是在另一个笔记本中准备的目录。（”sources_extend_01.cat”）
此目录也可以用于通用的光度红移/光谱能量分布拟合软件，如EAZY和gsf（见笔记本No.04）。
目前，我们使用一个输入目录。#
filts = ['f115w', 'f150w', 'f200w', 'f090w', 'f435w', 'f606w', 'f814w', 'f105w', 'f125w', 'f140w', 'f160w'] # 定义滤光片名称列表
eazy_filts = [309, 310, 311, 308, 1, 4, 6, 202, 203, 204, 205] # 定义EAZY滤光片对应的索引
magzp = 25.0 # 目录中的星等零点
# 读取目录;
fd_input = ascii.read(f'{DIR_DATA}sources_extend_01.cat') # 从指定路径读取源数据文件
ra_input = fd_input['x_or_RA'] # 获取输入数据中的RA（右升交角）
dec_input = fd_input['y_or_Dec'] # 获取输入数据中的Dec（天赤纬）
# 写入文件头;
fw = open(f'{DIR_OUT}l3_nis_flux.cat', 'w') # 打开输出文件以写入
fw.write('# id') # 写入文件头的ID
for ff in range(len(filts)): # 遍历所有滤光片
fw.write(f' F{eazy_filts[ff]} E{eazy_filts[ff]}') # 写入每个滤光片的EAZY索引
fw.write('\n') # 换行
# 写入内容;
for ii in range(len(fd['id'])): # 遍历每个源的ID
# 计算当前源与所有输入源之间的距离
rtmp = np.sqrt((fd['sky_centroid'].ra.value[ii] - ra_input[:])**2 + (fd['sky_centroid'].dec.value[ii] - dec_input[:])**2)
iix = np.argmin(rtmp) # 找到距离最近的输入源索引
for ff in range(len(filts)): # 遍历所有滤光片
if ff == 0: # 如果是第一个滤光片
fw.write(f"{fd['id'][ii]}") # 写入源的ID
mag = fd_input[f'niriss_{filts[ff]}_magnitude'][iix] # 获取最近源的星等
flux_nu = 10**((mag - magzp) / (-2.5)) # 计算对应的光通量
# 当前目录未提供适当的误差;
# 假设光通量的误差为5%.
flux_err_nu = flux_nu * 0.05 # 计算光通量误差
fw.write(f' {flux_nu:.5e} {flux_err_nu:.5e}') # 写入光通量及其误差
fw.write('\n') # 换行
fw.close() # 关闭输出文件
1. 加载 2D 光谱;#
# 选择滤光片、光栅和目标对象
filt = 'f200w' # 设置滤光片为 'f200w'
# grism = 'G150R' # 原本的光栅设置为 'G150R'
grism = 'G150C' # 设置光栅为 'G150C'
id = '00004' # 设置目标对象的ID
# 零索引的抖动编号 --- 测试数据有两个抖动位置，因此可以是0或1。
ndither = 0 # 设置抖动位置为0
# 构建文件路径，包含数据目录、滤光片、光栅和目标ID
file_2d = f'{DIR_DATA}l3_nis_{filt}_{grism}_s{id}_cal.fits'
# 打开FITS文件
hdu_2d = fits.open(file_2d)
# 对齐光栅方向
# - x方向 = 色散（波长）方向。
# - y方向 = 交叉色散。
# 在这个笔记本中。
if grism == 'G150C': # 如果光栅是 'G150C'
# 如果光谱是水平的；
data_2d = hdu_2d[ndither*7+1].data # 获取二维数据
dq_2d = hdu_2d[ndither*7+2].data # 获取数据质量数组
err_2d = hdu_2d[ndither*7+3].data # 获取误差数组
wave_2d = hdu_2d[ndither*7+4].data # 获取波长数组
else: # 如果光栅不是 'G150C'
data_2d = rotate(hdu_2d[ndither*7+1].data, 90) # 将数据旋转90度
dq_2d = rotate(hdu_2d[ndither*7+2].data, 90) # 将数据质量数组旋转90度
err_2d = rotate(hdu_2d[ndither*7+3].data, 90) # 将误差数组旋转90度
wave_2d = rotate(hdu_2d[ndither*7+4].data, 90) # 将波长数组旋转90度
# 获取观测的方位角；
hd_2d = hdu_2d[1].header # 获取FITS文件的头信息
PA_V3 = hd_2d['PA_V3'] # 从头信息中提取方位角PA_V3
PA_V3 # 输出方位角
plt.imshow(data_2d, vmin=0, vmax=1000) # 显示二维数据，设置最小值为0，最大值为1000
# 获取源的光谱轮廓；
# 再次，y是交叉色散方向，x是色散方向。
y2d, x2d = data_2d.shape[:] # 获取2D数据的形状，y和x的维度
# 切割出分割图；
iix = np.where(fd['id'] == int(id))[0][0] # 找到与给定id匹配的索引
# 从图像3目录中获取目标位置；
ycen = fd['ycentroid'][iix].value # 获取目标的y中心坐标
xcen = fd['xcentroid'][iix].value # 获取目标的x中心坐标
# 切割大小 = 2D光谱的y方向；
rsq = y2d # 设置切割大小为y方向的维度
# 从数据中切割出科学图像和分割图；
sci_cut = data[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)] # 切割科学图像
seg_cut = segdata[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)] # 切割分割图
# 根据Grism观测的PA旋转图像；
if grism == 'G150C': # 如果使用G150C光栅
sci_rot = rotate(sci_cut, PA_V3) # 旋转科学图像
else: # 否则
sci_rot = rotate(sci_cut, PA_V3+90) # 旋转科学图像，增加90度
WFSS光栅在下图中沿x轴方向分散。#
plt.imshow(sci_rot, vmin=0, vmax=1.0) # 显示旋转后的科学图像，设置最小值为0，最大值为1.0
plt.title('Direct image') # 设置图像标题为“直接图像”
plt.xlabel('Wavelength direction >>>', color='r', fontsize=18) # 设置x轴标签为“波长方向 >>>”，字体颜色为红色，字体大小为18
plt.ylabel('Cross-dispersion direction >>>', color='r', fontsize=18) # 设置y轴标签为“交叉色散方向 >>>”，字体颜色为红色，字体大小为18
2. 在不同的 x 位置获取光谱轮廓 — 这将用于最佳提取。#
for ii in range(sci_rot.shape[1]): # 遍历sci_rot的每一列
flux_tmp = sci_rot[:, ii] # 获取当前列的光谱数据
xx_tmp = np.arange(0, len(sci_rot[:, ii]), 1) # 创建一个与光谱数据长度相同的数组，用于x轴
plt.plot(xx_tmp, flux_tmp, label='x={}'.format(ii)) # 绘制当前列的光谱数据，并设置标签
plt.legend(loc=0, fontsize=8) # 显示图例，位置为0，字体大小为8
# 沿着 x 方向（分散方向）求和
flux_y = np.zeros(len(sci_rot[:, 0]), 'float') # 创建一个与 sci_rot 行数相同的零数组，用于存储 y 方向的 flux
for ii in range(sci_rot.shape[0]): # 遍历 sci_rot 的每一行
flux_y[ii] = np.sum(sci_rot[ii, :]) # 对每一行的所有列求和，并存储到 flux_y 中
# 如果需要，进行天空背景减法。
# sky = np.mean([flux_y[0], flux_y[-1]]) # 计算 flux_y 的首尾平均值，作为天空背景（注释掉）
# 归一化处理；
flux_y[:] /= flux_y.sum() # 将 flux_y 中的每个元素除以其总和，实现归一化
plt.plot(xx_tmp, flux_y) # 绘制 y 位置与源通量的关系图
plt.xlabel('y-position') # 设置 x 轴标签为 'y-position'
plt.ylabel('Source Flux') # 设置 y 轴标签为 'Source Flux'
3. 一维提取;#
显示管道 1D 提取作为示例；#
# 正常提取;
flux_disp1 = np.zeros(x2d, 'float') # 初始化光谱通量数组
err_disp1 = np.zeros(x2d, 'float') # 初始化误差数组
wave_disp1 = np.zeros(x2d, 'float') # 初始化波长数组
for ii in range(x2d): # 在波长方向上循环
mask_tmp = (dq_2d[:, ii] == 0) & (err_2d[:, ii] > 0) # 创建掩膜，选择有效数据
# 在一个盒子内求和;
flux_disp1[ii] = np.sum(data_2d[:, ii][mask_tmp]) # 计算有效光谱通量的总和
err_disp1[ii] = np.sqrt(np.sum(err_2d[:, ii][mask_tmp]**2)) # 计算有效误差的平方和的平方根
wave_disp1[ii] = wave_2d[0, ii] # 获取对应的波长值
plt.errorbar(wave_disp1, flux_disp1, yerr=err_disp1) # 绘制带误差条的光谱图
plt.xlim(1.7, 2.3) # 设置x轴范围
最优提取;#
# 根据 Horne(1986, PASP, 98, 609) 的方法进行处理
flux_disp = np.zeros(x2d, 'float') # 初始化光谱通量数组
err_disp = np.zeros(x2d, 'float') # 初始化误差数组
wave_disp = np.zeros(x2d, 'float') # 初始化波长数组
# Sigma clipping（标准差剪切）
sig = 5.0 # 设置标准差阈值
for ii in range(x2d): # ii : 波长元素索引
# 创建掩膜；
# 1. DQ 数组（数据质量）
# 2. 误差值
# 3. CR 检测（宇宙射线）
mask_tmp = (dq_2d[:, ii] == 0) & (err_2d[:, ii] > 0) & ((data_2d[:, ii] - flux_y[:] * flux_disp1[ii])**2 < sig**2 * err_2d[:, ii]**2)
ivar = 1. / err_2d[:, ii]**2 # 计算逆方差
num = flux_y[:] * data_2d[:, ii] * ivar # 计算分子
den = flux_y[:]**2 * ivar # 计算分母
# 计算光谱通量和误差
flux_disp[ii] = num[mask_tmp].sum(axis=0) / den[mask_tmp].sum(axis=0) # 计算有效光谱通量
err_disp[ii] = np.sqrt(1./den[mask_tmp].sum(axis=0)) # 计算有效误差
wave_disp[ii] = wave_2d[0, ii] # 获取对应的波长
# 绘制误差条形图
plt.errorbar(wave_disp, flux_disp, yerr=err_disp) # 绘制光谱通量与误差
plt.xlim(1.7, 2.3) # 设置x轴范围
# 比较不同光谱数据的绘图
plt.errorbar(wave_disp, flux_disp, yerr=err_disp, color='r', label='Optimal') # 绘制最佳光谱数据，红色线，带误差条
plt.errorbar(wave_disp1, flux_disp1, yerr=err_disp1, color='b', alpha=0.5, label='Box') # 绘制盒子光谱数据，蓝色线，带误差条，透明度为0.5
plt.ylim(-10, 20000) # 设置y轴的范围
plt.legend(loc=0) # 显示图例，位置自动选择
plt.xlabel('Wavelength') # 设置x轴标签为“波长”
4. 将一维光谱写入文件;#
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_1d_opt.fits' # 定义输出文件名，包含目录和文件格式信息
# Now make it into a Spectrum1D instance. # 创建一个Spectrum1D实例
obs = Spectrum1D(spectral_axis=wave_disp*u.um, # 设置光谱轴，单位为微米
flux=flux_disp*u.MJy, # 设置光谱的通量，单位为MJy
uncertainty=StdDevUncertainty(err_disp), unit='MJy') # 设置不确定性，使用标准偏差不确定性
obs.write(file_1d, format='tabular-fits', overwrite=True) # 将Spectrum1D实例写入文件，格式为tabular-fits，允许覆盖
5. 沿x轴的光谱轮廓 = 分散光谱的分辨率；#
由于WFSS没有光阑，任何分散光谱都受到源形态的影响。我们将在接下来的笔记本中需要对有效光谱分辨率的估计。在这里，我们尝试提前进行估算；
for ii in range(sci_rot.shape[0]): # 遍历sci_rot数组的每一行
flux_tmp = sci_rot[ii, :] # 获取当前行的光谱数据
xx_tmp = np.arange(0, len(sci_rot[ii, :]), 1) # 创建一个与当前行长度相同的数组，用于x轴
plt.plot(xx_tmp, flux_tmp, label=f'y={ii}') # 绘制当前行的光谱数据，并添加标签
plt.legend(loc=1, fontsize=8) # 显示图例，位置在右上角，字体大小为8
plt.xlabel('Wavelength direction') # 设置x轴标签为“波长方向”
plt.title('Source light profile along dispersed direction', fontsize=14) # 设置图表标题，字体大小为14
*除非您对空间分辨光谱感兴趣，否则您可以叠加数据并获得光谱轮廓作为一个良好的近似。#
# 沿着交叉色散方向求和
flux_x = np.zeros(len(sci_rot[0, :]), 'float') # 创建一个与sci_rot的列数相同的零数组，用于存储每列的总和
for ii in range(sci_rot.shape[0]): # 遍历sci_rot的每一行
flux_x[ii] = np.sum(sci_rot[:, ii]) # 计算每一列的总和并存储到flux_x中
# 归一化;
flux_x[:] /= flux_x.sum() # 将flux_x中的每个元素除以总和，以进行归一化
plt.plot(xx_tmp, flux_x, label='Convolution kernel') # 绘制归一化后的flux_x与xx_tmp的关系图，并添加标签
plt.legend(loc=2) # 显示图例，位置设置为2
使用Moffat函数拟合；#
# Fitting function with Moffat
# Moffat 函数定义
def moffat(xx, A, x0, gamma, alp):
# 计算 Moffat 函数值
yy = A * (1. + (xx - x0)**2 / gamma**2)**(-alp)
return yy # 返回计算结果
def fit_mof(xx, lsf):
# 拟合 Moffat 函数
# xx = lsf * 0 # 初始化 xx 为零（注释掉的代码）
# for ii in range(len(lsf)): # 注释掉的代码，原本用于计算 xx
# xx[ii] = ii - len(lsf) / 2.
popt, pcov = curve_fit(moffat, xx, lsf) # 使用 curve_fit 拟合 Moffat 函数
return popt # 返回拟合参数
def LSF_mof(xsf, lsf, f_plot=True):
'''
输入参数:
=======
xsf : 轮廓的 x 轴数据.
lsf : 光强轮廓.
'''
# for ii in range(len(sci[0, :])): # 注释掉的代码，原本用于计算 lsf
# lsf[ii] = np.mean(sci_rot[int(height/2.-5):int(height/2.+5), ii])
# xsf[ii] = ii - len(lsf) / 2.
try:
A, xm, gamma, alpha = fit_mof(xsf, lsf) # 尝试拟合 Moffat 函数
except RuntimeError: # 捕获拟合失败的异常
print('Fitting failed.') # 输出拟合失败信息
A, xm, gamma, alpha = -1, -1, -1, -1 # 设置参数为 -1 表示失败
pass # 继续执行
if A > 0: # 如果拟合参数 A 大于 0
lsf_mod = moffat(xsf, A, 0, gamma, alpha) # 计算模型光强轮廓
if f_plot: # 如果需要绘图
yy = moffat(xsf, A, xm, gamma, alpha) # 计算拟合的 Moffat 函数值
plt.plot(xsf, yy, 'r.', ls='-', label='Data') # 绘制数据点
plt.plot(xsf, lsf_mod, 'b+', ls='-', label=f'Model: gamma={gamma:2f}\nalpha={alpha:2f}') # 绘制拟合模型
plt.legend() # 添加图例
plt.show() # 显示图形
return A, xm, gamma, alpha # 返回拟合参数
# LSF, 线扩散函数
iix_peak = np.argmax(flux_x) # 找到flux_x中最大值的索引
xx_tmp_shift = xx_tmp - xx_tmp[iix_peak] # 将xx_tmp中的值减去峰值位置的值，以便中心化
A, xm, gamma, alpha = LSF_mof(xx_tmp_shift, flux_x) # 调用LSF_mof函数计算参数A, xm, gamma和alpha
# 写入数据到文件
# 参数单位为像素
fm = open(f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_moffat.txt', 'w') # 打开文件以写入模式
fm.write('# A x0 gamma alp\n') # 写入文件头信息
fm.write('# Moffat function\n') # 写入函数类型信息
fm.write(f'{A:.3f} {xm:.3f} {gamma:.3f} {alpha:.3f}\n') # 写入Moffat函数的参数值，格式化为小数点后三位
fm.close() # 关闭文件
对其他滤光片和其他天体重复上述操作。#
以下大列执行相同的过程，适用于其他滤光片和抖动位置。#
grism = 'G150C' # 设置光栅类型为G150C
id = '00004' # 设置ID
DIR_OUT = './output/' # 设置输出目录
if not os.path.exists(DIR_OUT): # 如果输出目录不存在
os.mkdir(DIR_OUT) # 创建输出目录
filts = ['f115w', 'f150w', 'f200w'] # 定义滤光片列表
ndithers = np.arange(0, 2, 1) # 定义偏移量范围
sig = 5.0 # 设置信号阈值
for filt in filts: # 遍历每个滤光片
print(filt) # 打印当前滤光片名称
# 读取2D光谱数据
file_2d = f'{DIR_DATA}l3_nis_{filt}_{grism}_s{id}_cal.fits' # 构建2D光谱文件名
hdu_2d = fits.open(file_2d) # 打开2D光谱文件
for ndither in ndithers: # 遍历每个偏移量
print(ndither) # 打印当前偏移量
if grism == 'G150C': # 如果光栅类型为G150C
# 如果光谱是水平的
data_2d = hdu_2d[ndither*7+1].data # 获取2D数据
dq_2d = hdu_2d[ndither*7+2].data # 获取数据质量数组
err_2d = hdu_2d[ndither*7+3].data # 获取误差数组
wave_2d = hdu_2d[ndither*7+4].data # 获取波长数组
else: # 如果光栅类型不是G150C
data_2d = rotate(hdu_2d[ndither*7+1].data, 90) # 旋转数据
dq_2d = rotate(hdu_2d[ndither*7+2].data, 90) # 旋转数据质量数组
err_2d = rotate(hdu_2d[ndither*7+3].data, 90) # 旋转误差数组
wave_2d = rotate(hdu_2d[ndither*7+4].data, 90) # 旋转波长数组
y2d, x2d = data_2d.shape[:] # 获取数据的形状
plt.close() # 关闭当前图形
plt.imshow(data_2d, vmin=0, vmax=300) # 显示2D数据图像
plt.show() # 展示图像
# 重新提取2D图像
# if ndither == 0: # 如果是第一个偏移量
rsq = y2d # 设置方形区域的边长
# 提取科学数据和分段数据
sci_cut = data[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)]
seg_cut = segdata[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)]
# 不确定提取框的偏移是否是bug
if grism == 'G150C': # 如果光栅类型为G150C
sci_rot = rotate(sci_cut, PA_V3+0) # 旋转科学数据
else: # 如果光栅类型不是G150C
sci_rot = rotate(sci_cut, PA_V3+0+90) # 旋转科学数据
# 这是为了光谱分辨率
# 沿x轴获取光谱轮廓
for ii in range(sci_rot.shape[0]): # 遍历每一行
flux_tmp = sci_rot[ii, :] # 获取当前行的光谱数据
xx_tmp = np.arange(0, len(sci_rot[ii, :]), 1) # 创建x轴数组
# 沿交叉色散方向求和
flux_x = np.zeros(len(sci_rot[0, :]), 'float') # 初始化光谱数据数组
for ii in range(sci_rot.shape[0]): # 遍历每一行
flux_x[ii] = np.sum(sci_rot[ii, :]) # 计算每一行的总光谱
# 归一化
flux_x[:] /= flux_x.sum() # 归一化光谱数据
# LSF
iix_peak = np.argmax(flux_x) # 找到光谱峰值的索引
xx_tmp_shift = xx_tmp - xx_tmp[iix_peak] # 计算x轴的偏移量
A, xm, gamma, alpha = LSF_mof(xx_tmp_shift, flux_x) # 拟合Moffat函数
if ndither == 0: # 如果是第一个偏移量
# 写入文件
fm = open(f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_moffat.txt', 'w') # 打开输出文件
fm.write('# A x0 gamma alp\n') # 写入标题
fm.write(f'{A:.3f} {xm:.3f} {gamma:.3f} {alpha:.3f}\n') # 写入拟合参数
fm.close() # 关闭文件
# 这是为了最佳提取
# 沿x（色散）方向求和
flux_y = np.zeros(len(sci_rot[:, 0]), 'float') # 初始化光谱数据数组
for ii in range(sci_rot.shape[0]): # 遍历每一行
flux_y[ii] = np.sum(sci_rot[ii, :]) # 计算每一行的总光谱
# 归一化
flux_y[:] /= flux_y.sum() # 归一化光谱数据
# 根据Horne的方法
flux_disp = np.zeros(x2d, 'float') # 初始化光谱数据数组
err_disp = np.zeros(x2d, 'float') # 初始化误差数组
wave_disp = np.zeros(x2d, 'float') # 初始化波长数组
for ii in range(x2d): # 遍历每一列
# 掩膜
# 1. 数据质量数组
# 2. 误差值
# 3. CR检测
mask_tmp = (dq_2d[:, ii] == 0) & (err_2d[:, ii] > 0) # 创建掩膜
ivar = 1. / err_2d[:, ii]**2 # 计算逆方差
num = flux_y[:] * data_2d[:, ii] * ivar # 计算分子
den = flux_y[:]**2 * ivar # 计算分母
flux_disp[ii] = num[mask_tmp].sum(axis=0)/den[mask_tmp].sum(axis=0) # 计算光谱数据
err_disp[ii] = np.sqrt(1./den[mask_tmp].sum(axis=0)) # 计算误差
wave_disp[ii] = wave_2d[0, ii] # 获取波长
plt.close() # 关闭当前图形
con_plot = (wave_disp > 0) # 创建有效波长的掩膜
plt.errorbar(wave_disp[con_plot], flux_disp[con_plot], yerr=err_disp[con_plot]) # 绘制误差条
plt.ylim(-0, 3000) # 设置y轴范围
plt.show() # 展示图像
# 写入文件
# 现在将其转换为Spectrum1D实例
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{ndither}_1d_opt.fits' # 构建1D光谱文件名
if wave_disp[1] - wave_disp[0] < 0: # 如果波长数组是递减的
obs = Spectrum1D(spectral_axis=wave_disp[::-1]*u.um, # 创建光谱对象
flux=flux_disp[::-1]*u.MJy,
uncertainty=StdDevUncertainty(err_disp[::-1]), unit='MJy')
else: # 如果波长数组是递增的
obs = Spectrum1D(spectral_axis=wave_disp*u.um, # 创建光谱对象
flux=flux_disp*u.MJy,
uncertainty=StdDevUncertainty(err_disp), unit='MJy')
obs.write(file_1d, format='tabular-fits', overwrite=True) # 写入1D光谱文件
另一个天体;#
吸收线星系
grism = 'G150C' # 设置光栅类型为G150C
id = '00003' # 设置ID
DIR_OUT = './output/' # 设置输出目录
filts = ['f115w', 'f150w', 'f200w'] # 定义滤波器列表
ndithers = np.arange(0, 2, 1) # 定义dither的范围，这里有两个dither
sig = 5.0 # 设置信号阈值
for filt in filts: # 遍历每个滤波器
print(filt) # 打印当前滤波器
# 读取2D光谱数据
file_2d = f'{DIR_DATA}l3_nis_{filt}_{grism}_s{id}_cal.fits' # 构建2D数据文件名
hdu_2d = fits.open(file_2d) # 打开2D数据文件
for ndither in ndithers: # 遍历每个dither
print(ndither) # 打印当前dither
if grism == 'G150C': # 如果光栅是G150C
# 如果光谱是水平的
data_2d = hdu_2d[ndither*7+1].data # 获取2D数据
dq_2d = hdu_2d[ndither*7+2].data # 获取数据质量数组
err_2d = hdu_2d[ndither*7+3].data # 获取误差数组
wave_2d = hdu_2d[ndither*7+4].data # 获取波长数组
else: # 如果光栅不是G150C
data_2d = rotate(hdu_2d[ndither*7+1].data, 90) # 旋转2D数据
dq_2d = rotate(hdu_2d[ndither*7+2].data, 90) # 旋转数据质量数组
err_2d = rotate(hdu_2d[ndither*7+3].data, 90) # 旋转误差数组
wave_2d = rotate(hdu_2d[ndither*7+4].data, 90) # 旋转波长数组
y2d, x2d = data_2d.shape[:] # 获取2D数据的形状
plt.close() # 关闭当前图像
plt.imshow(data_2d, vmin=0, vmax=20) # 显示2D数据图像
plt.show() # 展示图像
# 重新提取2D图像
# if ndither == 0:
rsq = y2d # 设置提取区域的大小
sci_cut = data[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)] # 提取科学数据
seg_cut = segdata[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)] # 提取分段数据
# 不确定提取框的偏移是否是bug
if grism == 'G150C': # 如果光栅是G150C
sci_rot = rotate(sci_cut, PA_V3+0) # 旋转科学数据
else: # 如果光栅不是G150C
sci_rot = rotate(sci_cut, PA_V3+0+90) # 旋转科学数据
# 这是为了光谱分辨率
# 沿x轴获取光谱轮廓
for ii in range(sci_rot.shape[0]): # 遍历每一行
flux_tmp = sci_rot[ii, :] # 获取当前行的光谱数据
xx_tmp = np.arange(0, len(sci_rot[ii, :]), 1) # 创建x轴坐标
# 沿交叉色散方向求和
flux_x = np.zeros(len(sci_rot[0, :]), 'float') # 初始化光谱数据数组
for ii in range(sci_rot.shape[0]): # 遍历每一行
flux_x[ii] = np.sum(sci_rot[ii, :]) # 计算每一行的总光谱
# 归一化
flux_x[:] /= flux_x.sum() # 归一化光谱数据
# LSF
iix_peak = np.argmax(flux_x) # 找到光谱数据的峰值索引
xx_tmp_shift = xx_tmp - xx_tmp[iix_peak] # 计算x轴的偏移量
A, xm, gamma, alpha = LSF_mof(xx_tmp_shift, flux_x) # 计算Moffat函数参数
if ndither == 0: # 如果是第一个dither
# 写入文件
fm = open(f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_moffat.txt', 'w') # 打开输出文件
fm.write('# A x0 gamma alp\n') # 写入参数说明
fm.write(f'{A:.3f} {xm:.3f} {gamma:.3f} {alpha:.3f}\n') # 写入Moffat函数参数
fm.close() # 关闭文件
# 这是为了最优提取
# 沿x（色散）方向求和
flux_y = np.zeros(len(sci_rot[:, 0]), 'float') # 初始化光谱数据数组
for ii in range(sci_rot.shape[0]): # 遍历每一行
flux_y[ii] = np.sum(sci_rot[ii, :]) # 计算每一行的总光谱
# 归一化
flux_y[:] /= flux_y.sum() # 归一化光谱数据
# 按照Horne的方法
flux_disp = np.zeros(x2d, 'float') # 初始化光谱数据数组
err_disp = np.zeros(x2d, 'float') # 初始化误差数组
wave_disp = np.zeros(x2d, 'float') # 初始化波长数组
for ii in range(x2d): # 遍历每个像素
# 掩膜;
# 1. 数据质量数组
# 2. 误差值
# 3. CR检测
mask_tmp = (dq_2d[:, ii] == 0) & (err_2d[:, ii] > 0) # 创建掩膜
ivar = 1. / err_2d[:, ii]**2 # 计算逆方差
num = flux_y[:] * data_2d[:, ii] * ivar # 计算分子
den = flux_y[:]**2 * ivar # 计算分母
flux_disp[ii] = num[mask_tmp].sum(axis=0)/den[mask_tmp].sum(axis=0) # 计算光谱数据
err_disp[ii] = np.sqrt(1./den[mask_tmp].sum(axis=0)) # 计算误差
wave_disp[ii] = wave_2d[0, ii] # 获取波长
plt.close() # 关闭当前图像
con_plot = (wave_disp > 0) # 创建有效波长的掩膜
plt.errorbar(wave_disp[con_plot], flux_disp[con_plot], yerr=err_disp[con_plot]) # 绘制误差条图
plt.ylim(-20, 100) # 设置y轴范围
plt.show() # 展示图像
# 写入:
# 现在，将其转换为Spectrum1D实例。
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{ndither}_1d_opt.fits' # 构建1D数据文件名
if wave_disp[1] - wave_disp[0] < 0: # 如果波长是递减的
obs = Spectrum1D(spectral_axis=wave_disp[::-1]*u.um, # 创建光谱对象，波长反向
flux=flux_disp[::-1]*u.MJy, # 创建光谱对象，光谱反向
uncertainty=StdDevUncertainty(err_disp[::-1]), unit='MJy') # 创建光谱对象，误差反向
else: # 如果波长是递增的
obs = Spectrum1D(spectral_axis=wave_disp*u.um, # 创建光谱对象，波长正常
flux=flux_disp*u.MJy, # 创建光谱对象，光谱正常
uncertainty=StdDevUncertainty(err_disp), unit='MJy') # 创建光谱对象，误差正常
obs.write(file_1d, format='tabular-fits', overwrite=True) # 写入1D光谱数据文件