扩展孔径光度测量

扩展孔径光度测量#

用例： 在一个区域内测量星系光度。相关于 JDox 科学用例 #22。

数据： WST 模拟的 NIRCam 图像来自 JADES JWST GTO 外星系空白场。

（Williams et al. 2018）https://ui.adsabs.harvard.edu/abs/2018ApJS..236…33W。

工具： photutils, matplotlib, scipy, scikit。

跨仪器： 可能涉及 MIRI。

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

引言#

本笔记本使用 photutils 检测 NIRCam 深成像中的物体/星系。首先在 F200W 图像中进行检测，然后在所有 9 个滤光片（F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, F444W）中获得等光度光度测量。目录被重新加载，并对完整目录和单个星系进行一些简单分析。

该笔记本仅分析 JADES 模拟的中央 1000 x 1000 像素（30” x 30”）。这些切片已获得 STScI 的授权，并由作者（Williams et al.）提供。

注意： 光度测量是经过孔径匹配的，但未进行 PSF 校正。为了获得更准确的颜色测量，应该实施 PSF 校正，因为波长范围很大（因此 PSF FWHM 跨越 >4 的因子）。

注意： 模拟的 JADES 图像的单位（e-/s）与 JWST 管道产品（MJy/sr）不同。

注意： 缺少曝光图，但计算通量不确定性是必需的。

待办事项#

  • PSF（点扩散函数）校正

  • 检查光度测量的准确性与模拟的JADES（James Webb Space Telescope Advanced Deep Extragalactic Survey）目录的对比

  • 需要曝光图以用于误差计算

  • AB（天文带）磁单位无法写入ecsv文件（即将发布astropy更新）

  • 带文本标签的图表看起来很糟糕（我希望光标悬停时能显示ID号码）

  • 修复图表的次要轴：磁量与通量

  • requirements.txt文件——但我不知道需要哪些版本

  • Robel的其他评论：https://github.com/spacetelescope/dat_pyinthesky/pull/82#pullrequestreview-355206337

下载WEBBPSF数据#

import os  # 导入操作系统模块

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

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

boxlink = 'https://stsci.box.com/shared/static/34o0keicz2iujyilg4uz617va46ks6u9.gz'  # 设置文件下载链接

boxfile = './webbpsf-data/webbpsf-data-1.0.0.tar.gz'  # 设置下载后保存的文件路径

webbpsf_folder = './webbpsf-data'  # 设置存放解压文件的文件夹路径

# 检查指定路径是否存在

isExist = os.path.exists(webbpsf_folder)

if not isExist:  # 如果路径不存在

  # 创建一个新目录，因为它不存在 

  os.makedirs(webbpsf_folder)  # 创建目录
    

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

gzf = tarfile.open(boxfile)  # 打开下载的tar.gz文件

gzf.extractall(webbpsf_folder, filter='data')  # 解压所有文件到指定文件夹，使用filter='data'跳过绝对路径或相对路径

导入包#

import os  # 导入操作系统模块，用于文件和目录操作

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

from astropy.convolution import convolve  # 从Astropy库导入卷积函数
from astropy.io import fits  # 从Astropy库导入FITS文件处理模块
from astropy.stats import gaussian_fwhm_to_sigma  # 导入高斯全宽半高转换函数
from astropy.table import QTable  # 从Astropy库导入QTable类，用于表格数据处理
import astropy.units as u  # 导入Astropy单位模块
from astropy.visualization import make_lupton_rgb, SqrtStretch, ImageNormalize, simple_norm  # 导入图像可视化相关函数
import astropy.wcs as wcs  # 导入Astropy的WCS模块，用于天文坐标转换

from photutils.background import Background2D, MedianBackground  # 导入背景处理模块
from photutils.segmentation import (detect_sources, deblend_sources, SourceCatalog,  # 导入源检测和分离模块
                                    make_2dgaussian_kernel, SegmentationImage)  # 导入2D高斯核和分割图像相关函数
from photutils.utils import calc_total_error  # 导入计算总误差的函数

Matplotlib 绘图设置#

有两个版本

  • notebook – 提供交互式图形，但会使整个笔记本的滚动变得稍微困难

  • inline – 提供非交互式图形，以便更好地滚动整体内容

import matplotlib.pyplot as plt  # 导入绘图库 matplotlib 的 pyplot 模块

import matplotlib as mpl  # 导入 matplotlib 库

