MOS 最优光谱提取

MOS 最优光谱提取#

用例： 最优光谱提取；方法参考 Horne (1986)。

数据： JWST 模拟的 NIRSpec MOS 数据；点源。

工具： jwst, stpsf, matplotlib, scipy, 自定义函数。

跨仪器： 任何光谱仪。

文档： 本笔记本是 STScI 更大后处理数据分析工具生态系统的一部分。

引言#

JWST 流水线为每种光谱模式生成来自组合曝光的 1-D 和 2-D 矫正光谱。目前，1D 产品是通过光圈提取生成的，计划通过 PSF 加权或拟合来实现最优提取。然而，在许多情况下，输出结果不一定是“最优”的，因此需要微调参数以改善结果。本笔记本旨在提供最优提取过程的逐步指导，并使用示例 JWST 数据进行说明。

定义术语#

最优提取： 一种首次由 Horne (1986) 定义的光圈提取方法。

S/N： 信噪比，衡量光谱噪声程度的指标。

WCS： 世界坐标系统，用于在不同参考框架之间转换。

导入库#

我们将使用以下库来执行最优光谱提取。

glob glob 用于收集文件名
numpy 用于处理数组函数以及其他各种活动
jwst.datamodels ImageModel, MultiSpecModel 用于访问示例数据的数模
astropy.io fits 用于低级 FITS 文件输入输出
astropy.modeling models, fitting 用于许多拟合任务
astropy.visualization astropy_mpl_style, simple_norm 用于显示美观的图像
scipy.interpolate interp1d, RegularGridInterpolator 用于所有插值需求
matplotlib.pyplot 用于绘制数据
matplotlib.patches Rectangle 用于在数据上绘制矩形
stpsf NIRSpec 用于生成和可视化来自仪器模型的 PSF（见附录 B）

%matplotlib inline  # 在Jupyter Notebook中内联显示matplotlib图形

import os  # 导入os模块，用于处理文件和目录

import tarfile  # 导入tarfile模块，用于处理tar文件

import zipfile  # 导入zipfile模块，用于处理zip文件

import urllib.request  # 导入urllib.request模块，用于处理URL请求

from urllib.parse import urlparse  # 从urllib.parse导入urlparse，用于解析URL

import requests  # 导入requests库，用于发送HTTP请求

from glob import glob  # 从glob导入glob，用于文件路径匹配

import numpy as np  # 导入numpy库，用于数值计算

from stdatamodels.jwst.datamodels import ImageModel  # 从stdatamodels导入JWST图像数据模型

from astropy.io import fits  # 从astropy.io导入fits模块，用于处理FITS文件

from astropy.modeling import models, fitting  # 从astropy.modeling导入模型和拟合工具

from astropy.visualization import astropy_mpl_style, simple_norm  # 导入可视化样式和简单归一化工具

from scipy.interpolate import interp1d, RegularGridInterpolator  # 从scipy.interpolate导入插值工具

import matplotlib.pyplot as plt  # 导入matplotlib.pyplot，用于绘图

from matplotlib.patches import Rectangle  # 从matplotlib.patches导入Rectangle，用于绘制矩形

plt.style.use(astropy_mpl_style)  # 使用导入的样式设置matplotlib的显示风格

import os  # 导入os模块，用于操作系统相关功能
import requests  # 导入requests模块，用于发送HTTP请求
import tarfile  # 导入tarfile模块，用于处理tar文件
from urllib.parse import urlparse  # 导入urlparse，用于解析URL

# 设置环境变量
os.environ["STPSF_PATH"] = "./data/stpsf-data"  # 设置STPSF数据的路径

# STPSF数据的URL和文件路径
stpsf_url = 'https://stsci.box.com/shared/static/kqfolg2bfzqc4mjkgmujo06d3iaymahv.gz'  # STPSF数据的下载链接
stpsf_file = './stpsf-data-LATEST.tar.gz'  # 下载后的文件名
stpsf_folder = "./data"  # 存放解压后文件的文件夹

# 定义下载文件的函数
def download_file(url, dest_path, timeout=60):
    parsed_url = urlparse(url)  # 解析URL

    # 检查URL的协议是否支持
    if parsed_url.scheme not in ["http", "https"]:
        raise ValueError(f"Unsupported URL scheme: {parsed_url.scheme}")  # 抛出不支持的协议错误

    # 发送GET请求下载文件
    response = requests.get(url, stream=True, timeout=timeout)  # 以流的方式下载文件
    response.raise_for_status()  # 检查请求是否成功

    # 将下载的内容写入目标文件
    with open(dest_path, "wb") as f:  # 以二进制写入模式打开文件
        for chunk in response.iter_content(chunk_size=8192):  # 分块读取内容
            f.write(chunk)  # 写入文件

# 检查STPSF文件夹是否存在
stpsfExist = os.path.exists(stpsf_folder)  # 检查文件夹是否存在

if not stpsfExist:  # 如果文件夹不存在
    os.makedirs(stpsf_folder)  # 创建文件夹
    download_file(stpsf_url, stpsf_file)  # 下载STPSF数据文件
    gzf = tarfile.open(stpsf_file)  # 打开下载的tar文件
    gzf.extractall(stpsf_folder)  # 解压文件到指定文件夹

加载数据#

我们将使用James Muzerolle提供的模拟级别3 MOS数据。这些文件来自一次包含多个点源的模拟访问，我们将从resample步骤的产品开始，这些文件的扩展名为s2d.fits。我们还将比较我们的最佳提取结果与extract1d步骤的产品，后者的扩展名为x1d.fits。有关这些文件的结构和格式的详细信息，请参见科学数据产品规范及其中的链接。

下面列出的最佳提取程序可以针对每个 s2d 文件中的每个 'SCI' 扩展重复进行。为了本笔记本的目的，我们将假设 resample 步骤已产生最佳输出，因此我们只需要访问这些扩展。（校正和组合输入光谱本身是一个复杂的过程，远超出本笔记本的范围！）

# 如果示例数据集已经下载，请注释掉这些行：

boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/optimal_extraction/optimal_extraction.zip'  # 数据集的下载链接

boxfile = './optimal_extraction.zip'  # 下载后保存的文件名

urllib.request.urlretrieve(boxlink, boxfile)  # 从链接下载文件并保存

zf = zipfile.ZipFile(boxfile, 'r')  # 打开下载的zip文件

list_names = zf.namelist()  # 获取zip文件内所有文件的名称列表

for item in list_names:  # 遍历文件名称列表

    zf.extract(member=item)  # 提取zip文件中的每个文件，使用extract而不是extractall更安全，避免绝对路径或相对路径问题

example_file = 'F170LP-G235M_MOS_observation-6_mod_correctedWCS_noflat_nooutlierdet_combined_s30263_'  # 示例文件的基本名称

s2d_file = os.path.join('s2d_files', example_file+'s2d.fits')  # 构建s2d文件的完整路径

x1d_file = os.path.join('x1d_files', example_file+'x1d.fits')  # 构建x1d文件的完整路径

data_model = ImageModel(s2d_file)  # 创建一个ImageModel对象，加载s2d_file数据

resampled_2d_image = data_model.data  # 获取重采样的2D图像数据；如果有多个SCI扩展，也需要指定EXTVER

weights_2d_image = data_model.wht  # 获取权重图像，用于后续估计每个像素的方差

image_shape = resampled_2d_image.shape  # 获取图像的形状

print(image_shape)  # 打印图像形状；注意x和y的顺序被交换

当我们想要查看二维光谱时，通常需要将像素在垂直方向上拉伸，以获得有用的图像。我们可以通过明确设置绘图的纵横比来实现这一点（我们会尽量保持矩形的比例）。

norm = simple_norm(resampled_2d_image, stretch='power')  # 使用简单归一化方法对重采样的2D图像进行归一化，采用'power'拉伸方式

aspect_ratio = image_shape[1] / (2 * image_shape[0])  # 计算图像的宽高比，宽度除以高度的两倍

fig1 = plt.figure()  # 创建一个新的图形对象，保存到dummy变量中以避免Jupyter Notebook的多余输出

