PSF 光度测量

用例： PSF光度测量，创建PSF，推导颜色-星等图。

数据： 使用MIRAGE获得的关于大麦哲伦云（LMC）天体测量校准场的NIRCam模拟图像，并通过JWST管道处理。这些模拟图像是通过对三个宽带滤光片对（F070W、F115W和F200W用于SW通道；F277W、F356W和F444W用于LW通道）进行4点亚像素抖动获得的。我们仅模拟了1个NIRCam SW探测器（即“ NRCB1”）。在本示例中，我们使用两个SW滤光片（F115W和F200W）的Level-2图像（.cal，已校准但未整形），并在每个滤光片中推导光度。其他滤光片的图像也可用，并可用于测试笔记本和/或不同的滤光片组合。

工具： photutils。

跨仪器： NIRSpec、NIRISS、MIRI。

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

PSF光度测量可以通过以下方式获得：

• 从STPSF获得的单一模型

• 来自STPSF的PSF模型网格

• 单一有效PSF（ePSF）

进行中的工作：

• 创建ePSF网格并使用ePSF网格进行减缩

• 使用ePSF网格来扰动STPSF模型

该笔记本展示了：

• 如何从STPSF获得PSF模型（或构建ePSF）

• 如何对图像进行PSF光度测量

• 如何交叉匹配不同图像的目录

• 如何推导和应用光度零点

最终图表显示：

• 4幅图像的仪器颜色-星等图

• 仪器颜色-星等图及误差

• 星等零点

• 校准的颜色-星等图（与输入颜色-星等图进行比较）

• 输入和输出光度的比较

关于pysynphot的说明：pysynphot的数据文件由校准参考数据系统（Calibration Reference Data System，CRDS）单独分发。它们预计会遵循特定的目录结构，该结构位于根目录下，由必须在使用此软件包之前设置的PYSYN_CDBS环境变量标识。在下面的示例中，根目录被任意命名为/my/local/dir/trds/。

export PYSYN_CDBS=/my/local/dir/trds/

有关数据文件的配置和下载，请参见文档。

导入

import glob as glob  # 导入glob模块，用于文件路径操作

import os  # 导入os模块，用于操作系统功能

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

import requests  # 导入requests模块，用于HTTP请求

import time  # 导入time模块，用于时间操作

import warnings  # 导入warnings模块，用于发出警告

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

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

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

import pandas as pd  # 导入pandas库，用于数据处理

import stpsf  # 导入stpsf模块，用于点扩散函数处理

from astropy import units as u  # 从astropy导入单位模块

from astropy.coordinates import SkyCoord, match_coordinates_sky  # 导入天文坐标相关功能

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

from astropy.modeling.fitting import LevMarLSQFitter  # 导入Levenberg-Marquardt最小二乘拟合器

from astropy.nddata import NDData  # 导入NDData类，用于处理多维数据

from astropy.stats import sigma_clipped_stats  # 导入sigma裁剪统计功能

from astropy.table import QTable, Table  # 导入QTable和Table类，用于表格数据处理

from astropy.utils.exceptions import AstropyUserWarning  # 导入Astropy用户警告类

from astropy.visualization import simple_norm  # 导入简单归一化功能

from jwst.datamodels import ImageModel  # 从jwst.datamodels导入图像模型类

from photutils.aperture import (CircularAnnulus, CircularAperture,  # 导入光度测量相关类
                                aperture_photometry)

from photutils.background import MADStdBackgroundRMS, MMMBackground  # 导入背景噪声估计类

from photutils.detection import DAOStarFinder  # 导入DAO星点查找器

from photutils.psf import (EPSFBuilder, GriddedPSFModel, IterativePSFPhotometry,  # 导入PSF相关功能
                           SourceGrouper, extract_stars)

导入绘图函数

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

from matplotlib import pyplot as plt  # 导入pyplot模块用于绘图

import matplotlib.ticker as ticker  # 导入ticker模块用于设置坐标轴刻度

# 设置图像的颜色映射为'viridis'
plt.rcParams['image.cmap'] = 'viridis'

# 设置图像的原点位置为左下角
plt.rcParams['image.origin'] = 'lower'

# 设置坐标轴标题和标签的字体大小为14
plt.rcParams['axes.titlesize'] = plt.rcParams['axes.labelsize'] = 14

# 设置x轴和y轴刻度标签的字体大小为14
plt.rcParams['xtick.labelsize'] = plt.rcParams['ytick.labelsize'] = 14

# 定义字体样式1
font1 = {'family': 'helvetica', 'color': 'black', 'weight': 'normal', 'size': '12'}

# 定义字体样式2
font2 = {'family': 'helvetica', 'color': 'black', 'weight': 'normal', 'size': '20'}

下载 STPSF 和 Synphot 数据

import os  # 导入操作系统模块
import requests  # 导入请求模块
import tarfile  # 导入tarfile模块
from urllib.parse import 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}")  # 抛出不支持的协议错误

    response = requests.get(url, stream=True, timeout=timeout)  # 发送GET请求下载文件
    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.gz文件
    gzf.extractall(stpsf_folder)  # 解压文件到指定文件夹

import os  # 导入操作系统模块
import tarfile  # 导入tarfile模块以处理tar文件
from urllib import request  # 从urllib导入request模块以处理URL请求

# 设置环境变量
os.environ["PYSYN_CDBS"] = "./grp/redcat/trds/"  # 设置PYSYN_CDBS环境变量，指向数据目录

# Synphot数据
synphot_url = 'http://ssb.stsci.edu/trds/tarfiles/synphot5.tar.gz'  # Synphot数据的URL
synphot_file = './synphot5.tar.gz'  # 下载的tar.gz文件的本地路径
synphot_folder = './grp'  # 存放解压后文件的文件夹路径

# 收集synphot文件
if not os.path.exists(synphot_folder):  # 如果文件夹不存在
    os.makedirs(synphot_folder)  # 创建文件夹
    request.urlretrieve(synphot_url, synphot_file)  # 从URL下载文件到本地
    gzf = tarfile.open(synphot_file)  # 打开下载的tar.gz文件
    gzf.extractall('./', filter='data')  # 解压文件到当前目录

加载图像并创建一些有用的字典

我们加载所有图像，并创建一个包含所有图像的字典，按探测器（detectors）和滤光片（filters）分类。这对于检查可用的探测器和滤光片非常有用，并帮助我们决定是否要对所有图像进行光度测量（photometry），还是仅对子集（例如，仅对SW滤光片）进行测量。

我们还创建了一个包含一些分析所需参数的字典。该字典包含光度零点（photometric zeropoints）（来自MIRAGE配置文件）和NIRCam点扩散函数（PSF）全宽半高（FWHM），这些数据来自NIRCam点扩散函数的JDox页面。FWHM是通过分析使用STPSF模拟的预期NIRCam PSF计算得出的。

注意：一旦在调试后每个探测器的零点和FWHM值可用，该字典将进行更新。

因此，我们有两个字典：

• 单个Level-2校准图像的字典

• 包含其他一些有用参数的字典

# 初始化一个字典，用于存储不同探测器的图像数据
dict_images = {'NRCA1': {}, 'NRCA2': {}, 'NRCA3': {}, 'NRCA4': {}, 'NRCA5': {},
               'NRCB1': {}, 'NRCB2': {}, 'NRCB3': {}, 'NRCB4': {}, 'NRCB5': {}}

# 初始化短波和长波滤光器的字典
dict_filter_short = {}
dict_filter_long = {}

# 初始化短波和长波的滤光器和探测器列表
ff_short = []
det_short = []
det_long = []
ff_long = []
detlist_short = []
detlist_long = []
filtlist_short = []
filtlist_long = []

# 检查当前目录下是否存在处理过的图像文件
if not glob.glob('./*cal*fits'):
    print("Downloading images")  # 如果没有，打印下载信息

    # 定义图像数据的下载链接和保存路径
    boxlink_images_lev2 = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/images_level2.tar.gz'
    boxfile_images_lev2 = './images_level2.tar.gz'

    # 下载图像数据
    request.urlretrieve(boxlink_images_lev2, boxfile_images_lev2)

    # 解压下载的tar文件
    tar = tarfile.open(boxfile_images_lev2, 'r')
    tar.extractall(filter='data')

    # 设置图像目录并获取所有处理过的图像文件
    images_dir = './'
    images = sorted(glob.glob(os.path.join(images_dir, "*cal.fits")))

else:
    # 如果文件已存在，直接获取图像文件
    images_dir = './'
    images = sorted(glob.glob(os.path.join(images_dir, "*cal.fits")))

# 遍历所有图像文件
for image in images:
    im = fits.open(image)  # 打开图像文件
    f = im[0].header['FILTER']  # 获取滤光器信息
    d = im[0].header['DETECTOR']  # 获取探测器信息

    # 根据探测器类型进行调整
    if d == 'NRCBLONG':
        d = 'NRCB5'  # 将NRCBLONG映射到NRCB5
    elif d == 'NRCALONG':
        d = 'NRCA5'  # 将NRCALONG映射到NRCA5
    else:
        d = d  # 保持原样

    wv = float(f[1:3])  # 提取波长信息

    # 根据波长将图像分类为短波或长波
    if wv > 24:         
        ff_long.append(f)  # 添加到长波滤光器列表
        det_long.append(d)  # 添加到长波探测器列表
    else:
        ff_short.append(f)  # 添加到短波滤光器列表
        det_short.append(d)  # 添加到短波探测器列表

    # 获取短波和长波探测器的唯一列表
    detlist_short = sorted(list(dict.fromkeys(det_short)))
    detlist_long = sorted(list(dict.fromkeys(det_long)))

    # 初始化唯一滤光器列表
    unique_list_filters_short = []
    unique_list_filters_long = []

    # 为短波滤光器创建字典
    for x in ff_short:
        if x not in unique_list_filters_short:
            dict_filter_short.setdefault(x, {})  # 如果滤光器不在列表中，则添加

    # 为长波滤光器创建字典
    for x in ff_long:
        if x not in unique_list_filters_long:
            dict_filter_long.setdefault(x, {})  # 如果滤光器不在列表中，则添加

    # 将短波探测器与其对应的滤光器字典关联
    for d_s in detlist_short:
        dict_images[d_s] = dict_filter_short

    # 将长波探测器与其对应的滤光器字典关联
    for d_l in detlist_long:
        dict_images[d_l] = dict_filter_long

    # 获取短波和长波滤光器的唯一列表
    filtlist_short = sorted(list(dict.fromkeys(dict_filter_short)))
    filtlist_long = sorted(list(dict.fromkeys(dict_filter_long)))

    # 将图像文件添加到相应的探测器和滤光器字典中
    if len(dict_images[d][f]) == 0:
        dict_images[d][f] = {'images': [image]}  # 如果没有图像，初始化列表
    else:
        dict_images[d][f]['images'].append(image)  # 否则，添加图像

# 打印可用的探测器和滤光器信息
print("Available Detectors for SW channel:", detlist_short)
print("Available Detectors for LW channel:", detlist_long)
print("Available SW Filters:", filtlist_short)
print("Available LW Filters:", filtlist_long)

# 定义JWST望远镜的滤光片列表
filters = ['F070W', 'F090W', 'F115W', 'F140M', 'F150W2', 'F150W', 'F162M', 'F164N', 'F182M',
           'F187N', 'F200W', 'F210M', 'F212N', 'F250M', 'F277W', 'F300M', 'F322W2', 'F323N',
           'F335M', 'F356W', 'F360M', 'F405N', 'F410M', 'F430M', 'F444W', 'F460M', 'F466N', 'F470N', 'F480M']

# 定义每个滤光片的点扩散函数全宽半最大值（FWHM）
psf_fwhm = [0.987, 1.103, 1.298, 1.553, 1.628, 1.770, 1.801, 1.494, 1.990, 2.060, 2.141, 2.304, 2.341, 1.340,
            1.444, 1.585, 1.547, 1.711, 1.760, 1.830, 1.901, 2.165, 2.179, 2.300, 2.302, 2.459, 2.507, 2.535, 2.574]

# 定义每个滤光片的零点磁度（ZP）对于模型A
zp_modA = [25.7977, 25.9686, 25.8419, 24.8878, 27.0048, 25.6536, 24.6957, 22.3073, 24.8258, 22.1775, 25.3677, 24.3296,
           22.1036, 22.7850, 23.5964, 24.8239, 23.6452, 25.3648, 20.8604, 23.5873, 24.3778, 23.4778, 20.5588,
           23.2749, 22.3584, 23.9731, 21.9502, 20.0428, 19.8869, 21.9002]

# 定义每个滤光片的零点磁度（ZP）对于模型B
zp_modB = [25.7568, 25.9771, 25.8041, 24.8738, 26.9821, 25.6279, 24.6767, 22.2903, 24.8042, 22.1499, 25.3391, 24.2909,
           22.0574, 22.7596, 23.5011, 24.6792, 23.5769, 25.3455, 20.8631, 23.4885, 24.3883, 23.4555, 20.7007,
           23.2763, 22.4677, 24.1562, 22.0422, 20.1430, 20.0173, 22.4086]

# 创建一个字典，将滤光片与其对应的PSF FWHM和零点磁度（模型A和模型B）关联起来
dict_utils = {filters[i]: {'psf fwhm': psf_fwhm[i], 'VegaMAG zp modA': zp_modA[i],
                           'VegaMAG zp modB': zp_modB[i]} for i in range(len(filters))}

选择用于分析的探测器和/或滤光片

如果我们只对分析中的某些滤光片（和/或某些探测器）感兴趣，如本例所示，我们可以从字典中选择这些滤光片（探测器）的二级（Level-2）校准图像，并仅分析这些图像。

在这个特定的例子中，我们分析探测器 NRCB1 的滤光片 F115W 和 F200W 的图像。

dets_short = ['NRCB1']  # 本例中感兴趣的探测器

filts_short = ['F115W', 'F200W']  # 本例中感兴趣的滤光片

显示图像

检查我们的图像是否存在伪影，并确认可以用于分析。

fig, ax = plt.subplots(ncols=2, figsize=(14, 14))  # 创建一个包含2列的子图，图形大小为14x14

for det in dets_short:  # 遍历短名称的探测器列表

    for i, filt in enumerate(filts_short):  # 遍历短名称的滤镜列表，并获取索引

        image = fits.open(dict_images[det][filt]['images'][0])  # 打开指定探测器和滤镜的图像文件

        data_sb = image[1].data  # 获取图像数据

        norm = simple_norm(data_sb, 'sqrt', percent=99.)  # 计算数据的归一化，以平方根方式和99%的百分位数

        ax[i].imshow(data_sb, norm=norm, cmap='Greys')  # 在子图中显示图像数据，使用灰度色图

        ax[i].set_xlabel("X [px]", fontdict=font2)  # 设置X轴标签

        ax[i].set_ylabel("Y [px]", fontdict=font2)  # 设置Y轴标签

        ax[i].set_title(filt, fontdict=font2)  # 设置子图标题为滤镜名称

plt.tight_layout()  # 调整子图布局以避免重叠

创建点扩散函数（PSF）模型

I. 使用 STPSF 创建 PSF 模型

我们创建一个字典，该字典将包含使用 STPSF 为上述选择的探测器和滤光片创建的 PSF。

dict_psfs_stpsf = {}  # 初始化一个空字典，用于存储PSF模型数据