# Use this version for non-interactive plots (easier scrolling of the notebook)
%matplotlib inline  # 设置为非交互式绘图模式，方便在笔记本中滚动查看

# Use this version if you want interactive plots
# %matplotlib notebook  # 如果需要交互式绘图，可以取消注释这一行

# These gymnastics are needed to make the sizes of the figures
# be the same in both the inline and notebook versions
%config InlineBackend.print_figure_kwargs = {'bbox_inches': None}  # 配置图形输出的边界框参数，使得图形在两种模式下大小一致

mpl.rcParams['savefig.dpi'] = 80  # 设置保存图形时的分辨率为 80 dpi
mpl.rcParams['figure.dpi'] = 80  # 设置图形显示时的分辨率为 80 dpi

创建待加载和分析的图像列表#

# 基础URL，用于获取JWST数据分析工具的NIRCam光度数据
baseurl = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/nircam_photometry/'

# 定义所需的滤光片
filters = 'F090W F115W F150W F200W F277W F335M F356W F410M F444W'.split()

# 数据图像 [电子/秒]
imagefiles = {}  # 创建一个空字典以存储图像文件路径

# 遍历每个滤光片
for filt in filters:
    # 根据滤光片名称生成文件名
    filename = f'jades_jwst_nircam_goods_s_crop_{filt}.fits'
    # 将完整的文件路径存储到字典中
    imagefiles[filt] = os.path.join(baseurl, filename)

# 权重图像（逆方差图；IVM）
weightfiles = {}  # 创建一个空字典以存储权重文件路径

# 遍历每个滤光片
for filt in filters:
    # 根据滤光片名称生成权重文件名
    filename = f'jades_jwst_nircam_goods_s_crop_{filt}_wht.fits'
    # 将完整的权重文件路径存储到字典中
    weightfiles[filt] = os.path.join(baseurl, filename)

加载探测图像：F200W#

filt = 'F200W'  # 定义滤波器类型为'F200W'

infile = imagefiles[filt]  # 从图像文件字典中获取对应滤波器的输入文件

hdu = fits.open(infile)  # 打开输入文件，返回一个HDU列表

data = hdu[0].data  # 获取HDU列表中第一个HDU的图像数据

imwcs = wcs.WCS(hdu[0].header, hdu)  # 创建WCS对象以处理图像的世界坐标系统

weight = fits.open(weightfiles[filt])[0].data  # 打开权重文件并获取其数据

报告图像大小和视场#

ny, nx = data.shape  # 获取数据的行数和列数

pixscale = wcs.utils.proj_plane_pixel_scales(imwcs)[0]  # 计算投影平面像素尺度

pixscale *= imwcs.wcs.cunit[0].to('arcsec')  # 将像素尺度转换为角秒

outline = '%d x %d pixels' % (ny, nx)  # 创建输出字符串，包含像素的行列数

outline += ' = %g" x %g"' % (ny * pixscale, nx * pixscale)  # 添加输出字符串，包含实际尺寸

outline += ' (%.2f" / pixel)' % pixscale  # 添加输出字符串，包含每个像素的角度

print(outline)  # 打印最终的输出字符串

创建彩色图像（可选）#

# 3 NIRCam短波长通道图像

r = fits.open(imagefiles['F200W'])[0].data  # 打开F200W图像文件并获取数据

g = fits.open(imagefiles['F150W'])[0].data  # 打开F150W图像文件并获取数据

b = fits.open(imagefiles['F090W'])[0].data  # 打开F090W图像文件并获取数据

rgb = make_lupton_rgb(r, g, b, Q=5, stretch=0.02)  # 创建RGB图像，使用Lupton算法，Q和stretch参数控制图像显示效果

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

ax = plt.subplot(projection=imwcs)  # 创建一个子图并设置投影为imwcs

plt.imshow(rgb, origin='lower')  # 显示RGB图像，设置原点在左下角

plt.xlabel('Right Ascension')  # 设置x轴标签为“右升角”

plt.ylabel('Declination')  # 设置y轴标签为“赤纬”

fig.tight_layout()  # 调整图形布局以避免重叠

plt.subplots_adjust(left=0.15)  # 调整子图左侧的边距

使用 photutils 检测源和去混叠#

https://photutils.readthedocs.io/en/latest/segmentation.html

# For detection, requiring 5 connected pixels 2-sigma above background
# 检测要求：5个相连的像素点，信号强度需高于背景的2倍标准差