img1 = plt.imshow(resampled_2d_image, cmap='gray', aspect=aspect_ratio,  # 显示重采样的2D图像，使用灰度色图和计算出的宽高比
                  norm=norm, interpolation='none')  # 应用归一化和无插值显示

clb1 = plt.colorbar()  # 添加颜色条以显示图像的颜色映射

抱歉，您似乎没有提供任何内容需要翻译。请提供需要翻译的Markdown内容，我将很乐意为您翻译。

最优提取算法#

以下是我们将要遵循的步骤大纲：

在2D图像上定义提取区域
识别高信噪比（S/N）交叉色散（已分箱和叠加）切片，以用于初始核拟合
定义提取核
1. 单一或复合点扩散函数（PSF）
2. 背景的多项式拟合
将提取核拟合到初始切片
跳过: 拟合几何失真
1. 确定用于轨迹拟合的交叉色散分箱
2. 对每个分箱进行核的第一次拟合，以找到轨迹中心
3. 轨迹中心的多项式拟合
将复合模型（核 | 轨迹）与2D图像结合，创建输出的1D光谱
将输出光谱与目录光度进行比较，以进行通量校准（尚不确定如何进行）

附录：

附录A：批处理
附录B：STPSF

定义提取区域#

我们首先确定在2D重采样图像中包含我们想要提取的光谱轨迹的区域。对于只有单一源的简单情况，我们理论上可以使用整个图像。然而，我们可能仍然希望排除背景中较大的系统性波动，这可能会使拟合变得复杂，或者排除几乎没有信号的轨迹部分，这将使得拟合轨迹中心变得困难。此外，在处理背景已进行节点减法的数据时，图像中将包含负轨迹，我们也希望将其排除。

fig2 = plt.figure(figsize=(9, 9)) # 创建一个9x9英寸的图形，尽可能大以适应笔记本

img2 = plt.imshow(resampled_2d_image, cmap='gray', aspect=aspect_ratio, 
                  norm=norm, interpolation='none') # 显示重采样的2D图像，使用灰度色图和之前定义的归一化

# 创建区域框和滑块

region_x = region_y = 0 # 初始化区域框的左上角坐标为(0, 0)

region_h, region_w = image_shape # 获取图像的高度和宽度

region_rectangle = Rectangle((region_x, region_y), region_w, region_h, 
                             facecolor='none', edgecolor='b', linestyle='--') # 创建一个矩形区域框，边框为蓝色虚线

current_axis = plt.gca() # 获取当前坐标轴

rect_patch = current_axis.add_patch(region_rectangle) # 将矩形区域框添加到当前坐标轴

我们将区域坐标设置为 x1=51, y1=3, x2=1268, y2=9，然后创建一个新数组，仅包含我们的提取区域（这样我们就不需要不断地索引原始数组）。

x1, x2 = 51, 1268  # 定义提取区域的x坐标范围

y1, y2 = 3, 9  # 定义提取区域的y坐标范围

# 创建网格以提取指定区域的y和x坐标
er_y, er_x = np.mgrid[y1:y2+1, x1:x2+1]

# 从重采样的2D图像中提取指定区域的像素值
extraction_region = resampled_2d_image[er_y, er_x]

# 从权重2D图像中提取指定区域的权重值
weights_region = weights_2d_image[er_y, er_x]

# 获取提取区域的形状
er_ny, er_nx = extraction_region.shape

# 计算长宽比
aspect_ratio = er_nx / (3. * er_ny)

# 对提取区域进行简单归一化处理
er_norm = simple_norm(extraction_region, stretch='power')

# 创建一个新的图形
fig3 = plt.figure()

# 显示提取区域的图像，设置颜色映射、长宽比、归一化和插值方式
img3 = plt.imshow(extraction_region, cmap='gray', aspect=aspect_ratio, 
                  norm=er_norm, interpolation='none')

# 添加颜色条以显示颜色映射
clb3 = plt.colorbar()

创建核切片#

我们现在定义一个交叉色散切片（cross-dispersion slice）作为我们的提取区域，以便拟合初始提取核（extraction kernel）。作为初始猜测，我们将对位于轨迹中间的30列进行叠加（coadd）。

slice_width = 30  # 定义切片的宽度

initial_column = er_nx // 2  # 计算初始列索引，取图像宽度的一半

def kernel_slice_coadd(width, column_idx):
    """
    对提取区域的多个列（= width）进行叠加，
    以column_idx为中心。
    """
    
    half_width = width // 2  # 计算切片的半宽度

    # 生成要叠加的列索引范围，确保不超出边界
    to_coadd = np.arange(max(0, column_idx - half_width), 
                         min(er_nx-1, column_idx + half_width))

    # 返回叠加后的结果，按列求和并归一化
    return extraction_region[:, to_coadd].sum(axis=1) / width

slice_0 = kernel_slice_coadd(slice_width, initial_column)  # 调用函数进行切片叠加

接下来，我们将绘制结果切片，并查看是否需要调整合成区域的宽度和中心。

# 创建一个包含两个子图的图形，纵向排列
fig4, (iax4, pax4) = plt.subplots(nrows=2, ncols=1, figsize=(8, 12))

# 调整子图之间的间距和图形的顶部和底部边距
plt.subplots_adjust(hspace=0.15, top=0.95, bottom=0.05)

# 在第一个子图中显示提取区域的图像，使用灰度色图
img4 = iax4.imshow(extraction_region, cmap='gray', aspect=aspect_ratio, 
                   norm=er_norm, interpolation='none')

# 创建切片框的函数
def make_slice(width, column_idx):
    sy, sh, sw = 0, er_ny, width  # 定义切片的起始y坐标、高度和宽度
    sx = column_idx - width // 2  # 计算切片的起始x坐标
    return sx, sy, sw, sh  # 返回切片的坐标和尺寸

# 使用切片宽度和初始列索引生成切片的坐标和尺寸
*sxy, sw, sh = make_slice(slice_width, initial_column)

# 创建一个矩形表示切片区域，设置边框颜色和样式
slice_rectangle = Rectangle(sxy, sw, sh, facecolor='none', 
                            edgecolor='b', linestyle='--')

# 将切片矩形添加到第一个子图中
iax4.add_patch(slice_rectangle)

# 绘制合并切片的信号
xd_pixels = np.arange(er_ny)  # 创建x轴像素数组

# 在第二个子图中绘制合并切片的信号曲线
lin4, = pax4.plot(xd_pixels, slice_0, 'k-')

# 设置第二个子图的x轴和y轴标签
xlbl4 = pax4.set_xlabel('Cross-dispersion pixel')  # x轴标签
ylbl4 = pax4.set_ylabel('Coadded signal')  # y轴标签

列索引为 670，宽度为 50 对于这个文件似乎效果不错，因此我们现在可以生成用于核拟合的最终切片。

# 调用函数 kernel_slice_coadd，传入参数 50 和 670，生成一个核切片
kernel_slice = kernel_slice_coadd(50, 670)

定义提取核#

现在我们将定义一个提取核（extraction kernel），该核将用于在色散方向上拟合每个像素的轨迹。这个核由两个部分组成：

一个点扩散函数（PSF）模板（或多个PSF的组合，用于去混叠）
一个用于背景拟合的多项式

选择一个PSF模板#

我们可以考虑多种PSF（点扩散函数）模板作为我们的内核，但全面比较超出了本笔记本的范围。我们将仅演示高斯（Gaussian）和莫法特（Moffat）轮廓。

需要注意两点：

这里展示的方法仅适用于真实的点源（true point source）。扩展源（extended sources）需要不同的方法论。
STPSF包可以直接从仪器模型构建复合PSF；然而，这一过程比使用astropy.modeling工具拟合1D轮廓要复杂得多，因此被放置在附录B中。

我们首先绘制两个轮廓与核切片的关系，采用简单的归一化方法，以便暂时忽略缩放，中心位于核的最大值所在的像素。（稍后我们会进行真正的拟合，别担心！）

max_pixel = np.argmax(kernel_slice)  # 找到kernel_slice中最大值的索引