for det in dets_short:  # 遍历短名称的探测器列表

    dict_psfs_stpsf.setdefault(det, {})  # 为每个探测器设置一个空字典，如果不存在则创建

    for j, filt in enumerate(filts_short):  # 遍历短名称的滤光片列表，并获取索引

        dict_psfs_stpsf[det].setdefault(filt, {})  # 为每个滤光片设置一个空字典，如果不存在则创建

        dict_psfs_stpsf[det][filt]['psf model grid'] = None  # 初始化PSF模型网格为None

        dict_psfs_stpsf[det][filt]['psf model single'] = None  # 初始化PSF单模型为None

下面的函数创建一个单一的点扩散函数（PSF）或一个PSF网格，并允许用户将PSF保存为FITS文件。默认情况下，模型PSF存储在psf字典中。对于PSF网格，用户可以选择要创建的PSF数量。PSF可以以探测器采样或过采样的方式创建（过采样因子可以在函数内部更改）。

注意：默认的源光谱是，如果安装了pysynphot，则为Castelli & Kurucz（2004）提供的G2V星光谱。如果没有安装pysynphot，则默认使用一个简单的平坦光谱，使得在每个波长上检测到的光子数量相同。

def create_psf_model(det='NRCB1', fov=11, create_grid=False, num=9, save_psf=False, detsampled=False):
    # 创建PSF模型的函数，参数包括探测器类型、视场大小、是否创建网格等

    nrc = stpsf.NIRCam()  # 初始化NIRCam对象

    nrc.detector = det  # 设置探测器类型

    nrc.filter = filt  # 设置滤光片类型

    src = stpsf.specFromSpectralType('G5V', catalog='phoenix')  # 从光谱类型生成源对象

    if detsampled:
        # 如果选择了探测器采样
        print("Creating a detector sampled PSF")  # 打印信息
        fov = 21  # 设置视场大小为21
    else:
        # 否则创建过采样的PSF
        print("Creating an oversampled PSF")  # 打印信息
        fov = fov  # 保持视场大小不变

    print(f"Using a {fov} px fov")  # 打印使用的视场大小

    if create_grid:
        # 如果选择创建PSF网格
        print("")  # 打印空行
        print(f"Creating a grid of PSF for filter {filt} and detector {det}")  # 打印信息
        print("")  # 打印空行

        outname = f'nircam_{det}_{filt}_fovp{fov}_samp4_npsf{num}.fits'  # 设置输出文件名

        if os.path.exists(outname):
            # 如果输出文件已存在
            grid_psf = GriddedPSFModel.read(outname)  # 从文件读取网格PSF
        else:
            # 否则创建新的网格PSF
            grid_psf = nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False,
                                    fov_pixels=fov, use_detsampled_psf=detsampled,
                                    save=save_psf)  # 生成网格PSF

        dict_psfs_stpsf[det][filt]['psf model grid'] = grid_psf  # 将网格PSF存入字典

    else:
        # 如果选择创建单个PSF
        print("")  # 打印空行
        print(f"Creating a single PSF for filter {filt} and detector {det}")  # 打印信息
        print("")  # 打印空行

        outname = f'nircam_{det}_{filt}_fovp{fov}_samp4_npsf{num}.fits'  # 设置输出文件名

        if os.path.exists(outname):
            # 如果输出文件已存在
            single_psf = GriddedPSFModel.read(outname)  # 从文件读取单个PSF
        else:
            # 否则创建新的单个PSF
            single_psf = nrc.psf_grid(num_psfs=1, oversample=4, source=src, all_detectors=False,
                                      fov_pixels=fov, use_detsampled_psf=detsampled,
                                      save=save_psf)  # 生成单个PSF

        dict_psfs_stpsf[det][filt]['psf model single'] = single_psf  # 将单个PSF存入字典

    return dict_psfs_stpsf  # 返回包含PSF模型的字典

单一点扩散函数模型 (Single PSF model)

for det in dets_short:  # 遍历短列表中的每个探测器

    for filt in filts_short:  # 遍历短列表中的每个滤光片

        create_psf_model(fov=11, num=1, create_grid=False, save_psf=True, detsampled=False)  # 创建点扩散函数模型

显示单个点扩散函数（PSF）模型

fig, ax = plt.subplots(ncols=2, figsize=(14, 14))  # 创建一个包含2列的子图，图像大小为14x14英寸

for det in dets_short:  # 遍历短名称的探测器列表

    for i, filt in enumerate(filts_short):  # 遍历短名称的滤镜列表，并获取索引

        img = dict_psfs_stpsf[det][filt]['psf model single'].data[0]  # 获取指定探测器和滤镜的PSF模型图像数据

        norm_epsf = simple_norm(img, 'log', percent=99.)  # 使用对数归一化方法对图像进行归一化处理，99%为最大值

        ax[i].imshow(img, norm=norm_epsf)  # 在对应的子图中显示图像，应用归一化

        ax[i].set_xlabel('X [px]', fontdict=font2)  # 设置X轴标签

        ax[i].set_ylabel('Y [px]', fontdict=font2)  # 设置Y轴标签

        ax[i].set_title(filt, fontdict=font2)  # 设置子图标题为当前滤镜名称

plt.tight_layout()  # 调整子图布局以避免重叠

PSF 网格

for det in dets_short:  # 遍历短列表中的每个探测器

    for filt in filts_short:  # 遍历短列表中的每个滤光片

        create_psf_model(fov=11, num=25, create_grid=True, save_psf=True, detsampled=False)  
        # 创建点扩散函数模型，参数包括视场(fov)、样本数量(num)、是否创建网格(create_grid)、是否保存PSF(save_psf)以及探测器是否采样(detsampled)

显示PSF网格

我们展示了1个滤光片（F115W）的PSF网格及其与均值的差异。

# 从字典中获取指定的PSF模型网格
griddedpsfmodel = dict_psfs_stpsf[dets_short[0]][filts_short[0]]['psf model grid']

# 绘制PSF模型网格，并设置图形大小为10x10英寸
fig = griddedpsfmodel.plot_grid(figsize=(10, 10))

# 使用griddedpsfmodel对象绘制网格图，设置图像大小为10x10
fig = griddedpsfmodel.plot_grid(figsize=(10, 10), deltas=True, cmap='viridis', vmax_scale=0.3)

II. 创建有效点扩散函数模型 (Effective PSF, ePSF)

有关 photutils 有效点扩散函数的更多信息，请参见 这里。

• 从我们想要用于构建点扩散函数 (PSF) 的图像中选择恒星。我们使用 DAOStarFinder 函数在图像中查找亮星（设置较高的检测阈值）。 DAOStarFinder 使用 DAOFIND (Stetson 1987) 算法在图像中检测恒星。DAOFIND 在图像中搜索局部密度最大值，这些最大值的峰值幅度大于 threshold（大约；阈值应用于卷积图像），并且具有与定义的二维高斯核相似的大小和形状。

注意：阈值和与最近邻的最大距离取决于用户的科学案例（即视场中的恒星数量、拥挤程度、亮源数量、构建 ePSF 所需的最小恒星数量等），必须相应地进行修改。

• 使用 EPSBuilder 函数构建有效点扩散函数（排除边界框超出探测器边缘的物体）。

我们创建一个字典，包含上述选择的探测器和滤光片的有效点扩散函数（PSF）。

dict_psfs_epsf = {}  # 初始化一个空字典，用于存储PSF和ePSF数据

for det in dets_short:  # 遍历短名称的探测器列表

    dict_psfs_epsf.setdefault(det, {})  # 为每个探测器设置一个空字典，如果不存在则创建

    for j, filt in enumerate(filts_short):  # 遍历短名称的滤光片列表，并获取索引

        dict_psfs_epsf[det].setdefault(filt, {})  # 为每个滤光片设置一个空字典，如果不存在则创建

        dict_psfs_epsf[det][filt]['table psf stars'] = {}  # 初始化表格PSF星体的字典

        dict_psfs_epsf[det][filt]['epsf single'] = {}  # 初始化单个ePSF的字典

        dict_psfs_epsf[det][filt]['epsf grid'] = {}  # 初始化ePSF网格的字典

        for i in np.arange(0, len(dict_images[det][filt]['images']), 1):  # 遍历图像列表的索引

            dict_psfs_epsf[det][filt]['table psf stars'][i + 1] = None  # 为表格PSF星体分配None值

            dict_psfs_epsf[det][filt]['epsf single'][i + 1] = None  # 为单个ePSF分配None值

            dict_psfs_epsf[det][filt]['epsf grid'][i + 1] = None  # 为ePSF网格分配None值

请注意，来自处理管道的二级（Level-2）和三级（Level-3）图像的单位是 MJy/sr（因此是表面亮度）。图像的实际单位可以通过头文件关键字 BUNIT 进行检查。标量转换常数被复制到头文件关键字 PHOTMJSR，该关键字提供了从 DN/s 到 MJy/立体角的转换。对于我们的分析，我们将数据转换回 DN/s。

def find_stars_epsf(img_num, filt_num, det='NRCA1', filt='F070W', dist_sel=False):
    # 定义函数，查找图像中的PSF星星

    bkgrms = MADStdBackgroundRMS()  # 初始化背景噪声计算类
    mmm_bkg = MMMBackground()  # 初始化背景计算类

    image = fits.open(dict_images[det][filt]['images'][img_num])  # 打开指定图像文件
    data_sb = image[1].data  # 获取图像数据
    imh = image[1].header  # 获取图像头信息

    print(f"Finding PSF stars on image {img_num + 1} of filter {filt}, detector {det}")  # 打印当前处理的图像信息

    data = data_sb / imh['PHOTMJSR']  # 将数据转换为DN/s
    units = imh['BUNIT']  # 获取数据单位
    print(f"Conversion factor from {units} to DN/s for filter {filt}: {imh['PHOTMJSR']}")  # 打印转换因子

    sigma_psf = dict_utils[filt]['psf fwhm']  # 获取PSF的FWHM值
    print(f"FWHM for the filter {filt}: {sigma_psf} px")  # 打印FWHM值

    std = bkgrms(data)  # 计算背景噪声的标准差
    bkg = mmm_bkg(data)  # 计算背景值

    # 初始化DAOStarFinder以查找星星
    daofind = DAOStarFinder(threshold=th[filt_num] * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
                            sharplo=0.30, sharphi=1.40)

    psf_stars = daofind(data)  # 查找星星并返回结果
    dict_psfs_epsf[det][filt]['table psf stars'][img_num + 1] = psf_stars  # 将结果存入字典

    if dist_sel:  # 如果选择计算距离
        print("")  # 打印空行
        print("Calculating closest neighbour distance")  # 打印正在计算最近邻距离

        d = []  # 初始化距离列表

        # 初始化DAOStarFinder以查找所有星星
        daofind_tot = DAOStarFinder(threshold=10 * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
                                    sharplo=0.30, sharphi=1.40)

        stars_tot = daofind_tot(data)  # 查找所有星星

        x_tot = stars_tot['xcentroid']  # 获取所有星星的x坐标
        y_tot = stars_tot['ycentroid']  # 获取所有星星的y坐标

        for xx, yy in zip(psf_stars['xcentroid'], psf_stars['ycentroid']):  # 遍历PSF星星
            sep = []  # 初始化分离距离列表
            dist = np.sqrt((x_tot - xx)**2 + (y_tot - yy)**2)  # 计算与所有星星的距离
            sep = np.sort(dist)[1:2][0]  # 找到最近的星星距离
            d.append(sep)  # 将距离添加到列表中

        psf_stars['min distance'] = d  # 将最小距离添加到PSF星星表中
        mask_dist = (psf_stars['min distance'] > min_sep[filt_num])  # 创建掩码以筛选星星

        psf_stars = psf_stars[mask_dist]  # 应用掩码筛选星星
        dict_psfs_epsf[det][filt]['table psf stars'][img_num + 1] = psf_stars  # 更新字典中的结果

        print("Minimum distance required:", min_sep[filt_num], "px")  # 打印所需的最小距离
        print("")  # 打印空行
        print(f"Number of isolated sources found in the image used to build ePSF for {filt}: {len(psf_stars)}")  # 打印找到的孤立源数量
        print("-----------------------------------------------------")  # 打印分隔线
        print("")  # 打印空行
    else:  # 如果不计算距离
        print("")  # 打印空行
        print(f"Number of sources used to build ePSF for {filt}: {len(psf_stars)}")  # 打印用于构建ePSF的源数量
        print("--------------------------------------------")  # 打印分隔线
        print("")  # 打印空行

tic = time.perf_counter()  # 记录开始时间

th = [700, 500]  # 两个滤波器的阈值水平（长度必须与分析的滤波器数量匹配）

min_sep = [10, 10]  # ePSF星星与最近邻星星之间可接受的最小间隔

for det in dets_short:  # 遍历短时间检测器列表

    for j, filt in enumerate(filts_short):  # 遍历短时间滤波器列表，并获取索引

        for i in np.arange(0, len(dict_images[det][filt]['images']), 1):  # 遍历当前检测器和滤波器的图像列表

            find_stars_epsf(i, j, det=det, filt=filt, dist_sel=False)  # 调用函数查找ePSF星星

toc = time.perf_counter()  # 记录结束时间

print("Elapsed Time for finding stars:", toc - tic)  # 输出查找星星所用的时间

II. 构建有效的点扩散函数 (PSF)

def build_epsf(det='NRCA1', filt='F070W'):
    # 创建一个函数来构建电子点扩散函数（ePSF）

    mmm_bkg = MMMBackground()  # 初始化背景估计器

    image = fits.open(dict_images[det][filt]['images'][i])  # 打开指定探测器和滤光片的图像文件
    data_sb = image[1].data  # 获取图像数据
    imh = image[1].header  # 获取图像头信息

    data = data_sb / imh['PHOTMJSR']  # 将数据标准化为每平方弧秒的光度

    hsize = (sizes[j] - 1) / 2  # 计算半尺寸

    x = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['xcentroid']  # 获取星星的x坐标
    y = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['ycentroid']  # 获取星星的y坐标

    # 创建掩码，筛选出在图像有效范围内的星星
    mask = ((x > hsize) & (x < (data.shape[1] - 1 - hsize)) & 
            (y > hsize) & (y < (data.shape[0] - 1 - hsize)))

    stars_tbl = Table()  # 创建一个表格来存储星星信息
    stars_tbl['x'] = x[mask]  # 将符合条件的x坐标存入表格
    stars_tbl['y'] = y[mask]  # 将符合条件的y坐标存入表格

    bkg = mmm_bkg(data)  # 计算背景

    data_bkgsub = data.copy()  # 复制数据以进行背景减法
    data_bkgsub -= bkg  # 从数据中减去背景

    nddata = NDData(data=data_bkgsub)  # 创建NDData对象以存储背景减去后的数据
    stars = extract_stars(nddata, stars_tbl, size=sizes[j])  # 提取星星

    print(f"Creating ePSF for image {i + 1} of filter {filt}, detector {det}")  # 打印当前处理的信息

    epsf_builder = EPSFBuilder(oversampling=oversample, maxiters=3, progress_bar=False)  # 初始化ePSF构建器

    epsf, fitted_stars = epsf_builder(stars)  # 构建ePSF并获取拟合的星星
    dict_psfs_epsf[det][filt]['epsf single'][i + 1] = epsf  # 将构建的ePSF存入字典