# Measure background and set detection threshold
# 测量背景并设置检测阈值

bkg_estimator = MedianBackground()  # 创建中位数背景估计器

bkg = Background2D(data, (50, 50), filter_size=(3, 3), bkg_estimator=bkg_estimator)  # 计算二维背景

threshold = bkg.background + (2. * bkg.background_rms)  # 设置阈值为背景加上2倍的背景标准差

# Before detection, smooth image with Gaussian FWHM = 3 pixels
# 在检测之前，用FWHM为3像素的高斯函数平滑图像

smooth_kernel = make_2dgaussian_kernel(3.0, size=3)  # 创建高斯平滑核

convolved_data = convolve(data, smooth_kernel)  # 对数据进行卷积处理

# Detect and deblend
# 检测并进行去混叠

segm_detect = detect_sources(convolved_data, threshold, npixels=5)  # 检测源，要求5个相连像素

segm_deblend = deblend_sources(convolved_data, segm_detect, npixels=5, nlevels=32, contrast=0.001)  # 去混叠处理

# Save segmentation map of detected objects
# 保存检测到的物体的分割图

segm_hdu = fits.PrimaryHDU(segm_deblend.data.astype(np.uint32), header=imwcs.to_header())  # 创建分割图的HDU

segm_hdu.writeto('JADES_detections_segm.fits', overwrite=True)  # 将分割图写入FITS文件，覆盖同名文件

在检测图像中测量光度（及更多）#

https://photutils.readthedocs.io/en/latest/segmentation.html#centroids-photometry-and-morphological-properties

#error = bkg.background_rms  # 背景噪声的均方根值

# Input weight should be exposure map. Fudging for now.  # 输入权重应为曝光图。暂时使用简单处理。

error = calc_total_error(data, bkg.background_rms, weight/500)  # 计算总误差，使用数据、背景均方根和权重

#cat = source_properties(data-bkg.background, segm_deblend, wcs=imwcs, background=bkg.background, error=error)  # 原代码：获取源属性

cat = SourceCatalog(data-bkg.background, segm_deblend, wcs=imwcs, background=bkg.background, error=error)  # 创建源目录，包含数据、分割图、WCS、背景和误差

显示检测结果与图像（可选）#

fig, ax = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(9.5, 6))  # 创建一个2行3列的子图，设置共享x和y轴，并指定图形大小

# For RA,Dec axes instead of pixels, add: , subplot_kw={'projection': imwcs})
# 如果需要使用RA和Dec坐标轴而不是像素，可以添加: , subplot_kw={'projection': imwcs})

# Color image
ax[0, 0].imshow(rgb, origin='lower')  # 显示RGB彩色图像，原点在下方
ax[0, 0].set_title('Color Image')  # 设置子图标题为“彩色图像”

# Data
norm = simple_norm(data, 'sqrt', percent=99.)  # 使用平方根归一化数据，99%分位数
ax[0, 1].imshow(data, origin='lower', cmap='Greys_r', norm=norm)  # 显示数据图像，使用灰度反转色图
ax[0, 1].set_title('Detection Image F200W')  # 设置子图标题为“检测图像 F200W”

# Segmentation map
cmap = segm_deblend.make_cmap(seed=12345)  # 创建分割图的颜色映射，使用随机种子
ax[0, 2].imshow(segm_deblend, origin='lower', cmap=cmap, interpolation='nearest')  # 显示分割图像，使用最近邻插值
ax[0, 2].set_title('Detections (Segmentation Image)')  # 设置子图标题为“检测（分割图像）”

# Weight
ax[1, 0].imshow(weight, origin='lower', cmap='Greys_r', vmin=0)  # 显示权重图像，使用灰度反转色图，最小值为0
ax[1, 0].set_title('Weight Image F200W')  # 设置子图标题为“权重图像 F200W”

# RMS
ax[1, 1].imshow(bkg.background_rms, origin='lower', norm=None)  # 显示背景RMS图像
ax[1, 1].set_title('Background RMS')  # 设置子图标题为“背景RMS”

# Total error including Poisson noise
norm = simple_norm(error, 'sqrt', percent=99.)  # 使用平方根归一化误差数据，99%分位数
ax[1, 2].imshow(error, origin='lower', norm=norm)  # 显示误差图像
ax[1, 2].set_title('RMS + Poisson noise')  # 设置子图标题为“RMS + 泊松噪声”

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

查看探测图像中的所有测量量（可选）#