fwhm = 1.  # 设置全宽半高（FWHM）参数

# 创建Moffat型一维模型，设置幅度、gamma、中心位置和alpha参数
moffat_profile = models.Moffat1D(amplitude=1, gamma=fwhm, x_0=max_pixel, alpha=1)

# 创建高斯型一维模型，设置幅度、均值和标准差
gauss_profile = models.Gaussian1D(amplitude=1, mean=max_pixel, stddev=fwhm)

fig5 = plt.figure()  # 创建一个新的图形窗口

# 绘制归一化的kernel_slice
kern5 = plt.plot(xd_pixels, kernel_slice / kernel_slice[max_pixel], label='Kernel Slice')

# 绘制Moffat型曲线
moff5 = plt.plot(xd_pixels, moffat_profile(xd_pixels), label='Moffat Profile')

# 绘制高斯型曲线
gaus5 = plt.plot(xd_pixels, gauss_profile(xd_pixels), label='Gaussian Profile')

lgd5 = plt.legend()  # 添加图例

高斯轮廓看起来是一个更好的近似，因此我们将使用这个轮廓来处理该光谱。在下面的单元格中，我们可以使用模型操作添加更多的点扩散函数（PSF）模板；这留给读者作为练习。

我们需要对幅度进行去归一化，因此我们将其设置为切片的最大像素值。

psf_template = gauss_profile  # 将高斯轮廓赋值给psf_template

psf_template.amplitude = kernel_slice[max_pixel]  # 设置psf_template的幅度为kernel_slice中最大像素的值

print(psf_template)  # 打印psf_template以查看其内容

# 如果需要对多个源进行解混合，可以在这里添加更多的PSF模板：

多项式背景#

我们将使用多项式来拟合背景。建议进行一些实验，以找到最适合数据的多项式阶数；在这个例子中，我们将使用二次多项式（2nd-degree polynomial）。

对于经过点位减法（nod-subtracted）处理的数据，提取区域内可能没有足够的像素来准确拟合残差。在这种情况下，可以使用零阶多项式（0th-order polynomial）或 Const1D 模型来表示背景；如果希望完全避免拟合背景，可以将参数设置为 fixed = True。

# 创建一个二次多项式模型
background_poly = models.Polynomial1D(2)

# 打印多项式模型的详细信息
print(background_poly)

最后一步是将点扩散函数（PSF）和背景结合起来，以创建我们的复合模型。

# 定义提取核，包含点扩散函数模板和背景多项式
extraction_kernel = psf_template + background_poly

# 打印提取核的内容
print(extraction_kernel)

拟合提取核#

现在我们已经有了一个提取核（extraction kernel），我们希望将其拟合到我们的核切片（kernel slice），以便在下一步中拥有最佳的工具来拟合轨迹中心（trace centers）。我们还绘制了拟合组件（fit components）以及拟合与核切片的对比图（fit vs the kernel slice），作为视觉检查；如果结果不可接受，我们可以返回到前一部分，调整参数，然后再试一次。

fitter = fitting.LevMarLSQFitter()  # 创建一个Levenberg-Marquardt最小二乘拟合器

fit_extraction_kernel = fitter(extraction_kernel, xd_pixels, kernel_slice)  # 拟合提取核

print(fit_extraction_kernel)  # 打印拟合结果

fit_line = fit_extraction_kernel(xd_pixels)  # 计算拟合线

fig6, (fax6, fln6) = plt.subplots(nrows=2, ncols=1, figsize=(8, 12))  # 创建一个包含两个子图的图形

plt.subplots_adjust(hspace=0.15, top=0.95, bottom=0.05)  # 调整子图之间的间距和图形的边距

psf6 = fax6.plot(xd_pixels, fit_extraction_kernel[0](xd_pixels), label="PSF")  # 绘制点扩散函数（PSF）

poly6 = fax6.plot(xd_pixels, fit_extraction_kernel[1](xd_pixels), label="Background")  # 绘制背景

sum6 = fax6.plot(xd_pixels, fit_line, label="Composite Kernel")  # 绘制复合核

lgd6a = fax6.legend()  # 添加图例到第一个子图

lin6 = fln6.plot(xd_pixels, kernel_slice, label='Kernel Slice')  # 绘制核切片

fit6 = fln6.plot(xd_pixels, fit_line, 'o', label='Extraction Kernel')  # 绘制提取核的拟合结果

lgd6b = fln6.legend()  # 添加图例到第二个子图

波长变化的全宽半最大值 (FWHM)（已跳过）#

NIRSpec 的点扩散函数（PSF）宽度会随波长变化，因此对于科学数据，沿光谱轨迹的多个位置进行拟合可能是有益的。下面是该过程的演示；然而，请注意，对于这个示例数据集，未经过优化的重采样和合并抖动输入光谱引入了宽度变化伪影，因此我们实际上不会使用此步骤的结果进行提取。

如果我们希望考虑变化的全宽半最大值（FWHM），我们可以在色散方向上对二维光谱进行分箱，并对每个箱进行拟合。我们上面定义的核可以作为我们的初始估计，这在光谱的非常微弱区域中可能会很有帮助，因为astropy.modeling的拟合例程对初始估计非常敏感。

（在计算并绘制了分箱的核FWHM后，下一步将是找到合适的模型，并将FWHM作为箱中心的函数进行拟合。拟合模型将被包含在下面的最终一维提取中。）

n_bin = 100  # 设置要分的bin的数量

bin_width = er_nx // n_bin  # 计算每个bin的宽度

# 计算每个bin的中心位置
bin_centers = np.arange(0, er_nx, bin_width+1, dtype=float) + bin_width // 2

# 将提取区域中的光谱数据按bin宽度进行求和，形成binned_spectrum
binned_spectrum = np.hstack([extraction_region[:, i:i+bin_width+1].sum(axis=1)[:, None] 
                            for i in range(0, er_nx, bin_width+1)])

# 初始化一个与bin_centers相同形状的数组，用于存储每个bin的FWHM
bin_fwhms = np.zeros_like(bin_centers, dtype=float)

# 遍历每个bin中心，进行拟合以获取FWHM
for y in range(bin_centers.size):
    bin_fit = fitter(fit_extraction_kernel, xd_pixels, binned_spectrum[:, y])  # 拟合当前bin的光谱
    bin_fwhms[y] = bin_fit.stddev_0.value  # 存储拟合得到的FWHM值

# 获取binned_spectrum的形状
bin_ny, bin_nx = binned_spectrum.shape

# 计算图像的宽高比
bin_ar = bin_nx / (3 * bin_ny)

# 创建一个包含两个子图的图形
fig_fwhm, ax_fwhm = plt.subplots(nrows=2, ncols=1, figsize=(6, 10))

plt.subplots_adjust(hspace=0.05)  # 调整子图之间的间距

# 在第一个子图中显示binned_spectrum的图像
fwhm_img = ax_fwhm[0].imshow(binned_spectrum, aspect=bin_ar, interpolation='none',
                             cmap='gray')

# 在第二个子图中绘制bin中心与FWHM的关系
fwhm_plot = ax_fwhm[1].plot(bin_centers, bin_fwhms)

# 设置第二个子图的x轴标签
xlbl_fwhm = ax_fwhm[1].set_xlabel("Bin center (px)")

# 设置第二个子图的y轴标签
ylbl_fwhm = ax_fwhm[1].set_ylabel("FWHM (arcsec)")

拟合几何畸变 (跳过)#

管道 resample 步骤将所有输入的二维光谱滴水（drizzle）到一个矩形网格上，因此我们最优提取过程中的这一特定步骤通常不是必需的。这里简要讨论该过程，以作为提取未矩形光谱（后缀为 _cal.fits）的指南，其中光谱的轨迹可能具有显著的曲率，并且轨迹的色散并不与列对齐。

定义轨迹拟合的区间#

根据2D重采样光谱的噪声水平，定义色散方向上的区间可能是有益的。这些区间可以是均匀或不均匀分布的，一旦定义完成，就可以对每个区间中的列进行合并（可能使用s2d文件中的WHT扩展），并创建一个区间中心位置的数组。