注意：在这里我们将最大迭代次数限制为3次（以限制运行时间），但实际上应该使用大约10次或更多的迭代。

tic = time.perf_counter()  # 记录开始时间

sizes = [11, 11]  # 每个PSF星体的切割区域大小 - 必须与分析的滤光片数量匹配

oversample = 4  # 过采样因子

for det in dets_short:  # 遍历短时间内的探测器

    for j, filt in enumerate(filts_short):  # 遍历短时间内的滤光片

        for i in np.arange(0, len(dict_images[det][filt]['images']), 1):  # 遍历每个图像

            with warnings.catch_warnings():  # 捕获警告

                # 忽略关于星体靠近图像边缘的警告
                warnings.simplefilter("ignore", category=AstropyUserWarning) 

                build_epsf(det=det, filt=filt)  # 构建有效的PSF

toc = time.perf_counter()  # 记录结束时间

print("Time to build the Effective PSF:", toc - tic)  # 输出构建有效PSF所需的时间

显示ePSFs

我们仅为每个滤光片显示1个ePSF。

fig, ax = plt.subplots(ncols=2, figsize=(14, 14))  # 创建一个包含2列的子图，图像大小为14x14

for det in dets_short:  # 遍历短名称的探测器列表

    for i, filt in enumerate(filts_short):  # 遍历短名称的滤光片列表，并获取索引i

        img = dict_psfs_epsf[det][filt]['epsf single'][i + 1].data  # 从字典中获取对应探测器和滤光片的图像数据

        norm_epsf = simple_norm(img, 'log', percent=99.)  # 使用对数归一化方法对图像进行归一化，99%分位数

        ax[i].imshow(img, norm=norm_epsf)  # 在第i个子图中显示归一化后的图像

        ax[i].set_title(filt, fontdict=font2)  # 设置子图标题为当前滤光片名称，使用font2字体字典

进行中的工作 - 构建有效的PSF网格

两个函数：

• 统计网格中的PSF星星数量

• 创建网格化的ePSF

第一个函数的目的是统计在用户定义的网格编号N所划分的每个子区域中有多少个良好的PSF星星。该函数应从用户提供的数字开始，并迭代直到最小网格大小为2x2。根据用户希望在网格每个单元中包含的PSF星星数量，他们可以选择合适的网格大小或修改星星检测的阈值，这些阈值是在创建单个ePSF时选择的（见上面的寻找星星单元）。

第二个函数使用EPSFBuilder创建PSF网格。该函数将返回一个GriddedEPSFModel对象，其中包含一个N × n × n的3D数组。该3D数组表示创建的N个2D n × n ePSF。它应包括一个grid_xypos键，该键将指明每个PSF在探测器上的位置。grid_xypos中的元组顺序对应于3D数组中PSF的编号。

I. 统计网格中每个区域的PSF星星数量

def count_PSFstars_grid(grid_points=5, size=15, min_numpsf=40):
    # 定义一个函数，用于计算给定网格中PSF星星的数量

    num_grid_calc = np.arange(2, grid_points + 1, 1)  # 创建一个从2到grid_points的数组
    num_grid_calc = num_grid_calc[::-1]  # 将数组反转

    for num in num_grid_calc:  # 遍历每个网格大小
        print(f"Calculating the number of PSF stars in a {num} x {num} grid")  # 打印当前计算的网格大小
        print("")  # 打印空行

        image = fits.open(dict_images[det][filt]['images'][i])  # 打开指定的图像文件
        data_sb = image[1].data  # 获取图像数据

        points = np.int16((data_sb.shape[0] / num) / 2)  # 计算每个网格的点数
        x_center = np.arange(points, 2 * points * (num), 2 * points)  # 计算x轴中心点
        y_center = np.arange(points, 2 * points * (num), 2 * points)  # 计算y轴中心点

        centers = np.array(np.meshgrid(x_center, y_center)).T.reshape(-1, 2)  # 创建中心点的网格坐标

        for n, val in enumerate(centers):  # 遍历每个中心点
            x = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['xcentroid']  # 获取PSF星星的x坐标
            y = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['ycentroid']  # 获取PSF星星的y坐标
            # flux = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['flux']  # 获取PSF星星的flux（注释掉）

            half_size = (size - 1) / 2  # 计算网格的一半大小

            lim1 = val[0] - points + half_size  # 计算x轴的下限
            lim2 = val[0] + points - half_size  # 计算x轴的上限
            lim3 = val[1] - points + half_size  # 计算y轴的下限
            lim4 = val[1] + points - half_size  # 计算y轴的上限

            test = (x > lim1) & (x < lim2) & (y > lim3) & (y < lim4)  # 检查PSF星星是否在当前网格内

            if np.count_nonzero(test) < min_numpsf:  # 如果当前网格内的PSF星星数量少于最小要求
                print(f"Center Coordinates of grid cell {i + 1} are ({val[0]}, {val[1]}) --- Not enough PSF stars in the cell (number of PSF stars < {min_numpsf})")  # 打印不足的提示
            else:  # 如果当前网格内的PSF星星数量满足要求
                print(f"Center Coordinate of grid cell {n + 1} are ({val[0]}, {val[1]}) --- Number of PSF stars: {np.count_nonzero(test)}")  # 打印满足的提示
        print("")  # 打印空行

for det in dets_short:  # 遍历短列表中的每个探测器

    for j, filt in enumerate(filts_short):  # 遍历短列表中的每个滤光片，并获取其索引

        for i in np.arange(0, len(dict_images[det][filt]['images']), 1):  # 遍历当前探测器和滤光片下的每个图像

            print(f"Analyzing image {i + 1} of filter {filt}, detector {det}")  # 打印当前正在分析的图像信息

            print("")  # 打印空行以便于阅读

            count_PSFstars_grid(grid_points=5, size=15, min_numpsf=40)  # 计算PSF星星的网格

TODO - 创建ePSF网格

以下是创建 ePSF 网格的函数，可以将其保存到 epsf 字典中。

TODO - 使用 ePSF 网格来扰动 STPSF 模型

以下是创建一个PSF模型网格的函数，该网格是通过对上述创建的ePSF网格进行扰动得到的STPSF PSF模型。

执行PSF光度测量

我们对图像执行PSF光度测量，默认情况下将输出目录和残差图像保存在下面创建的字典中。也可以使用参数save_output和save_residuals将输出目录（pickle pandas对象）和残差图像（fits文件）保存在当前目录中。

dict_phot = {}  # 初始化一个空字典，用于存储光度数据

for det in dets_short:  # 遍历短名称的探测器列表

    dict_phot.setdefault(det, {})  # 为每个探测器设置一个空字典，如果不存在则创建

    for j, filt in enumerate(filts_short):  # 遍历短名称的滤光片列表，并获取索引

        dict_phot[det].setdefault(filt, {})  # 为每个滤光片设置一个空字典，如果不存在则创建

        dict_phot[det][filt]['residual images'] = {}  # 初始化残差图像字典

        dict_phot[det][filt]['output photometry tables'] = {}  # 初始化输出光度表字典

        for i in np.arange(0, len(dict_images[det][filt]['images']), 1):  # 遍历当前探测器和滤光片的图像列表

            dict_phot[det][filt]['residual images'][i + 1] = None  # 为每个图像索引初始化残差图像为None

            dict_phot[det][filt]['output photometry tables'][i + 1] = None  # 为每个图像索引初始化输出光度表为None

注意：为了加快笔记本的运行速度，我们在查找算法中使用了一个较高的阈值（阈值 ~ 2000），并且在下面的分析中，我们将使用从之前的减缩运行中获得的sigma阈值 = 10的目录。为了进行有意义的数据减缩，用户应相应地修改阈值。

在这里，我们使用STPSF PSF的网格作为PSF模型，但用户可以更改模型并使用其他可用的模型（即，单个STPSF PSF，单个ePSF），通过修改函数中的psf参数。

def psf_phot(det='NRCA1', filt='F070W', th=2000, psf='grid_stpsf', save_residuals=False, save_output=False):
    # 定义函数 psf_phot，接受多个参数，包括探测器、滤光片、阈值、PSF类型等

    bkgrms = MADStdBackgroundRMS()  # 创建背景噪声标准差计算对象
    mmm_bkg = MMMBackground()  # 创建背景计算对象
    fitter = LevMarLSQFitter()  # 创建最小二乘拟合器

    im = fits.open(dict_images[det][filt]['images'][i])  # 打开指定探测器和滤光片的图像
    imh = im[1].header  # 获取图像头信息
    data_sb = im[1].data  # 获取图像数据

    d = im[0].header['DETECTOR']  # 获取探测器信息
    prim_dith_pos = im[0].header['PATT_NUM']  # 获取主偏移位置
    prim_dith_num = im[0].header['NUMDTHPT']  # 获取主偏移点数量
    subpx_dith_pos = im[0].header['SUBPXNUM']  # 获取子像素偏移位置
    subpx_dith_num = im[0].header['SUBPXPNS']  # 获取子像素偏移点数量

    data = data_sb / imh['PHOTMJSR']  # 将数据转换为DN/s单位
    units = imh['BUNIT']  # 获取数据单位
    print(f"Conversion factor from {units} to DN/s for filter {filt}: {imh['PHOTMJSR']}")  # 打印转换因子
    print("Applying conversion to the data")  # 打印正在应用转换

    sigma_psf = dict_utils[filt]['psf fwhm']  # 获取PSF的全宽半最大值
    print(f"FWHM for the filter {filt}: {sigma_psf}")  # 打印FWHM值

    std = bkgrms(data)  # 计算数据的背景噪声标准差
    bkg = mmm_bkg(data)  # 计算数据的背景值

    daofind = DAOStarFinder(threshold=th * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
                            sharplo=0.30, sharphi=1.40)  # 创建星源查找器

    grouper = SourceGrouper(5.0 * sigma_psf)  # 创建源分组器

    # grid PSF
    if psf == 'grid_stpsf':  # 如果选择网格PSF
        print("Using as PSF model STPSF PSFs grid")  # 打印使用网格PSF模型
        psf_model = dict_psfs_stpsf[det][filt]['psf model grid'].copy()  # 复制网格PSF模型

    # single psf:
    if psf == 'single_stpsf':  # 如果选择单个PSF
        print("Using as PSF model STPSF single PSF")  # 打印使用单个PSF模型
        psf_model = dict_psfs_stpsf[det][filt]['psf model single'].copy()  # 复制单个PSF模型

    # epsf:
    if psf == 'single_epsf':  # 如果选择单个ePSF
        print("Using as PSF model single ePSF")  # 打印使用单个ePSF模型
        psf_model = dict_psfs_epsf[det][filt]['epsf single'][i + 1].copy()  # 复制单个ePSF模型

    print(f"Performing the photometry on image {i + 1} of filter {filt}, detector {det}")  # 打印正在进行光度测量的信息

    tic = time.perf_counter()  # 记录开始时间

    data_sub = data - mmm_bkg(data)  # 从数据中减去背景
    psf_shape = (11, 11)  # 定义PSF形状

    phot = IterativePSFPhotometry(psf_model, psf_shape, daofind,
                                  grouper=grouper, fitter=fitter,
                                  fitter_maxiters=500,
                                  maxiters=2, aperture_radius=ap_radius[j])  # 创建迭代PSF光度测量对象
    result = phot(data_sub)  # 执行光度测量

    toc = time.perf_counter()  # 记录结束时间

    dtime = (toc - tic)  # 计算光度测量所需时间
    print(f"Time needed to perform photometry on image {i + 1}: {dtime:.2f} sec")  # 打印所需时间
    print(f"Number of sources detected in image {i + 1} for filter {filt}: {len(result)}")  # 打印检测到的源数量

    residual_image = phot.make_residual_image(data_sub, psf_shape=psf_shape)  # 创建残差图像

    dict_phot[det][filt]['residual images'][i + 1] = residual_image  # 保存残差图像
    dict_phot[det][filt]['output photometry tables'][i + 1] = result  # 保存光度测量结果

    # save the residual images as fits file:
    if save_residuals:  # 如果选择保存残差图像
        hdu = fits.PrimaryHDU(residual_image)  # 创建FITS主HDU
        hdul = fits.HDUList([hdu])  # 创建HDU列表
        residual_outname = f'residual_{d}_{filt}_STPSF_gridPSF_{prim_dith_pos}of{prim_dith_num}_{subpx_dith_pos}of{subpx_dith_num}.fits'  # 定义输出文件名

        dir_output_phot = './'  # 定义输出目录

        hdul.writeto(os.path.join(dir_output_phot, residual_outname))  # 写入FITS文件
        outname = 'phot_{d}_{filt}_STPSF_gridPSF_level2_{prim_dith_pos}of{prim_dith_num}_{subpx_dith_pos}of{subpx_dith_num}.pkl'  # 定义光度测量结果文件名

    # save the output photometry Tables
    if save_output:  # 如果选择保存光度测量结果
        tab = result.to_pandas()  # 将结果转换为Pandas DataFrame
        tab.to_pickle(os.path.join(dir_output_phot, outname))  # 保存为pickle文件

tic_tot = time.perf_counter()  # 记录开始时间

ap_radius = [3.0, 3.5]  # 光圈半径，必须与分析的滤波器数量匹配

# 检查当前目录是否存在残差图像文件
if glob.glob('./*residual*.fits'):
    print("Deleting Residual images from directory")  # 打印删除信息
    files = glob.glob('./residual*.fits')  # 获取所有残差图像文件
    for file in files:
        os.remove(file)  # 删除每个残差图像文件

# 遍历短名称探测器
for det in dets_short:
    # 遍历短名称滤波器
    for j, filt in enumerate(filts_short):
        # 遍历当前探测器和滤波器的图像数量
        for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
            # 调用psf_phot函数进行光度测量
            psf_phot(det=det, filt=filt, th=2000, psf='grid_stpsf', save_residuals=True, save_output=False) 

toc_tot = time.perf_counter()  # 记录结束时间

number = len(filts_short) * len(dict_images[det][filt]['images'])  # 计算处理的图像总数
dtime = (toc_tot - tic_tot)  # 计算总耗时

# 打印处理所有图像所耗费的时间
print(f"Time elapsed to perform the photometry of the {number} images: {dtime:.2f} sec")

输出光度表

# 访问字典 dict_phot 中的 NRCB1 关键字
dict_phot['NRCB1']  # 获取 NRCB1 的数据

# 在 NRCB1 中访问 F115W 关键字
dict_phot['NRCB1']['F115W']  # 获取 F115W 的数据

# 访问 output photometry tables 关键字
dict_phot['NRCB1']['F115W']['output photometry tables']  # 获取输出光度表的数据

# 获取 output photometry tables 中的第二个元素（索引为1）
dict_phot['NRCB1']['F115W']['output photometry tables'][1]  # 返回第二个光度表

显示减法图像

作为一个示例，我们展示了一个科学图像与经过数据处理后的残差图像之间的比较，适用于两个滤光器。请注意，残差图像是通过在上面的单元中进行光度测量（photometry）时使用非常高的检测阈值（detection threshold）获得的。

fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(14, 14))  # 创建一个2x2的子图，图像大小为14x14英寸

