交叉滤波 PSF 匹配#
用例： 在一个星系外的“空白”区域测量星系光度。相关于 JDox 科学用例 #22。
数据： JWST 模拟的 NIRCam 图像来自 JADES JWST GTO 星系外空白区域。
（Williams et al. 2018）https://ui.adsabs.harvard.edu/abs/2018ApJS..236…33W。
工具： photutils, matplotlib。
跨仪器： 可能涉及 NIRISS, MIRI。
文档： 本笔记本是 STScI 更大 后处理数据分析工具生态系统 的一部分。
引言#
本笔记本使用 photutils 在 NIRCam 深成像中检测天体/星系。首先在 F200W 图像中进行检测，然后在所有 9 个滤光片（F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, F444W）中获取等光度光度测量。使用 PSF 匹配来校正在长波长图像（红色于 F200W 之上）中测量的光度。
在保存光度目录后，笔记本演示了如何在新会话中加载该笔记本。
它展示了对完整目录和单个星系的一些简单分析。
通过将测量的颜色与模拟输入颜色进行比较，
它显示经过 PSF 校正后，测量更为准确，尽管仍可进一步改进。
该笔记本仅分析完整 JADES 模拟的中心 1000 x 1000 像素（30” x 30”）。这些切片已获得作者（Williams et al.）的许可，并在 STScI 进行存档。
注意事项：
这是一项正在进行中的工作。可能获得更准确的光度测量。
这些 JADES 图像是使用 Guitarra 模拟的，而不是 MIRAGE。它们具有不同的属性和单位（e-/s），与 JWST 管道产品（MJy/sr）不同。
所有图像在分析前都已对齐至相同的 0.03” 像素网格。如果需要，可以使用
reproject进行此对齐。通过考虑相关噪声和/或更直接地测量每个图像中的噪声，可以进一步改善通量不确定性计算。
开发者笔记：
以下是报告的问题摘要：
单位在
plot中可以正常工作，但与errorbar和text不兼容。流量单位可以自动转换为 AB 亮度，但流量 不确定性 不能自动转换为亮度不确定性。
次要轴应该自动处理流量和 AB 亮度单位之间的转换。
sharex和sharey在 WCSprojection中无法正常工作。
我也无法弄明白如何：
使绘图坐标轴自动缩放到仅 部分 绘制的数据。
在交互式图中使工具提示悬停在光标下的数据点上。
# 检查PEP 8风格格式
# %load_ext pycodestyle_magic # 加载pycodestyle_magic扩展以检查代码风格
# %flake8_on --ignore E261,E501,W291,W2,E302,E305 # 启用flake8检查，忽略特定的错误代码
导入包#
Numpy 用于一般数组计算
Photutils 用于光度计算、PSF 匹配
Astropy 用于处理 FITS、WCS、表格、单位、彩色图像、绘图、卷积
os 和 glob 用于文件管理
copy 用于表格修改
Matplotlib 用于绘图
Watermark 用于检查所有导入包的版本（可选）
import numpy as np # 导入NumPy库，用于数值计算
import os # 导入os库，用于操作系统相关功能
from glob import glob # 从glob模块导入glob函数，用于文件路径匹配
from copy import deepcopy # 从copy模块导入deepcopy函数，用于深拷贝对象
import astropy # 导入Astropy库，版本4.2用于将亮度写入ecsv文件
from astropy.io import fits # 从astropy.io模块导入fits，用于处理FITS文件
import astropy.wcs as wcs # 导入astropy.wcs模块，用于处理世界坐标系统（WCS）
from astropy.table import QTable, Table # 从astropy.table模块导入QTable和Table，用于表格数据处理
import astropy.units as u # 导入astropy.units模块，用于处理物理单位
from astropy.visualization import make_lupton_rgb, SqrtStretch, LogStretch, hist # 导入可视化相关函数
from astropy.visualization.mpl_normalize import ImageNormalize # 导入图像归一化工具
from astropy.convolution import Gaussian2DKernel # 导入高斯2D卷积核
from astropy.stats import gaussian_fwhm_to_sigma # 导入将FWHM转换为sigma的函数
from astropy.coordinates import SkyCoord # 导入SkyCoord类，用于天球坐标处理
from photutils.background import Background2D, MedianBackground # 导入背景处理工具
from photutils.segmentation import (detect_sources, deblend_sources, SourceCatalog, # 导入源检测和分离工具
make_2dgaussian_kernel, SegmentationImage)
from photutils.utils import calc_total_error # 导入计算总误差的工具
from photutils.psf.matching import resize_psf, SplitCosineBellWindow, create_matching_kernel # 导入PSF匹配工具
from astropy.convolution import convolve, convolve_fft # 导入卷积函数
Matplotlib 设置用于绘图#
有两个版本
notebook– 提供交互式图形，但使整个笔记本的滚动变得稍微困难inline– 提供非交互式图形，以便更好地滚动整体内容
开发者笔记：
%matplotlib notebook 有时会生成过大的图形；需要重新运行单元格以使其恢复正常
import matplotlib.pyplot as plt # 导入绘图库，用于绘制图形
import matplotlib as mpl # 导入matplotlib库的其他功能
# 使用此版本进行非交互式绘图（更容易滚动笔记本）
%matplotlib inline # 设置为内联模式，图形将在笔记本中直接显示
# 如果您想要交互式图形，请使用此版本
# %matplotlib notebook # 设置为交互模式，允许与图形交互
import matplotlib.ticker as ticker # 导入刻度器，用于自定义坐标轴刻度
mpl_colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] # 获取当前绘图颜色循环
# 这些技巧是为了使图形的大小在内联和笔记本版本中保持一致
%config InlineBackend.print_figure_kwargs = {'bbox_inches': None} # 配置图形输出时的边界框设置
mpl.rcParams['savefig.dpi'] = 100 # 设置保存图形时的分辨率为100 DPI
mpl.rcParams['figure.dpi'] = 100 # 设置图形显示时的分辨率为100 DPI
import mplcursors # 可选：用于悬停在绘制的点上并显示ID号码
# 显示Python及导入库的版本信息
try:
import watermark # 导入watermark库，用于显示版本信息
%load_ext watermark # 加载watermark扩展
# %watermark -n -v -m -g -iv # 原始命令，显示更多信息（已注释）
%watermark -iv -v # 显示导入库的版本和Python的版本
except ImportError:
pass # 如果没有安装watermark库，则跳过
创建待加载和分析的图像列表#
所有数据和权重图像必须对齐到相同的像素网格。
（如有需要，请使用 reproject 来实现。）
input_file_url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/nircam_photometry/' # 输入文件的URL地址
filters = 'F090W F115W F150W F200W F277W F335M F356W F410M F444W'.split() # 定义滤光片列表并分割成单个滤光片
wavelengths = np.array([int(filt[1:4]) / 100 for filt in filters]) * u.um # 将滤光片的波长转换为微米，例如F115W = 1.15微米
# 数据图像 [e-/s]，与JWST管道图像的单位 [Myr/sr] 不同
image_files = {} # 创建一个字典来存储图像文件路径
for filt in filters: # 遍历每个滤光片
filename = f'jades_jwst_nircam_goods_s_crop_{filt}.fits' # 构建图像文件名
image_files[filt] = os.path.join(input_file_url, filename) # 将完整路径存入字典
# 权重图像（逆方差图；IVM）
weight_files = {} # 创建一个字典来存储权重文件路径
for filt in filters: # 遍历每个滤光片
filename = f'jades_jwst_nircam_goods_s_crop_{filt}_wht.fits' # 构建权重文件名
weight_files[filt] = os.path.join(input_file_url, filename) # 将完整路径存入字典
加载探测图像：F200W#
detection_filter = filt = 'F200W' # 设置检测滤光器为F200W
infile = image_files[filt] # 从图像文件中获取对应滤光器的文件名
hdu = fits.open(infile) # 打开图像文件，返回HDU对象
data = hdu[0].data # 获取HDU对象中第一个扩展的数据部分
imwcs = wcs.WCS(hdu[0].header, hdu) # 创建WCS对象以处理图像的世界坐标系统
weight = fits.open(weight_files[filt])[0].data # 打开权重文件并获取数据部分
报告图像大小和视场#
ny, nx = data.shape # 获取数据的行数和列数
# image_pixel_scale = np.abs(hdu[0].header['CD1_1']) * 3600 # 原始像素比例计算（已注释掉）
image_pixel_scale = wcs.utils.proj_plane_pixel_scales(imwcs)[0] # 计算图像的像素比例
image_pixel_scale *= imwcs.wcs.cunit[0].to('arcsec') # 将单位转换为角秒
outline = '%d x %d pixels' % (ny, nx) # 创建输出字符串，包含像素的行数和列数
outline += ' = %g" x %g"' % (ny * image_pixel_scale, nx * image_pixel_scale) # 添加图像的实际尺寸
outline += ' (%.2f" / pixel)' % image_pixel_scale # 添加每个像素的尺寸
print(outline) # 打印输出结果
创建彩色图像（可选）#
# 3 NIRCam短波长通道图像
r = fits.open(image_files['F200W'])[0].data # 打开F200W图像文件并获取数据，赋值给r
g = fits.open(image_files['F150W'])[0].data # 打开F150W图像文件并获取数据，赋值给g
b = fits.open(image_files['F090W'])[0].data # 打开F090W图像文件并获取数据，赋值给b
color_image_short_wavelength = make_lupton_rgb(r, g, b, Q=5, stretch=0.02) # 使用Lupton RGB方法生成短波长彩色图像
# 3 NIRCam长波长通道图像
r = fits.open(image_files['F444W'])[0].data # 打开F444W图像文件并获取数据
g = fits.open(image_files['F356W'])[0].data # 打开F356W图像文件并获取数据
b = fits.open(image_files['F277W'])[0].data # 打开F277W图像文件并获取数据
color_image_long_wavelength = make_lupton_rgb(r, g, b, Q=5, stretch=0.02) # 生成长波长的RGB图像
开发者注释：
sharex 和 sharey 在 projection=imwcs 下似乎存在一些兼容性问题
作为一种解决方法，我尝试过但未能关闭下面右侧图的 y 轴标签
fig = plt.figure(figsize=(9.5, 4)) # 创建一个9.5x4英寸的图形
ax_sw = fig.add_subplot(1, 2, 1, projection=imwcs) # 添加第一个子图，使用给定的WCS投影
# , sharex=True, sharey=True) # 注释掉的代码用于共享x和y轴
ax_sw.imshow(color_image_short_wavelength, origin='lower') # 显示短波长通道的图像，原点在下方
ax_sw.set_xlabel('Right Ascension') # 设置x轴标签为“赤经”
ax_sw.set_ylabel('Declination') # 设置y轴标签为“赤纬”
ax_sw.set_title('Short Wavelength Channel') # 设置子图标题为“短波长通道”
ax_lw = fig.add_subplot(1, 2, 2, projection=imwcs, sharex=ax_sw, sharey=ax_sw) # 添加第二个子图，使用相同的WCS投影，并共享x和y轴
ax_lw.imshow(color_image_long_wavelength, origin='lower') # 显示长波长通道的图像，原点在下方
ax_lw.set_xlabel('Right Ascension') # 设置x轴标签为“赤经”
ax_lw.set_title('Long Wavelength Channel') # 设置子图标题为“长波长通道”
ax_lw.set_ylabel('') # 设置y轴标签为空
# ax_lw.set(yticklabels=[]) # 注释掉的代码用于隐藏y轴刻度标签
# plt.subplots_adjust(left=0.15) # 注释掉的代码用于调整子图的左边距
print('Interactive zoom / pan controls both images simultaneously') # 打印信息，表示交互式缩放/平移同时控制两个图像
使用 astropy.photutils 检测源和去混叠#
https://photutils.readthedocs.io/en/latest/segmentation.html
# 在此定义所有检测和测量参数，以便对每个图像进行一致的测量
class Photutils_Catalog:
def __init__(self, filt, image_file=None, verbose=True):
# 初始化函数，接收滤波器、图像文件和详细输出参数
self.image_file = image_file or image_files[filt] # 如果未提供图像文件，则使用默认图像文件
self.hdu = fits.open(self.image_file) # 打开图像文件
self.data = self.hdu[0].data # 获取图像数据
self.imwcs = wcs.WCS(self.hdu[0].header, self.hdu) # 获取图像的WCS信息
self.zeropoint = self.hdu[0].header['ABMAG'] * u.ABmag # 获取零点星等
self.weight_file = weight_files[filt] # 获取权重文件名
self.weight = fits.open(weight_files[filt])[0].data # 打开权重文件并获取数据
if verbose:
# 如果verbose为真，输出滤波器、零点和图像文件信息
print(filt, ' zeropoint =', self.zeropoint)
print(self.image_file)
print(self.weight_file)
def measure_background_map(self, bkg_size=50, filter_size=3, verbose=True):
# 计算50x50像素单元中的sigma裁剪背景，然后在3x3单元上进行中值滤波
# 为获得最佳结果，图像应在两个维度上跨越整数个单元（这里，1000=20x50像素）
# https://photutils.readthedocs.io/en/stable/background.html
self.background_map = Background2D(self.data, bkg_size, filter_size=filter_size) # 计算背景图
def run_detect_sources(self, nsigma, npixels, smooth_fwhm=2, kernel_size=5,
deblend_levels=32, deblend_contrast=0.001, verbose=True):
# 设置检测阈值图为背景基线以上nsigma倍的RMS
detection_threshold = (nsigma * self.background_map.background_rms) + self.background_map.background
# 在检测之前，用高斯函数对数据进行卷积
smooth_kernel = make_2dgaussian_kernel(smooth_fwhm, size=kernel_size) # 创建高斯卷积核
convolved_data = convolve(self.data, smooth_kernel) # 对数据进行卷积处理
# 在经过卷积的图像中检测具有npixels个连接像素的源
# https://photutils.readthedocs.io/en/stable/segmentation.html
self.segm_detect = detect_sources(convolved_data, detection_threshold, npixels=npixels) # 检测源
# 去混叠：分离连接/重叠的源
# https://photutils.readthedocs.io/en/stable/segmentation.html#source-deblending
self.segm_deblend = deblend_sources(convolved_data, self.segm_detect, npixels=npixels,
nlevels=deblend_levels, contrast=deblend_contrast) # 去混叠处理
if verbose:
# 如果verbose为真，输出检测和去混叠的结果
output = 'Cataloged %d objects' % self.segm_deblend.nlabels # 输出去混叠后的对象数量
output += ', deblended from %d detections' % self.segm_detect.nlabels # 输出检测到的对象数量
median_threshold = (nsigma * self.background_map.background_rms_median) \
+ self.background_map.background_median # 计算中值阈值
output += ' with %d pixels above %g-sigma threshold' % (npixels, nsigma) # 输出阈值信息
# 背景输出等同于SourceExtractor报告的内容
output += '\n'
output += 'Background median %g' % self.background_map.background_median # 输出背景中值
output += ', RMS %g' % self.background_map.background_rms_median # 输出背景RMS
output += ', threshold median %g' % median_threshold # 输出阈值中值
print(output) # 打印输出信息
def measure_source_properties(self, exposure_time, local_background_width=24):
# "effective_gain" = 曝光时间图（从数据速率单位转换为计数）
# 权重 = 逆方差图 = 1 / sigma_background**2（没有源）
# https://photutils.readthedocs.io/en/latest/api/photutils.utils.calc_total_error.html
self.exposure_time_map = exposure_time * self.background_map.background_rms_median**2 * self.weight # 计算曝光时间图
# 数据RMS不确定性是背景RMS和源泊松不确定性的组合
background_rms = 1 / np.sqrt(self.weight) # 计算背景RMS
# 有效增益参数需要在每个地方都是正值（不能为零），因此添加小值1e-8
self.data_rms = calc_total_error(self.data, background_rms, self.exposure_time_map+1e-8) # 计算数据RMS
self.catalog = SourceCatalog(self.data-self.background_map.background, self.segm_deblend, wcs=self.imwcs,
error=self.data_rms, background=self.background_map.background,
localbkg_width=local_background_width) # 创建源目录
# 创建一个Photutils_Catalog对象，使用指定的检测滤波器
detection_catalog = Photutils_Catalog(detection_filter)
# 测量背景图，以便后续源检测使用
detection_catalog.measure_background_map()
# 运行源检测，设置检测阈值为2个标准差，最小连通像素数为5
detection_catalog.run_detect_sources(nsigma=2, npixels=5)
在检测图像中测量光度（及更多）#
https://photutils.readthedocs.io/en/latest/segmentation.html#centroids-photometry-and-morphological-properties
# exposure_time = hdu[0].header.get('EXPTIME') # seconds # 从图像头文件中获取曝光时间（秒）
exposure_time = 49500 # seconds; Adding by hand because it's missing from the image header # 手动添加曝光时间，因为图像头文件中缺失
detection_catalog.measure_source_properties(exposure_time) # 使用曝光时间测量源属性
# 保存检测到的物体的分割图
# 创建一个主HDU对象，包含分割图数据，并将WCS头信息添加到HDU中
segm_hdu = fits.PrimaryHDU(detection_catalog.segm_deblend.data.astype(np.uint32), header=imwcs.to_header())
# 将HDU写入FITS文件，文件名为'JADES_detections_segm.fits'，如果文件已存在则覆盖
segm_hdu.writeto('JADES_detections_segm.fits', overwrite=True)
显示检测结果与图像（可选）#
开发者注释：
交互式图表可以同时缩放/平移所有帧。（用户可以对单个星系进行缩放。）
fig, ax = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(9.5, 6)) # 创建一个2行3列的子图，x和y轴共享，设置图像大小
# fig, ax = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(9.5, 6), subplot_kw={'projection': imwcs})
# 如果需要使用RA和Dec坐标轴而不是像素坐标，可以添加: , subplot_kw={'projection': imwcs})
# 彩色图像
ax[0,0].imshow(color_image_short_wavelength, origin='lower') # 显示短波长的彩色图像，原点在下方
ax[0,0].set_title('Color Image') # 设置子图标题为“彩色图像”
# 数据
norm = ImageNormalize(stretch=SqrtStretch(), vmin=0) # 创建图像归一化对象，使用平方根拉伸，最小值为0
ax[0,1].imshow(data, origin='lower', cmap='Greys_r', norm=norm) # 显示检测图像，使用灰度反转色图
ax[0,1].set_title('Detection Image F200W') # 设置子图标题为“检测图像 F200W”
# 分割图
cmap = detection_catalog.segm_deblend.make_cmap(seed=12345) # 创建分割图的色图，使用随机种子
ax[0,2].imshow(detection_catalog.segm_deblend, origin='lower', cmap=cmap, interpolation='none') # 显示分割图，原点在下方
ax[0,2].set_title('Detections (Segmentation Image)') # 设置子图标题为“检测（分割图）”
# norm = ImageNormalize(stretch=SqrtStretch())
# ax[1,0].imshow(weight, origin='lower', cmap='Greys_r', vmin=0)
# ax[1,0].set_title('Weight Image F200W')
ax[1,0].imshow(detection_catalog.exposure_time_map, origin='lower', cmap='Greys_r', vmin=0) # 显示曝光时间图，原点在下方
ax[1,0].set_title('Exposure Time Map') # 设置子图标题为“曝光时间图”
# ax[1,0].imshow(background_map.background, origin='lower', cmap='Greys_r')
# ax[1,0].set_title('Background Pedestal')
# RMS
# norm = ImageNormalize()
vmin, vmax = 0.001, 0.0035 # 设置RMS图像的最小值和最大值
ax[1,1].imshow(detection_catalog.background_map.background_rms, origin='lower', vmin=vmin, vmax=vmax) # 显示背景RMS图，原点在下方
ax[1,1].set_title('Background RMS') # 设置子图标题为“背景RMS”
# 包括泊松噪声的总误差
ax[1,2].imshow(detection_catalog.data_rms, origin='lower', vmin=vmin, vmax=vmax) # 显示数据RMS图，原点在下方
ax[1,2].set_title('RMS + Poisson noise') # 设置子图标题为“RMS + 泊松噪声”
fig.tight_layout() # 调整子图布局，使其更紧凑
print('Interactive zoom / pan controls all images simultaneously') # 打印信息，表示所有图像的交互缩放/平移控制同时生效
查看 / 保存检测图像中的测量量（可选）#
仅保留某些量#
# 定义所需的列名
columns = 'label xcentroid ycentroid sky_centroid area semimajor_sigma semiminor_sigma'.split()
# 添加更多的列名
columns += 'ellipticity orientation gini'.split()
# 再添加一些列名
columns += 'kron_radius local_background segment_flux segment_fluxerr kron_flux kron_fluxerr'.split()
# 注释掉的代码，可能是为了添加更多的列名
# columns += 'source_sum source_sum_err kron_flux kron_fluxerr kron_radius local_background'.split()
# 将检测目录转换为表格格式，并选择指定的列
source_table = detection_catalog.catalog.to_table(columns=columns)
# 将 'semimajor_sigma' 列重命名为 'a'
source_table.rename_column('semimajor_sigma', 'a')
# 将 'semiminor_sigma' 列重命名为 'b'
source_table.rename_column('semiminor_sigma', 'b')
# 用 'ra' 和 'dec' 替换 'sky_centroid'
source_table['ra'] = source_table['sky_centroid'].ra.degree * u.degree # 获取天球坐标的右升角
source_table['dec'] = source_table['sky_centroid'].dec.degree * u.degree # 获取天球坐标的 declination
# 获取当前表格的所有列名
columns = list(source_table.columns)
# 重新排列列的顺序，保持 'ra' 和 'dec' 在前面
columns = columns[:3] + ['ra', 'dec'] + columns[4:-2]
# 根据新的列顺序更新表格
source_table = source_table[columns]
# 如果感兴趣，可以查看/保存输出，但其他滤光片的光度测量将很快添加！
# 将源表写入名为 'JADES_detections.ecsv' 的文件，允许覆盖
# source_table.write('JADES_detections.ecsv', overwrite=True)
# 将源表写入名为 'JADES_detections.cat' 的文件，格式为固定宽度的两行ASCII，使用空格作为分隔符，允许覆盖
# source_table.write('JADES_detections.cat', format='ascii.fixed_width_two_line', delimiter=' ', overwrite=True)
# 显示源表内容
# source_table
将测量的通量（数据单位）转换为星等#
有关单位的详细信息，请参见以下链接：
开发者注释：
流量单位可以自动转换为星等单位
但流量不确定性 *不能* 自动转换为星等不确定性
因此我编写了这个函数来实现转换
对于检测到的星等不确定性，应该使用 u.mag 而不是 u.ABmag
但对于未检测到的星等不确定性，则引用 u.ABmag 上限
它们需要保持一致，因此我们选择使用 u.ABmag
# 未检测到: mag = 99; magerr = 1-σ上限假设零通量
# 未观察到: mag = -99; magerr = 0
def fluxes2mags(flux, fluxerr):
nondet = flux < 0 # 如果通量为负，则为非检测
unobs = (fluxerr <= 0) + (fluxerr == np.inf) # 如果通量不确定性为负或无穷大，则为未观察到
mag = flux.to(u.ABmag) # 将通量转换为AB星等
magupperlimit = fluxerr.to(u.ABmag) # 如果通量=0，则为1-σ上限
mag = np.where(nondet, 99 * u.ABmag, mag) # 对于非检测，设置星等为99
mag = np.where(unobs, -99 * u.ABmag, mag) # 对于未观察到，设置星等为-99
magerr = 2.5 * np.log10(1 + fluxerr/flux) # 计算星等的不确定性
magerr = magerr.value * u.ABmag # 将不确定性转换为AB星等
magerr = np.where(nondet, magupperlimit, magerr) # 对于非检测，使用上限作为不确定性
magerr = np.where(unobs, 0 * u.ABmag, magerr) # 对于未观察到，设置不确定性为0
return mag, magerr # 返回星等和不确定性
# 包含我在astropy中找不到的特性:
# 对于非检测/未观察到，mag = 99 / -99
# 通量不确定性 -> 星等不确定性
使用在检测图像中定义的等亮度光圈进行多波段光度测量#
（类似于在双图像模式下运行SourceExtractor）
（尚未进行点扩散函数（PSF）校正）
# filters = 'F090W F200W F444W'.split() # 为测试而快速定义过滤器
for filt in filters: # 遍历每个过滤器
filter_catalog = Photutils_Catalog(filt) # 创建一个新的Photutils_Catalog对象，使用当前过滤器
filter_catalog.measure_background_map() # 测量背景图
# 在检测到的图像中测量该过滤器的天体光度
# 分割图将定义等光度光圈
filter_catalog.segm_deblend = detection_catalog.segm_deblend # 将分割去混叠信息赋值给当前过滤器的分割图
filter_catalog.measure_source_properties(exposure_time) # 测量源属性，使用曝光时间
# 将测量到的通量转换为nJy和AB星等
filter_table = filter_catalog.catalog.to_table() # 将过滤器目录转换为表格格式
source_table[filt+'_flux'] = flux = filter_table['segment_flux'] * filter_catalog.zeropoint.to(u.nJy) # 计算并存储通量
source_table[filt+'_fluxerr'] = fluxerr = filter_table['segment_fluxerr'] * filter_catalog.zeropoint.to(u.nJy) # 计算并存储通量误差
mag, magerr = fluxes2mags(flux, fluxerr) # 将通量和通量误差转换为星等和星等误差
source_table[filt+'_mag'] = mag # 存储星等
source_table[filt+'_magerr'] = magerr # 存储星等误差
光圈修正：等光度（isophotal）到总（Kron光圈）通量#
total_flux_table = deepcopy(source_table) # 复制源表到新表，以包含总光度校正
reference_flux_auto = total_flux_table['kron_flux'] # 获取自动克朗光流
reference_flux_iso = total_flux_table['segment_flux'] # 获取分段光流
kron_flux_corrections = reference_flux_auto / reference_flux_iso # 计算克朗光流校正因子
total_flux_table['total_flux_cor'] = kron_flux_corrections # 将校正因子添加到总光流表中
for filt in filters: # 遍历每个滤光片
total_flux_table[filt+'_flux'] *= kron_flux_corrections # 校正每个滤光片的光流
total_flux_table[filt+'_fluxerr'] *= kron_flux_corrections # 校正每个滤光片的光流误差
# total_flux_table[filt+'_mag'] = total_flux_table[filt+'_flux'].to(u.ABmag) # 不处理未检测到的情况
total_flux_table[filt+'_mag'] = fluxes2mags(total_flux_table[filt+'_flux'], total_flux_table[filt+'_fluxerr'])[0] # 将光流转换为星等
# 星等的不确定性 magerr 保持不变
重新格式化输出目录以提高可读性（可选）#
开发者注释：
'd' 格式与单位不兼容 (pix2)
作为解决方法，我将格式设置为 '.0f'（一个没有小数的浮点数）
ValueError: 无法解析格式字符串 "d" 对于列 "area" 中的条目 "101.0"
isophotal_table = deepcopy(total_flux_table) # 复制总通量表到新表，以包含PSF校正
old_columns = list(isophotal_table.columns) # 获取旧列名列表
# 重新排列列的顺序
i = old_columns.index('segment_flux') # 找到'segment_flux'的索引
j = old_columns.index(filters[0]+'_flux') # 找到第一个滤光片的通量索引
columns = old_columns[:i] # 获取检测目录（不包括source_sum和kron_flux）
columns += old_columns[-1:] # 添加总通量校正列
columns += old_columns[j:-1] # 添加所有滤光片的光度列
isophotal_table = isophotal_table[columns] # 根据新列顺序重新构建表
for column in columns: # 遍历每一列
isophotal_table[column].info.format = '.4f' # 设置每列的格式为小数点后4位
isophotal_table['ra'].info.format = '11.7f' # 设置'ra'列格式为小数点后7位，总宽度11位
isophotal_table['dec'].info.format = ' 11.7f' # 设置'dec'列格式为小数点后7位，总宽度11位
isophotal_table['label'].info.format = 'd' # 设置'label'列格式为整数
isophotal_table['area'].info.format = '.0f' # 设置'area'列格式为整数（不带小数），'d'会引发错误
无点扩散函数（PSF）校正的等光度测光#
我们建议在长波长通道中进行点扩散函数（PSF）校正的测光。但如果您有兴趣，您可以现在保存目录。
# 将等光度表写入ECSV格式文件，文件名为'JADES_isophotal_photometry.ecsv'，并允许覆盖
# isophotal_table.write('JADES_isophotal_photometry.ecsv', overwrite=True)
# 将等光度表写入ASCII格式文件，文件名为'JADES_isophotal_photometry.cat'，格式为固定宽度两行，分隔符为空格，并允许覆盖
# isophotal_table.write('JADES_isophotal_photometry.cat', format='ascii.fixed_width_two_line', delimiter=' ', overwrite=True)
# 输出等光度表的内容
# isophotal_table
PSF幅度校正#
颜色校正的执行方式如下所述：
https://www.stsci.edu/~dcoe/ColorPro/color
请注意，这与一种常见的方法不同，该方法是将每个图像降级到最宽的点扩散函数（PSF）。
光度校正仅在长波长图像中进行，波长大于2.4微米。F200W探测图像被卷积（模糊）到每个长波长图像的PSF。然后，我们根据在该孔径中损失的幅度来校正颜色：
PSF幅度校正 = 探测图像幅度 - 模糊探测图像幅度
在实践中，我们实际上校正的是通量：
PSF通量校正 = 探测图像通量 / 模糊探测图像通量
加载点扩散函数 (PSFs)#
本笔记本中使用的PSF文件是从JDox上的tar文件中提取的：
https://jwst-docs.stsci.edu/near-infrared-camera/nircam-predicted-performance/nircam-point-spread-functions#NIRCamPointSpreadFunctions-SimulatedNIRCamPSFs
PSFs_SW_filters（短波长通道）：https://stsci.box.com/s/s2lepxr9086gq4sogr3kwftp54n1c5vl
PSFs_LW_filters（长波长通道）：https://stsci.box.com/s/gzl7blxb1k3p4n66gs7jvt7xorfrotyb
每个FITS文件包含：
– hdu[0]：一个4倍过采样的PSF
– hdu[1]：在探测器像素尺度下的PSF（短波长和长波长通道分别为0.031”和0.063”）
本笔记本将使用后者：在探测器像素尺度下的PSF。我们认为前者对本笔记本没有优势。
PSF_inputs = {} # 初始化字典以存储PSF输入数据
PSF_images = {} # 初始化字典以存储PSF图像
detector_pixel_scales = {'SW': 0.031, 'LW': 0.063} # 定义探测器像素比例，SW和LW通道的像素比例
PSF_url = os.path.join(input_file_url, 'NIRCam_PSFs') # 构建PSF文件的URL路径
for i, filt in enumerate(filters): # 遍历每个滤光片及其索引
lam = wavelengths[i] # 获取对应的波长
if lam < 2.4 * u.um: # 判断波长是否小于2.4微米
channel = 'SW' # 短波长通道 < 2.4 微米
else:
channel = 'LW' # 长波长通道 > 2.4 微米
# 加载PSF文件
PSF_file = 'PSF_%scen_G5V_fov299px_ISIM41.fits' % filt # 构建PSF文件名
# PSF_file = os.path.join('NIRCam_PSFs_' + channel, PSF_file) # 原始路径构建（已注释）
PSF_file = os.path.join(PSF_url, PSF_file) # 组合成完整的PSF文件路径
print(PSF_file) # 打印PSF文件路径
PSF_hdu = fits.open(PSF_file) # 打开PSF文件
PSF_inputs[filt] = data = PSF_hdu[1].data # 将扩展[1]的数据存储到PSF输入字典中（像素比例未过采样）
ny, nx = data.shape # 获取PSF数据的形状
# 从探测器像素比例调整到图像像素比例（此处为0.03" / 像素）
detector_pixel_scale = detector_pixel_scales[channel] # 获取当前通道的探测器像素比例
ny_resize = ny * detector_pixel_scale / image_pixel_scale # 计算调整后的高度
ny_resize = np.round(ny_resize) # 四舍五入
ny_resize = int((ny_resize // 2) * 2 + 1) # 确保像素数为奇数，以确保PSF居中
PSF_pixel_scale = ny_resize / ny * image_pixel_scale # 计算PSF的像素比例
PSF_image = resize_psf(PSF_inputs[filt], PSF_pixel_scale, image_pixel_scale) # 调整PSF大小
r = (ny_resize - ny) // 2 # 计算需要裁剪的边界
PSF_images[filt] = PSF_image[r:-r, r:-r] # 裁剪PSF图像以与输入PSF大小相同
# 注意PSF现在未归一化，但将在卷积步骤中进行归一化
# print(filt, ny, ny_resize, PSF_image.shape, PSF_images[filt].shape, PSF_pixel_scale) # 打印调试信息（已注释）
np.sum(PSF_images[filt]) # 计算指定滤波器下PSF图像的总和
np.sum(PSF_image[r:-r, r:-r]) # 计算PSF_image数组中去掉边缘r个像素后的区域的所有元素的总和
# 显示点扩散函数（PSF）（可选）
fig, ax = plt.subplots(1, len(filters), figsize=(9.5, 1.5), sharex=True, sharey=True) # 创建子图，设置图形大小和共享坐标轴
r = 15 # PSF 将显示到半径 r（像素）
for i, filt in enumerate(filters): # 遍历每个滤光片
data = PSF_images[filt] # 获取当前滤光片的 PSF 数据
ny, nx = data.shape # 获取数据的形状（高度和宽度）
yc = ny // 2 # 计算图像中心的 y 坐标
xc = nx // 2 # 计算图像中心的 x 坐标
stamp = data[yc-r:yc+r+1, xc-r:xc+r+1] # 提取中心区域的 PSF 图像
norm = ImageNormalize(stretch=LogStretch()) # 为每个滤光片单独缩放
ax[i].imshow(stamp, cmap='Greys_r', norm=norm, origin='lower') # 显示 PSF 图像，使用灰度反转色图
ax[i].set_title(filt.upper()) # 设置子图标题为滤光片名称（大写）
ax[i].axis('off') # 关闭坐标轴
PSF 匹配#
https://photutils.readthedocs.io/en/stable/psf_matching.html
确定PSF卷积核#
PSF_kernels = {} # 创建一个空字典用于存储PSF核
reference_filter = 'F200W' # 定义参考滤波器为'F200W'
reference_PSF = PSF_images[reference_filter] # 获取参考滤波器的PSF图像
i_reference = filters.index(reference_filter) # 找到参考滤波器在滤波器列表中的索引
window = SplitCosineBellWindow(alpha=0.35, beta=0.3) # 创建一个分割余弦贝尔窗口
for filt in filters[i_reference+1:]: # 遍历参考滤波器之后的所有滤波器
PSF_kernels[filt] = create_matching_kernel(reference_PSF, PSF_images[filt], window=window) # 创建匹配核并存储在字典中
将F200W探测图像卷积到长波长点扩散函数（PSFs）#
# 打开参考图像文件
reference_image_hdu = fits.open(image_files[reference_filter])
# 获取参考图像的数据
reference_image_data = reference_image_hdu[0].data[:]
# 遍历输出滤波器，从参考滤波器的下一个开始
for output_filter in filters[i_reference+1:]:
# 定义输出图像的文件名
output_image = 'jades_convolved_%s_to_%s.fits' % (reference_filter, output_filter)
# 检查输出图像文件是否已存在
if os.path.exists(output_image):
print(output_image, 'EXISTS') # 如果存在，打印提示信息
else:
print(output_filter + '...') # 打印当前处理的输出滤波器
# 获取对应输出滤波器的PSF核
PSF_kernel = PSF_kernels[output_filter][yc-r:yc+r+1, xc-r:xc+r+1]
# 使用卷积进行图像处理，normalize_kernel=True表示归一化核
convolved_image = convolve(reference_image_data, PSF_kernel, normalize_kernel=True)
# 将卷积后的图像数据更新到参考图像中
reference_image_hdu[0].data = convolved_image
print('SAVING %s' % output_image) # 打印保存信息
# 将处理后的图像数据写入输出文件
reference_image_hdu.writeto(output_image)
卷积图像中的多波段光度测量#
# 测量并保存每个模糊图像中的通量
blurry_catalog = QTable() # 创建一个空的QTable以存储模糊图像的通量数据
for blurry_filter in filters[i_reference+1:]: # 遍历参考滤光片之后的所有模糊滤光片
image_file = 'jades_convolved_%s_to_%s.fits' % (reference_filter, blurry_filter) # 构建模糊图像文件名
filter_catalog = Photutils_Catalog(blurry_filter, image_file=image_file) # 创建Photutils_Catalog对象以处理当前滤光片的图像
filter_catalog.measure_background_map() # 测量背景图以便后续处理
# 在检测到的图像中测量该滤光片的光度
# 分割图将定义等光度孔径
filter_catalog.segm_deblend = detection_catalog.segm_deblend # 将检测到的图像的分割图应用于当前滤光片
filter_catalog.measure_source_properties(exposure_time) # 测量源属性，包括通量等
# 将测量的通量转换为nJy单位
filter_table = filter_catalog.catalog.to_table() # 将测量结果转换为表格格式
blurry_catalog[blurry_filter+'_flux'] = filter_table['segment_flux'] * filter_catalog.zeropoint.to(u.nJy) # 将通量转换为nJy并存储在模糊目录中
PSF 亮度修正#
抱歉，我无法直接访问外部链接。不过，如果您能提供该链接中的具体内容或文本，我将很乐意帮助您进行翻译和处理。请将需要翻译的内容粘贴在这里。
PSF_corrected_table = deepcopy(isophotal_table) # 深拷贝isophotal_table，创建PSF_corrected_table
reference_fluxes = PSF_corrected_table[reference_filter+'_flux'] # 获取参考滤光片的流量，det_flux_auto
for filt in filters[i_reference+1:]: # 遍历参考滤光片之后的所有滤光片
# 将等光度流量转换为总流量
blurry_total_fluxes = blurry_catalog[filt+'_flux'] * kron_flux_corrections # 计算模糊总流量
PSF_flux_corrections = reference_fluxes / blurry_total_fluxes # 计算PSF流量校正因子
PSF_corrected_table[filt+'_flux'] *= PSF_flux_corrections # 更新PSF校正后的流量
PSF_corrected_table[filt+'_fluxerr'] *= PSF_flux_corrections # 更新PSF校正后的流量误差
# PSF_corrected_table[filt+'_mag'] = PSF_corrected_fluxes.to(u.ABmag) # 不处理未探测到的情况
PSF_corrected_table[filt+'_mag'], PSF_corrected_table[filt+'_magerr'] = \ # 计算校正后的星等和误差
fluxes2mags(PSF_corrected_table[filt+'_flux'], PSF_corrected_table[filt+'_fluxerr'])
PSF_corrected_table[filt+'_PSF_flux_cor'] = PSF_flux_corrections # 将PSF流量校正因子存入表中
# 将PSF校正后的表格写入ECSV格式文件，文件名为'JADES_photometry.ecsv'，允许覆盖
PSF_corrected_table.write('JADES_photometry.ecsv', overwrite=True)
# 将PSF校正后的表格写入ASCII格式文件，文件名为'JADES_photometry.cat'，使用固定宽度的两行格式，分隔符为空格，允许覆盖
PSF_corrected_table.write('JADES_photometry.cat', format='ascii.fixed_width_two_line', delimiter=' ', overwrite=True)
# PSF_corrected_table
# 导入必要的库
import numpy as np # 导入NumPy库用于数值计算
import pandas as pd # 导入Pandas库用于数据处理
from astropy.io import fits # 导入Astropy库用于处理FITS文件
from jwst import datamodels # 导入JWST数据模型库
# 定义函数以读取FITS文件并返回数据
def read_fits_file(file_path):
# 使用Astropy读取FITS文件
with fits.open(file_path) as hdul:
# 返回数据部分
return hdul[1].data # 返回第二个HDU的数据
# 定义函数以处理PSF校正
def psf_correction(data):
# 在这里进行PSF校正的具体实现
# 假设我们进行简单的归一化处理
corrected_data = data / np.max(data) # 将数据归一化
return corrected_data # 返回校正后的数据
# 定义主函数
def main():
# 读取FITS文件
file_path = 'path/to/your/fits/file.fits' # 指定FITS文件路径
data = read_fits_file(file_path) # 调用函数读取数据
# 进行PSF校正
corrected_data = psf_correction(data) # 调用函数进行校正
# 将校正后的数据转换为DataFrame
df = pd.DataFrame(corrected_data) # 创建Pandas DataFrame
# 保存校正后的数据到CSV文件
df.to_csv('psf_corrected_data.csv', index=False) # 保存为CSV文件，不包含索引
# 如果该脚本是主程序，则执行主函数
if __name__ == '__main__':
main() # 调用主函数
# !cat JADES_photometry.cat # 显示名为JADES_photometry.cat的文件内容
开始新会话并分析结果#
只需运行上面的前几个代码块，包括导入库、定义文件列表和创建彩色图像。
加载目录和分割图#
# Catalog: ecsv格式保留单位，以便在Python笔记本中加载
output_catalog = QTable.read('JADES_photometry.ecsv') # 读取名为'JADES_photometry.ecsv'的文件并将其存储为QTable对象
# 重新构建滤光片列表
filters = [] # 初始化滤光片列表
for param in output_catalog.columns: # 遍历输出目录中的每一列
if param[-4:] == '_mag': # 检查列名是否以'_mag'结尾
filters.append(param[:-4]) # 将列名去掉'_mag'后添加到滤光片列表中
# 分割图像文件名
segmfile = 'JADES_detections_segm.fits'
# 打开FITS文件并读取数据
segm = fits.open(segmfile)[0].data
# 将数据转换为分割图像对象
segm = SegmentationImage(segm)
输入模拟 JADES JAGUAR 目录#
# full_catalog_file = 'JADES_SF_mock_r1_v1.1.fits.gz' # 原始数据文件在此演示中不可用
# 从302,515个模拟星系中裁剪出653个，适用于本演示的小图像区域
cropped_catalog_file = 'JADES_SF_mock_r1_v1.1_crop.fits.gz' # 裁剪后的星系目录文件名
input_catalog_file = os.path.join(input_file_url, cropped_catalog_file) # 生成完整的输入文件路径
simulated_catalog = Table.read(input_catalog_file) # 读取裁剪后的星系目录数据
将天体与 photutils 目录匹配#
# 将输入坐标转换为SkyCoord对象，使用RA和DEC列，并指定单位为度
input_coordinates = SkyCoord(ra=simulated_catalog['RA']*u.degree, dec=simulated_catalog['DEC']*u.degree)
# 如果在表中保存了output_catalog['sky_centroid']，可以使用它，但此示例保存了（ra, dec）坐标：
# 将输出目录中的RA和DEC列转换为SkyCoord对象
detected_coordinates = SkyCoord(ra=output_catalog['ra'], dec=output_catalog['dec'])
# 将photutils检测到的源与输入对象目录进行匹配：
# input_indices为匹配的索引，separation2d为二维分离距离，distance3d为三维距离
input_indices, separation2d, distance3d = detected_coordinates.match_to_catalog_sky(input_coordinates)
确定输入对象与匹配输出源之间的阈值最大距离
fig = plt.figure(figsize=(9.5, 4)) # 创建一个图形，设置大小为9.5x4英寸
plt.plot(np.sort(separation2d.to(u.arcsec)), zorder=10) # 绘制分离度的排序图，zorder设置图层顺序
separation_max = 0.036 * u.arcsec # 通过观察图形确定的最大分离度
plt.axhline(0, c='k', ls=':') # 绘制y=0的水平虚线，颜色为黑色
plt.axhline(separation_max.value, c='r', ls='--', label='"good" matches') # 绘制最大分离度的水平虚线，颜色为红色，线型为虚线，并添加标签
plt.title('Matches between output (detected) and input (simulated) catalogs') # 设置图形标题
plt.xlabel('Detected sources (sorted)') # 设置x轴标签
plt.ylabel('Separation (arcsec)') # 设置y轴标签
plt.legend() # 显示图例
# 判断哪些匹配的距离小于最大分离距离
good_matches = separation2d < separation_max
# 找到唯一的匹配项及其计数
unique_matches, index_counts = np.unique(input_indices[good_matches], return_counts=True)
# 打印匹配的数量和唯一匹配的数量
print('%d matches (%d unique) between input catalog (%d galaxies) and photutils catalog (%d detected sources)'
% (np.sum(good_matches), len(unique_matches), len(simulated_catalog), len(output_catalog)))
# 找到匹配次数大于1的多个匹配项
multiple_matches = unique_matches[index_counts > 1]
# 如果存在多个匹配项，则打印这些输入源的ID
if len(multiple_matches):
print('Input sources matched multiple times:', list(output_catalog['id'][multiple_matches]))
# 输入与输出颜色
# 定义两个滤波器
filt1, filt2 = 'F200W F444W'.split()
# 从模拟目录中获取输入光通量，使用滤波器1
input_flux1 = simulated_catalog['NRC_%s_fnu' % filt1][input_indices][good_matches]
# 从模拟目录中获取输入光通量，使用滤波器2
input_flux2 = simulated_catalog['NRC_%s_fnu' % filt2][input_indices][good_matches]
# 将输入光通量转换为AB星等，滤波器1
input_mag1 = (input_flux1 * u.nJy).to(u.ABmag).value
# 将输入光通量转换为AB星等，滤波器2
input_mag2 = (input_flux2 * u.nJy).to(u.ABmag).value
# 从输出目录中获取输出星等，滤波器1
output_mag1 = output_catalog[filt1 + '_mag'][good_matches].value
# 从输出目录中获取输出星等，滤波器2
output_mag2 = output_catalog[filt2 + '_mag'][good_matches].value
# 获取输出目录中的标签
output_ids = output_catalog['label'][good_matches].data.astype(int)
# 仅绘制检测到的对象
det1 = (0 < output_mag1) & (output_mag1 < 90) # 检测到的滤波器1的对象
det2 = (0 < output_mag2) & (output_mag2 < 90) # 检测到的滤波器2的对象
det = det1 * det2 # 仅保留同时在两个滤波器中检测到的对象
# 根据检测条件筛选输出星等
output_mag1 = output_mag1[det]
output_mag2 = output_mag2[det]
# 计算输入颜色
input_color = input_mag1 - input_mag2
# 计算输出颜色
output_color = output_mag1 - output_mag2
# 保存未校正的输出星等
output_mag1_uncor = output_mag1
# 获取滤波器2的PSF校正光通量
PSF_flux_cor = output_catalog[filt2+'_PSF_flux_cor']
# 计算PSF校正后的星等
PSF_mag_cor = -2.5 * np.log10(PSF_flux_cor)
# 计算未校正的输出星等
output_mag2_uncor = output_mag2 - PSF_mag_cor[good_matches][det].value
# 计算未校正的输出颜色
output_color_uncor = output_mag1_uncor - output_mag2_uncor
# 创建绘图
plt.figure(figsize=(8, 6))
# 绘制PSF校正后的数据点
plt.plot(input_mag1, output_color - input_color, 'bo', alpha=0.5, label='PSF corrected')
# 绘制未校正的数据点
plt.plot(input_mag1, output_color_uncor - input_color, 'r.', alpha=0.5, label='uncorrected')
# 添加水平线表示正确的颜色差
plt.axhline(0, c='g', label='correct')
# 设置x轴标签
plt.xlabel('Input ' + filt1 + ' (mag)')
# 设置y轴标签
plt.ylabel('Measured $-$ Input colors [%s $-$ %s]' % (filt1, filt2))
# 设置图表标题
plt.title('Accuracy of measured JADES photometry colors')
# 显示图例
plt.legend()
# 保存图像（注释掉以避免自动保存）
# plt.savefig('JADES_color_deviation_%s-%s.png' % (filt1, filt2))
绘制 F200W 与 F090W 亮度图并寻找掉落星系#
mag1 = output_catalog['F090W_mag'] # 获取F090W波段的星等数据
mag2 = output_catalog['F200W_mag'] # 获取F200W波段的星等数据
# 仅绘制F200W波段的检测结果
det2 = (0*u.ABmag < mag2) & (mag2 < 90*u.ABmag) # 筛选出F200W波段星等在0到90之间的检测结果
mag1 = mag1[det2] # 根据检测结果筛选F090W波段星等
mag2 = mag2[det2] # 根据检测结果筛选F200W波段星等
ids = output_catalog['label'][det2].data.astype(int) # 获取对应的标签并转换为整数类型
plt.figure(figsize=(8, 4)) # 创建一个8x4英寸的图形
plt.scatter(mag1, mag2, c=ids) # 绘制散点图，x轴为F090W星等，y轴为F200W星等，点的颜色由ids决定
plt.xlabel('F090W AB magnitude') # 设置x轴标签
plt.ylabel('F200W AB magnitude') # 设置y轴标签
plt.xlim(plt.xlim()[::-1]) # 反转x轴，使亮度较高的星等在右侧
plt.ylim(plt.ylim()[::-1]) # 反转y轴，使亮度较高的星等在顶部
mplcursors.cursor(hover=mplcursors.HoverMode.Transient) # 添加悬停光标功能
print('Hover cursor over data point to reveal magnitudes and catalog ID number') # 提示用户悬停以查看星等和目录ID
# 选择最亮的 F090W 脱落源
dropouts = output_catalog['F090W_mag'] > 90 * u.ABmag # 筛选出 F090W 星等大于 90 的脱落源
i_brightest_dropout = output_catalog[dropouts]['F200W_mag'].argmin() # 在脱落源中找到 F200W 星等最小的索引
output_id = output_catalog[dropouts][i_brightest_dropout]['label'] # 获取最亮脱落源的标签
output_id # 输出最亮脱落源的标签
# 从目录中选择具有特定ID的对象
output_index = segm.get_index(output_id) # 获取输出ID在分段中的索引
output_obj = output_catalog[output_index] # 根据索引从输出目录中获取对象
segmobj = segm.segments[segm.get_index(output_id)] # 获取与输出ID对应的分段对象
print(output_id, output_obj['F090W_mag'], output_obj['F200W_mag']) # 打印输出ID及其对应的F090W和F200W波段的星等
# 可选择通过位置选择一个对象
# x, y = 905, 276 # 定义对象的 x 和 y 坐标
# id = segm.data[y,x] # 从分段数据中获取指定位置的对象 ID
在所有滤光片中显示选定对象#
fig, ax = plt.subplots(2, len(filters)+1, figsize=(9.5, 3.5), sharex=True, sharey=True) # 创建一个2行(len(filters)+1)列的子图
ax[0, 0].imshow(color_image_short_wavelength[segmobj.slices], origin='lower', extent=segmobj.bbox.extent) # 显示短波长的彩色图像
ax[0, 0].set_title('Color') # 设置标题为'Color'
cmap = segm.make_cmap(seed=12345) # 创建一个颜色映射（cmap），使用随机种子12345
ax[1, 0].imshow(segm.data[segmobj.slices], origin='lower', extent=segmobj.bbox.extent, cmap=cmap, # 显示分段数据
interpolation='nearest') # 使用最近邻插值
ax[1, 0].set_title('Segment') # 设置标题为'Segment'
for i in range(1, len(filters)+1): # 遍历每个滤波器
filt = filters[i-1] # 获取当前滤波器
# Show data on top row
data = fits.open(image_files[filt])[0].data # 打开当前滤波器对应的图像文件
stamp = data[segmobj.slices] # 提取感兴趣区域的数据
norm = ImageNormalize(stretch=SqrtStretch()) # 创建一个归一化对象，使用平方根拉伸
ax[0, i].imshow(stamp, extent=segmobj.bbox.extent, cmap='Greys_r', norm=norm, origin='lower') # 显示图像数据
ax[0, i].set_title(filt.upper()) # 设置标题为滤波器名称（大写）
# Show weights on bottom row
weight = fits.open(weight_files[filt])[0].data # 打开当前滤波器对应的权重文件
stamp = weight[segmobj.slices] # 提取感兴趣区域的权重数据
# set black to zero weight (no exposure time / bad pixel)
ax[1, i].imshow(stamp, extent=segmobj.bbox.extent, vmin=0, cmap='Greys_r', origin='lower') # 显示权重数据，黑色表示权重为零
ax[0, 0].set_ylabel('Data') # 设置左侧y轴标签为'Data'
ax[1, 0].set_ylabel('Weight'); # 设置左侧y轴标签为'Weight'
绘制单个天体的光谱能量分布 (SED)#
# output_obj = output_catalog[index] # 已在上面完成
# 从模拟目录中获取输入对象，使用输入索引和输出索引
input_obj = simulated_catalog[input_indices][output_index]
# 从输入对象中提取每个滤波器的光通量，并将其存储为NumPy数组
input_fluxes = np.array([input_obj['NRC_%s_fnu' % filt] for filt in filters])
# 从输出对象中提取每个滤波器的光通量，并将其存储为NumPy数组
output_fluxes = np.array([output_obj[filt+'_flux'].value for filt in filters])
# 从输出对象中提取每个滤波器的光通量误差，并将其存储为NumPy数组
output_flux_errs = np.array([output_obj[filt+'_fluxerr'].value for filt in filters])
# 测量的通量未能恢复总输入通量
# 给定已知的模拟输入，确定恢复的通量占比
# 使用此比例将测量的SED缩放到输入SED以进行比较，下面将绘制结果
# 仅使用F200W缩放输出到输入通量（其他滤波器可能有不正确的通量校正）
filt = 'F200W' # 定义滤波器为F200W
# 计算输出通量与输入通量的比例
flux_scale_factor = output_obj[filt+'_flux'].value / input_obj['NRC_%s_fnu' % filt]
# 计算通量因子，即输入通量与输出通量的比值
flux_factor = 1 / flux_scale_factor # 输入通量 / 输出通量
# 打印通过photutils恢复的输入通量的百分比
print('%d%% of input flux recovered by photutils' % (100 * flux_scale_factor))
if 0:
# 将输出缩放到输入通量，考虑所有测量的通量和不确定性
# Benitez+00 方程 8 和 9
FOT = np.sum(input_fluxes * output_fluxes / output_flux_errs**2) # 计算加权输入和输出通量的总和
FTT = np.sum(input_fluxes**2 / output_flux_errs**2) # 计算加权输入通量平方的总和
flux_scale_factor = FOT / FTT # a_m: 观察到的 / 理论的 (输出 / 输入)
flux_factor = 1 / flux_scale_factor # 输入 / 输出
print('%d%% of input flux recovered by photutils' % (100 * flux_scale_factor)) # 打印通过 photutils 恢复的输入通量的百分比
PSF_flux_corrections = np.ones(len(filters)) # 初始化一个与filters长度相同的数组，所有元素为1
for i, filt in enumerate(filters): # 遍历filters列表，i为索引，filt为当前过滤器名称
PSFcor_column = filt+'_PSF_flux_cor' # 生成对应的PSF_flux_cor列名
if PSFcor_column in list(output_obj.columns): # 检查该列名是否在output_obj的列中
PSF_flux_corrections[i] = output_obj[PSFcor_column] # 如果存在，则将该列的值赋给PSF_flux_corrections数组的对应位置
PSF_flux_corrections # 输出最终的PSF_flux_corrections数组
开发者注释：
当通量扩展到负值时，自动次轴的幅度无法正常工作
所以我编写了自己的代码来实现这一功能
def add_magnitude_axis(ax, flux_units=u.nJy, plothuge=True):
# 获取当前y轴的上下限并转换为指定的流量单位
ylo, yhi = plt.ylim() * flux_units
# 将y轴的上限转换为AB星等
maghi = yhi.to(u.ABmag).value
# 计算y轴标签的第一个值，向上取整到小数点后一位
ytx1 = np.ceil(maghi * 10) / 10. # 例如：24.101 -> 24.2
# 计算y轴标签的第二个值，向上取整到整数
ytx2 = np.ceil(maghi) # 例如：24.1 -> 25
# 计算ytx1的小数部分
fpart = ytx1 - int(ytx1) # 例如：0.2
# 根据小数部分的值决定ytx1的内容
if np.isclose(fpart, 0) or np.isclose(fpart, 0.9):
ytx1 = [] # 如果小数部分接近0或0.9，则ytx1为空
elif np.isclose(fpart, 0.1) or np.isclose(fpart, 0.2):
# 如果小数部分接近0.1或0.2，生成多个y轴标签
ytx1 = np.array([ytx1, ytx2-0.7, ytx2-0.5, ytx2-0.3]) # 例如：24.1, 24.3, 24.5, 24.7
elif np.isclose(fpart, 0.3) or np.isclose(fpart, 0.4):
# 如果小数部分接近0.3或0.4，生成多个y轴标签
ytx1 = np.array([ytx1, ytx2-0.5, ytx2-0.3]) # 例如：24.3, 24.5, 24.7
elif np.isclose(fpart, 0.5):
# 如果小数部分接近0.5，生成两个y轴标签
ytx1 = np.array([ytx1, ytx2-0.3]) # 例如：24.5, 24.7
elif np.isclose(fpart, 0.6):
# 如果小数部分接近0.6，生成一个y轴标签
ytx1 = np.array([ytx1, ytx2-0.2]) # 例如：24.6, 24.8
# 如果ytx1是浮点数，则将其转换为数组
if isinstance(ytx1, float):
ytx1 = np.array([ytx1])
# 根据plothuge的值决定ytx3的内容
if plothuge:
ytx3 = ytx2 + np.array([0, 0.2, 0.5, 1, 2]) # 大范围的y轴标签
else:
ytx3 = ytx2 + np.array([0, 0.2, 0.5, 1, 1.5, 2, 3]) # 小范围的y轴标签
# 将ytx2转换为数组
ytx2 = np.array([ytx2])
# 合并ytx1和ytx3
ytx = np.concatenate([ytx1, ytx3])
# 将y轴标签格式化为字符串
yts = ['%g' % mag for mag in ytx]
# 将y轴标签从AB星等转换为指定的流量单位
ytx = (ytx * u.ABmag).to(flux_units).value
# 创建一个共享y轴的第二个坐标轴
ax2 = ax.twinx()
# 设置y轴标签位置
ax.yaxis.set_label_position('left')
ax2.yaxis.set_label_position('right')
# 设置y轴刻度方向
ax.yaxis.tick_left()
ax2.yaxis.tick_right()
# 设置第二个坐标轴的y轴刻度和标签
ax2.set_yticks(ytx)
ax2.set_yticklabels(yts)
ax2.set_ylabel('Magnitude (AB)') # 设置y轴标签为“Magnitude (AB)”
ax2.set_ylim(ylo.value, yhi.value) # 设置y轴的上下限
ax2.set_zorder(-100) # 设置z轴顺序，确保交互光标输出左侧坐标轴ax
开发者注释：
errorbar 也不支持单位
fig, ax = plt.subplots(figsize=(8, 6)) # 创建一个8x6英寸的图形和坐标轴
plt.plot(wavelengths, input_fluxes, 'o-', label='Input fluxes', zorder=10) # 绘制输入通量的折线图，点标记为圆圈
label = 'Measured PSF-corrected fluxes $\\times$ %.1f' % flux_factor # 创建标签，表示经过PSF校正的通量
plt.errorbar(wavelengths.value, output_fluxes * flux_factor, output_flux_errs * flux_factor, # 绘制经过PSF校正的通量的误差条
ms=8, marker='s', mfc=mpl_colors[1], c='k', lw=3, alpha=0.5, ls='none', label=label) # 设置标记样式和颜色
label = 'Measured uncorrected fluxes $\\times$ %.1f' % flux_factor # 创建标签，表示未经校正的通量
plt.errorbar(wavelengths.value, output_fluxes * flux_factor / PSF_flux_corrections, output_flux_errs * flux_factor, # 绘制未经校正的通量的误差条
ms=8, marker='s', mfc='r', c='k', lw=3, alpha=0.5, ls='none', label=label, zorder=-10) # 设置标记样式和颜色
plt.legend() # 显示图例
plt.axhline(0, c='k', ls=':') # 绘制一条水平的虚线，y=0
plt.xlim(0, 5) # 设置x轴的范围为0到5
plt.xlabel('Wavelength ($\\mu$m)') # 设置x轴标签为波长
plt.ylabel('Flux (nJy)') # 设置y轴标签为通量
add_magnitude_axis(ax) # 添加幅度轴
plt.savefig('JADES photutils SED linear.png') # 保存图形为PNG文件
开发者笔记：
理想情况下，辅助轴应该能够知道如何在单位之间转换。
作为一种解决方法，我编写了函数并将其输入。
即便如此，这仅在通量（flux）为正值时有效。
将所有通量裁剪为正值也无法解决问题。
我希望能够自动缩放y轴的限制，仅针对某些绘制的数据，类似于：
https://stackoverflow.com/questions/7386872/make-matplotlib-autoscaling-ignore-some-of-the-plots
但那个旧的解决方案并没有奏效。
所以目前，我只是硬编码了一个适用于该对象的范围：y=10-70。
fig, ax = plt.subplots(figsize=(8, 6)) # 创建一个8x6英寸的图形和坐标轴
if 0: # 这段代码没有生效
ax.autoscale(True) # 自动缩放坐标轴
detections = output_fluxes > 0 # 检测输出通量大于0的值
ax.plot(wavelengths[detections], output_fluxes[detections] * flux_factor, visible=False) # 绘制可见的输出通量
ax.autoscale(False) # 关闭自动缩放
plt.plot(wavelengths, input_fluxes, 'o-', label='Input fluxes', zorder=10, scaley=False) # 绘制输入通量
label = 'Measured PSF-corrected fluxes $\\times$ %.1f' % flux_factor # 创建标签，显示PSF校正后的通量
plt.errorbar(wavelengths.value, output_fluxes * flux_factor, output_flux_errs * flux_factor, # 绘制带误差条的输出通量
ms=8, marker='s', mfc=mpl_colors[1], c='k', lw=3, alpha=0.5, ls='none', label=label)
label = 'Measured uncorrected fluxes $\\times$ %.1f' % flux_factor # 创建标签，显示未校正的通量
plt.errorbar(wavelengths.value, output_fluxes * flux_factor / PSF_flux_corrections, output_flux_errs * flux_factor, # 绘制带误差条的未校正通量
ms=8, marker='s', mfc='r', c='k', lw=3, alpha=0.5, ls='none', label=label, zorder=-10)
plt.legend() # 显示图例
plt.axhline(0, c='k', ls=':') # 绘制y=0的水平线
plt.xlim(0, 5) # 设置x轴范围
plt.ylim(15, 70) # 设置y轴范围
plt.xlabel('Wavelength ($\\mu$m)') # 设置x轴标签
plt.ylabel('Flux (nJy)') # 设置y轴标签
plt.semilogy() # 使用对数y轴
ax.yaxis.set_major_formatter(ticker.FormatStrFormatter("%g")) # 设置y轴主刻度格式
ax.yaxis.set_minor_formatter(ticker.FormatStrFormatter("%g")) # 设置y轴次刻度格式
# 添加AB星等作为右侧的次坐标轴
# （注意：如果任何通量为负，则此功能会失效）
# https://matplotlib.org/gallery/subplots_axes_and_figures/secondary_axis.html#sphx-glr-gallery-subplots-axes-and-figures-secondary-axis-py
def AB2nJy(mAB): # 定义AB星等到nJy的转换函数
m = mAB * u.ABmag # 将AB星等转换为量纲
f = m.to(u.nJy) # 转换为nJy
return f.value # 返回nJy值
def nJy2AB(F_nJy): # 定义nJy到AB星等的转换函数
f = F_nJy * u.nJy # 将nJy值转换为量纲
m = f.to(u.ABmag) # 转换为AB星等
return m.value # 返回AB星等值
# secondary_axis = add_magnitude_axis(ax, flux_units) # 这行代码被注释掉
secax = ax.secondary_yaxis('right', functions=(nJy2AB, AB2nJy)) # 创建右侧次坐标轴
secax.set_ylabel('magnitude (AB)') # 设置次坐标轴的y轴标签
secax.yaxis.set_major_formatter(ticker.FormatStrFormatter("%g")) # 设置次坐标轴主刻度格式
secax.yaxis.set_minor_formatter(ticker.FormatStrFormatter("%g")) # 设置次坐标轴次刻度格式
plt.title('JADES SED PSF-corrected') # 设置图形标题
# plt.savefig('JADES photutils SED.png') # 保存图形（这行代码被注释掉）