# 将cat对象转换为表格格式
table = cat.to_table()  # 将cat对象转换为表格

仅保留某些量#

# 定义要提取的列名
columns = 'label xcentroid ycentroid sky_centroid area semimajor_sigma semiminor_sigma ellipticity orientation gini'.split()

# 从cat对象中提取指定列并转换为表格格式
tbl = cat.to_table(columns=columns)

# 将列名'semimajor_sigma'重命名为'a'
tbl.rename_column('semimajor_sigma', 'a')

# 将列名'semiminor_sigma'重命名为'b'
tbl.rename_column('semiminor_sigma', 'b')

看起来您提供的内容不完整，似乎是一个表格的标题或标识符。如果您有特定的Python代码或JWST数据处理的代码需要我添加中文注释，请提供完整的代码片段。我将很乐意帮助您添加注释并保持代码结构不变。

将测量的通量（数据单位）转换为星等#

https://docs.astropy.org/en/stable/units/

https://docs.astropy.org/en/stable/units/equivalencies.html#photometric-zero-point-equivalency

https://docs.astropy.org/en/stable/units/logarithmic_units.html#logarithmic-units

# 未探测：mag = 99；magerr = 1-sigma 上限假设零通量

# 未观察：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-sigma上限

    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）

filters = 'F090W F115W F150W F200W F277W F335M F356W F410M F444W'.split()  # 定义过滤器列表

for filt in filters:  # 遍历每个过滤器

    infile = imagefiles[filt]  # 获取对应过滤器的输入文件

    print(filt)  # 打印当前过滤器名称

    print(infile)  # 打印当前输入文件路径

    print(weightfiles[filt])  # 打印当前过滤器的权重文件路径

    hdu = fits.open(infile)  # 打开输入文件

    data = hdu[0].data  # 获取数据部分

    zp = hdu[0].header['ABMAG'] * u.ABmag  # 获取零点（AB magnitude）

    weight = fits.open(weightfiles[filt])[0].data  # 打开权重文件并获取数据

    # Measure background  # 测量背景

    bkg = Background2D(data, (50, 50), filter_size=(3, 3), bkg_estimator=bkg_estimator)  # 计算二维背景

    #error = bkg.background_rms  # 计算背景的均方根误差（注释掉）

    error = calc_total_error(data, bkg.background_rms, weight/500)  # 计算总误差

    # Measure properties in each image of previously detected objects  # 测量先前检测到的物体在每幅图像中的属性

    filtcat = SourceCatalog(data-bkg.background, segm_deblend, wcs=imwcs, background=bkg.background, error=error)  # 创建源目录

    # Convert measured fluxes to fluxes in nJy and to AB magnitudes  # 将测量的通量转换为nJy和AB星等

    filttbl = filtcat.to_table()  # 将源目录转换为表格

    tbl[filt+'_flux']    = flux    = filttbl['segment_flux']     * zp.to(u.nJy)  # 计算并存储通量

    tbl[filt+'_fluxerr'] = fluxerr = filttbl['segment_fluxerr'] * zp.to(u.nJy)  # 计算并存储通量误差

    mag, magerr = fluxes2mags(flux, fluxerr)  # 将通量和通量误差转换为星等和星等误差

    #mag = mag * u.ABmag  # incompatible with file writing  # 注释掉的行，表示与文件写入不兼容

    tbl[filt+'_mag']    = mag.value  # 存储星等

    tbl[filt+'_magerr'] = magerr.value  # 存储星等误差

查看完整结果（可选）#

看起来您提供的信息不完整，您提到的“tbl”似乎是一个表格或数据结构的名称，但没有具体的代码或上下文。如果您能提供更多的代码或具体的内容，我将能够更好地帮助您添加中文注释并保持代码的结构和功能。请提供相关的Python代码或数据处理的示例。谢谢！

将光度数据保存为输出目录#

tbl.write('JADESphotometry.ecsv', overwrite=True)  # 将表格数据写入名为'JADESphotometry.ecsv'的文件中，设置overwrite=True以允许覆盖现有文件

!head -175 JADESphotometry.ecsv  # 显示前175行
 

这段代码使用了shell命令`head`来显示文件`JADESphotometry.ecsv`的前175行。

重新格式化输出目录以提高可读性（可选）#

# 移除面积的单位（像素）
tbl['area'] = tbl['area'].value.astype(int)  # 将'area'列的值转换为整数类型