for det in dets_short:  # 遍历短名称的探测器

    for i, filt in enumerate(filts_short):  # 遍历短名称的滤镜，并获取索引

        image = fits.open(dict_images[det][filt]['images'][0])  # 打开对应探测器和滤镜的图像文件

        data_sb = image[1].data  # 获取图像数据

        norm = simple_norm(data_sb, 'sqrt', percent=99.)  # 计算归一化参数，使用平方根法和99%的百分位数

        ax[0, i].imshow(data_sb, norm=norm, cmap='Greys')  # 在第一个子图中显示图像数据，使用灰度色图

        ax[0, i].set_xlabel("X [px]", fontdict=font2)  # 设置X轴标签

        ax[0, i].set_ylabel("Y [px]", fontdict=font2)  # 设置Y轴标签

        ax[0, i].set_title(filt, fontdict=font2)  # 设置子图标题为滤镜名称

for det in dets_short:  # 遍历短名称的探测器

    for i, filt in enumerate(filts_short):  # 遍历短名称的滤镜，并获取索引

        res = dict_phot[det][filt]['residual images'][1]  # 获取对应探测器和滤镜的残差图像数据

        norm = simple_norm(res, 'sqrt', percent=99.)  # 计算归一化参数，使用平方根法和99%的百分位数

        ax[1, i].imshow(res, norm=norm, cmap='Greys')  # 在第二个子图中显示残差图像数据，使用灰度色图

        ax[1, i].set_xlabel("X [px]", fontdict=font2)  # 设置X轴标签

        ax[1, i].set_ylabel("Y [px]", fontdict=font2)  # 设置Y轴标签

plt.tight_layout()  # 调整子图布局以避免重叠

第二部分 - 数据分析

注意：在这里，我们使用通过STPSF PSF网格获得的降噪结果作为PSF模型。用户可以使用不同的PSF模型（单一PSF模型、PSF网格等）进行数据分析，并比较结果。

加载具有点扩散函数（PSF）光度的表格

# 如果当前目录下没有符合条件的文件
if not glob.glob('./*phot*gridPSF*.pkl'):
    
    print("Downloading Photometry Output")  # 打印下载提示

    # F115W波段的光度数据链接
    boxlink_cat_f115w = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/phot_cat_F115W.tar.gz'
    boxfile_cat_f115w = './phot_cat_F115W.tar.gz'  # 下载后的文件名

    # 下载F115W数据
    request.urlretrieve(boxlink_cat_f115w, boxfile_cat_f115w)

    # 打开下载的tar文件
    tar = tarfile.open(boxfile_cat_f115w, 'r')
    tar.extractall(filter='data')  # 解压到指定目录

    # F200W波段的光度数据链接
    boxlink_cat_f200w = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/phot_cat_F200W.tar.gz'
    boxfile_cat_f200w = './phot_cat_F200W.tar.gz'  # 下载后的文件名

    # 下载F200W数据
    request.urlretrieve(boxlink_cat_f200w, boxfile_cat_f200w)

    # 打开下载的tar文件
    tar = tarfile.open(boxfile_cat_f200w, 'r')
    tar.extractall(filter='data')  # 解压到指定目录

    cat_dir = './'  # 数据目录

    # 获取F115W和F200W波段的光度数据文件列表
    phots_pkl_f115w = sorted(glob.glob(os.path.join(cat_dir, '*F115W*gridPSF*.pkl')))
    phots_pkl_f200w = sorted(glob.glob(os.path.join(cat_dir, '*F200W*gridPSF*.pkl')))                       

else:
    cat_dir = './'  # 数据目录

    # 获取F115W和F200W波段的光度数据文件列表
    phots_pkl_f115w = sorted(glob.glob(os.path.join(cat_dir, '*F115W*gridPSF*.pkl')))
    phots_pkl_f200w = sorted(glob.glob(os.path.join(cat_dir, '*F200W*gridPSF*.pkl')))                      

results_f115w = []  # 存储F115W结果的列表
results_f200w = []  # 存储F200W结果的列表

# 遍历F115W和F200W的光度数据文件
for phot_pkl_f115w, phot_pkl_f200w in zip(phots_pkl_f115w, phots_pkl_f200w):
    
    ph_f115w = pd.read_pickle(phot_pkl_f115w)  # 读取F115W数据
    ph_f200w = pd.read_pickle(phot_pkl_f200w)  # 读取F200W数据

    result_f115w = QTable.from_pandas(ph_f115w)  # 将F115W数据转换为QTable格式
    result_f200w = QTable.from_pandas(ph_f200w)  # 将F200W数据转换为QTable格式

    results_f115w.append(result_f115w)  # 将结果添加到列表中
    results_f200w.append(result_f200w)  # 将结果添加到列表中

将图像转换为数据模型（DataModel）

为了分配WCS坐标并进行图像的交叉匹配，我们需要将图像转换为数据模型（DataModel）。坐标是在JWST管道的步骤assign_wcs中分配的，这使我们能够对每个滤光片获得的不同目录进行交叉匹配。

images_f115w = []  # 初始化F115W图像列表

images_f200w = []  # 初始化F200W图像列表

# 遍历F115W图像字典中的每个图像
for i in np.arange(0, len(dict_images['NRCB1']['F115W']['images']), 1):
    
    # 创建F115W图像模型
    image_f115w = ImageModel(dict_images['NRCB1']['F115W']['images'][i])
    
    # 将图像模型添加到F115W图像列表中
    images_f115w.append(image_f115w)
        

# 遍历F200W图像字典中的每个图像
for i in np.arange(0, len(dict_images['NRCB1']['F200W']['images']), 1):
    
    # 创建F200W图像模型
    image_f200w = ImageModel(dict_images['NRCB1']['F200W']['images'][i])
    
    # 将图像模型添加到F200W图像列表中
    images_f200w.append(image_f200w)

对两个滤镜的四幅图像进行目录交叉匹配

我们对目录进行交叉匹配，以获得单一的颜色-星等图（color-magnitude diagrams）。

如果匹配之间的距离小于 0.5 像素（px），则将两个滤镜中的星体关联起来。

results_clean_f115w = []  # 初始化F115W清理结果列表

results_clean_f200w = []  # 初始化F200W清理结果列表

# 遍历F115W和F200W图像的索引
for i in np.arange(0, len(images_f115w), 1):

    # 创建F115W的掩码，筛选有效的x_fit, y_fit和flux_fit
    mask_f115w = ((results_f115w[i]['x_fit'] > 0) & (results_f115w[i]['x_fit'] < 2048) &
                  (results_f115w[i]['y_fit'] > 0) & (results_f115w[i]['y_fit'] < 2048) &
                  (results_f115w[i]['flux_fit'] > 0))

    # 根据掩码过滤F115W结果
    result_clean_f115w = results_f115w[i][mask_f115w]

    # 使用WCS转换获取F115W的RA和Dec
    ra_f115w, dec_f115w = images_f115w[i].meta.wcs(result_clean_f115w['x_fit'], result_clean_f115w['y_fit'])

    # 创建SkyCoord对象以存储RA和Dec
    radec_f115w = SkyCoord(ra_f115w, dec_f115w, unit='deg')

    # 将RA和Dec添加到清理结果中
    result_clean_f115w['radec'] = radec_f115w

    # 将清理后的结果添加到F115W结果列表中
    results_clean_f115w.append(result_clean_f115w)

    # 创建F200W的掩码，筛选有效的x_fit, y_fit和flux_fit
    mask_f200w = ((results_f200w[i]['x_fit'] > 0) & (results_f200w[i]['x_fit'] < 2048) &
                  (results_f200w[i]['y_fit'] > 0) & (results_f200w[i]['y_fit'] < 2048) &
                  (results_f200w[i]['flux_fit'] > 0))

    # 根据掩码过滤F200W结果
    result_clean_f200w = results_f200w[i][mask_f200w]

    # 使用WCS转换获取F200W的RA和Dec
    ra_f200w, dec_f200w = images_f200w[i].meta.wcs(result_clean_f200w['x_fit'], result_clean_f200w['y_fit'])

    # 创建SkyCoord对象以存储RA和Dec
    radec_f200w = SkyCoord(ra_f200w, dec_f200w, unit='deg')

    # 将RA和Dec添加到清理结果中
    result_clean_f200w['radec'] = radec_f200w

    # 将清理后的结果添加到F200W结果列表中
    results_clean_f200w.append(result_clean_f200w)

max_sep = 0.015 * u.arcsec  # 定义最大分离角度为0.015角秒

matches_phot_single = []  # 初始化匹配结果列表

filt1 = 'F115W'  # 定义第一个滤光片
filt2 = 'F200W'  # 定义第二个滤光片

# 遍历两个结果集，分别为F115W和F200W
for res1, res2 in zip(results_clean_f115w, results_clean_f200w):

    idx, d2d, _ = match_coordinates_sky(res1['radec'], res2['radec'])  # 匹配两个结果集的坐标

    sep_constraint = d2d < max_sep  # 根据最大分离角度筛选匹配结果

    match_phot_single = Table()  # 初始化单个匹配结果表

    # 从F115W结果中提取符合条件的参数
    x_0_f115w = res1['x_0'][sep_constraint]  # 提取x_0坐标
    y_0_f115w = res1['y_0'][sep_constraint]  # 提取y_0坐标
    x_fit_f115w = res1['x_fit'][sep_constraint]  # 提取拟合x坐标
    y_fit_f115w = res1['y_fit'][sep_constraint]  # 提取拟合y坐标
    radec_f115w = res1['radec'][sep_constraint]  # 提取天球坐标
    mag_f115w = (-2.5 * np.log10(res1['flux_fit']))[sep_constraint]  # 计算F115W的星等
    emag_f115w = (1.086 * (res1['flux_unc'] / res1['flux_fit']))[sep_constraint]  # 计算F115W的误差

    # 从F200W结果中提取符合条件的参数
    x_0_f200w = res2['x_0'][idx[sep_constraint]]  # 提取x_0坐标
    y_0_f200w = res2['y_0'][idx[sep_constraint]]  # 提取y_0坐标
    x_fit_f200w = res2['x_fit'][idx[sep_constraint]]  # 提取拟合x坐标
    y_fit_f200w = res2['y_fit'][idx[sep_constraint]]  # 提取拟合y坐标
    radec_f200w = res2['radec'][idx][sep_constraint]  # 提取天球坐标
    mag_f200w = (-2.5 * np.log10(res2['flux_fit']))[idx[sep_constraint]]  # 计算F200W的星等
    emag_f200w = (1.086 * (res2['flux_unc'] / res2['flux_fit']))[idx[sep_constraint]]  # 计算F200W的误差

    # 将F115W的参数添加到匹配结果表中
    match_phot_single['x_0_' + filt1] = x_0_f115w
    match_phot_single['y_0_' + filt1] = y_0_f115w
    match_phot_single['x_fit_' + filt1] = x_fit_f115w
    match_phot_single['y_fit_' + filt1] = y_fit_f115w
    match_phot_single['radec_' + filt1] = radec_f115w
    match_phot_single['mag_' + filt1] = mag_f115w
    match_phot_single['emag_' + filt1] = emag_f115w

    # 将F200W的参数添加到匹配结果表中
    match_phot_single['x_0_' + filt2] = x_0_f200w
    match_phot_single['y_0_' + filt2] = y_0_f200w
    match_phot_single['x_fit_' + filt2] = x_fit_f200w
    match_phot_single['y_fit_' + filt2] = y_fit_f200w
    match_phot_single['radec_' + filt2] = radec_f200w
    match_phot_single['mag_' + filt2] = mag_f200w
    match_phot_single['emag_' + filt2] = emag_f200w

    matches_phot_single.append(match_phot_single)  # 将当前匹配结果添加到总列表中

四幅图像的颜色-星等图（仪器星等）

plt.figure(figsize=(12, 16))  # 创建一个12x16英寸的图形

plt.clf()  # 清空当前图形

for i in np.arange(0, len(matches_phot_single), 1):  # 遍历所有匹配的单个光度数据

    ax = plt.subplot(2, 2, i + 1)  # 创建一个2x2的子图，选择第i+1个子图

    j = str(i + 1)  # 将索引i转换为字符串形式

    xlim0 = -0.5  # 设置x轴下限
    xlim1 = 0.8   # 设置x轴上限
    ylim0 = -1    # 设置y轴下限
    ylim1 = -9    # 设置y轴上限

    ax.set_xlim(xlim0, xlim1)  # 应用x轴限制
    ax.set_ylim(ylim0, ylim1)  # 应用y轴限制

    ax.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
    ax.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
    ax.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
    ax.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

    f115w_single = matches_phot_single[i]['mag_' + filt1]  # 获取filt1滤光片的光度数据
    f200w_single = matches_phot_single[i]['mag_' + filt2]  # 获取filt2滤光片的光度数据

    ax.scatter(f115w_single - f200w_single, f115w_single, s=1, color='k')  # 绘制散点图，x轴为f115w与f200w的差值，y轴为f115w的光度

    ax.set_xlabel(filt1 + '-' + filt2, fontdict=font2)  # 设置x轴标签
    ax.set_ylabel(filt1, fontdict=font2)  # 设置y轴标签
    ax.text(xlim0 + 0.1, -8.65, f"Image {j}", fontdict=font2)  # 在图中添加文本，标识图像序号
    

plt.tight_layout()  # 调整子图参数，使之填充整个图像区域

daofind与PSF算法获取位置的差异（以像素为单位）

我们展示了通过daofind和PSF拟合算法获得的星体位置之间的差异。我们还展示了差异 \(\Delta\)X 和 \(\Delta\)Y 与仪器星等（instrumental magnitudes）之间的关系。

plt.figure(figsize=(12, 6))  # 创建一个12x6英寸的图形

ax1 = plt.subplot(1, 2, 1)  # 创建1行2列的子图，选择第一个子图

xlim0 = -1  # x轴的最小值
xlim1 = 1   # x轴的最大值
ylim0 = -1  # y轴的最小值
ylim1 = 1   # y轴的最大值

ax1.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax1.set_ylim(ylim0, ylim1)  # 设置y轴范围

ax1.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax1.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

x_find_f115w = results_clean_f115w[0]['x_0']  # 获取f115w结果中的x坐标
y_find_f115w = results_clean_f115w[0]['y_0']  # 获取f115w结果中的y坐标

x_psf_f115w = results_clean_f115w[0]['x_fit']  # 获取f115w结果中的拟合x坐标
y_psf_f115w = results_clean_f115w[0]['y_fit']  # 获取f115w结果中的拟合y坐标

delta_x_f115w = x_find_f115w - x_psf_f115w  # 计算x坐标的偏差
delta_y_f115w = y_find_f115w - y_psf_f115w  # 计算y坐标的偏差

_, d_x_f115w, sigma_d_x_f115w = sigma_clipped_stats(delta_x_f115w)  # 计算x偏差的统计量
_, d_y_f115w, sigma_d_y_f115w = sigma_clipped_stats(delta_y_f115w)  # 计算y偏差的统计量

ax1.scatter(delta_x_f115w, delta_y_f115w, s=1, color='gray')  # 在子图中绘制散点图

ax1.set_xlabel(r'$\Delta$ X (px)', fontdict=font2)  # 设置x轴标签
ax1.set_ylabel(r'$\Delta$ Y (px)', fontdict=font2)  # 设置y轴标签
ax1.set_title(filt1, fontdict=font2)  # 设置子图标题

# 在子图中添加文本，显示x和y的偏差及其标准差
ax1.text(xlim0 + 0.05, ylim1 - 0.15, rf'$\Delta$ X = {d_x_f115w:5.3f} $\pm$ {sigma_d_x_f115w:5.3f}', 
         color='k', fontdict=font2)