如果2D光谱的信噪比（S/N）很高，这一步可能就不必要了，每个交叉色散列可以在下一步中单独拟合。

使用修改后的提取核拟合每个箱体#

我们希望用我们的提取核拟合每个定义的箱体，但我们不希望其他伪影或噪声干扰轨迹。因此，我们复制提取核，然后将除轮廓中心（上面示例中的 mean_0）以外的每个参数设置为 fixed = True。从轨迹的一端开始，遍历每个箱体，使用提取核拟合切片，并将结果轨迹中心存储在一个数组中。

使用一维多项式拟合轨迹中心#

这一步很简单：创建一个 Polynomial1D 模型，然后将其拟合到上一步骤得到的轨迹中心。

由于我们不会进行拟合，而是创建一个占位符轨迹中心模型：一个零阶多项式。

trace_center_model = models.Polynomial1D(0)  # 我们使用常数，因为光谱已经被校正
trace_center_model.c0 = fit_extraction_kernel.mean_0.value  # 使用PSF轮廓中心的参数

print(trace_center_model)  # 打印模型信息

构建最终的1D光谱#

我们计算最终的1D光谱作为2D光谱在交叉色散方向上的加权和，使用我们基于轨迹中心的复合模型（提取核）作为权重。我们还需要纳入每个像素的方差，这将通过重采样步骤输出的WHT扩展来估计。

创建方差图像#

Horne 的算法需要每个像素的方差。目前，误差并未在重采样步骤中传播；然而，根据 DrizzlePac 手册，我们可以通过雨滴权重图像来估计方差：\( Var \approx 1 / (W \times s^4) \)，其中 \(s\) 是像素尺度。目前，NIRSpec 的雨滴参数设置为 PIXFRAC = 1.0。

scale = 1.0  # 调整此值以应对NIRSpec PIXFRAC的变化

# 我们希望在下一步中排除任何权重为0的像素
# 因此，我们将使用掩码数组操作。

bad_pixels = weights_region == 0  # 标记权重为0的像素

masked_wht = np.ma.array(weights_region, mask=bad_pixels)  # 创建掩码数组，掩盖坏像素

variance_image = np.ma.divide(1., weights_region * scale**4)  # 计算方差图像，使用掩码数组

我们可以显示方差图像，以查看提取区域中是否有任何区域不会包含在光谱中（如下图所示，红色区域表示）。对于这个特定的光谱示例，每个像素都有非零权重。

from copy import copy  # 从copy模块导入copy函数

fig_var = plt.figure()  # 创建一个新的图形窗口

palette = copy(plt.cm.gray)  # 复制灰度调色板

palette.set_bad('r', alpha=0.7)  # 设置无效值的颜色为红色，透明度为0.7

var_norm = simple_norm(variance_image, stretch='log', vmin=0.006, vmax=0.1)  # 对方差图像进行归一化处理，使用对数伸缩

img_var = plt.imshow(variance_image, interpolation='none', aspect=aspect_ratio, norm=var_norm, cmap=palette)  # 显示方差图像，设置插值方式、纵横比、归一化和调色板

生成一维光谱#

现在，我们最终计算我们的1D光谱，通过对交叉色散列进行求和：

\[S_x = \frac{1}{G_x}\sum_{y} \frac{I_{xy}\cdot K_y(x)}{V_{xy}}\]

其中 \(I\) 是二维重采样图像中的像素值，\(K\) 是设置为列的轨迹中心的提取核，\(V\) 是方差图像中的像素值，而 \(G\) 是由以下公式给出的核归一化：

\[G_x = \sum_y \frac{K_y^2(x)}{V_{xy}}\]

spectrum = np.zeros(er_nx, dtype=float) # 初始化我们的光谱，填充为零

column_pixels = np.arange(er_nx) # 创建一个包含列像素的数组

trace_centers = trace_center_model(column_pixels) # 计算我们的轨迹中心数组

# 遍历每一列

for x in column_pixels:

    # 为这一列创建内核，使用拟合的轨迹中心

    kernel_column = fit_extraction_kernel.copy() # 复制提取内核

    kernel_column.mean_0 = trace_centers[x] # 设置内核的均值为轨迹中心

    # kernel_column.stddev_0 = fwhm_fit(x) # 如果考虑变化的FWHM，取消注释这一行。

    kernel_values = kernel_column(xd_pixels) # 计算当前列的内核值

    # 根据Horne1986方程5进行单位归一化，P_x = P/sum(P_x)。 

    kernel_values = np.ma.masked_outside(kernel_values, 0, np.inf) # 遮蔽无效值

    kernel_values = kernel_values/np.ma.sum(kernel_values) # 归一化内核值

    # 隔离光谱和方差图像中相关的列

    variance_column = variance_image[:, x] # 记住，numpy数组是行，列的顺序

    image_pixels = extraction_region[:, x] # 获取提取区域的像素值

    # 计算内核归一化

    g_x = np.ma.sum(kernel_values**2 / variance_column) # 计算归一化因子

    if np.ma.is_masked(g_x): # 如果这一列无效，则跳过

        continue

    # 现在对加权列进行求和

    weighted_column = np.ma.divide(image_pixels * kernel_values, variance_column) # 计算加权列

    spectrum[x] = np.ma.sum(weighted_column) / g_x # 将加权列的和除以归一化因子，得到光谱值

我们需要一个波长数组来显示光谱，这可以通过从数据模型的元数据中存储的WCS对象创建。

wcs = data_model.meta.wcs  # 获取数据模型中的WCS（世界坐标系统）信息

print(wcs.__repr__())  # 打印WCS对象的字符串表示

alpha_C, delta_C, y = wcs(er_x, er_y)  # 使用WCS将像素坐标（er_x, er_y）转换为天球坐标（alpha_C, delta_C）和波长（y）

wavelength = y[0]  # 提取波长信息，y数组的第一个元素

fig7 = plt.figure()  # 创建一个新的图形窗口

spec7 = plt.plot(wavelength, spectrum)  # 绘制波长与光谱的关系图

# 将提取的光谱写入文件

# 这部分留给读者练习

# 假设我们有一个光谱数据的数组和文件名
spectrum_data = [...]  # 光谱数据数组
output_filename = 'spectrum.txt'  # 输出文件名

# 打开文件以写入数据
with open(output_filename, 'w') as file:  # 以写入模式打开文件
    for value in spectrum_data:  # 遍历光谱数据中的每个值
        file.write(f"{value}\n")  # 将每个值写入文件，并换行

在这个代码中，我们假设有一个光谱数据的数组`spectrum_data`和一个输出文件名`output_filename`。通过打开文件并遍历光谱数据，我们将每个值写入文件中。

我们还想将我们优化提取的光谱与 x1d 管道产品进行比较。我们将对光谱进行归一化，以便能够在同一坐标轴上绘制它们。

（请注意，x1d 光谱包含来自背景减除步骤的负值，这通常会导致负的通量计算。在与我们优化提取的版本进行比较时，我们需要对此进行修正。）

x1d_model = fits.open(x1d_file)  # 打开指定的x1d文件，返回一个FITS模型对象

# For a file with multiple spectra, the index to .spec is EXTVAR  # 对于包含多个光谱的文件，.spec的索引为EXTVAR

tmp = x1d_model[1].data  # 获取FITS模型中第二个扩展的数据

x1d_wave = tmp['WAVELENGTH']  # 提取波长数据

x1d_flux = tmp['FLUX']  # 提取光谱数据

if x1d_flux.sum() <= 0:  # 如果光谱数据的总和小于等于0
    x1d_flux = -x1d_flux  # 将光谱数据取反

fig8 = plt.figure()  # 创建一个新的图形窗口

x1d8 = plt.plot(x1d_wave, x1d_flux / x1d_flux.max(), label="Pipeline")  # 绘制归一化后的光谱数据，标签为"Pipeline"

opt8 = plt.plot(wavelength, spectrum / spectrum.max(), label="Optimal", alpha=0.7)  # 绘制优化光谱数据，标签为"Optimal"，透明度为0.7