# 用ra和dec替换sky_centroid
tbl['ra'] = tbl['sky_centroid'].ra.degree  # 提取sky_centroid的RA值并转换为度
tbl['dec'] = tbl['sky_centroid'].dec.degree  # 提取sky_centroid的Dec值并转换为度

# 获取当前表格的列名
columns = list(tbl.columns)  # 将表格的列名转换为列表

# 重新排列列名，将ra和dec插入到适当的位置
columns = columns[:3] + ['ra', 'dec'] + columns[4:-2]  # 保留前3列，插入ra和dec，然后保留剩余列

# 根据新的列顺序重新构建表格
tbl = tbl[columns]  # 重新排列表格的列

for column in columns:  # 遍历所有列
    tbl[column].info.format = '.4f'  # 设置每列的格式为浮点数，保留四位小数

tbl['ra'].info.format = '11.7f'  # 设置'ra'列的格式为浮点数，总宽度11，保留7位小数
tbl['dec'].info.format = '11.7f'  # 设置'dec'列的格式为浮点数，总宽度11，保留7位小数

tbl['label'].info.format = 'd'  # 设置'label'列的格式为整数
tbl['area'].info.format = 'd'  # 设置'area'列的格式为整数

tbl.write('JADESphotometry.cat', format='ascii.fixed_width_two_line', delimiter=' ', overwrite=True)  # 将表格数据写入文件，格式为固定宽度两行，使用空格作为分隔符，允许覆盖已有文件

!head -10 JADESphotometry.cat  # 显示文件 JADESphotometry.cat 的前10行

开始新会话并分析结果#

加载目录和分割图#

# Catalog: ecsv格式保留单位以便在Python笔记本中加载
tbl = QTable.read('JADESphotometry.ecsv')  # 读取JADESphotometry.ecsv文件并存储为表格

# 重新构建滤光片列表
filters = []  # 初始化滤光片列表

for param in tbl.columns:  # 遍历表格中的每一列
    if param[-4:] == '_mag':  # 检查列名是否以'_mag'结尾
        filters.append(param[:-4])  # 将列名去掉'_mag'后添加到滤光片列表中

# 分割图像文件
segmfile = 'JADES_detections_segm.fits'  # 定义分割图像文件名
segm = fits.open(segmfile)[0].data  # 打开FITS文件并获取数据部分
segm = SegmentationImage(segm)  # 将数据转换为分割图像对象

绘制数量计数与亮度的关系#

fig = plt.figure(figsize=(8, 4))  # 创建一个8x4英寸的图形

filt = 'F200W'  # 设置滤光片为F200W

mag1 = tbl[filt + '_mag']  # 从表格中提取对应滤光片的星等数据

mag1 = mag1[(0 < mag1) & (mag1 < 90)]  # 仅保留有效的观测数据（星等在0到90之间）

n = plt.hist(mag1, histtype='step', label=filt)  # 绘制星等的直方图，使用线条样式并添加标签

plt.xlabel('AB magnitude')  # 设置x轴标签为AB星等

plt.ylabel('Number counts')  # 设置y轴标签为数量计数

plt.legend()  # 显示图例

绘制 F200W 与 F090W 亮度并寻找缺失值#

#import mplcursors  # 导入mplcursors库（注释掉，可能是为了后续使用）

# Would love a better solution here!  # 这里希望有更好的解决方案！

mag1 = tbl['F090W_mag']  # 获取F090W波段的星等数据

mag2 = tbl['F200W_mag']  # 获取F200W波段的星等数据

# Only plot detections in F200W  # 仅绘制F200W波段的检测数据

det2 = (0 < mag2) & (mag2 < 90)  # 筛选出F200W波段星等在0到90之间的检测

mag1 = mag1[det2]  # 根据det2筛选F090W波段的星等
mag2 = mag2[det2]  # 根据det2筛选F200W波段的星等
ids = tbl['label'][det2]  # 获取对应的标签

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

plt.plot(mag1, mag2, '.')  # 绘制F090W与F200W星等的散点图

for i in range(len(mag1)):  # 遍历每个点
    plt.text(mag1[i], mag2[i], ids[i])  # 在每个点旁边添加对应的标签

plt.xlabel('F090W AB magnitude')  # 设置x轴标签
plt.ylabel('F200W AB magnitude')  # 设置y轴标签

查看一个天体#

# 可以通过位置选择对象

# x, y = 905, 276  # 定义x和y坐标