ax1.text(xlim0 + 0.05, ylim1 - 0.30, rf'$\Delta$ Y = {d_y_f115w:5.3f} $\pm$ {sigma_d_y_f115w:5.3f}', 
         color='k', fontdict=font2)

ax1.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')  # 绘制x=0的虚线
ax1.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')  # 绘制y=0的虚线

ax2 = plt.subplot(1, 2, 2)  # 创建1行2列的子图，选择第二个子图

ax2.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax2.set_ylim(ylim0, ylim1)  # 设置y轴范围

ax2.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax2.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

x_find_f200w = results_clean_f200w[0]['x_0']  # 获取f200w结果中的x坐标
y_find_f200w = results_clean_f200w[0]['y_0']  # 获取f200w结果中的y坐标

x_psf_f200w = results_clean_f200w[0]['x_fit']  # 获取f200w结果中的拟合x坐标
y_psf_f200w = results_clean_f200w[0]['y_fit']  # 获取f200w结果中的拟合y坐标

delta_x_f200w = x_find_f200w - x_psf_f200w  # 计算x坐标的偏差
delta_y_f200w = y_find_f200w - y_psf_f200w  # 计算y坐标的偏差

_, d_x_f200w, sigma_d_x_f200w = sigma_clipped_stats(delta_x_f200w)  # 计算x偏差的统计量
_, d_y_f200w, sigma_d_y_f200w = sigma_clipped_stats(delta_y_f200w)  # 计算y偏差的统计量

ax2.scatter(delta_x_f200w, delta_y_f200w, s=1, color='gray')  # 在子图中绘制散点图

# 在子图中添加文本，显示x和y的偏差及其标准差
ax2.text(xlim0 + 0.05, ylim1 - 0.15, rf'$\Delta$ X = {d_x_f200w:5.3f} $\pm$ {sigma_d_x_f200w:5.3f}', 
         color='k', fontdict=font2)
ax2.text(xlim0 + 0.05, ylim1 - 0.30, rf'$\Delta$ Y = {d_y_f200w:5.3f} $\pm$ {sigma_d_y_f200w:5.3f}', 
         color='k', fontdict=font2)

ax2.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')  # 绘制x=0的虚线
ax2.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')  # 绘制y=0的虚线

ax2.set_xlabel(r'$\Delta$ X (px)', fontdict=font2)  # 设置x轴标签
ax2.set_ylabel(r'$\Delta$ Y (px)', fontdict=font2)  # 设置y轴标签
ax2.set_title(filt2, fontdict=font2)  # 设置子图标题

plt.tight_layout()  # 调整子图布局以避免重叠

plt.figure(figsize=(12, 8))  # 创建一个12x8英寸的图形

ax1 = plt.subplot(2, 2, 1)  # 创建一个2x2的子图，选择第一个子图

xlim0 = -9  # 设置x轴的最小值
xlim1 = -1  # 设置x轴的最大值
ylim0 = -1  # 设置y轴的最小值
ylim1 = 1   # 设置y轴的最大值

ax1.set_xlim(xlim0, xlim1)  # 设置第一个子图的x轴范围
ax1.set_ylim(ylim0, ylim1)  # 设置第一个子图的y轴范围

ax1.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax1.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

mag_inst_f115w = -2.5 * np.log10(results_clean_f115w[0]['flux_fit'])  # 计算F115W滤光片的瞬时星等

ax1.scatter(mag_inst_f115w, delta_x_f115w, s=1, color='gray')  # 在第一个子图中绘制散点图
ax1.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')  # 绘制y=0的虚线

ax1.set_xlabel(filt1 + '_inst', fontdict=font2)  # 设置x轴标签
ax1.set_ylabel(r'$\Delta$ X (px)', fontdict=font2)  # 设置y轴标签

ax2 = plt.subplot(2, 2, 2)  # 创建第二个子图

ax2.set_xlim(xlim0, xlim1)  # 设置第二个子图的x轴范围
ax2.set_ylim(ylim0, ylim1)  # 设置第二个子图的y轴范围

ax2.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax2.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

ax2.scatter(mag_inst_f115w, delta_y_f115w, s=1, color='gray')  # 在第二个子图中绘制散点图
ax2.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')  # 绘制y=0的虚线

ax2.set_xlabel(filt1 + '_inst', fontdict=font2)  # 设置x轴标签
ax2.set_ylabel(r'$\Delta$ Y (px)', fontdict=font2)  # 设置y轴标签

ax3 = plt.subplot(2, 2, 3)  # 创建第三个子图

ax3.set_xlim(xlim0, xlim1)  # 设置第三个子图的x轴范围
ax3.set_ylim(ylim0, ylim1)  # 设置第三个子图的y轴范围

ax3.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax3.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax3.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax3.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

mag_inst_f200w = -2.5 * np.log10(results_clean_f200w[0]['flux_fit'])  # 计算F200W滤光片的瞬时星等

ax3.scatter(mag_inst_f200w, delta_x_f200w, s=1, color='gray')  # 在第三个子图中绘制散点图
ax3.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')  # 绘制y=0的虚线

ax3.set_xlabel(filt2 + '_inst', fontdict=font2)  # 设置x轴标签
ax3.set_ylabel(r'$\Delta$ X (px)', fontdict=font2)  # 设置y轴标签

ax4 = plt.subplot(2, 2, 4)  # 创建第四个子图

ax4.set_xlim(xlim0, xlim1)  # 设置第四个子图的x轴范围
ax4.set_ylim(ylim0, ylim1)  # 设置第四个子图的y轴范围

ax4.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax4.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax4.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax4.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

ax4.scatter(mag_inst_f200w, delta_y_f200w, s=1, color='gray')  # 在第四个子图中绘制散点图
ax4.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')  # 绘制y=0的虚线

ax4.set_xlabel(filt2 + '_inst', fontdict=font2)  # 设置x轴标签
ax4.set_ylabel(r'$\Delta$ Y (px)', fontdict=font2)  # 设置y轴标签

plt.tight_layout()  # 自动调整子图参数，使之填充整个图像区域

对每个滤光片交叉匹配4个目录

为了获得最终的颜色-亮度图，我们需要对每个滤光片的所有目录进行交叉匹配，然后再交叉匹配得到的最终目录。

注意：这是最保守的方法，因为我们要求一颗星星必须在所有4个目录中都能找到。

开发者注意事项：

我找不到更简单的方法来编写这个函数，你需要匹配前两个目录，得出一个仅包含匹配项的子目录，然后对所有其他可用目录进行迭代。我们还应该考虑如何创建一个函数，以便在X个目录中保留星星，即使它们在Y个目录中可用（例如，如果由于某种原因，某个测量在1个图像中不可用，但在其他3个图像中测量良好，那么丢弃该对象是没有意义的）。

def crossmatch_filter(table=None):
    # 初始化计数器
    num = 0

    # 创建一个字符串数组，表示每个分类的编号
    num_cat = np.char.mod('%d', np.arange(1, len(table) + 1, 1))

    # 匹配第一个和第二个表的坐标
    idx_12, d2d_12, _ = match_coordinates_sky(table[num]['radec'], table[num + 1]['radec'])

    # 设置分离约束条件
    sep_constraint_12 = d2d_12 < max_sep

    # 创建一个表来存储匹配结果
    matches_12 = Table()

    # 存储第一个表的匹配结果
    matches_12['radec_' + num_cat[num]] = table[num]['radec'][sep_constraint_12]
    matches_12['mag_' + num_cat[num]] = (-2.5 * np.log10(table[num]['flux_fit']))[sep_constraint_12]
    matches_12['emag_' + num_cat[num]] = (1.086 * (table[num]['flux_unc'] / 
                                                   table[num]['flux_fit']))[sep_constraint_12]

    # 存储第二个表的匹配结果
    matches_12['radec_' + num_cat[num + 1]] = table[num + 1]['radec'][idx_12[sep_constraint_12]]
    matches_12['mag_' + num_cat[num + 1]] = (-2.5 * np.log10(table[num + 1]['flux_fit']))[idx_12[sep_constraint_12]]
    matches_12['emag_' + num_cat[num + 1]] = (1.086 * (table[num + 1]['flux_unc'] /
                                                       table[num + 1]['flux_fit']))[idx_12[sep_constraint_12]]

    # 匹配结果与第三个表的坐标
    idx_123, d2d_123, _ = match_coordinates_sky(matches_12['radec_' + num_cat[num]], table[num + 2]['radec'])

    # 设置分离约束条件
    sep_constraint_123 = d2d_123 < max_sep

    # 创建一个表来存储匹配结果
    matches_123 = Table()

    # 存储第一个和第二个表的匹配结果
    matches_123['radec_' + num_cat[num]] = matches_12['radec_' + num_cat[num]][sep_constraint_123]
    matches_123['mag_' + num_cat[num]] = matches_12['mag_' + num_cat[num]][sep_constraint_123]
    matches_123['emag_' + num_cat[num]] = matches_12['emag_' + num_cat[num]][sep_constraint_123]
    matches_123['radec_' + num_cat[num + 1]] = matches_12['radec_' + num_cat[num + 1]][sep_constraint_123]
    matches_123['mag_' + num_cat[num + 1]] = matches_12['mag_' + num_cat[num + 1]][sep_constraint_123]
    matches_123['emag_' + num_cat[num + 1]] = matches_12['emag_' + num_cat[num + 1]][sep_constraint_123]
    matches_123['radec_' + num_cat[num + 2]] = table[num + 2]['radec'][idx_123[sep_constraint_123]]
    matches_123['mag_' + num_cat[num + 2]] = (-2.5 * np.log10(table[num + 2]['flux_fit']))[idx_123[sep_constraint_123]]
    matches_123['emag_' + num_cat[num + 2]] = (1.086 * (table[num + 2]['flux_unc'] /
                                                        table[num + 2]['flux_fit']))[idx_123[sep_constraint_123]]

    # 匹配结果与第四个表的坐标
    idx_1234, d2d_1234, _ = match_coordinates_sky(matches_123['radec_' + num_cat[num]], table[num + 3]['radec'])

    # 设置分离约束条件
    sep_constraint_1234 = d2d_1234 < max_sep

    # 创建一个表来存储匹配结果
    matches_1234 = Table()

    # 存储第一个、第二个和第三个表的匹配结果
    matches_1234['radec_' + num_cat[num]] = matches_123['radec_' + num_cat[num]][sep_constraint_1234]
    matches_1234['mag_' + num_cat[num]] = matches_123['mag_' + num_cat[num]][sep_constraint_1234]
    matches_1234['emag_' + num_cat[num]] = matches_123['emag_' + num_cat[num]][sep_constraint_1234]
    matches_1234['radec_' + num_cat[num + 1]] = matches_123['radec_' + num_cat[num + 1]][sep_constraint_1234]
    matches_1234['mag_' + num_cat[num + 1]] = matches_123['mag_' + num_cat[num + 1]][sep_constraint_1234]
    matches_1234['emag_' + num_cat[num + 1]] = matches_123['emag_' + num_cat[num + 1]][sep_constraint_1234]
    matches_1234['radec_' + num_cat[num + 2]] = matches_123['radec_' + num_cat[num + 2]][sep_constraint_1234]
    matches_1234['mag_' + num_cat[num + 2]] = matches_123['mag_' + num_cat[num + 2]][sep_constraint_1234]
    matches_1234['emag_' + num_cat[num + 2]] = matches_123['emag_' + num_cat[num + 2]][sep_constraint_1234]
    matches_1234['radec_' + num_cat[num + 3]] = table[num + 3]['radec'][idx_1234[sep_constraint_1234]]
    matches_1234['mag_' + num_cat[num + 3]] = (-2.5 * np.log10(table[num + 3]['flux_fit']))[idx_1234[sep_constraint_1234]]
    matches_1234['emag_' + num_cat[num + 3]] = (1.086 * (table[num + 3]['flux_unc'] /
                                                         table[num + 3]['flux_fit']))[idx_1234[sep_constraint_1234]]

    # 返回最终的匹配结果
    return matches_1234

# 使用 crossmatch_filter 函数对 F115W 结果进行交叉匹配
matches_f115w = crossmatch_filter(table=results_clean_f115w)

# 使用 crossmatch_filter 函数对 F200W 结果进行交叉匹配
matches_f200w = crossmatch_filter(table=results_clean_f200w)

对于最终目录，我们假设亮度（magnitude）是4个测量值的平均值，而亮度的误差是其标准差。

为了方便在表格上执行这个算术操作，我们将表格转换为pandas数据框（dataframe）。

# 将 matches_f115w 转换为 pandas DataFrame
df_f115w = matches_f115w.to_pandas()

# 将 matches_f200w 转换为 pandas DataFrame
df_f200w = matches_f200w.to_pandas()

# 计算 RA 坐标的均值并存储在新列中
df_f115w['RA_' + filt1] = df_f115w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].mean(axis=1)

# 计算 RA 坐标的标准差并存储在新列中
df_f115w['e_RA_' + filt1] = df_f115w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].std(axis=1)

# 计算 Dec 坐标的均值并存储在新列中
df_f115w['Dec_' + filt1] = df_f115w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].mean(axis=1)

# 计算 Dec 坐标的标准差并存储在新列中
df_f115w['e_Dec_' + filt1] = df_f115w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].std(axis=1)

# 计算滤光片 f115w 的仪器测量值均值并存储在新列中
df_f115w[filt1 + '_inst'] = df_f115w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].mean(axis=1)

# 计算滤光片 f115w 的仪器测量值标准差并存储在新列中
df_f115w['e' + filt1 + '_inst'] = df_f115w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].std(axis=1)

# 计算 RA 坐标的均值并存储在新列中
df_f200w['RA_' + filt2] = df_f200w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].mean(axis=1)

# 计算 RA 坐标的标准差并存储在新列中
df_f200w['e_RA_' + filt2] = df_f200w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].std(axis=1)

# 计算 Dec 坐标的均值并存储在新列中
df_f200w['Dec_' + filt2] = df_f200w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].mean(axis=1)

# 计算 Dec 坐标的标准差并存储在新列中
df_f200w['e_Dec_' + filt2] = df_f200w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].std(axis=1)

# 计算滤光片 f200w 的仪器测量值均值并存储在新列中
df_f200w[filt2 + '_inst'] = df_f200w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].mean(axis=1)

# 计算滤光片 f200w 的仪器测量值标准差并存储在新列中
df_f200w['e' + filt2 + '_inst'] = df_f200w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].std(axis=1)

最终颜色-星等图（仪器星等）

plt.figure(figsize=(12, 14))  # 创建一个12x14英寸的图形

plt.clf()  # 清空当前图形

ax1 = plt.subplot(1, 2, 1)  # 创建一个1行2列的子图，选择第一个子图

ax1.set_xlabel(filt1 + '_inst -' + filt2 + '_inst', fontdict=font2)  # 设置x轴标签
ax1.set_ylabel(filt1 + '_inst', fontdict=font2)  # 设置y轴标签

xlim0 = -0.5  # x轴下限
xlim1 = 0.8   # x轴上限
ylim0 = -1.5  # y轴下限
ylim1 = -9    # y轴上限

ax1.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax1.set_ylim(ylim0, ylim1)  # 设置y轴范围