lgd8 = plt.legend()  # 显示图例

看起来您没有提供任何Markdown内容。请提供需要翻译的文本，我将帮助您进行翻译并保持格式。

附录 A：批处理#

当需要对大量光谱进行最佳提取时，逐步按照上述过程处理每个光谱可能不够实际。在这种情况下，我们可以最初对一两个光谱使用这些交互式方法，以便决定一些提取参数（例如，使用什么点扩散函数（PSF）模板轮廓，或者用什么阶数的多项式来拟合背景），然后使用这些参数以非交互方式处理所有光谱。之后，我们可以检查每个提取光谱的输出，并重新审视任何需要更个性化处理的光谱。

我们可以通过为上述每个步骤定义函数，以及一个单一的主函数来迭代单个目录中的所有光谱，从而非交互式地提取大量光谱。

定义提取区域#

这个步骤无法以非交互方式执行，因此我们在这里跳过。不过，对于真实数据集，有两种好的方法（和一种不好的方法）来处理这个问题：

在批处理之前为每个二维光谱定义一个提取区域。您可以将区域边界框保存到一个 Python 字典中（或者将它们写入文件，然后在迭代时读取）。
直观地检查二维光谱，仅对那些不需要定义特定提取区域（即小于完整二维光谱）的光谱进行批处理。其余光谱可以单独提取。
跳过此步骤，并假设任何需要定义特定提取区域的光谱都需要单独重新处理。虽然不推荐这种方法，但我们将在这里使用它。

创建核切片#

def batch_kernel_slice(extraction_region, slice_width=30, column_idx=None):
    """
    在提取区域 `extraction_region` 中创建一个切片， 
    该切片在交叉色散方向上，以 `column_idx` 为中心， 
    宽度为 `slice_width` 像素。如果未给定 `column_idx`， 
    则使用总信号最大的列。
    """

    # 如果没有提供列索引，则找到信号总和最大的列索引
    if column_idx is None:
        column_idx = np.argmax(extraction_region.sum(axis=0))

    # 获取提取区域的形状
    ny, nx = extraction_region.shape
    half_width = slice_width // 2  # 计算切片的一半宽度

    # 确保切片不会超出提取区域的边界
    to_coadd = np.arange(max(0, column_idx - half_width), 
                         min(nx-1, column_idx + half_width))

    # 返回切片在交叉方向上的平均值
    return extraction_region[:, to_coadd].sum(axis=1) / slice_width

创建并拟合提取核#

def batch_fit_extraction_kernel(xd_slice, psf_profile=models.Gaussian1D, 

                                height_param_name='amplitude', height_param_value=None,

                                width_param_name='stddev', width_param_value=1.,

                                center_param_name='mean', center_param_value=None,

                                other_psf_args=[], other_psf_kw={},

                                bg_model=models.Polynomial1D,

                                bg_args=[3], bg_kw={}):
    """
    初始化一个复合提取核，然后将其拟合到
    一维数组 `xd_slice`，该数组通常是通过
    上面定义的 `kernel_slice` 函数生成的。

    为了允许具有不同参数名称的 PSF 模板模型，
    我们使用 `height_param_*`、`width_param_*` 和
    `center_param_*` 关键字参数。我们将任何其他
    位置或关键字参数收集到 PSF 模型中，
    存储在 `other_psf_*` 中。如果高度或中心值为 `None`，
    则将从数据中计算它们。

    同样，任何希望传递给背景拟合模型（默认为 `Polynomial1D`）的
    位置或关键字参数都可以通过 `bg_args` 和 `bg_kw` 接受。

    请注意，此函数无法处理涉及多个 PSF 进行解混合的情况。
    建议单独处理此类光谱，使用上述交互程序。
    """

    xd_pixels = np.arange(xd_slice.size)  # 创建一个像素索引数组

    if center_param_value is None:  # 如果中心参数值为 None
        center_param_value = np.argmax(xd_slice)  # 计算最大值的索引作为中心

    if height_param_value is None:  # 如果高度参数值为 None
        # 如果通过 center_param_value 传递了非整数值，
        # 我们需要进行插值。
        slice_interp = interp1d(xd_pixels, xd_slice)  # 创建插值函数
        height_param_value = slice_interp(center_param_value)  # 计算高度值

    # 创建 PSF 和背景模型
    psf_kw = dict([(height_param_name, height_param_value),  # 设置高度参数
                   (width_param_name, width_param_value),  # 设置宽度参数
                   (center_param_name, center_param_value)])  # 设置中心参数
    psf_kw.update(other_psf_kw)  # 更新其他 PSF 参数
    psf = psf_profile(*other_psf_args, **psf_kw)  # 创建 PSF 实例

    bg = bg_model(*bg_args, **bg_kw)  # 创建背景模型实例

    composite_kernel = psf + bg  # 组合 PSF 和背景模型

    fitter = fitting.LevMarLSQFitter()  # 创建拟合器实例
    return fitter(composite_kernel, xd_pixels, xd_slice)  # 返回拟合结果

考虑变化的全宽半最大值（FWHM）#

根据这里所示的过程，这部分留给用户自行练习。请注意，如果需要，下面的 batch_extract_spectrum 和 batch_optimal_extraction 也需要进行修改，以纳入此功能。

def batch_vary_fwhm(extraction_region, kernel):
    # 定义一个函数，拟合波长变化的全宽半最大值（FWHM）

    pass  # 实现该函数的具体逻辑

拟合轨迹中心#

如果需要，请将此替换为一个实际的执行拟合的函数。

def batch_fit_trace_centers(extraction_region, kernel,
                            trace_model=models.Polynomial1D,  # 使用多项式模型进行轨迹拟合
                            trace_args=[0], trace_kw={}):  # 模型参数和关键字参数

    """
    使用模型拟合轨迹的几何失真。
    目前这是一个占位符函数，因为几何失真通常在`resample`步骤中被移除。
    但是，如果需要此功能，请使用此函数签名以保持与本附录其余部分的兼容性。
    """

    trace_centers = trace_model(*trace_args, **trace_kw)  # 创建轨迹中心模型
    trace_centers.c0 = kernel.mean_0  # 设置轨迹中心的初始值为核的均值

    return trace_centers  # 返回拟合后的轨迹中心

生成一维光谱#