# id = segm.data[y,x]  # 根据坐标获取对象的ID

# 通过ID号选择对象
id = 261  # F090W dropout，定义要选择的对象ID

obj = tbl[id-1]  # 从表格中获取对应ID的对象

看起来您提供的信息不完整。请提供您希望我添加中文注释的Python代码，以便我可以帮助您。

# 获取对象的椭圆率属性
ellipticity = obj['ellipticity']  # 从对象中提取椭圆率

# 从segm对象中获取指定id的段对象
segmobj = segm.segments[segm.get_index(id)]

# 输出获取的段对象
segmobj

在所有图像中显示对象#

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(rgb[segmobj.slices], origin='lower', extent=segmobj.bbox.extent)  # 在第一个子图中显示RGB图像
ax[0, 0].set_title('Color')  # 设置第一个子图的标题为'Color'

cmap = segm.make_cmap(seed=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(imagefiles[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(weightfiles[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（光谱能量分布）#

fig, ax = plt.subplots(figsize=(8, 6))  # 创建一个8x6英寸的图形和坐标轴

for filt in filters:  # 遍历所有过滤器

    lam = int(filt[1:4]) / 100  # 提取波长并转换为微米

    plt.errorbar(lam, obj[filt+'_flux'].value, obj[filt+'_fluxerr'].value, marker='.', c='b')  # 绘制带误差条的光谱数据点

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轴标签为通量（nJy）

mlim = 31.4  # 设置AB星等的限制值
flim = mlim * u.ABmag  # 将限制值转换为AB星等
flim = flim.to(u.nJy).value  # 将限制值转换为nJy单位

# 添加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星等转换为ABmag单位
    f = m.to(u.nJy)  # 转换为nJy单位
    f = f.value  # 获取数值部分
    f = np.where(f > flim, f, flim)  # 如果小于flim，则使用flim
    return f  # 返回nJy值

def nJy2AB(F_nJy):  # 定义nJy到AB星等的转换函数
    f = F_nJy * u.nJy  # 将输入的nJy值转换为nJy单位
    m = f.to(u.ABmag)  # 转换为AB星等
    m = m.value  # 获取数值部分
    m = np.where(m < mlim, m, mlim)  # 如果大于mlim，则使用mlim
    return m  # 返回AB星等值

plt.ylim(flim, plt.ylim()[1])  # 设置y轴下限为flim，上限保持不变

secax = ax.secondary_yaxis('right', functions=(nJy2AB, AB2nJy))  # 创建右侧次坐标轴并设置转换函数
secax.set_ylabel('magnitude (AB)')  # 设置右侧坐标轴标签为AB星等

当通量 <= 0 时，星等转换失败#

fig, ax = plt.subplots(figsize=(8, 6))  # 创建一个8x6英寸的图形和坐标轴

for filt in filters:  # 遍历每个过滤器

    lam = int(filt[1:4]) / 100  # 提取波长并转换为微米

    plt.errorbar(lam, obj[filt+'_flux'].value, obj[filt+'_fluxerr'].value, marker='.', c='b')  # 绘制带误差条的光谱数据

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轴标签为通量（nJy）

f0 = 10**(0.4 * 31.4)  # 零星等时的通量（nJy）

b0 = 1.e-12  # 这个值应该依赖于过滤器

# 添加AB星等作为右侧的次坐标轴
# https://matplotlib.org/gallery/subplots_axes_and_figures/secondary_axis.html#sphx-glr-gallery-subplots-axes-and-figures-secondary-axis-py
def AB2nJy(m):  # 定义从AB星等到nJy的转换函数
    f = np.sinh(-0.4 * m * np.log(10) - np.log(b0)) * 2 * b0 * f0  # 计算nJy
    return f  # 返回计算结果

# Luptitudes
# https://www.sdss.org/dr12/algorithms/magnitudes/
def nJy2AB(f):  # 定义从nJy到AB星等的转换函数
    m = -2.5 / np.log(10) * (np.arcsinh((f / f0) / (2 * b0)) + np.log(b0))  # 计算AB星等
    return m  # 返回计算结果

#plt.ylim(flim, plt.ylim()[1])  # 可选：设置y轴范围

secax = ax.secondary_yaxis('right', functions=(nJy2AB, AB2nJy))  # 创建右侧的次坐标轴并设置转换函数
secax.set_ylabel('asinh magnitude (AB)')  # 设置次坐标轴的标签为asinh星等（AB）
Contents