ax1.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax1.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

# 创建SkyCoord对象，表示F115W和F200W的天球坐标
radec_f115w_inst = SkyCoord(df_f115w['RA_' + filt1], df_f115w['Dec_' + filt1], unit='deg')
radec_f200w_inst = SkyCoord(df_f200w['RA_' + filt2], df_f200w['Dec_' + filt2], unit='deg')

# 匹配F115W和F200W的坐标
idx_inst, d2d_inst, _ = match_coordinates_sky(radec_f115w_inst, radec_f200w_inst)

sep_constraint_inst = d2d_inst < max_sep  # 根据最大分离限制筛选匹配

# 提取满足条件的F115W数据
f115w_inst = np.array(df_f115w[filt1 + '_inst'][sep_constraint_inst])
ef115w_inst = np.array(df_f115w['e' + filt1 + '_inst'][sep_constraint_inst])
radec_f115w = radec_f115w_inst[sep_constraint_inst]

# 提取满足条件的F200W数据
f200w_inst = np.array(df_f200w[filt2 + '_inst'][idx_inst[sep_constraint_inst]])
ef200w_inst = np.array(df_f200w['e' + filt2 + '_inst'][idx_inst[sep_constraint_inst]])
radec_f200w = radec_f200w_inst[idx_inst[sep_constraint_inst]]

ax1.scatter(f115w_inst - f200w_inst, f115w_inst, s=1, color='k')  # 在ax1上绘制散点图

ax2 = plt.subplot(2, 2, 2)  # 创建第二个子图

ax2.set_xlabel(filt1 + '_inst', fontdict=font2)  # 设置x轴标签
ax2.set_ylabel(r'$\sigma$' + filt1, fontdict=font2)  # 设置y轴标签

xlim0 = -9  # x轴下限
xlim1 = -1.5  # x轴上限
ylim0 = -0.01  # y轴下限
ylim1 = 1  # y轴上限

ax2.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax2.set_ylim(ylim0, ylim1)  # 设置y轴范围

ax2.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax2.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

ax2.scatter(df_f115w[filt1 + '_inst'], df_f115w['e' + filt1 + '_inst'], s=1, color='k')  # 在ax2上绘制散点图

ax3 = plt.subplot(2, 2, 4)  # 创建第四个子图

ax3.set_xlabel(filt2 + '_inst', fontdict=font2)  # 设置x轴标签
ax3.set_ylabel(r'$\sigma$' + filt2, fontdict=font2)  # 设置y轴标签

ax3.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax3.set_ylim(ylim0, ylim1)  # 设置y轴范围

ax3.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax3.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax3.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax3.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

ax3.scatter(df_f200w[filt2 + '_inst'], df_f200w['e' + filt2 + '_inst'], s=1, color='k')  # 在ax3上绘制散点图

plt.tight_layout()  # 调整子图布局以避免重叠

光度零点 (Photometric Zeropoints)

为了获得最终校准的颜色-星等图，我们需要计算光度零点。因此，我们需要对校准后的图像（Level-3）进行光圈光度测量 (aperture photometry)，应用适当的光圈修正（上面字典中提供的值是针对无限光圈的），然后将其与点扩散函数 (PSF) 光度测量进行比较。因此，我们可以将步骤总结如下：

• 进行光圈光度测量

• 应用适当的光圈修正

• 应用表格化的零点

• 与PSF光度测量进行交叉匹配

加载已校准和矫正的图像（三级成像处理管道）

# 初始化一个字典，用于存储不同探测器的图像数据
dict_images_combined = {'NRCA1': {}, 'NRCA2': {}, 'NRCA3': {}, 'NRCA4': {}, 'NRCA5': {},
                        'NRCB1': {}, 'NRCB2': {}, 'NRCB3': {}, 'NRCB4': {}, 'NRCB5': {}}

# 初始化短波和长波滤光器的字典
dict_filter_short = {}
dict_filter_long = {}

# 初始化短波和长波的滤光器、探测器列表
ff_short = []  # 短波滤光器列表
det_short = []  # 短波探测器列表
det_long = []  # 长波探测器列表
ff_long = []  # 长波滤光器列表
detlist_short = []  # 短波探测器唯一列表
detlist_long = []  # 长波探测器唯一列表
filtlist_short = []  # 短波滤光器唯一列表
filtlist_long = []  # 长波滤光器唯一列表

# 检查是否存在合并后的fits文件
if not glob.glob('./*combined*fits'):
    print("Downloading images")  # 打印下载图像的提示

    # 定义图像数据的下载链接和本地文件名
    boxlink_images_lev3 = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/images_level3.tar.gz'
    boxfile_images_lev3 = './images_level3.tar.gz'

    # 下载图像数据
    request.urlretrieve(boxlink_images_lev3, boxfile_images_lev3)

    # 解压下载的tar文件
    tar = tarfile.open(boxfile_images_lev3, 'r')
    tar.extractall(filter='data')

    # 设置图像目录
    images_dir = './'
    # 获取合并后的fits文件列表
    files_singles = sorted(glob.glob(os.path.join(images_dir, "*combined*fits")))
else:
    # 如果文件已存在，直接获取合并后的fits文件列表
    images_dir = './'
    files_singles = sorted(glob.glob(os.path.join(images_dir, "*combined*fits")))

# 遍历每个合并后的fits文件
for file in files_singles:
    im = fits.open(file)  # 打开fits文件
    f = im[0].header['FILTER']  # 获取滤光器信息
    d = im[0].header['DETECTOR']  # 获取探测器信息

    # 根据探测器类型进行调整
    if d == 'NRCBLONG':
        d = 'NRCB5'  # 将NRCBLONG映射为NRCB5
    elif d == 'NRCALONG':
        d = 'NRCA5'  # 将NRCALONG映射为NRCA5
    else:
        d = d  # 保持原探测器名称

    wv = float(f[1:3])  # 提取波长信息

    # 根据波长将滤光器分类
    if wv > 24:
        ff_long.append(f)  # 添加到长波滤光器列表
        det_long.append(d)  # 添加到长波探测器列表
    else:
        ff_short.append(f)  # 添加到短波滤光器列表
        det_short.append(d)  # 添加到短波探测器列表

    # 获取短波和长波探测器的唯一列表
    detlist_short = sorted(list(dict.fromkeys(det_short)))
    detlist_long = sorted(list(dict.fromkeys(det_long)))

    unique_list_filters_short = []  # 短波滤光器唯一列表
    unique_list_filters_long = []  # 长波滤光器唯一列表

    # 更新短波滤光器字典
    for x in ff_short:
        if x not in unique_list_filters_short:
            dict_filter_short.setdefault(x, {})  # 如果滤光器不在列表中，则添加

    # 更新长波滤光器字典
    for x in ff_long:
        if x not in unique_list_filters_long:
            dict_filter_long.setdefault(x, {})  # 如果滤光器不在列表中，则添加

    # 将短波探测器与滤光器字典关联
    for d_s in detlist_short:
        dict_images_combined[d_s] = dict_filter_short

    # 将长波探测器与滤光器字典关联
    for d_l in detlist_long:
        dict_images_combined[d_l] = dict_filter_long

    # 获取短波和长波滤光器的唯一列表
    filtlist_short = sorted(list(dict.fromkeys(dict_filter_short)))
    filtlist_long = sorted(list(dict.fromkeys(dict_filter_long)))

    # 将图像文件添加到对应的探测器和滤光器字典中
    if len(dict_images_combined[d][f]) == 0:
        dict_images_combined[d][f] = {'images': [file]}  # 如果没有图像，初始化列表
    else:
        dict_images_combined[d][f]['images'].append(file)  # 否则，添加图像文件

# 打印可用的探测器和滤光器信息
print("Available Detectors for SW channel:", detlist_short)  # 打印短波探测器
print("Available Detectors for LW channel:", detlist_long)  # 打印长波探测器
print("Available SW Filters:", filtlist_short)  # 打印短波滤光器
print("Available LW Filters:", filtlist_long)  # 打印长波滤光器

显示图像

plt.figure(figsize=(14, 14))  # 创建一个14x14英寸的图形

for det in dets_short:  # 遍历短列表中的每个探测器

    for i, filt in enumerate(filts_short):  # 遍历短列表中的每个滤光片，并获取索引

        image = fits.open(dict_images_combined[det][filt]['images'][0])  # 打开对应探测器和滤光片的图像文件

        data_sb = image[1].data  # 获取图像数据

        ax = plt.subplot(1, len(filts_short), i + 1)  # 创建子图，行数为1，列数为滤光片数量，当前索引为i+1

        norm = simple_norm(data_sb, 'sqrt', percent=99.)  # 使用平方根归一化数据，99%分位数

        plt.xlabel("X [px]", fontdict=font2)  # 设置x轴标签

        plt.ylabel("Y [px]", fontdict=font2)  # 设置y轴标签

        plt.title(filt, fontdict=font2)  # 设置子图标题为当前滤光片名称

        ax.imshow(data_sb, norm=norm, cmap='Greys')  # 显示图像数据，使用灰度色图和归一化

plt.tight_layout()  # 调整子图布局以避免重叠

光圈光度法 (Aperture Photometry)

如我们之前所做的，我们创建一个字典，其中包含每个图像的衍生光圈光度（aperture photometry）表格。

dict_aper = {}  # 初始化一个空字典，用于存储光圈数据

for det in dets_short:  # 遍历短列表中的每个探测器

    dict_aper.setdefault(det, {})  # 如果字典中没有该探测器，则初始化一个空字典

    for j, filt in enumerate(filts_short):  # 遍历短列表中的每个滤光片，并获取其索引

        dict_aper[det].setdefault(filt, {})  # 如果该探测器下没有该滤光片，则初始化一个空字典

        dict_aper[det][filt]['stars for ap phot'] = None  # 初始化用于光圈光度测量的星星数据为None

        dict_aper[det][filt]['stars for ap phot matched'] = None  # 初始化匹配的光圈光度测量星星数据为None

        dict_aper[det][filt]['aperture phot table'] = None  # 初始化光圈光度表为None

寻找明亮的孤立恒星

def find_bright_stars(det='NRCA1', filt='F070W', dist_sel=False):
    # 初始化背景噪声计算方法
    bkgrms = MADStdBackgroundRMS()
    # 初始化背景均值计算方法
    mmm_bkg = MMMBackground()

    # 打开指定探测器和滤光片的图像文件
    image = fits.open(dict_images_combined[det][filt]['images'][i])
    # 获取图像数据
    data_sb = image[1].data
    # 获取图像头信息
    imh = image[1].header

    # 打印当前处理的图像信息
    print(f"Selecting stars for aperture photometry on image {i + 1} of filter {filt}, detector {det}")

    # 将数据转换为DN/s单位
    data = data_sb / imh['PHOTMJSR']
    units = imh['BUNIT']
    # 打印单位转换因子
    print(f"Conversion factor from {units} to DN/s for filter {filt}: {imh['PHOTMJSR']}")

    # 获取滤光片的点扩散函数（PSF）全宽半最大值（FWHM）
    sigma_psf = dict_utils[filt]['psf fwhm']
    # 打印FWHM值
    print(f"FWHM for the filter {filt}: {sigma_psf} px")

    # 计算背景标准差
    std = bkgrms(data)
    # 计算背景值
    bkg = mmm_bkg(data)

    # 使用DAOStarFinder寻找星点
    daofind = DAOStarFinder(threshold=th[j] * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
                            sharplo=0.30, sharphi=1.40)
    # 获取经过校正的星点
    apcorr_stars = daofind(data)
    # 将找到的星点存储在字典中
    dict_aper[det][filt]['stars for ap phot'] = apcorr_stars

    if dist_sel:
        # 如果需要计算距离选择
        print("")
        print("Calculating closest neighbour distance")

        d = []  # 初始化距离列表

        # 使用DAOStarFinder寻找所有星点
        daofind_tot = DAOStarFinder(threshold=10 * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
                                    sharplo=0.30, sharphi=1.40)
        stars_tot = daofind_tot(data)

        # 获取所有星点的坐标
        x_tot = stars_tot['xcentroid']
        y_tot = stars_tot['ycentroid']

        # 计算每个经过校正的星点与所有其他星点的距离
        for xx, yy in zip(apcorr_stars['xcentroid'], apcorr_stars['ycentroid']):
            sep = []  # 初始化距离列表
            dist = np.sqrt((x_tot - xx)**2 + (y_tot - yy)**2)  # 计算距离
            sep = np.sort(dist)[1:2][0]  # 获取最近邻距离
            d.append(sep)  # 将距离添加到列表中

        # 将最小距离添加到校正星点数据中
        apcorr_stars['min distance'] = d
        mask_dist = (apcorr_stars['min distance'] > min_sep[j])  # 创建距离掩码

        # 应用掩码过滤星点
        apcorr_stars = apcorr_stars[mask_dist]

        # 更新字典中的星点数据
        dict_aper[det][filt]['stars for ap phot'] = apcorr_stars

        # 打印最小距离和找到的星点数量
        print("Minimum distance required:", min_sep[j], "px")
        print("")
        print(f"Number of bright isolated sources found in the image for {filt}: {len(apcorr_stars)}")
        print("-----------------------------------------------------")
        print("")
    else:
        # 如果不需要距离选择，直接打印找到的星点数量
        print("")
        print(f"Number of bright sources found in the image for {filt}: {len(apcorr_stars)}")
        print("--------------------------------------------")
        print("")

    return  # 返回函数结束

tic = time.perf_counter()  # 记录开始时间

th = [700, 500]  # 两个滤光器的阈值水平（长度必须与分析的滤光器数量匹配）

min_sep = [10, 10]  # zp 星星与最近邻星星之间可接受的最小距离

# 遍历短时间检测列表
for det in dets_short:

    # 遍历短时间滤光器列表
    for j, filt in enumerate(filts_short):

        # 遍历合并图像的索引
        for i in np.arange(0, len(dict_images_combined[det][filt]['images']), 1):

            # 查找明亮的星星，dist_sel参数设为False
            find_bright_stars(det=det, filt=filt, dist_sel=False)

toc = time.perf_counter()  # 记录结束时间

# 打印用于光度测量的星星查找所用的时间
print("Elapsed Time for finding stars for Aperture Photometry:", toc - tic)

为了进一步获得高质量的样本，我们对两个滤光器的目录进行交叉匹配，仅保留共同的恒星。

for det in dets_short:  # 遍历短时间序列中的每个探测器

    for j, filt in enumerate(filts_short):  # 遍历短时间序列中的每个滤光片

        for i in np.arange(0, len(dict_images_combined[det][filt]['images']), 1):  # 遍历当前探测器和滤光片的所有图像

            image = ImageModel(dict_images_combined[det][filt]['images'][i])  # 创建图像模型对象

            ra, dec = image.meta.wcs(dict_aper[det][filt]['stars for ap phot']['xcentroid'],  # 获取天体的RA和Dec
                                     dict_aper[det][filt]['stars for ap phot']['ycentroid'])

            radec = SkyCoord(ra, dec, unit='deg')  # 创建SkyCoord对象以存储RA和Dec坐标

            dict_aper[det][filt]['stars for ap phot']['radec'] = radec  # 将RA和Dec坐标存储到字典中

# 使用 match_coordinates_sky 函数匹配两个滤波器中的星体坐标
idx_ap, d2d_ap, _ = match_coordinates_sky(dict_aper[det][filt1]['stars for ap phot']['radec'],
                                          dict_aper[det][filt2]['stars for ap phot']['radec'])