def batch_extract_spectrum(extraction_region, trace, kernel, 

                           weights_image, 

                           trace_center_param='mean',

                           scale=1.0):

    """

    从提取区域优化提取1D光谱。

    

    从`weights_image`（应与`extraction_region`具有相同的维度）创建方差图像。

    然后，对于光谱的每一列，我们根据上述定义的方程求和光圈，

    遮蔽权重为零的像素。 

    

    注意，与交互式逐步方法不同，这里我们将进行向量化以提高速度。

    这需要使用模型集来处理核，但这是允许的，因为

    我们不需要拟合任何东西。

    

    `trace_center_param`是定义轨迹中心的参数名称，

    *不带模型编号下标*（因为我们将单独处理组件）。

    

    `scale`是输入像素与输出像素的大小比率，

    在滴水时等同于`resample`步骤中的滴水参数中的PIXFRAC。

    """

    

    bad_pixels = weights_image == 0.  # 标记权重为零的坏像素
    masked_wht = np.ma.array(weights_image, mask=bad_pixels)  # 创建带掩码的权重数组
    variance_image = np.ma.divide(1., masked_wht * scale**4)  # 计算方差图像
    
    ny, nx = extraction_region.shape  # 获取提取区域的形状
    trace_pixels = np.arange(nx)  # 创建轨迹像素数组
    xd_pixels = np.arange(ny)  # 创建y方向像素数组
    trace_centers = trace(trace_pixels)  # 计算轨迹中心数组
    
    # 创建用于向量化的核图像，这需要一些技巧...
    # ******************************************************************
    # * 重要提示:                                                     *
    # * ----------                                                     *
    # * 注意，由于模型集的实现方式，修改现有模型实例以使用它们是不可行的。  *
    # * 相反，我们将创建一个新的核实例，使用原始核的拟合参数。                           *
    # *                                                                *
    # * 警告：这假设PSF是第一个元素，背景是第二个。如果您在创建复合核时更改了这一点，*
    # * 请确保相应地更新此部分，否则将无法正常工作！                                *
    # ******************************************************************
    
    psf0, bg0 = kernel  # 解包核中的PSF和背景
    psf_params = {}  # 初始化PSF参数字典

    for pname, pvalue in zip(psf0.param_names, psf0.parameters):  # 遍历PSF参数
        if pname == trace_center_param:  # 如果参数是轨迹中心参数
            psf_params[pname] = trace_centers  # 使用轨迹中心
        else:
            psf_params[pname] = np.full(nx, pvalue)  # 其他参数填充为常数

    psf_set = psf0.__class__(n_models=nx, **psf_params)  # 创建新的PSF模型集

    # 如果不使用Polynomial1D作为背景模型，请编辑此部分：
    bg_set = bg0.__class__(len(bg0.param_names)-1, n_models=nx)  # 创建背景模型集

    for pname, pvalue in zip(bg0.param_names, bg0.parameters):  # 遍历背景参数
        setattr(bg_set, pname, np.full(nx, pvalue))  # 将背景参数填充为常数

    kernel_set = psf_set + bg_set  # 合并PSF和背景模型集

    # 我们传递model_set_axis=False，以便模型集中每个模型使用相同的输入，
    # 并转置结果以修正方向。
    kernel_image = kernel_set(xd_pixels, model_set_axis=False).T  # 计算核图像
    
    # 现在我们使用numpy.ma例程执行加权求和，以保留我们的掩码
    # 单位归一化核，遵循Horne1986公式5，P_x = P/sum(P_x)。 
    kernel_image = np.ma.masked_outside(kernel_image, 0, np.inf)  # 遮蔽非正值
    kernel_image = kernel_image / np.ma.sum(kernel_image)  # 归一化核图像
    
    g = np.ma.sum(kernel_image**2 / variance_image, axis=0)  # 计算归一化因子
    weighted_spectrum = np.ma.divide(kernel_image * extraction_region, variance_image)  # 计算加权光谱
    spectrum1d = np.ma.sum(weighted_spectrum, axis=0) / g  # 计算1D光谱
    
    # 将任何被遮蔽的值设置为0。
    return spectrum1d.filled(0.)  # 返回填充后的光谱

便利函数#

def batch_wavelength_from_wcs(datamodel, pix_x, pix_y):
    """
    从数据模型的元数据中获取WCS对象， 
    根据给定的像素坐标生成世界坐标，并返回1D波长。
    """
    
    wcs = datamodel.meta.wcs  # 获取WCS对象
    aC, dC, y = wcs(pix_x, pix_y)  # 将像素坐标转换为世界坐标
    return y[0]  # 返回第一个波长值

def batch_save_extracted_spectrum(filename, wavelength, spectrum):
    """
    提取光谱的快速FITS输出。
    替换为您首选的输出格式和函数。
    """
    
    wcol = fits.Column(name='wavelength', format='E',  # 创建波长列
                       array=wavelength)  # 设置波长数组
    scol = fits.Column(name='spectrum', format='E',  # 创建光谱列
                       array=spectrum)  # 设置光谱数组
    cols = fits.ColDefs([wcol, scol])  # 定义列
    hdu = fits.BinTableHDU.from_columns(cols)  # 创建二进制表HDU
    hdu.writeto(filename, overwrite=True)  # 将HDU写入文件，覆盖同名文件

def batch_plot_output(resampled_image, extraction_bbox, 
                      kernel_slice, kernel_model,
                      wavelength, spectrum, filename):
    """
    方便的汇总输出图形， 
    允许对每个处理文件的结果进行可视检查。
    """
    
    fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, ncols=1, 
                                        figsize=(8, 12))  # 创建3行1列的子图
    fig.suptitle(filename)  # 设置图形标题
    
    ny, nx = resampled_image.shape  # 获取重采样图像的形状
    aspect = nx / (2 * ny)  # 计算纵横比
    
    # 子图1：提取区域
    power_norm = simple_norm(resampled_image, 'power')  # 归一化图像
    _ = ax1.imshow(resampled_image, interpolation='none',
                   aspect=aspect, norm=power_norm, cmap='gray')  # 显示重采样图像
    rx, ry, rw, rh = extraction_bbox  # 解包提取边界框
    region = Rectangle((rx, ry), rw, rh, facecolor='none', 
                       edgecolor='b', linestyle='--')  # 创建提取区域矩形
    _ = ax1.add_patch(region)  # 将矩形添加到子图1
    
    # 子图2：核拟合
    xd_pixels = np.arange(kernel_slice.size)  # 创建像素索引
    fit_line = kernel_model(xd_pixels)  # 计算拟合线
    _ = ax2.plot(xd_pixels, kernel_slice, label='Kernel Slice')  # 绘制核切片
    _ = ax2.plot(xd_pixels, fit_line, 'o', label='Extraction Kernel')  # 绘制提取核
    _ = ax2.legend()  # 添加图例
    
    # 子图3：提取的光谱
    _ = ax3.plot(wavelength, spectrum)  # 绘制光谱
    
    fig.savefig(filename, bbox_inches='tight')  # 保存图形
    plt.close(fig)  # 关闭图形

遍历所需文件#

def batch_optimal_extraction(file_list):
    """
    遍历一个fits文件路径的列表，优化提取每个文件中的SCI扩展，
    生成输出摘要图像，然后保存结果光谱。

    注意，在示例数据集中，每个文件中只有一个SCI扩展。
    对于具有多个SCI扩展的数据，需要对这些扩展进行第二次循环。
    """

    # 对于这个示例数据，我们将使用所有函数的默认值
    for i, fitsfile in enumerate(file_list):
        # 打印当前处理的文件信息
        print("Processing file {} of {}: {}".format(i+1, len(file_list), fitsfile))

        # 创建图像模型对象，读取fits文件
        dmodel = ImageModel(fitsfile)

        # 获取数据和权重
        spec2d = dmodel.data  # 获取二维光谱数据
        wht2d = dmodel.wht    # 获取权重数据

        # 批量获取光谱的核切片
        k_slice = batch_kernel_slice(spec2d)

        # 批量拟合提取核
        k_model = batch_fit_extraction_kernel(k_slice)

        # 批量拟合轨迹中心
        trace = batch_fit_trace_centers(spec2d, k_model)

        # 批量提取光谱
        spectrum = batch_extract_spectrum(spec2d, trace, k_model, wht2d)

        # 获取光谱数据的形状
        ny, nx = spec2d.shape

        # 创建二维网格
        y2d, x2d = np.mgrid[:ny, :nx]

        # 从WCS获取波长
        wavelength = batch_wavelength_from_wcs(dmodel, x2d, y2d)

        # 定义边界框
        bbox = [0, 0, nx-1, ny-1]

        # 生成输出文件名
        outfile = fitsfile.replace('s2d.fits', 'x1d_optimal')

        # 绘制输出图像并保存
        batch_plot_output(spec2d, bbox, k_slice, k_model,
                          wavelength, spectrum, 
                          outfile+'.png')

        # 保存提取的光谱到fits文件
        batch_save_extracted_spectrum(outfile+'.fits', wavelength, spectrum)

在示例数据集上运行#

请特别注意在拟合过程中产生警告的任何光谱 - 这些光谱很可能是交互式重处理的良好候选者。

plt.ioff()  # 如果不关闭这个选项，matplotlib会尝试为每个光谱显示一个（不可见的）图

s2d_files = glob(os.path.join('s2d_files', '*s2d.fits'))  # 获取所有s2d.fits文件的路径

batch_optimal_extraction(s2d_files)  # 对获取的s2d文件进行批量最优提取

plt.ion()  # 现在重新打开这个选项，以便后续的图形能够正常显示！

看起来您没有提供任何Markdown内容。如果您有需要翻译的文本，请将其粘贴在这里，我将很乐意为您翻译。

附录 B: STPSF#

