点源光度测量#
使用案例： 光度测量、交叉匹配目录、推导并应用光度零点。
数据： 使用 MIRAGE 模拟的 NIRCam 图像，并通过 JWST 流水线 处理大麦哲伦云 (LMC) 天体测量校准场。模拟使用了 4 点子像素抖动，针对三对宽滤光片：SW 通道的 F070W、F115W 和 F200W，以及 LW 通道的 F277W、F356W 和 F444W。我们仅模拟了 1 个 NIRCam SW 探测器（即 “NRCB1”）。在此示例中，我们使用两个 SW 滤光片（即 F115W 和 F200W）的 Level-3 图像（已校准和整形），并在每个图像中推导光度测量。其他滤光片的图像也可用，并可用于测试笔记本和/或不同滤光片组合。
工具： astropy, photutils。
跨仪器： NIRSpec, NIRISS, MIRI。
文档： 此笔记本是 STScI 更大 后处理数据分析工具生态系统 的一部分。
最终图表显示：
仪器颜色-星等图及误差
星等零点
校准颜色-星等图（与输入颜色-星等图比较）
输入与输出光度测量的比较
导入函数#
import os # 导入操作系统模块
import sys # 导入系统模块
import time # 导入时间模块
import numpy as np # 导入NumPy库，用于数值计算
import pandas as pd # 导入Pandas库，用于数据处理
import glob as glob # 导入glob模块，用于文件路径匹配
import tarfile # 导入tarfile模块，用于处理tar文件
import urllib.request # 导入urllib.request模块，用于处理URL请求
import jwst # 导入JWST相关模块
from jwst.datamodels import ImageModel # 从JWST数据模型中导入图像模型
from astropy.io import fits # 从Astropy库导入FITS文件处理模块
from astropy.table import Table # 从Astropy库导入表格处理模块
from astropy.visualization import (ZScaleInterval, SqrtStretch, ImageNormalize) # 导入图像可视化相关模块
from astropy.visualization import simple_norm # 导入简单归一化函数
from astropy.stats import sigma_clipped_stats, SigmaClip # 导入统计分析相关模块
from astropy.coordinates import SkyCoord, match_coordinates_sky # 导入天文坐标相关模块
from astropy import units as u # 导入单位模块
from photutils.detection import DAOStarFinder # 从Photutils库导入星点检测器
from photutils.background import MMMBackground, MADStdBackgroundRMS, Background2D # 导入背景处理相关模块
from photutils.aperture import CircularAperture, CircularAnnulus, aperture_photometry # 导入光度测量相关模块
导入绘图函数#
%matplotlib inline # 在Jupyter Notebook中内联显示图形
from matplotlib import style, pyplot as plt # 导入matplotlib的样式和绘图模块
import matplotlib.patches as patches # 导入绘图补丁模块，用于添加形状
import matplotlib.ticker as ticker # 导入刻度模块，用于自定义坐标轴刻度
# 设置图像的颜色映射为'viridis'
plt.rcParams['image.cmap'] = 'viridis'
# 设置图像的原点位置为左下角
plt.rcParams['image.origin'] = 'lower'
# 设置坐标轴标题和标签的字体大小为30
plt.rcParams['axes.titlesize'] = plt.rcParams['axes.labelsize'] = 30
# 设置x轴和y轴刻度标签的字体大小为30
plt.rcParams['xtick.labelsize'] = plt.rcParams['ytick.labelsize'] = 30
# 定义字体样式1
font1 = {'family': 'helvetica', 'color': 'black', 'weight': 'normal', 'size': '12'}
# 定义字体样式2
font2 = {'family': 'helvetica', 'color': 'black', 'weight': 'normal', 'size': '20'}
加载图像并创建一些有用的字典#
我们加载所有图像，并创建一个包含所有图像的字典，按探测器（detectors）和滤光片（filters）进行分类。这对于检查可用的探测器和滤光片非常有用，并帮助我们决定是否对所有图像进行光度测量（photometry），还是仅对子集进行（例如，仅对SW滤光片进行）。
我们还创建了一个包含一些分析所需参数的字典。该字典包含光度零点（photometric zeropoints，来自MIRAGE配置文件）和NIRCam点扩散函数（PSF）全宽半最大值（FWHM），这些数据来自NIRCam点扩散函数的JDox页面。FWHM是通过分析使用WebbPSF模拟的预期NIRCam PSF计算得出的。
注意：一旦每个探测器/滤光片的零点和FWHM值在调试完成后可用，该字典将会更新。
因此，我们有两个字典：
Level-3图像的字典
其他一些有用参数的字典
# 初始化一个字典，用于存储不同探测器的图像数据
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 = []
# 检查当前目录下是否存在已合并的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'
# 下载图像数据
urllib.request.urlretrieve(boxlink_images_lev3, boxfile_images_lev3)
# 解压下载的tar文件
tar = tarfile.open(boxfile_images_lev3, 'r')
tar.extractall(filter='data') # 添加filter='data'以跳过绝对路径或相对路径
# 设置图像目录
images_dir = './'
images = sorted(glob.glob(os.path.join(images_dir, "*combined*fits"))) # 获取所有合并的fits文件
else:
# 如果已存在合并的fits文件，直接获取文件列表
images_dir = './'
images = sorted(glob.glob(os.path.join(images_dir, "*combined*fits")))
# 遍历所有图像文件
for image in images:
im = fits.open(image) # 打开fits文件
f = im[0].header['FILTER'] # 获取过滤器信息
d = im[0].header['DETECTOR'] # 获取探测器信息
# 根据探测器类型进行重命名
if d == 'NRCBLONG':
d = 'NRCB5'
elif d == 'NRCALONG':
d = '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) # 打印可用的长波过滤器
# 定义过滤器列表
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]
# 定义模型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]
# 定义模型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和零点磁量关联起来
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))}
选择用于分析的探测器和/或滤光片#
如果我们只对分析中的某些滤光片（和/或某些探测器）感兴趣，如本例所示，我们可以从字典中选择这些滤光片（探测器）的图像，并仅分析这些图像。
在这个特定的例子中，我们分析探测器 NRCB1 的滤光片 F115W 和 F200W 的图像。
dets_short = ['NRCB1'] # 本例中感兴趣的探测器
filts_short = ['F115W', 'F200W'] # 本例中感兴趣的滤光片
显示图像#
为了检查我们的图像是否存在伪影并且可以用于分析，我们将其显示出来。
plt.figure(figsize=(14, 14)) # 创建一个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 # 获取图像数据
ax = plt.subplot(1, len(filts_short), i + 1) # 创建子图，行数为1，列数为滤镜数量，当前索引为i+1
plt.xlabel("X [px]", fontdict=font2) # 设置X轴标签
plt.ylabel("Y [px]", fontdict=font2) # 设置Y轴标签
plt.title(filt, fontdict=font2) # 设置子图标题为当前滤镜名称
norm = simple_norm(data_sb, 'sqrt', percent=99.) # 使用平方根归一化数据，99%分位数
ax.imshow(data_sb, norm=norm, cmap='Greys') # 显示图像数据，使用灰度色图
plt.tight_layout() # 调整子图布局以避免重叠
光圈光度法 (Aperture Photometry)#
有关使用 Photutils 进行光圈光度法的更多信息，请访问这里：光圈光度法
首先，我们创建一个有用的字典，其中包含：使用我们的查找算法检测到的源（DAOStarFinder），包含光圈光度结果的表格，以及一个最终的表格，包含检测到的源的位置（包括 x,y 和 RA,Dec）、仪器星等（instrumental magnitudes）及其误差。
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]['sources found'] = None # 初始化“找到的源”键，值为None
dict_aper[det][filt]['aperture phot table'] = None # 初始化“光圈光度表”键，值为None
dict_aper[det][filt]['final aperture phot table'] = None # 初始化“最终光圈光度表”键，值为None
1. 寻找源#
DAOStarFinder 使用 DAOFIND (Stetson 1987) 算法在图像中检测星星。DAOFIND 在图像中搜索局部密度最大值，这些最大值的峰值幅度大于 threshold（大致上；阈值应用于卷积图像），并且具有与定义的二维高斯核相似的大小和形状。
注意：我们采用了作为背景估计器的函数 MMMBackground，该函数使用 DAOPHOT MMM 算法计算整个图像的背景。
在处理可变背景和/或需要屏蔽没有数据的区域（例如，如果我们正在分析包含所有 4 个 NIRCam SW 探测器的图像，即包含芯片间隙的图像）时，我们可以设置 var_bkg = True，并使用一种更复杂的算法来考虑这些问题。
请注意，来自管道的 Level-2 和 Level-3 图像的单位是 MJy/sr（因此是表面亮度）。图像的实际单位可以通过头文件关键字 BUNIT 检查。标量转换常数被复制到头文件关键字 PHOTMJSR，该关键字提供从 DN/s 到 megaJy/steradian 的转换。对于我们的分析，我们回退到 DN/s。
def find_stars(det='NRCA1', filt='F070W', var_bkg=False):
# 定义查找星星的函数，参数包括探测器、滤光片和是否使用变背景
bkgrms = MADStdBackgroundRMS() # 初始化MAD标准背景RMS对象
mmm_bkg = MMMBackground() # 初始化MMM背景估计对象
image = fits.open(dict_images[det][filt]['images'][i]) # 打开指定探测器和滤光片的图像文件
data_sb = image[1].data # 获取图像数据
imh = image[1].header # 获取图像头信息
print("Finding sources on image {number}, filter {f}, detector {d}:".format(number=i + 1, f=filt, d=det)) # 打印当前处理的图像信息
print("") # 打印空行
data = data_sb / imh['PHOTMJSR'] # 将图像数据转换为DN/S单位
print("Conversion factor from {units} to DN/S for filter {f}:".format(units=imh['BUNIT'], f=filt), imh['PHOTMJSR']) # 打印转换因子
sigma_psf = dict_utils[filt]['psf fwhm'] # 获取滤光片的PSF FWHM值
print("FWHM for the filter {f}:".format(f=filt), sigma_psf, "px") # 打印FWHM值
if var_bkg: # 如果使用变背景
print("Using 2D Background") # 打印使用2D背景的信息
sigma_clip = SigmaClip(sigma=3.) # 初始化SigmaClip对象
coverage_mask = (data == 0) # 创建覆盖掩码，标记数据为0的区域
# 计算2D背景
bkg = Background2D(data, (25, 25), filter_size=(3, 3), sigma_clip=sigma_clip, bkg_estimator=mmm_bkg,
coverage_mask=coverage_mask, fill_value=0.0)
data_bkgsub = data.copy() # 复制数据
data_bkgsub = data_bkgsub - bkg.background # 从数据中减去背景
_, _, std = sigma_clipped_stats(data_bkgsub) # 计算经过sigma裁剪的统计量
else: # 如果不使用变背景
std = bkgrms(data) # 计算背景的标准差
bkg = mmm_bkg(data) # 估计背景
data_bkgsub = data.copy() # 复制数据
data_bkgsub -= bkg # 从数据中减去背景
# 初始化DAOStarFinder以查找星源
daofind = DAOStarFinder(threshold=th[j] * std, fwhm=sigma_psf, roundhi=1.0,
roundlo=-1.0, sharplo=0.30, sharphi=1.40)
found_stars = daofind(data_bkgsub) # 在背景减去后的数据中查找星源
dict_aper[det][filt]['sources found'] = found_stars # 将找到的星源存入字典
print("") # 打印空行
print("Number of sources found in the image:", len(found_stars)) # 打印找到的星源数量
print("-------------------------------------") # 打印分隔线
print("") # 打印空行
tic = time.perf_counter() # 记录开始时间
th = [10, 10] # 两个滤波器的阈值水平（长度必须与分析的滤波器数量匹配）
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(det=det, filt=filt, var_bkg=False) # 调用函数查找星星，背景变量设为False
toc = time.perf_counter() # 记录结束时间
print("Elapsed Time for finding stars:", toc - tic) # 输出查找星星的耗时
2. 光圈光度法 (Aperture Photometry)#
在下面的函数中，我们需要指定4个参数，以执行光圈光度法：用于分析的滤光片（即图像）、光圈半径以及用于确定局部背景的环形区域的内半径和外半径。光圈半径和环形区域的值取决于所采用的滤光片和用户的科学案例。
请注意，光圈对象支持多个位置，允许在每个位置执行多个光圈的光度测量。对于多个光圈，输出表的列名会附加位置索引。
def aperture_phot(det=det, filt=filt, radius=[3.5], sky_in=7, sky_out=10):
# 定义一个函数，用于进行光度测量，输入参数包括探测器、滤光片、半径、内外天区半径
positions = np.transpose((dict_aper[det][filt]['sources found']['xcentroid',
dict_aper[det][filt]['sources found']['ycentroid'])))
# 获取源的位置坐标，x和y中心坐标
image = fits.open(dict_images[det][filt]['images'][0])
# 打开指定的图像文件
data_sb = image[1].data
# 获取图像数据
imh = image[1].header
# 获取图像头信息
data = data_sb / imh['PHOTMJSR']
# 将数据归一化，除以光度单位
aa = np.argwhere(data < 0)
# 找到数据中小于0的元素的索引
for i in np.arange(0, len(aa), 1):
data[aa[i][0], aa[i][1]] = 0
# 将小于0的值置为0
error = np.sqrt(data)
# 计算数据的误差，使用平方根
tic = time.perf_counter()
# 记录开始时间
table_aper = Table()
# 创建一个空的表格，用于存储光度测量结果
for rad in radius:
print("Performing aperture photometry for filter {1}; radius r = {0} px".format(rad, filt))
# 打印当前正在处理的滤光片和半径信息
rr = str(rad)
# 将半径转换为字符串
aperture = CircularAperture(positions, r=rad)
# 创建一个圆形光圈对象
annulus_aperture = CircularAnnulus(positions, r_in=sky_in, r_out=sky_out)
# 创建一个圆环天区对象
annulus_mask = annulus_aperture.to_mask(method='center')
# 将圆环天区转换为掩模
bkg_median = []
# 初始化背景中值列表
bkg_stdev = []
# 初始化背景标准差列表
for mask in annulus_mask:
annulus_data = mask.multiply(data)
# 使用掩模提取圆环区域的数据
annulus_data_1d = annulus_data[mask.data > 0]
# 将掩模区域的数据展平为一维数组
_, median_sigclip, stdev_sigclip = sigma_clipped_stats(annulus_data_1d)
# 计算背景的中值和标准差，使用sigma剪切法
bkg_median.append(median_sigclip)
# 将中值添加到列表中
bkg_stdev.append(stdev_sigclip)
# 将标准差添加到列表中
bkg_median = np.array(bkg_median)
# 将背景中值列表转换为数组
bkg_stdev = np.array(bkg_stdev)
# 将背景标准差列表转换为数组
phot = aperture_photometry(data, aperture, method='exact', error=error)
# 进行光度测量，使用精确方法
phot['annulus_median'] = bkg_median
# 将背景中值添加到光度测量结果中
phot['aper_bkg'] = bkg_median * aperture.area
# 计算并添加光圈背景值
phot['aper_sum_bkgsub'] = phot['aperture_sum'] - phot['aper_bkg']
# 计算并添加背景校正后的光圈总和
table_aper.add_column(phot['aperture_sum'], name='aper_sum_' + rr + 'px')
# 将光圈总和添加到结果表中
table_aper.add_column(phot['annulus_median'], name='annulus_median_' + rr + 'px')
# 将背景中值添加到结果表中
table_aper.add_column(phot['aper_bkg'], name='aper_bkg_' + rr + 'px')
# 将光圈背景值添加到结果表中
table_aper.add_column(phot['aper_sum_bkgsub'], name='aper_sum_bkgsub_' + rr + 'px')
# 将背景校正后的光圈总和添加到结果表中
error_poisson = phot['aperture_sum_err']
# 获取光圈总和的泊松误差
error_scatter_sky = aperture.area * bkg_stdev**2
# 计算由于背景散射引起的误差
error_mean_sky = bkg_stdev**2 * aperture.area**2 / annulus_aperture.area
# 计算平均背景引起的误差
fluxerr = np.sqrt(error_poisson + error_scatter_sky + error_mean_sky)
# 计算总误差
table_aper.add_column(fluxerr, name='flux_err_' + rr + 'px')
# 将总误差添加到结果表中
dict_aper[det][filt]['aperture phot table'] = table_aper
# 将结果表存储到字典中
toc = time.perf_counter()
# 记录结束时间
print("Time Elapsed:", toc - tic)
# 打印耗时信息
return
# 返回函数结束
det = 'NRCB1' # 设置探测器为NRCB1
filt1 = 'F115W' # 设置第一个滤光片为F115W
filt2 = 'F200W' # 设置第二个滤光片为F200W
# 对于第一个滤光片进行光度测量，使用半径为3.5的光圈，天空区域内外半径分别为7和10
aperture_phot(det=det, filt=filt1, radius=[3.5], sky_in=7, sky_out=10)
# 对于第二个滤光片进行光度测量，使用半径为4.0的光圈，天空区域内外半径分别为6和9
aperture_phot(det=det, filt=filt2, radius=[4.0], sky_in=6, sky_out=9)
3. 包含星等和位置的表格#
我们创建最终的表格，包含每个目录的位置信息、星等和误差。
radii = [3.5, 4.0] # 定义半径列表
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): # 遍历图像索引
radius = str(radii[j]) # 将当前半径转换为字符串
image = fits.open(dict_images[det][filt]['images'][0]) # 打开图像文件
data = image[1].data # 获取图像数据
image_model = ImageModel(dict_images[det][filt]['images'][0]) # 创建图像模型
# 创建掩膜，筛选出有效的源
mask = ((dict_aper[det][filt]['sources found']['xcentroid'] > 0) &
(dict_aper[det][filt]['sources found']['xcentroid'] < data.shape[1]) &
(dict_aper[det][filt]['sources found']['ycentroid'] > 0) &
(dict_aper[det][filt]['sources found']['ycentroid'] < data.shape[0]) &
(dict_aper[det][filt]['aperture phot table']['aper_sum_bkgsub_' + radius + 'px'] > 0))
table_phot = Table() # 创建一个新的表格对象
# 将掩膜筛选后的x和y坐标存入表格
table_phot['x'] = dict_aper[det][filt]['sources found']['xcentroid'][mask]
table_phot['y'] = dict_aper[det][filt]['sources found']['ycentroid'][mask]
# 使用WCS将x和y坐标转换为RA和Dec
ra, dec = image_model.meta.wcs(table_phot['x'], table_phot['y'])
table_phot['radec'] = SkyCoord(ra, dec, unit='deg') # 将RA和Dec存入表格
# 计算并存储光度
table_phot[filt + '_inst'] = -2.5 * np.log10(dict_aper[det][filt]['aperture phot table']['aper_sum_bkgsub_' + radius + 'px'][mask])
# 计算并存储光度误差
table_phot['e' + filt + '_inst'] = 1.086 * (dict_aper[det][filt]['aperture phot table']['flux_err_' + radius + 'px'][mask] /
dict_aper[det][filt]['aperture phot table']['aper_sum_bkgsub_' + radius + 'px'][mask])
dict_aper[det][filt]['final aperture phot table'] = table_phot # 将最终的光度表存入字典
匹配源表#
我们交叉匹配目录以获得单一的颜色-亮度图（color-magnitude diagram）。
如果匹配之间的距离小于 0.5 像素（px），则来自两个滤光片的星体被关联。
max_sep = 0.015 * u.arcsec # 定义最大分离角度为0.015角秒
# 匹配两个滤波器的天体坐标
idx_inst, d2d_inst, _ = match_coordinates_sky(dict_aper[det][filt1]['final aperture phot table']['radec'],
dict_aper[det][filt2]['final aperture phot table']['radec'])
sep_constraint_inst = d2d_inst < max_sep # 计算分离角度是否小于最大分离角度
matched_sources = Table() # 创建一个新的表格以存储匹配的源
# 将匹配的源的坐标和其他信息添加到表格中
matched_sources['radec'] = dict_aper[det][filt1]['final aperture phot table']['radec'][sep_constraint_inst] # 添加滤波器1的坐标
matched_sources['x_' + filt1] = dict_aper[det][filt1]['final aperture phot table']['x'][sep_constraint_inst] # 添加滤波器1的x坐标
matched_sources['y_' + filt1] = dict_aper[det][filt1]['final aperture phot table']['y'][sep_constraint_inst] # 添加滤波器1的y坐标
matched_sources['x_' + filt2] = dict_aper[det][filt2]['final aperture phot table']['x'][idx_inst[sep_constraint_inst]] # 添加滤波器2的x坐标
matched_sources['y_' + filt2] = dict_aper[det][filt2]['final aperture phot table']['y'][idx_inst[sep_constraint_inst]] # 添加滤波器2的y坐标
matched_sources[filt1 + '_inst'] = dict_aper[det][filt1]['final aperture phot table'][filt1 + '_inst'][sep_constraint_inst] # 添加滤波器1的测量值
matched_sources['e' + filt1 + '_inst'] = dict_aper[det][filt1]['final aperture phot table']['e' + filt1 + '_inst'][sep_constraint_inst] # 添加滤波器1的误差
matched_sources[filt2 + '_inst'] = dict_aper[det][filt2]['final aperture phot table'][filt2 + '_inst'][idx_inst[sep_constraint_inst]] # 添加滤波器2的测量值
matched_sources['e' + filt2 + '_inst'] = dict_aper[det][filt2]['final aperture phot table']['e' + filt2 + '_inst'][idx_inst[sep_constraint_inst]] # 添加滤波器2的误差
仪器颜色-星等图及误差#
在误差图中，灰色点表示在每个滤光片中所有检测的测量值，而黑色点表示在两个目录之间匹配的检测值。
plt.figure(figsize=(12, 14)) # 创建一个12x14英寸的图形
plt.clf() # 清空当前图形
ax1 = plt.subplot(1, 2, 1) # 创建1行2列的子图，选择第1个子图
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 = -0.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轴次刻度定位器
ax1.scatter(matched_sources[filt1 + '_inst'] - matched_sources[filt2 + '_inst'], matched_sources[filt1 + '_inst'], # 绘制散点图
s=1, color='k') # 设置点的大小和颜色
ax2 = plt.subplot(2, 2, 2) # 创建2行2列的子图，选择第2个子图
ax2.set_xlabel(filt1 + '_inst', fontdict=font2) # 设置x轴标签
ax2.set_ylabel('$\sigma$ ' + filt1, fontdict=font2) # 设置y轴标签
xlim0 = -11 # x轴下限
xlim1 = -0.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(dict_aper[det][filt1]['final aperture phot table'][filt1 + '_inst'], # 绘制散点图
dict_aper[det][filt1]['final aperture phot table']['e' + filt1 + '_inst'], s=1, c='gray') # 设置点的大小和颜色
ax2.scatter(matched_sources[filt1 + '_inst'], matched_sources['e' + filt1 + '_inst'], s=1, c='k') # 绘制散点图
ax3 = plt.subplot(2, 2, 4) # 创建2行2列的子图，选择第4个子图
ax3.set_xlabel(filt2 + '_inst', fontdict=font2) # 设置x轴标签
ax3.set_ylabel('$\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(dict_aper[det][filt2]['final aperture phot table'][filt2 + '_inst'], # 绘制散点图
dict_aper[det][filt2]['final aperture phot table']['e' + filt2 + '_inst'], s=1, c='gray') # 设置点的大小和颜色
ax3.scatter(matched_sources[filt2 + '_inst'], matched_sources['e' + filt2 + '_inst'], s=1, c='k') # 绘制散点图
plt.tight_layout() # 调整子图间距以避免重叠
光度零点 (Photometric Zeropoints)#
为了获得最终校准的颜色-星等图 (color-magnitude diagram)，我们需要计算光度零点 (photometric zeropoints)。因此，我们需要对图像进行光圈光度测量 (aperture photometry)，使用与包围能量 (encircled energy, EE) 值相对应的光圈半径，并应用适当的光圈修正 (aperture correction)，因为所采用的光圈是有限的（上面字典中提供的值是针对无限光圈的）。然后，将其与上述确定的光圈光度测量进行比较。因此，我们可以将步骤总结如下：
使用适当的半径（相当于某个EE值）进行光圈光度测量
应用适当的光圈修正
应用表格化的零点
与光圈光度测量进行交叉匹配
注意：为了选择用于零点计算的恒星，我们可以使用CMD作为参考选择一个子样本，或者对一组明亮的孤立恒星进行光圈光度测量。要更改所使用的方法，我们可以分别设置 sample = cmd 或 sample = found。
加载孔径校正表#
注意：这些值是通过对合成的WebbPSF点扩散函数（PSFs）进行研究获得的。一旦我们获得了飞行中的测量数据，这些值将会更新。
import os # 导入os模块，用于文件和目录操作
import urllib.request # 导入urllib.request模块，用于下载文件
import pandas as pd # 导入pandas库，用于数据处理
# 检查是否存在名为'aperture_correction_table.txt'的文件
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' # 定义下载后保存的文件路径
# 下载文件并保存到指定路径
urllib.request.urlretrieve(boxlink_apcorr_table, boxfile_apcorr_table)
ap_tab = './aperture_correction_table.txt' # 设置文件路径
# 读取'aperture_correction_table.txt'文件，使用空格作为分隔符
aper_table = pd.read_csv(ap_tab, header=None, sep='\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()
1a. 使用CMD作为参考进行子采样#
注意：如果用户希望施加更严格的条件，可以在子样本中选择没有邻近星星距离小于X像素的星星，其中X可以适当固定，这可以通过设置dist_sel = True来实现。
限制星等和最近邻星的最大距离取决于用户的科学案例（即视场中的星星数量、拥挤程度、亮源数量等），必须相应地进行修改。
最终表格xy_cmd包含在两个滤光片中所选星星的位置（X,Y）。
xy_cmd = Table() # 创建一个空的表格用于存储选中的星体坐标
mag_lim_bright = -7.0 # 明亮星体的亮度限制
mag_lim_faint = -4.5 # 微弱星体的亮度限制
# 创建一个掩码，筛选出亮度在限制范围内的星体
mask_xy = ((matched_sources[filt1 + '_inst'] < mag_lim_faint) & (matched_sources[filt1 + '_inst'] > mag_lim_bright))
# 将符合条件的星体坐标添加到xy_cmd表中
xy_cmd['x_' + filt1] = matched_sources['x_' + filt1][mask_xy] # 添加filt1的x坐标
xy_cmd['y_' + filt1] = matched_sources['y_' + filt1][mask_xy] # 添加filt1的y坐标
xy_cmd['x_' + filt2] = matched_sources['x_' + filt2][mask_xy] # 添加filt2的x坐标
xy_cmd['y_' + filt2] = matched_sources['y_' + filt2][mask_xy] # 添加filt2的y坐标
print("Number of stars selected from CMD:", len(xy_cmd)) # 打印选中的星体数量
dist_sel = False # 是否计算距离的标志
if dist_sel: # 如果选择计算距离
print("") # 打印空行
print("Calculating closest neighbour distance") # 打印计算最近邻距离的提示
d_filt1 = [] # 存储filt1的最小距离
d_filt2 = [] # 存储filt2的最小距离
min_sep_filt1 = 10 # filt1的最小距离阈值
min_sep_filt2 = 10 # filt2的最小距离阈值
# 计算filt1的最小距离
for i, j in zip(xy_cmd['x_' + filt1], xy_cmd['y_' + filt1]):
sep = [] # 存储距离
# 计算当前星体到所有其他星体的距离
dist = np.sqrt((dict_aper[det][filt1]['sources found']['xcentroid'] - i)**2 + (dict_aper[det][filt1]['sources found']['ycentroid'] - j)**2)
sep = np.sort(dist)[1:2][0] # 找到最近的邻居距离
d_filt1.append(sep) # 将距离添加到列表中
# 计算filt2的最小距离
for i, j in zip(xy_cmd['x_' + filt2], xy_cmd['y_' + filt2]):
sep = [] # 存储距离
# 计算当前星体到所有其他星体的距离
dist = np.sqrt((dict_aper[det][filt2]['sources found']['xcentroid'] - i)**2 + (dict_aper[det][filt2]['sources found']['ycentroid'] - j)**2)
sep = np.sort(dist)[1:2][0] # 找到最近的邻居距离
d_filt2.append(sep) # 将距离添加到列表中
# 将最小距离添加到xy_cmd表中
xy_cmd['min distance ' + filt1] = d_filt1
xy_cmd['min distance ' + filt2] = d_filt2
# 创建一个掩码，筛选出符合最小距离要求的星体
mask_dist = ((xy_cmd['min distance ' + filt1] > min_sep_filt1) & (xy_cmd['min distance ' + filt2] > min_sep_filt2))
xy_cmd = xy_cmd[mask_dist] # 应用掩码，更新xy_cmd表
print("") # 打印空行
print("Minimum distance required for filter {f}:".format(f=filt1), min_sep_filt1, "px") # 打印filt1的最小距离要求
print("Minimum distance required for filter {f}:".format(f=filt2), min_sep_filt2, "px") # 打印filt2的最小距离要求
print("") # 打印空行
print("Number of stars selected from CMD including distance constraints:", len(xy_cmd)) # 打印考虑距离约束后选中的星体数量
如我们之前所做，我们创建一个字典，其中包含每个图像的衍生光圈光度（aperture photometry，针对明亮恒星）的表格。
dict_zp = {} # 初始化一个空字典，用于存储每个探测器和滤光片的相关信息
for det in dets_short: # 遍历每个短探测器
dict_zp.setdefault(det, {}) # 如果字典中没有该探测器，则初始化一个空字典
for j, filt in enumerate(filts_short): # 遍历每个短滤光片
dict_zp[det].setdefault(filt, {}) # 如果该探测器下没有该滤光片，则初始化一个空字典
dict_zp[det][filt]['sources found'] = None # 初始化'找到的源'为None
dict_zp[det][filt]['aperture phot table'] = None # 初始化'光圈光度表'为None
dict_zp[det][filt]['final aperture phot table'] = None # 初始化'最终光圈光度表'为None
def find_bright_stars(det='NRCA1', filt='F070W', var_bkg=False, dist_sel=False):
# 定义寻找亮星的函数，参数包括探测器、滤光片、是否使用变背景和是否选择距离
bkgrms = MADStdBackgroundRMS() # 初始化MAD标准背景RMS对象
mmm_bkg = MMMBackground() # 初始化MMM背景估计对象
image = fits.open(dict_images[det][filt]['images'][i]) # 打开指定探测器和滤光片的图像文件
data_sb = image[1].data # 获取图像数据
imh = image[1].header # 获取图像头信息
print("Finding sources on image {number}, filter {f}, detector {d}:".format(number=i + 1, f=filt, d=det)) # 打印当前处理的图像信息
print("") # 打印空行
data = data_sb / imh['PHOTMJSR'] # 将图像数据转换为DN/s
print("Conversion factor from {units} to DN/s for filter {f}:".format(units=imh['BUNIT'], f=filt), imh['PHOTMJSR']) # 打印转换因子
sigma_psf = dict_utils[filt]['psf fwhm'] # 获取滤光片的PSF FWHM
print("FWHM for the filter {f}:".format(f=filt), sigma_psf, "px") # 打印FWHM信息
if var_bkg: # 如果使用变背景
print("Using 2D Background") # 打印使用2D背景的信息
sigma_clip = SigmaClip(sigma=3.) # 初始化SigmaClip对象
coverage_mask = (data == 0) # 创建覆盖掩码，标记数据为0的区域
# 计算2D背景
bkg = Background2D(data, (25, 25), filter_size=(3, 3), sigma_clip=sigma_clip,
bkg_estimator=mmm_bkg, coverage_mask=coverage_mask, fill_value=0.0)
data_bkgsub = data.copy() # 复制数据
data_bkgsub = data_bkgsub - bkg.background # 从数据中减去背景
_, _, std = sigma_clipped_stats(data_bkgsub) # 计算背景减去后的标准差
else: # 如果不使用变背景
std = bkgrms(data) # 计算背景RMS
bkg = mmm_bkg(data) # 估计背景
data_bkgsub = data.copy() # 复制数据
data_bkgsub -= bkg # 从数据中减去背景
# 初始化DAOStarFinder以寻找星源
daofind = DAOStarFinder(threshold=th[j] * std, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
zp_stars = daofind(data_bkgsub) # 在背景减去后的数据中寻找星源
dict_zp[det][filt]['sources found'] = zp_stars # 将找到的星源存入字典
if dist_sel: # 如果选择距离
print("") # 打印空行
print("Calculating closest neighbour distance") # 打印计算最近邻距离的信息
d = [] # 初始化距离列表
# 遍历找到的星源
for xx, yy in zip(zp_stars['xcentroid'], zp_stars['ycentroid']):
sep = [] # 初始化分离距离列表
dist = np.sqrt((dict_aper[det][filt]['sources found']['xcentroid'] - xx)**2 +
(dict_aper[det][filt]['sources found']['ycentroid'] - yy)**2) # 计算距离
sep = np.sort(dist)[1:2][0] # 获取最近的距离
d.append(sep) # 将距离添加到列表中
zp_stars['min distance'] = d # 将最小距离添加到星源数据中
mask_dist = (zp_stars['min distance'] > min_sep[j]) # 创建距离掩码
zp_stars = zp_stars[mask_dist] # 应用距离掩码
dict_zp[det][filt]['sources found'] = zp_stars # 更新找到的星源字典
print("Minimum distance required:", min_sep[j], "px") # 打印所需的最小距离
print("") # 打印空行
print("Number of bright isolated sources found in the image:", len(zp_stars)) # 打印找到的孤立亮星数量
print("-----------------------------------------------------") # 打印分隔线
print("") # 打印空行
else: # 如果不选择距离
print("") # 打印空行
print("Number of bright sources found in the image:", len(zp_stars)) # 打印找到的亮星数量
print("--------------------------------------------") # 打印分隔线
print("") # 打印空行
阈值（threshold）和最近邻（closest neighbour）的最大距离（maximum distance）取决于用户的科学案例（science case）（即视场内的星星数量、拥挤程度、亮源数量等），必须相应地进行调整。
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[det][filt]['images']), 1):
# 查找明亮的星星，参数包括检测、滤波器、是否考虑背景变化、是否选择距离
find_bright_stars(det=det, filt=filt, var_bkg=False, dist_sel=False)
toc = time.perf_counter() # 记录结束时间
# 打印查找星星所用的时间
print("Elapsed Time for finding stars:", toc - tic)
明亮孤立星的光圈光度测量#
要选择我们希望在光圈光度测量中使用的星星样本，请修改 sample 关键字，使用 cmd 或 found。
def aperture_phot_zp(det=det, filt=filt, sample='found'):
# 定义一个函数，用于进行光度测量，输入参数包括探测器、滤光片和样本类型
radii = [aper_table.loc[filt]['r70']] # 获取对应滤光片的70%半径
ees = '70'.split() # 定义有效半径的列表
ee_radii = dict(zip(ees, radii)) # 将有效半径与其对应的值组合成字典
if sample == 'found':
# 如果样本类型为'found'，则获取已找到的源的位置
positions = np.transpose((dict_zp[det][filt]['sources found']['xcentroid'],
dict_zp[det][filt]['sources found']['ycentroid']))
if sample == 'cmd':
# 如果样本类型为'cmd'，则获取命令中的源位置
positions = np.transpose((xy_cmd['x_' + filt], xy_cmd['y_' + filt]))
image = fits.open(dict_images[det][filt]['images'][0]) # 打开对应的图像文件
data_sb = image[1].data # 获取图像数据
imh = image[1].header # 获取图像头信息
data = data_sb / imh['PHOTMJSR'] # 进行数据的光度校正
aa = np.argwhere(data < 0) # 查找数据中小于0的值的索引
for i in np.arange(0, len(aa), 1):
# 将小于0的值替换为0
data[aa[i][0], aa[i][1]] = 0
# 从光圈校正表中获取天空背景的内外半径
sky = {"sky_in": aper_table.loc[filt]['r80'], "sky_out": aper_table.loc[filt]['r85']}
tic = time.perf_counter() # 记录开始时间
table_zp = Table() # 创建一个新的表格用于存储结果
for ee, radius in ee_radii.items():
# 对每个有效半径进行光度测量
print("Performing aperture photometry for radius equivalent to EE = {0}% for filter {1}".format(ee, 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) # 计算背景中值
bkg_median.append(median_sigclip) # 将中值添加到列表中
bkg_median = np.array(bkg_median) # 转换为数组
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 # 计算校正后的光度
# 将结果添加到表格中
table_zp.add_column(phot['aperture_sum'], name='aper_sum_' + ee)
table_zp.add_column(phot['annulus_median'], name='annulus_median_' + ee)
table_zp.add_column(phot['aper_bkg'], name='aper_bkg_ee_' + ee)
table_zp.add_column(phot['aper_sum_bkgsub'], name='aper_sum_bkgsub_' + ee)
table_zp.add_column(phot['aper_sum_corrected'], name='aper_sum_corrected_' + ee + '_' + filt)
dict_zp[det][filt]['aperture phot table'] = table_zp # 将结果表格存储到字典中
toc = time.perf_counter() # 记录结束时间
print("Time Elapsed:", toc - tic) # 输出耗时
return # 返回函数结束
# 使用指定的检测器和滤光片进行光度测量，样本为'cmd'
aperture_phot_zp(det=det, filt=filt1, sample='cmd') # 对第一个滤光片进行光度测量
# 使用指定的检测器和第二个滤光片进行光度测量，样本为'cmd'
aperture_phot_zp(det=det, filt=filt2, sample='cmd') # 对第二个滤光片进行光度测量
ee = '70' # 设置等效光度的值
sample = 'cmd' # 设置样本类型为'cmd'
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): # 遍历图像列表
image = fits.open(dict_images[det][filt]['images'][0]) # 打开图像文件
data = image[1].data # 获取图像数据
image_model = ImageModel(dict_images[det][filt]['images'][0]) # 创建图像模型
table_zp_mag = Table() # 创建一个新的表格用于存储零点光度数据
if sample == 'found': # 如果样本类型为'found'
mask = ((dict_zp[det][filt]['sources found']['xcentroid'] > 0) & # 创建掩码，筛选有效的x坐标
(dict_zp[det][filt]['sources found']['xcentroid'] < data.shape[1]) & # x坐标小于图像宽度
(dict_zp[det][filt]['sources found']['ycentroid'] > 0) & # 筛选有效的y坐标
(dict_zp[det][filt]['sources found']['ycentroid'] < data.shape[0]) & # y坐标小于图像高度
(dict_zp[det][filt]['aperture phot table']['aper_sum_corrected_' + ee + '_' + filt] > 0)) # 确保光度值大于0
table_zp_mag['x'] = dict_zp[det][filt]['sources found']['xcentroid'][mask] # 存储x坐标
table_zp_mag['y'] = dict_zp[det][filt]['sources found']['ycentroid'][mask] # 存储y坐标
ra, dec = image_model.meta.wcs(table_zp_mag['x'], table_zp_mag['y']) # 将像素坐标转换为天球坐标
table_zp_mag['radec'] = SkyCoord(ra, dec, unit='deg') # 存储天球坐标
table_zp_mag[filt + '_inst'] = -2.5 * np.log10(dict_zp[det][filt]['aperture phot table']['aper_sum_corrected_' + ee + '_' + filt][mask]) # 计算并存储瞬时光度
table_zp_mag[filt + '_zp'] = table_zp_mag[filt + '_inst'] + dict_utils[filt]['VegaMAG zp modB'] # 计算并存储零点光度
dict_zp[det][filt]['final aperture phot table'] = table_zp_mag # 更新字典中的最终光度表
if sample == 'cmd': # 如果样本类型为'cmd'
mask = ((dict_zp[det][filt]['aperture phot table']['aper_sum_corrected_' + ee + '_' + filt] > 0)) # 创建掩码，筛选有效的光度值
table_zp_mag['x'] = xy_cmd['x_' + filt][mask] # 存储x坐标
table_zp_mag['y'] = xy_cmd['y_' + filt][mask] # 存储y坐标
ra, dec = image_model.meta.wcs(table_zp_mag['x'], table_zp_mag['y']) # 将像素坐标转换为天球坐标
table_zp_mag['radec'] = SkyCoord(ra, dec, unit='deg') # 存储天球坐标
table_zp_mag[filt + '_inst'] = -2.5 * np.log10(dict_zp[det][filt]['aperture phot table']['aper_sum_corrected_' + ee + '_' + filt][mask]) # 计算并存储瞬时光度
table_zp_mag[filt + '_zp'] = table_zp_mag[filt + '_inst'] + dict_utils[filt]['VegaMAG zp modB'] # 计算并存储零点光度
dict_zp[det][filt]['final aperture phot table'] = table_zp_mag # 更新字典中的最终光度表
推导零点 (Zeropoints)#
plt.figure(figsize=(14, 8)) # 创建一个14x8英寸的图形
plt.clf() # 清除当前图形中的所有内容
ax1 = plt.subplot(2, 1, 1) # 创建一个2行1列的子图，选择第一个子图
ax1.set_xlabel(filt1 + '_inst', fontdict=font2) # 设置x轴标签
ax1.set_ylabel(filt1 + '_ZP', fontdict=font2) # 设置y轴标签
# 匹配天文坐标，获取索引和距离
idx_zp_1, d2d_zp_1, _ = match_coordinates_sky(dict_zp[det][filt1]['final aperture phot table']['radec'],
dict_aper[det][filt1]['final aperture phot table']['radec'])
sep_constraint_zp_1 = d2d_zp_1 < max_sep # 设置距离约束条件
# 根据约束条件提取匹配的零点和观测值
f115w_apZP_matched = np.array(dict_zp[det][filt1]['final aperture phot table'][filt1 + '_zp'][sep_constraint_zp_1])
f115w_ap_matched = np.array(dict_aper[det][filt1]['final aperture phot table'][filt1 + '_inst'][idx_zp_1[sep_constraint_zp_1]])
diff_f115w = f115w_apZP_matched - f115w_ap_matched # 计算差值
_, zp_f115w, zp_sigma_f115w = sigma_clipped_stats(diff_f115w) # 计算零点和标准差
# 设置x轴和y轴的限制
xlim0 = np.min(f115w_ap_matched) - 0.25
xlim1 = np.max(f115w_ap_matched) + 0.25
ylim0 = np.mean(diff_f115w) - 0.35
ylim1 = np.mean(diff_f115w) + 0.35
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_ap_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.1, filt1 + ' Zeropoint = %5.3f $\pm $ %5.3f' % (zp_f115w, zp_sigma_f115w),
color='k', fontdict=font2) # 添加文本注释
ax2 = plt.subplot(2, 1, 2) # 创建第二个子图
ax2.set_xlabel(filt2 + '_inst', fontdict=font2) # 设置x轴标签
ax2.set_ylabel(filt2 + '_ZP', fontdict=font2) # 设置y轴标签
# 匹配天文坐标，获取索引和距离
idx_zp_2, d2d_zp_2, _ = match_coordinates_sky(dict_zp[det][filt2]['final aperture phot table']['radec'],
dict_aper[det][filt2]['final aperture phot table']['radec'])
sep_constraint_zp_2 = d2d_zp_2 < max_sep # 设置距离约束条件
# 根据约束条件提取匹配的零点和观测值
f200w_apZP_matched = np.array(dict_zp[det][filt2]['final aperture phot table'][filt2 + '_zp'][sep_constraint_zp_2])
f200w_ap_matched = np.array(dict_aper[det][filt2]['final aperture phot table'][filt2 + '_inst'][idx_zp_2[sep_constraint_zp_2]])
diff_f200w = f200w_apZP_matched - f200w_ap_matched # 计算差值
_, zp_f200w, zp_sigma_f200w = sigma_clipped_stats(diff_f200w) # 计算零点和标准差
# 设置x轴和y轴的限制
xlim0 = np.min(f200w_ap_matched) - 0.25
xlim1 = np.max(f200w_ap_matched) + 0.25
ylim0 = np.mean(diff_f200w) - 0.35
ylim1 = np.mean(diff_f200w) + 0.35
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_ap_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.1, filt2 + ' Zeropoint = %5.3f $\pm$ %5.3f' % (zp_f200w, zp_sigma_f200w),
color='k', fontdict=font2) # 添加文本注释
plt.tight_layout() # 调整子图间距
导入输入目录#
import os # 导入os模块，用于文件和目录操作
import urllib.request # 导入urllib.request模块，用于下载文件
import pandas as pd # 导入pandas库，用于数据处理
# 检查当前目录下是否存在'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' # 定义下载后保存的文件名
# 下载文件并保存到指定路径
urllib.request.urlretrieve(boxlink_input_cat, boxfile_input_cat)
input_cat = './pointsource.cat' # 设置输入目录为下载的文件
# 读取输入的点源目录文件，指定列名和数据格式
cat = pd.read_csv(input_cat, header=None, sep='\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() # 输出数据框的前5行
# 计算给定探测器和滤波器的最小和最大RA值
lim_ra_min = np.min(dict_aper[det][filt1]['final aperture phot table']['radec'].ra) # 最小RA值
lim_ra_max = np.max(dict_aper[det][filt1]['final aperture phot table']['radec'].ra) # 最大RA值
# 计算给定探测器和滤波器的最小和最大DEC值
lim_dec_min = np.min(dict_aper[det][filt1]['final aperture phot table']['radec'].dec) # 最小DEC值
lim_dec_max = np.max(dict_aper[det][filt1]['final aperture phot table']['radec'].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) # RA范围条件
& (cat['dec_in'] < lim_dec_max)] # DEC范围条件
校准的颜色-星等图#
plt.figure(figsize=(12, 14)) # 创建一个12x14英寸的图形
plt.clf() # 清空当前图形
font2 = {'size': '30'} # 设置字体大小为30
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') # 绘制散点图，x为差值，y为f115w的值
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 = matched_sources[filt1 + '_inst'] + zp_f115w # 计算f115w的值
f200w = matched_sources[filt2 + '_inst'] + zp_f200w # 计算f200w的值
maglim = np.arange(18, 25, 1) # 创建一个从18到25的范围，步长为1
mags = [] # 初始化mags列表
errs_mag = [] # 初始化errs_mag列表
errs_col = [] # 初始化errs_col列表
for i in np.arange(0, len(maglim) - 1, 1): # 遍历maglim的每个区间
mag = (maglim[i] + maglim[i + 1]) / 2 # 计算当前区间的中值
err_mag1 = matched_sources['e' + filt1 + '_inst'][(f115w > maglim[i]) & (f115w < maglim[i + 1])] # 获取f115w在当前区间的误差
err_mag2 = matched_sources['e' + filt2 + '_inst'][(f115w > maglim[i]) & (f115w < maglim[i + 1])] # 获取f200w在当前区间的误差
err_mag = np.mean(err_mag1[i]) # 计算f115w的平均误差
err_temp = np.sqrt(err_mag1**2 + err_mag2**2) # 计算误差的平方和的平方根
err_col = np.mean(err_temp[i]) # 计算颜色的平均误差
errs_mag.append(err_mag) # 将当前误差添加到errs_mag列表
errs_col.append(err_col) # 将当前颜色误差添加到errs_col列表
mags.append(mag) # 将当前中值添加到mags列表
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') # 绘制散点图，x为f115w和f200w的差值，y为f115w的值
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
ax1.set_ylabel('$\Delta$ Mag', fontdict=font2) # 设置y轴标签为Δ Mag
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, matched_sources['radec']) # 匹配输入坐标与已匹配源的坐标
sep_f115w_cfr = d2d_f115w_cfr < max_sep # 计算匹配源与输入源之间的距离，筛选出小于最大分离的源
f115w_inp_cfr = np.array(cat_sel['f115w_in'][sep_f115w_cfr]) # 获取输入的f115w数据
f115w_ap_cfr = np.array(f115w[idx_f115w_cfr[sep_f115w_cfr]]) # 获取匹配的f115w数据
diff_f115w_cfr = f115w_inp_cfr - f115w_ap_cfr # 计算输入和匹配的f115w数据之间的差异
diff_f115w_cfr_sel = (f115w_inp_cfr - f115w_ap_cfr)[(f115w_ap_cfr > 17.75) & (f115w_ap_cfr < 21.5)] # 筛选出特定范围内的差异
_, med_diff_f115w_cfr, sig_diff_f115w_cfr = sigma_clipped_stats(diff_f115w_cfr_sel) # 计算中位数和标准差
xlim0 = 16 # 设置x轴下限
xlim1 = 24 # 设置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_ap_cfr, diff_f115w_cfr, s=5, color='k') # 绘制散点图，显示f115w匹配数据与差异
ax1.plot([xlim0, xlim1], [0, 0], color='r', lw=5, ls='--') # 绘制红色虚线，表示0差异线
ax1.plot([17.75, 17.75], [ylim0, ylim1], color='k', lw=1, ls='--') # 绘制f115w下限虚线
ax1.plot([21.5, 21.5], [ylim0, ylim1], color='k', lw=1, ls='--') # 绘制f115w上限虚线
ax1.text(xlim0 + 0.05, ylim1 - 0.15, filt1 + ' $\Delta$ Mag = %5.3f $\pm$ %5.3f' # 在图中添加文本，显示filt1的差异统计
% (med_diff_f115w_cfr, sig_diff_f115w_cfr), color='k', fontdict=font2)
ax2 = plt.subplot(2, 1, 2) # 创建第二个子图
ax2.set_xlabel(filt2, fontdict=font2) # 设置x轴标签为filt2
ax2.set_ylabel('$\Delta$ Mag', fontdict=font2) # 设置y轴标签为Δ Mag
idx_f200w_cfr, d2d_f200w_cfr, _ = match_coordinates_sky(radec_input, matched_sources['radec']) # 匹配输入坐标与已匹配源的坐标
sep_f200w_cfr = d2d_f200w_cfr < max_sep # 计算匹配源与输入源之间的距离，筛选出小于最大分离的源
f200w_inp_cfr = np.array(cat_sel['f200w_in'][sep_f200w_cfr]) # 获取输入的f200w数据
f200w_ap_cfr = np.array(f200w[idx_f200w_cfr[sep_f200w_cfr]]) # 获取匹配的f200w数据
diff_f200w_cfr = f200w_inp_cfr - f200w_ap_cfr # 计算输入和匹配的f200w数据之间的差异
diff_f200w_cfr_sel = (f200w_inp_cfr - f200w_ap_cfr)[(f200w_ap_cfr > 16.5) & (f200w_ap_cfr < 20.5)] # 筛选出特定范围内的差异
_, med_diff_f200w_cfr, sig_diff_f200w_cfr = sigma_clipped_stats(diff_f200w_cfr_sel) # 计算中位数和标准差
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_ap_cfr, diff_f200w_cfr, s=5, color='k') # 绘制散点图，显示f200w匹配数据与差异
ax2.plot([xlim0, xlim1], [0, 0], color='r', lw=5, ls='--') # 绘制红色虚线，表示0差异线
ax2.plot([16.5, 16.5], [ylim0, ylim1], color='k', lw=1, ls='--') # 绘制f200w下限虚线
ax2.plot([20.5, 20.5], [ylim0, ylim1], color='k', lw=1, ls='--') # 绘制f200w上限虚线
ax2.text(xlim0 + 0.05, ylim1 - 0.15, filt2 + ' $\Delta$ Mag = %5.3f $\pm$ %5.3f' # 在图中添加文本，显示filt2的差异统计
% (med_diff_f200w_cfr, sig_diff_f200w_cfr), color='k', fontdict=font2)
plt.tight_layout() # 调整子图间距，使其更美观
关于此笔记本#
作者: Matteo Correnti, JWST/NIRCam STScI 科学家 II \
更新时间: 2021-01-15