# 设置分离约束条件，筛选出距离小于最大分离的匹配
sep_constraint_ap = d2d_ap < max_sep

# 创建两个空的表格用于存储匹配结果
matched_apcorr_f115w = Table()  # 用于存储滤波器 f115w 的匹配结果
matched_apcorr_f200w = Table()   # 用于存储滤波器 f200w 的匹配结果

# 根据分离约束条件从 f115w 中提取匹配的星体
matched_apcorr_f115w = dict_aper[det][filt1]['stars for ap phot'][sep_constraint_ap]

# 根据索引从 f200w 中提取匹配的星体
matched_apcorr_f200w = dict_aper[det][filt2]['stars for ap phot'][idx_ap[sep_constraint_ap]]

# 将匹配结果存储回字典中
dict_aper[det][filt1]['stars for ap phot matched'] = matched_apcorr_f115w
dict_aper[det][filt2]['stars for ap phot matched'] = matched_apcorr_f200w

加载孔径校正表

注意：这些值是通过对合成STPSF点扩散函数（PSFs）的研究获得的。一旦我们获得了飞行中的测量数据，这些值将会更新。

if os.path.isfile('./aperture_correction_table.txt'):  # 检查文件是否存在

    ap_tab = './aperture_correction_table.txt'  # 如果存在，设置文件路径

else:  # 如果文件不存在

    print("Downloading the aperture correction table")  # 打印下载提示

    boxlink_apcorr_table = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/aperture_correction_table.txt'  # 定义下载链接

    boxfile_apcorr_table = './aperture_correction_table.txt'  # 定义保存路径

    request.urlretrieve(boxlink_apcorr_table, boxfile_apcorr_table)  # 下载文件并保存

    ap_tab = './aperture_correction_table.txt'  # 设置文件路径

# 读取光圈校正表，指定无表头，分隔符为任意空白字符，第一列作为索引
aper_table = pd.read_csv(ap_tab, header=None, sep=r'\s+', index_col=0,
                         names=['filter', 'pupil', 'wave', 'r10', 'r20', 'r30', 'r40', 'r50', 'r60', 'r70', 'r80',
                                'r85', 'r90', 'sky_flux_px', 'apcorr10', 'apcorr20', 'apcorr30', 'apcorr40',
                                'apcorr50', 'apcorr60', 'apcorr70', 'apcorr80', 'apcorr85', 'apcorr90', 'sky_in',
                                'sky_out'], comment='#', skiprows=0, usecols=range(0, 26))

aper_table.head()  # 显示数据表的前几行

执行光圈光度法

def aperture_phot(det=det, filt='F070W'):
    # 定义一个函数，用于进行光度测量，输入参数为探测器和滤光片

    radii = [aper_table.loc[filt]['r70']]
    # 从光度表中获取对应滤光片的半径

    ees = '70'.split()
    # 定义等效半径的列表

    ee_radii = dict(zip(ees, radii))
    # 创建一个字典，将等效半径与其对应的值关联

    positions = np.transpose((dict_aper[det][filt]['stars for ap phot matched']['xcentroid'],
                              dict_aper[det][filt]['stars for ap phot matched']['ycentroid']))
    # 获取匹配的星体的中心坐标

    image = fits.open(dict_images_combined[det][filt]['images'][0])
    # 打开合成图像文件

    data_sb = image[1].data
    # 获取图像数据

    imh = image[1].header
    # 获取图像头信息

    data = data_sb / imh['PHOTMJSR']
    # 进行数据的光度校正

    # 从光度校正表中获取天空背景的参数：
    sky = {"sky_in": aper_table.loc[filt]['r80'], "sky_out": aper_table.loc[filt]['r85']}
    # 定义天空背景的内外半径

    tic = time.perf_counter()
    # 记录开始时间

    table_aper = Table()
    # 创建一个表格用于存储光度测量结果

    for ee, radius in ee_radii.items():
        # 遍历每个等效半径
        print(f"Performing aperture photometry for radius equivalent to EE = {ee}% for filter {filt}")
        # 打印当前处理的信息

        aperture = CircularAperture(positions, r=radius)
        # 创建一个圆形光圈用于光度测量

        annulus_aperture = CircularAnnulus(positions, r_in=sky["sky_in"], r_out=sky["sky_out"])
        # 创建一个圆环用于背景测量

        annulus_mask = annulus_aperture.to_mask(method='center')
        # 将圆环转换为掩膜

        bkg_median = []
        # 初始化背景中值列表

        for mask in annulus_mask:
            # 遍历每个掩膜
            annulus_data = mask.multiply(data)
            # 使用掩膜提取背景数据

            annulus_data_1d = annulus_data[mask.data > 0]
            # 将有效数据转换为一维数组

            _, median_sigclip, _ = sigma_clipped_stats(annulus_data_1d)
            # 计算背景数据的中值，使用sigma剪切方法

            bkg_median.append(median_sigclip)
            # 将中值添加到列表中

        bkg_median = np.array(bkg_median)
        # 将背景中值列表转换为NumPy数组

        phot = aperture_photometry(data, aperture, method='exact')
        # 进行光度测量

        phot['annulus_median'] = bkg_median
        # 将背景中值添加到光度测量结果中

        phot['aper_bkg'] = bkg_median * aperture.area
        # 计算光圈背景

        phot['aper_sum_bkgsub'] = phot['aperture_sum'] - phot['aper_bkg']
        # 计算背景校正后的光度总和

        apcorr = [aper_table.loc[filt]['apcorr70']]
        # 获取光度校正因子

        phot['aper_sum_corrected'] = phot['aper_sum_bkgsub'] * apcorr
        # 应用校正因子

        phot['mag_corrected'] = -2.5 * np.log10(phot['aper_sum_corrected']) + dict_utils[filt]['VegaMAG zp modB']
        # 计算校正后的星等

        table_aper.add_column(phot['aperture_sum'], name='aper_sum_' + ee)
        # 将光度总和添加到结果表中

        table_aper.add_column(phot['annulus_median'], name='annulus_median_' + ee)
        # 将背景中值添加到结果表中

        table_aper.add_column(phot['aper_bkg'], name='aper_bkg_ee_' + ee)
        # 将光圈背景添加到结果表中

        table_aper.add_column(phot['aper_sum_bkgsub'], name='aper_sum_bkgsub_' + ee)
        # 将背景校正后的光度总和添加到结果表中

        table_aper.add_column(phot['aper_sum_corrected'], name='aper_sum_corrected_' + filt) 
        # 将校正后的光度总和添加到结果表中

        table_aper.add_column(phot['mag_corrected'], name='mag_corrected_' + filt)
        # 将校正后的星等添加到结果表中

        dict_aper[det][filt]['aperture phot table'] = table_aper
        # 将结果表存储到字典中

    toc = time.perf_counter()
    # 记录结束时间

    print("Time Elapsed:", toc - tic)
    # 打印耗时

    return
    # 返回函数结束

# 使用指定的滤镜对检测数据进行光度测量
aperture_phot(det=det, filt=filt1)  # 对det数据应用filt1滤镜进行光度测量

aperture_phot(det=det, filt=filt2)  # 对det数据应用filt2滤镜进行光度测量

推导零点（Zeropoints）

plt.figure(figsize=(14, 8))  # 创建一个14x8英寸的图形

plt.clf()  # 清除当前图形中的所有内容

ax1 = plt.subplot(2, 1, 1)  # 创建一个2行1列的子图，选择第一个子图

ax1.set_xlabel(filt1, fontdict=font2)  # 设置x轴标签为filt1
ax1.set_ylabel('Zeropoint', fontdict=font2)  # 设置y轴标签为'Zeropoint'

# 匹配星体坐标，获取索引和距离
idx_zp_1, d2d_zp_1, _ = match_coordinates_sky(dict_aper[det][filt1]['stars for ap phot matched']['radec'], radec_f115w_inst)

sep_constraint_zp_1 = d2d_zp_1 < max_sep  # 设置距离约束条件

# 获取匹配的光度数据
f115w_ap_matched = np.array(dict_aper[det][filt1]['aperture phot table']['mag_corrected_' + filt1][sep_constraint_zp_1])
f115w_psf_matched = np.array(df_f115w[filt1 + '_inst'][idx_zp_1[sep_constraint_zp_1]])

# 计算差异并进行统计
diff_f115w = f115w_ap_matched - f115w_psf_matched
_, zp_f115w, zp_sigma_f115w = sigma_clipped_stats(diff_f115w)  # 计算零点和标准差

# 设置x和y轴的限制
xlim0 = -9
xlim1 = -5
ylim0 = np.mean(diff_f115w) - 0.5
ylim1 = np.mean(diff_f115w) + 0.5

ax1.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax1.set_ylim(ylim0, ylim1)  # 设置y轴范围

# 设置坐标轴的刻度
ax1.xaxis.set_major_locator(ticker.AutoLocator())  # 主刻度自动定位
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 次刻度自动定位
ax1.yaxis.set_major_locator(ticker.AutoLocator())  # 主刻度自动定位
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 次刻度自动定位

# 绘制散点图和零点线
ax1.scatter(f115w_psf_matched, diff_f115w, s=50, color='k')  # 绘制散点图
ax1.plot([xlim0, xlim1], [zp_f115w, zp_f115w], color='r', lw=5, ls='--')  # 绘制零点线
ax1.text(xlim0 + 0.05, ylim1 - 0.15, filt1 + rf' Zeropoint = {zp_f115w:5.3f} $\pm$ {zp_sigma_f115w:5.3f}', color='k', fontdict=font2)  # 添加文本注释

ax2 = plt.subplot(2, 1, 2)  # 创建第二个子图

ax2.set_xlabel(filt2, fontdict=font2)  # 设置x轴标签为filt2
ax2.set_ylabel('Zeropoint', fontdict=font2)  # 设置y轴标签为'Zeropoint'

# 匹配星体坐标，获取索引和距离
idx_zp_2, d2d_zp_2, _ = match_coordinates_sky(dict_aper[det][filt2]['stars for ap phot matched']['radec'], radec_f200w_inst)

sep_constraint_zp_2 = d2d_zp_2 < max_sep  # 设置距离约束条件

# 获取匹配的光度数据
f200w_ap_matched = np.array(dict_aper[det][filt2]['aperture phot table']['mag_corrected_' + filt2][sep_constraint_zp_2])
f200w_psf_matched = np.array(df_f200w[filt2 + '_inst'][idx_zp_2[sep_constraint_zp_2]])

# 计算差异并进行统计
diff_f200w = f200w_ap_matched - f200w_psf_matched
_, zp_f200w, zp_sigma_f200w = sigma_clipped_stats(diff_f200w)  # 计算零点和标准差

# 设置x和y轴的限制
xlim0 = -9
xlim1 = -5
ylim0 = np.mean(diff_f200w) - 0.5
ylim1 = np.mean(diff_f200w) + 0.5

ax2.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax2.set_ylim(ylim0, ylim1)  # 设置y轴范围

# 设置坐标轴的刻度
ax2.xaxis.set_major_locator(ticker.AutoLocator())  # 主刻度自动定位
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 次刻度自动定位
ax2.yaxis.set_major_locator(ticker.AutoLocator())  # 主刻度自动定位
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 次刻度自动定位

# 绘制散点图和零点线
ax2.scatter(f200w_psf_matched, diff_f200w, s=50, color='k')  # 绘制散点图
ax2.plot([xlim0, xlim1], [zp_f200w, zp_f200w], color='r', lw=5, ls='--')  # 绘制零点线
ax2.text(xlim0 + 0.05, ylim1 - 0.15, filt2 + rf' Zeropoint = {zp_f200w:5.3f} $\pm$ {zp_sigma_f200w:5.3f}', color='k', fontdict=font2)  # 添加文本注释

plt.tight_layout()  # 调整子图间距

导入输入光度数据

import os  # 导入os模块，用于文件和目录操作
import pandas as pd  # 导入pandas库，用于数据处理
from urllib import request  # 从urllib库导入request模块，用于下载文件

# 检查当前目录下是否存在'pointsource.cat'文件
if os.path.isfile('./pointsource.cat'):
    input_cat = './pointsource.cat'  # 如果文件存在，设置输入目录为该文件
else:
    print("Downloading input pointsource catalog")  # 如果文件不存在，打印下载提示

    # 设置要下载的文件链接
    boxlink_input_cat = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/pointsource.cat'
    boxfile_input_cat = './pointsource.cat'  # 设置下载后保存的文件名

    # 下载文件并保存到指定路径
    request.urlretrieve(boxlink_input_cat, boxfile_input_cat)
    input_cat = './pointsource.cat'  # 设置输入目录为下载的文件

# 读取输入的点源目录文件，指定列名和数据格式
cat = pd.read_csv(input_cat, header=None, sep=r'\s+', names=['ra_in', 'dec_in', 'f070w_in', 'f115w_in',
                                                             'f200w_in', 'f277w_in', 'f356w_in', 'f444w_in'],
                  comment='#', skiprows=7, usecols=range(0, 8))

cat.head()  # 显示数据框的前几行

从输入目录中提取与所分析区域相同的星星

lim_ra_min = np.min(radec_f115w.ra)  # 计算ra的最小值

lim_ra_max = np.max(radec_f115w.ra)  # 计算ra的最大值

lim_dec_min = np.min(radec_f115w.dec)  # 计算dec的最小值

lim_dec_max = np.max(radec_f115w.dec)  # 计算dec的最大值

# 选择在指定ra和dec范围内的天体数据
cat_sel = cat[(cat['ra_in'] > lim_ra_min) & (cat['ra_in'] < lim_ra_max) & (cat['dec_in'] > lim_dec_min)

              & (cat['dec_in'] < lim_dec_max)]

校准的颜色-星等图

plt.figure(figsize=(12, 14))  # 创建一个12x14英寸的图形

plt.clf()  # 清除当前图形中的所有内容

ax1 = plt.subplot(1, 2, 1)  # 创建一个1行2列的子图，选择第一个子图

mag1_in = np.array(cat_sel['f115w_in'])  # 从cat_sel中提取f115w_in数据并转换为numpy数组
mag2_in = np.array(cat_sel['f200w_in'])  # 从cat_sel中提取f200w_in数据并转换为numpy数组
diff_in = mag1_in - mag2_in  # 计算f115w和f200w之间的差值

xlim0 = -0.25  # 设置x轴的最小值
xlim1 = 1.75   # 设置x轴的最大值
ylim0 = 25     # 设置y轴的最小值
ylim1 = 15     # 设置y轴的最大值

ax1.set_xlim(xlim0, xlim1)  # 设置第一个子图的x轴范围
ax1.set_ylim(ylim0, ylim1)   # 设置第一个子图的y轴范围

ax1.xaxis.set_major_locator(ticker.AutoLocator())  # 自动选择x轴主刻度
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 自动选择x轴次刻度
ax1.yaxis.set_major_locator(ticker.AutoLocator())  # 自动选择y轴主刻度
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 自动选择y轴次刻度

ax1.scatter(mag1_in - mag2_in, mag1_in, s=1, color='k')  # 在第一个子图中绘制散点图