我们可以直接通过仪器模型使用 STPSF 生成点扩散函数（PSF），而不是使用PSF模板。

STPSF的主要功能是生成成像PSF；然而，它可以生成一组单色PSF，我们可以将这些PSF进行组合。

stpsf 仅在此处需要，因此我们在本附录的开头导入它：

from stpsf import NIRSpec, display_psf  # 从stpsf模块导入NIRSpec类和display_psf函数

仪器属性#

有关仪器设置的完整列表，请参见 STPSF 文档。

instrument = NIRSpec()  # 创建NIRSpec仪器的实例

print(instrument.filter_list)  # 打印仪器可用的滤光片列表

# 参考：允许的掩模列表
allowed_masks = ('S200A1', 'S200A2', 'S400A1', 'S1600A1', 'S200B1', 
                 'MSA all open', 'Single MSA open shutter', 
                 'Three adjacent MSA open shutters')  # 定义允许的掩模选项

# 根据需要编辑这些参数
instrument.filter = 'F110W'  # 设置仪器的滤光片为F110W
instrument.image_mask = 'Three adjacent MSA open shutters'  # 设置图像掩模为三个相邻的MSA开放快门

单色点扩散函数（PSFs）#

我们可以使用的最严格的方法是为2D光谱中的每个波长生成一个点扩散函数（PSF），并将它们全部结合。然而，这种方法所需的计算时间和内存通常非常大，除非光谱在色散方向上相当短。一个更合理的方法（这正是我们在这里要做的）是创建一组均匀间隔的单色点扩散函数，并在它们之间进行插值。

psf_wavelengths = np.linspace(wavelength[0], wavelength[-1], num=10) * 1.0e-6  # 生成从最小波长到最大波长的10个波长值，单位转换为米

cube_hdul = instrument.calc_datacube(psf_wavelengths)  # 计算数据立方体，输出为HDUList对象

psf_cube = cube_hdul[1].data  # 从HDUList中提取数据立方体的实际数据

psf_cube.shape  # 获取数据立方体的形状

# 显示数据立方体的内容

fig9, ax9 = plt.subplots(nrows=5, ncols=2, figsize=(8, 12))  # 创建一个5行2列的子图，图形大小为8x12英寸

plt.subplots_adjust(hspace=0.15, wspace=0.01, left=0.06,  # 调整子图之间的间距和边距
                    right=0.94, bottom=0.05, top=0.95)

for row in range(5):  # 遍历行

    for col in range(2):  # 遍历列

        ax = ax9[row, col]  # 获取当前子图的坐标轴

        w = row * 2 + col  # 计算当前的波长索引

        wl = psf_wavelengths[w]  # 获取对应的波长
        
        # 显示当前波长的点扩散函数（PSF）
        display_psf(cube_hdul, ax=ax, cube_slice=w,
                    title="$\lambda$ = {:.3f} $\mu$m".format(wl*1e6),  # 设置标题，显示波长
                    vmax=.2, vmin=1e-4, ext=1, colorbar=False)  # 设置显示的最大值和最小值，不显示颜色条

        ax.xaxis.set_visible(False)  # 隐藏x轴
        ax.yaxis.set_visible(False)  # 隐藏y轴

插值方法#

我们选择的插值方法在很大程度上取决于点扩散函数（PSF）如何随波长变化。为了评估不同的方法，我们将创建另一个单色点扩散函数（PSF）进行比较。

# 计算给定波长（3.0微米）的点扩散函数（PSF）
reference_psf_hdul = instrument.calc_psf(monochromatic=3.0e-6)

# 从PSF数据中提取第二个HDU（Header Data Unit）的数据
reference_psf = reference_psf_hdul[1].data

# 使用简单归一化方法对PSF进行归一化，采用对数伸缩，设置最小值和最大值
ref_norm = simple_norm(reference_psf, stretch='log', vmin=1e-4, vmax=0.2)

最简单的方法是进行三维线性插值，接下来我们来看看效果如何。在下图中，左上角的图像是参考点扩散函数（PSF），右上角是线性插值后的点扩散函数，左下角是差异图像，右下角是参考（X）和插值（Y）点扩散函数的对数-对数图。

ref_pix = reference_psf >= 1e-4  # 创建一个布尔数组，标记参考PSF中大于等于1e-4的像素

psf_x = psf_y = np.arange(48)  # 创建一个包含48个元素的数组，用于PSF的x和y坐标

out_x, out_y = np.meshgrid(psf_x, psf_y, indexing='ij')  # 创建网格坐标，用于插值

interpolator = RegularGridInterpolator((psf_wavelengths, psf_x, psf_y), psf_cube, method='linear')  # 创建线性插值器

linear_psf = interpolator((3.0e-6, out_x, out_y))  # 使用插值器计算给定波长下的线性PSF

diff_lin_psf = reference_psf - linear_psf  # 计算参考PSF与线性PSF之间的差异

# 打印参考PSF的最小值和最大值
print("Reference: min {:.3e}, max {:.3e}".format(reference_psf.min(), reference_psf.max()))

# 打印线性PSF的最小值和最大值
print("Linear: min {:.3e}, max {:.3e}".format(linear_psf.min(), linear_psf.max()))

# 打印差异PSF的最小值和最大值
print("Diff: min {:.3e}, max {:.3e}".format(diff_lin_psf.min(), diff_lin_psf.max()))

# 打印总误差，计算差异PSF的平方和的平方根
print("Total error: {:.5e}".format(np.sqrt((diff_lin_psf**2).sum())))

# 创建一个2x2的子图
figA, axA = plt.subplots(nrows=2, ncols=2, figsize=(10, 10))

plt.subplots_adjust(wspace=0.01, left=0.05, right=0.95)  # 调整子图之间的间距

axA[0, 0].imshow(reference_psf, interpolation='none', norm=ref_norm)  # 显示参考PSF图像
axA[0, 0].xaxis.set_visible(False)  # 隐藏x轴
axA[0, 0].yaxis.set_visible(False)  # 隐藏y轴

axA[0, 1].imshow(linear_psf, interpolation='none', norm=ref_norm)  # 显示线性PSF图像
axA[0, 1].xaxis.set_visible(False)  # 隐藏x轴
axA[0, 1].yaxis.set_visible(False)  # 隐藏y轴

axA[1, 0].imshow(diff_lin_psf, interpolation='none', vmin=-5e-4, vmax=5e-4)  # 显示差异PSF图像
axA[1, 0].xaxis.set_visible(False)  # 隐藏x轴
axA[1, 0].yaxis.set_visible(False)  # 隐藏y轴

axA[1, 1].loglog(reference_psf[ref_pix], linear_psf[ref_pix], 'k+')  # 在对数坐标系中绘制参考PSF与线性PSF的关系
axA[1, 1].set_aspect('equal', 'box')  # 设置坐标轴比例相等

下一个方法计算量更大，但可能更准确。我们逐像素遍历PSF立方体，并沿着波长轴使用一维三次样条插值。

cubic_psf = np.zeros_like(psf_cube[0])  # 创建一个与psf_cube第一层相同形状的零数组，用于存储立方插值结果

for row in np.arange(48):  # 遍历行

    for col in np.arange(48):  # 遍历列

        spline = interp1d(psf_wavelengths, psf_cube[:, row, col], kind='cubic')  # 对每个像素的光谱数据进行立方插值

        cubic_psf[row, col] = spline(3.0e-6)  # 在指定波长处计算插值并存储结果

diff_cub_psf = reference_psf - cubic_psf  # 计算参考PSF与立方插值PSF之间的差异

print("Reference: min {:.3e}, max {:.3e}".format(reference_psf.min(), reference_psf.max()))  # 打印参考PSF的最小值和最大值

print("Cubic: min {:.3e}, max {:.3e}".format(cubic_psf.min(), cubic_psf.max()))  # 打印立方插值PSF的最小值和最大值

print("Diff: min {:.3e}, max {:.3e}".format(diff_cub_psf.min(), diff_cub_psf.max()))  # 打印差异PSF的最小值和最大值