ax1.set_xlabel(filt1 + ' - ' + filt2, fontdict=font2)  # 设置x轴标签
ax1.set_ylabel(filt1, fontdict=font2)  # 设置y轴标签
ax1.text(xlim0 + 0.15, 15.5, "Input", fontdict=font2)  # 在图中添加文本“Input”

ax2 = plt.subplot(1, 2, 2)  # 创建一个1行2列的子图，选择第二个子图

ax2.set_xlim(xlim0, xlim1)  # 设置第二个子图的x轴范围
ax2.set_ylim(ylim0, ylim1)   # 设置第二个子图的y轴范围

ax2.xaxis.set_major_locator(ticker.AutoLocator())  # 自动选择x轴主刻度
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 自动选择x轴次刻度
ax2.yaxis.set_major_locator(ticker.AutoLocator())  # 自动选择y轴主刻度
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 自动选择y轴次刻度

f115w = f115w_inst + zp_f115w  # 计算f115w的实际值
f200w = f200w_inst + zp_f200w  # 计算f200w的实际值

maglim = np.arange(18, 25, 1)  # 创建一个从18到25的范围，步长为1
mags = []  # 初始化一个空列表用于存储星等
errs_mag = []  # 初始化一个空列表用于存储星等误差
errs_col = []  # 初始化一个空列表用于存储颜色误差

for i in np.arange(0, len(maglim) - 1, 1):  # 遍历maglim的每个区间

    mag = (maglim[i] + maglim[i + 1]) / 2  # 计算当前区间的中间值
    err_mag1 = ef115w_inst[(f115w > maglim[i]) & (f115w < maglim[i + 1])]  # 计算f115w在当前区间的误差
    err_mag2 = ef200w_inst[(f115w > maglim[i]) & (f115w < maglim[i + 1])]  # 计算f200w在当前区间的误差
    err_mag = np.mean(err_mag1[i])  # 计算当前区间的平均误差
    err_temp = np.sqrt(err_mag1**2 + err_mag2**2)  # 计算误差的平方和的平方根
    err_col = np.mean(err_temp[i])  # 计算颜色误差的平均值

    errs_mag.append(err_mag)  # 将当前星等误差添加到列表中
    errs_col.append(err_col)  # 将当前颜色误差添加到列表中
    mags.append(mag)  # 将当前星等添加到列表中

col = [0] * (len(maglim) - 1)  # 初始化颜色列表，长度为maglim的长度减1

ax2.errorbar(col, mags, yerr=errs_mag, xerr=errs_col, fmt='o', color='k')  # 在第二个子图中绘制误差条

ax2.scatter(f115w - f200w, f115w, s=1, color='k')  # 在第二个子图中绘制散点图
ax2.text(xlim0 + 0.15, 15.5, "Output", fontdict=font2)  # 在图中添加文本“Output”

ax2.set_xlabel(filt1 + ' - ' + filt2, fontdict=font2)  # 设置第二个子图的x轴标签
ax2.set_ylabel(filt1, fontdict=font2)  # 设置第二个子图的y轴标签

plt.tight_layout()  # 调整子图参数，使之填充整个图像区域

输入与输出光度的比较

plt.figure(figsize=(14, 8))  # 创建一个14x8英寸的图形

plt.clf()  # 清除当前图形中的所有内容

ax1 = plt.subplot(2, 1, 1)  # 创建一个2行1列的子图，选择第一个子图

ax1.set_xlabel(filt1, fontdict=font2)  # 设置x轴标签为filt1，使用font2字体字典
ax1.set_ylabel(r'$\Delta$ Mag', fontdict=font2)  # 设置y轴标签为Δ Mag，使用font2字体字典

radec_input = SkyCoord(cat_sel['ra_in'], cat_sel['dec_in'], unit='deg')  # 创建SkyCoord对象，使用输入的RA和Dec

idx_f115w_cfr, d2d_f115w_cfr, _ = match_coordinates_sky(radec_input, radec_f115w)  # 匹配输入坐标与f115w坐标

sep_f115w_cfr = d2d_f115w_cfr < max_sep  # 计算匹配的分离度，判断是否小于最大分离度

f115w_inp_cfr = np.array(cat_sel['f115w_in'][sep_f115w_cfr])  # 获取输入的f115w数据
f115w_psf_cfr = np.array(f115w[idx_f115w_cfr[sep_f115w_cfr]])  # 获取PSF的f115w数据

diff_f115w_cfr = f115w_inp_cfr - f115w_psf_cfr  # 计算输入与PSF之间的差异
_, med_diff_f115w_cfr, sig_diff_f115w_cfr = sigma_clipped_stats(diff_f115w_cfr)  # 计算差异的中位数和标准差

xlim0 = 16  # 设置x轴下限
xlim1 = 24.5  # 设置x轴上限
ylim0 = np.mean(diff_f115w_cfr) - 0.5  # 设置y轴下限
ylim1 = np.mean(diff_f115w_cfr) + 0.5  # 设置y轴上限

ax1.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax1.set_ylim(ylim0, ylim1)  # 设置y轴范围

ax1.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax1.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

ax1.scatter(f115w_psf_cfr, diff_f115w_cfr, s=5, color='k')  # 绘制散点图，x为PSF数据，y为差异
ax1.plot([xlim0, xlim1], [0, 0], color='r', lw=5, ls='--')  # 绘制水平参考线

text = rf'{filt1} $\Delta$ Mag = {med_diff_f115w_cfr:5.3f} $\pm$ {sig_diff_f115w_cfr:5.3f}'  # 创建文本字符串，显示中位数和标准差
ax1.text(xlim0 + 0.05, ylim1 - 0.15, text, color='k', fontdict=font2)  # 在图中添加文本

ax2 = plt.subplot(2, 1, 2)  # 创建第二个子图

ax2.set_xlabel(filt2, fontdict=font2)  # 设置x轴标签为filt2，使用font2字体字典
ax2.set_ylabel(r'$\Delta$ Mag', fontdict=font2)  # 设置y轴标签为Δ Mag，使用font2字体字典

idx_f200w_cfr, d2d_f200w_cfr, _ = match_coordinates_sky(radec_input, radec_f200w)  # 匹配输入坐标与f200w坐标

sep_f200w_cfr = d2d_f200w_cfr < max_sep  # 计算匹配的分离度，判断是否小于最大分离度

f200w_inp_cfr = np.array(cat_sel['f200w_in'][sep_f200w_cfr])  # 获取输入的f200w数据
f200w_psf_cfr = np.array(f200w[idx_f200w_cfr[sep_f200w_cfr]])  # 获取PSF的f200w数据

diff_f200w_cfr = f200w_inp_cfr - f200w_psf_cfr  # 计算输入与PSF之间的差异
_, med_diff_f200w_cfr, sig_diff_f200w_cfr = sigma_clipped_stats(diff_f200w_cfr)  # 计算差异的中位数和标准差

xlim0 = 16  # 设置x轴下限
xlim1 = 24  # 设置x轴上限
ylim0 = np.mean(diff_f200w_cfr) - 0.5  # 设置y轴下限
ylim1 = np.mean(diff_f200w_cfr) + 0.5  # 设置y轴上限

ax2.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax2.set_ylim(ylim0, ylim1)  # 设置y轴范围

ax2.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax2.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

ax2.scatter(f200w_psf_cfr, diff_f200w_cfr, s=5, color='k')  # 绘制散点图，x为PSF数据，y为差异
ax2.plot([xlim0, xlim1], [0, 0], color='r', lw=5, ls='--')  # 绘制水平参考线

text = rf'{filt2} $\Delta$ Mag = {med_diff_f200w_cfr:5.3f} $\pm$ {sig_diff_f200w_cfr:5.3f}'  # 创建文本字符串，显示中位数和标准差

ax2.text(xlim0 + 0.05, ylim1 - 0.15, text, color='k', fontdict=font2)  # 在图中添加文本

plt.tight_layout()  # 调整子图参数，使之填充整个图像区域

plt.figure(figsize=(12, 6))  # 创建一个12x6英寸的图形

ax1 = plt.subplot(1, 2, 1)  # 创建1行2列的第一个子图

xlim0 = -10  # x轴下限
xlim1 = 10   # x轴上限
ylim0 = -10  # y轴下限
ylim1 = 10   # y轴上限

ax1.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax1.set_ylim(ylim0, ylim1)   # 设置y轴范围

ax1.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax1.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

ax1.set_xlabel(r'$\Delta$ RA (mas)', fontdict=font2)  # 设置x轴标签
ax1.set_ylabel(r'$\Delta$ Dec (mas)', fontdict=font2)  # 设置y轴标签
ax1.set_title(filt1, fontdict=font2)  # 设置子图标题

# 获取输入和PSF的RA值
ra_f115w_inp_cfr = np.array(cat_sel['ra_in'][sep_f115w_cfr])  # 输入RA值
ra_f115w_psf_cfr = np.array(radec_f115w.ra[idx_f115w_cfr[sep_f115w_cfr]])  # PSF RA值

# 获取输入和PSF的Dec值
dec_f115w_inp_cfr = np.array(cat_sel['dec_in'][sep_f115w_cfr])  # 输入Dec值
dec_f115w_psf_cfr = np.array(radec_f115w.dec[idx_f115w_cfr[sep_f115w_cfr]])  # PSF Dec值

dec_rad_f115w = np.radians(dec_f115w_psf_cfr)  # 将Dec值转换为弧度

# 计算RA差异
diffra_f115w_cfr = ((((ra_f115w_inp_cfr - ra_f115w_psf_cfr) * np.cos(dec_rad_f115w)) * u.deg).to(u.mas) / (1 * u.mas))  # RA差异转换为毫角秒

_, med_diffra_f115w_cfr, sig_diffra_f115w_cfr = sigma_clipped_stats(diffra_f115w_cfr)  # 计算RA差异的中位数和标准差

# 计算Dec差异
diffdec_f115w_cfr = (((dec_f115w_inp_cfr - dec_f115w_psf_cfr) * u.deg).to(u.mas) / (1 * u.mas))  # Dec差异转换为毫角秒

_, med_diffdec_f115w_cfr, sig_diffdec_f115w_cfr = sigma_clipped_stats(diffdec_f115w_cfr)  # 计算Dec差异的中位数和标准差

ax1.scatter(diffra_f115w_cfr, diffdec_f115w_cfr, s=1, color='k')  # 绘制RA和Dec差异的散点图
ax1.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')  # 绘制x=0的虚线
ax1.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')  # 绘制y=0的虚线

# 添加RA差异的文本
text = rf'$\Delta$ RA (mas) = {med_diffra_f115w_cfr:5.3f} $\pm$ {sig_diffra_f115w_cfr:5.3f}'
ax1.text(xlim0 + 0.05, ylim1 - 1.50, text, color='k', fontdict=font2)  # 在图中添加文本

# 添加Dec差异的文本
text = rf'$\Delta$ Dec (mas) = {med_diffdec_f115w_cfr:5.3f} $\pm$ {sig_diffdec_f115w_cfr:5.3f}'
ax1.text(xlim0 + 0.05, ylim1 - 3.0, text, color='k', fontdict=font2)  # 在图中添加文本

ax2 = plt.subplot(1, 2, 2)  # 创建1行2列的第二个子图

xlim0 = -10  # x轴下限
xlim1 = 10   # x轴上限
ylim0 = -10  # y轴下限
ylim1 = 10   # y轴上限

ax2.set_xlim(xlim0, xlim1)  # 设置x轴范围
ax2.set_ylim(ylim0, ylim1)   # 设置y轴范围
ax2.set_title(filt2, fontdict=font2)  # 设置子图标题

ax2.xaxis.set_major_locator(ticker.AutoLocator())  # 设置x轴主刻度定位器
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置x轴次刻度定位器
ax2.yaxis.set_major_locator(ticker.AutoLocator())  # 设置y轴主刻度定位器
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())  # 设置y轴次刻度定位器

ax2.set_xlabel(r'$\Delta$ RA (mas)', fontdict=font2)  # 设置x轴标签
ax2.set_ylabel(r'$\Delta$ Dec (mas)', fontdict=font2)  # 设置y轴标签

# 获取输入和PSF的RA值
ra_f200w_inp_cfr = np.array(cat_sel['ra_in'][sep_f200w_cfr])  # 输入RA值
ra_f200w_psf_cfr = np.array(radec_f200w.ra[idx_f200w_cfr[sep_f200w_cfr]])  # PSF RA值

# 获取输入和PSF的Dec值
dec_f200w_inp_cfr = np.array(cat_sel['dec_in'][sep_f200w_cfr])  # 输入Dec值
dec_f200w_psf_cfr = np.array(radec_f200w.dec[idx_f200w_cfr[sep_f200w_cfr]])  # PSF Dec值

dec_rad_f200w = np.radians(dec_f200w_psf_cfr)  # 将Dec值转换为弧度

# 计算RA差异
diffra_f200w_cfr = ((((ra_f200w_inp_cfr - ra_f200w_psf_cfr) * np.cos(dec_rad_f200w)) * u.deg).to(u.mas) / (1 * u.mas))  # RA差异转换为毫角秒

_, med_diffra_f200w_cfr, sig_diffra_f200w_cfr = sigma_clipped_stats(diffra_f200w_cfr)  # 计算RA差异的中位数和标准差

# 计算Dec差异
diffdec_f200w_cfr = (((dec_f200w_inp_cfr - dec_f200w_psf_cfr) * u.deg).to(u.mas) / (1 * u.mas))  # Dec差异转换为毫角秒

_, med_diffdec_f200w_cfr, sig_diffdec_f200w_cfr = sigma_clipped_stats(diffdec_f200w_cfr)  # 计算Dec差异的中位数和标准差

ax2.scatter(diffra_f200w_cfr, diffdec_f200w_cfr, s=1, color='k')  # 绘制RA和Dec差异的散点图
ax2.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')  # 绘制x=0的虚线
ax2.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')  # 绘制y=0的虚线

# 添加Dec差异的文本
text = rf'$\Delta$ Dec (mas) = {med_diffdec_f200w_cfr:5.3f} $\pm$ {sig_diffdec_f200w_cfr:5.3f}'
ax2.text(xlim0 + 0.05, ylim1 - 1.50, text, color='k', fontdict=font2)  # 在图中添加文本

# 添加RA差异的文本
text = rf'$\Delta$ RA (mas) = {med_diffra_f200w_cfr:5.3f} $\pm$ {sig_diffra_f200w_cfr:5.3f}'
ax2.text(xlim0 + 0.05, ylim1 - 3.0, text, color='k', fontdict=font2)  # 在图中添加文本

plt.tight_layout()  # 调整子图间距

最后说明

本笔记本提供了如何使用 photutils 包进行点扩散函数（PSF）光度测量的总体概述。在所有减缩步骤中采用的不同参数的选择以及PSF模型的选择取决于具体的用户科学案例。此外，只有在仪器调试后获得真实数据时，才能进行详细分析，以提供关于如何设置这些参数的建议，并概述在采用不同PSF模型（单一模型与PSF网格、网格中的PSF数量等）时输出光度测量的差异。在此背景下，我们注意到，选定的一个ERS项目（ERS 1334 - 解析恒星种群早期发布科学项目）将提供一个基本的测试基准，以探索上述不同选择如何影响拥挤恒星区域中PSF光度测量的质量。

关于此笔记本

作者: Matteo Correnti, JWST/NIRCam STScI 科学家 II & Larry Bradley, 数据分析工具分支副主任\

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

页面顶部