print("Total error: {:.5e}".format(np.sqrt((diff_cub_psf**2).sum())))  # 计算并打印总误差（均方根误差）

figB, axB = plt.subplots(nrows=2, ncols=2, figsize=(10, 10))  # 创建一个2x2的子图

plt.subplots_adjust(wspace=0.01, left=0.05, right=0.95)  # 调整子图之间的间距

axB[0, 0].imshow(reference_psf, interpolation='none', norm=ref_norm)  # 显示参考PSF图像

axB[0, 0].xaxis.set_visible(False)  # 隐藏x轴

axB[0, 0].yaxis.set_visible(False)  # 隐藏y轴

axB[0, 1].imshow(cubic_psf, interpolation='none', norm=ref_norm)  # 显示立方插值PSF图像

axB[0, 1].xaxis.set_visible(False)  # 隐藏x轴

axB[0, 1].yaxis.set_visible(False)  # 隐藏y轴

axB[1, 0].imshow(diff_cub_psf, interpolation='none', vmin=-5e-4, vmax=5e-4)  # 显示差异PSF图像，设置颜色范围

axB[1, 0].xaxis.set_visible(False)  # 隐藏x轴

axB[1, 0].yaxis.set_visible(False)  # 隐藏y轴

axB[1, 1].loglog(reference_psf[ref_pix], cubic_psf[ref_pix], 'k+')  # 在对数坐标系中绘制参考PSF与立方插值PSF的关系

axB[1, 1].set_aspect('equal', 'box')  # 设置坐标轴比例相等

虽然对数-对数（log-log）图看起来与线性情况几乎相同，但在样条（spline）情况下的差异图显示某些中心像素的误差略大。这与“总误差”（total error）统计量（差异图的平方和）一致，在第二种情况下更大。

我们可以在下面的图中看到，这两种方法之间的差异非常小，但线性插值（linearly-interpolated）点扩散函数（PSF）的总误差大约准确了~3倍。

figC = plt.figure()  # 创建一个新的图形对象

plt.loglog(linear_psf[ref_pix], cubic_psf[ref_pix], 'k+')  # 在对数坐标系中绘制线性插值和立方插值的点，使用黑色加号标记

plt.xlabel('Linear interpolation')  # 设置x轴标签为“线性插值”

plt.ylabel('Cubic interpolation')  # 设置y轴标签为“立方插值”

完整的轨迹点扩散函数 (PSF)#

现在我们可以为光谱轨迹生成一个完整的点扩散函数 (PSF)。请注意，在每个波长下，PSF将是重叠的相邻单色PSF的线性组合。如果存在几何畸变，在拟合轨迹中心之后生成这个PSF可能会更有利。

# 创建一个三维网格，分别对应波长、PSF的x和y坐标
cube_w, cube_x, cube_y = np.meshgrid(wavelength * 1e-6, psf_x, psf_y, indexing='ij')

# 使用插值器生成完整的PSF立方体
full_psf_cube = interpolator((cube_w, cube_x, cube_y))

# 获取PSF立方体的形状
nw, ny, nx = full_psf_cube.shape

# 计算y维度的一半，用于后续索引
half = ny // 2

# 初始化一个二维数组，用于存储追踪数据
trace = np.zeros((ny, nw), dtype=float)

# 遍历每个波长和对应的PSF
for wl, psf in enumerate(full_psf_cube):
    # 计算当前波长的下界索引
    lo = wl - half

    # 确保下界索引不小于0
    lo_w = max(lo, 0)

    # 计算下界的x坐标索引
    lo_x = lo_w - lo

    # 计算当前波长的上界索引
    hi = wl + half

    # 确保上界索引不超过nw
    hi_w = min(hi, nw)

    # 计算上界的x坐标索引
    hi_x = nx - (hi - hi_w)

    # 将当前PSF的值累加到trace数组的相应位置
    trace[:, lo_w:hi_w] += psf[:, lo_x:hi_x]

wpsf_aspect = nw / (2. * ny)  # 计算宽度与高度的比例，用于设置图像的纵横比

figD = plt.figure(figsize=(10, 8))  # 创建一个新的图形，设置图形大小为10x8英寸

trace_norm = simple_norm(trace, stretch='log', vmin=1e-4, vmax=0.2)  # 对数据进行归一化处理，使用对数伸缩，设置最小值和最大值

plt.imshow(trace, interpolation='none', aspect=wpsf_aspect, norm=trace_norm)  # 显示图像，设置插值方式、纵横比和归一化

plt.colorbar()  # 添加颜色条以显示数据值的颜色映射

重采样轨迹#

目前，我们的点扩散函数（PSF）数组的大小和位置与提取区域中的轨迹不一致。虽然我们可以通过移动和裁剪来调整到正确的大小，但光谱通常不会恰好位于一个像素中心，并且由于采样不足，PSF中的分数像素偏移可能会导致最终提取中出现显著误差。因此，我们将对光谱在提取区域中的位置进行最终重采样。为此，我们可以使用我们熟悉的 RegularGridInterpolator。我们将STPSF轨迹的中心（最初位于第23行）设置为我们的拟合轨迹中心，并进行适当的重采样。

trace_row = np.arange(ny)  # 创建一个从0到ny的数组，表示行索引

trace_interpolator = RegularGridInterpolator((trace_row, wavelength), trace)  # 创建一个正则网格插值器，用于插值trace数据

center_0 = 23  # 定义中心位置0

center_1 = fit_extraction_kernel.mean_0  # 从拟合提取核中获取中心位置1

out_lo = center_0 - center_1  # 计算输出的下限

out_hi = out_lo + er_ny  # 计算输出的上限

resample_row = np.linspace(out_lo, out_hi, er_ny)  # 在下限和上限之间生成均匀的重采样行

resample_y, resample_w = np.meshgrid(resample_row, wavelength, indexing='ij')  # 创建重采样的y和w网格

resampled_trace = trace_interpolator((resample_y, resample_w))  # 使用插值器对重采样的y和w进行插值，得到重采样的trace

figE, axE = plt.subplots(nrows=2, ncols=1, figsize=(10, 8))  # 创建一个包含2行1列的子图

plt.subplots_adjust(hspace=0.1)  # 调整子图之间的垂直间距

trace_renorm = simple_norm(resampled_trace, stretch='log')  # 对重采样的trace进行简单归一化，使用对数伸缩

axE[0].imshow(resampled_trace, interpolation='none', aspect=aspect_ratio, norm=trace_renorm)  # 在第一个子图中显示重采样的trace

axE[1].imshow(extraction_region, cmap='gray', aspect=aspect_ratio,  # 在第二个子图中显示提取区域
              norm=er_norm, interpolation='none')  # 使用灰度色图和指定的归一化

关于这个笔记本#

作者： Graham Kanarek，科学支持部门的工作人员科学家

更新时间： 2025-02-27，使用 STPSF 代替 WebbPSF。

最佳提取算法改编自 Horne (1986)。

看起来您没有提供任何内容需要翻译。如果您有特定的Markdown内容需要翻译，请将其粘贴在这里，我将很高兴为您提供帮助。

页面顶部

MOS 最优光谱提取

Contents

MOS 最优光谱提取#

引言#

定义术语#

导入库#

加载数据#

最优提取算法#

定义提取区域#

创建核切片#

定义提取核#

选择一个PSF模板#

多项式背景#

拟合提取核#

波长变化的全宽半最大值 (FWHM)（已跳过）#

拟合几何畸变 (跳过)#

定义轨迹拟合的区间#

使用修改后的提取核拟合每个箱体#

使用一维多项式拟合轨迹中心#

构建最终的1D光谱#

创建方差图像#

生成一维光谱#

附录 A：批处理#

定义提取区域#

创建核切片#

创建并拟合提取核#

考虑变化的全宽半最大值（FWHM）#

拟合轨迹中心#

生成一维光谱#

便利函数#

遍历所需文件#

在示例数据集上运行#

附录 B: STPSF#

仪器属性#

单色点扩散函数（PSFs）#

插值方法#

完整的轨迹点扩散函数 (PSF)#

重采样轨迹#

关于这个笔记本#