photutils

Photutils - 天文测光

Photutils 是从天文图像中提取定量数据的现代标准。它以模块化、符合 Python 风格的接口取代了 SExtractor 或 DAOPHOT 等旧版工具，并与 Astropy 生态系统集成。

适用场景

• 检测天文图像中的恒星或星系（DAOStarFinder、IRAFStarFinder）。
• 使用圆形、椭圆形或矩形孔径测量天体亮度。
• 执行 PSF 拟合测光（在拥挤的恒星场中尤为重要）。
• 估计并扣除复杂的二维背景（天空背景）。
• 计算天文源的质心（精确中心）。
• 执行“等光度分析”（测量星系形状）。
• 基于阈值对图像进行分割。

参考文档

官方文档：https://photutils.readthedocs.io/
GitHub：https://github.com/astropy/photutils
搜索模式：photutils.aperture、photutils.detection、photutils.background、photutils.psf

核心原则

先处理背景

如果没有正确扣除背景，就无法进行准确测光。Photutils 提供了估计整体和局部天空背景的工具。

孔径测光与 PSF 测光

• 孔径测光：累加固定形状内的像素值。最适合孤立天体或星系。
• PSF 测光：将恒星的数学模型拟合到数据中。最适合彼此重叠的恒星或高精度工作。

单位和 WCS

Photutils 全面支持单位。如果图像包含 WCS 元数据，孔径可以用天球坐标（度）定义，而不只是像素坐标。

快速参考

安装

pip install photutils

标准导入

import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u
from photutils.detection import DAOStarFinder
from photutils.aperture import CircularAperture, aperture_photometry
from photutils.background import Background2D, MedianBackground

基本模式 - 检测和孔径测光

import numpy as np
from photutils.detection import DAOStarFinder
from photutils.aperture import CircularAperture, aperture_photometry

# 1. Detect sources
# fwhm: Full Width at Half Maximum of the stars
finder = DAOStarFinder(fwhm=3.0, threshold=5.*std)
sources = finder(image - background)

# 2. Define apertures at detected positions
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))
apertures = CircularAperture(positions, r=4.0)

# 3. Perform photometry
phot_table = aperture_photometry(image - background, apertures)
print(phot_table)

关键规则

✅ 应当

• 先扣除背景 - 在已扣除背景的图像上进行检测和测光，效果总是更好。
• 使用 Background2D - 对于天空背景变化的图像，二维网格估计比单个中值更准确。
• 指定增益 - 要准确估计误差（不确定度），请向 aperture_photometry 提供探测器增益。
• 匹配 FWHM - 确保查找器中的 fwhm 参数与图像的实际 PSF 相符。
• 向量化孔径 - 一个 CircularAperture 对象可以包含数千个位置；不要在循环中创建对象。
• 使用 SkyAperture - 如果有 WCS，请用 RA/Dec 定义孔径，以确保图像旋转或重新投影后孔径仍对准目标。

❌ 不应

• 忽略掩膜像素 - 如果图像中有坏像素（宇宙射线、饱和恒星），请将掩膜传递给测光函数。
• 使用过小的孔径 - 孔径太小会损失光（需要进行孔径修正）；太大则会引入噪声。
• 盲目相信默认值 - 务必将孔径叠加到图像上进行可视化，确认它们与恒星对齐。
• 混用像素和天球单位 - 请留意传入的是 (x, y) 还是 (RA, Dec)。

反模式（切勿采用）

# ❌ BAD: Manual pixel summation
# light_sum = np.sum(image[y-r:y+r, x-r:x+r]) # No circular masking!

# ✅ GOOD: Proper aperture photometry
from photutils.aperture import CircularAperture
ap = CircularAperture((x, y), r=r)
phot = aperture_photometry(image, ap)

# ❌ BAD: Detecting sources on raw data with high background
# finder = DAOStarFinder(threshold=100)(raw_image)

# ✅ GOOD: Detection on noise-subtracted data
from photutils.background import Background2D, MedianBackground
bkg = Background2D(image, (50, 50), filter_size=(3, 3), bkg_estimator=MedianBackground())
sources = finder(image - bkg.background)

背景估计（photutils.background）

稳健的二维背景映射

from photutils.background import Background2D, MedianBackground
from astropy.stats import SigmaClip

sigma_clip = SigmaClip(sigma=3.0)
bkg_estimator = MedianBackground()

# box_size: size of the mesh grid
bkg = Background2D(image, box_size=(50, 50), 
                   filter_size=(3, 3), 
                   sigma_clip=sigma_clip, 
                   bkg_estimator=bkg_estimator)

print(f"Global median background: {bkg.background_median}")
clean_image = image - bkg.background

源检测（photutils.detection）

恒星查找算法

from photutils.detection import DAOStarFinder, IRAFStarFinder

# DAOStarFinder: Optimized for stars
daofind = DAOStarFinder(fwhm=3.0, threshold=5.0*std)
sources = daofind(data)

# IRAFStarFinder: Alternative implementation
iraffind = IRAFStarFinder(threshold=5.0*std, fwhm=3.0)
sources_iraf = iraffind(data)

# Extract coordinates for photometry
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))

孔径测光（photutils.aperture）

处理多种形状

from photutils.aperture import CircularAperture, CircularAnnulus, aperture_photometry

# 1. Main aperture and Annulus (for local background)
positions = [(100.5, 120.3), (250.1, 88.4)]
aperture = CircularAperture(positions, r=5.0)
annulus = CircularAnnulus(positions, r_in=10.0, r_out=15.0)

# 2. Local background subtraction
# Calculate median in the annulus
annulus_masks = annulus.to_mask(method='center')
bkg_median = []
for mask in annulus_masks:
    annulus_data = mask.multiply(image)
    annulus_data_1d = annulus_data[mask.data > 0]
    bkg_median.append(np.median(annulus_data_1d))

# 3. Perform photometry and subtract local background
phot = aperture_photometry(image, aperture)
phot['local_bkg'] = bkg_median
phot['aper_sum_corrected'] = phot['aperture_sum'] - (bkg_median * aperture.area)

PSF 测光（photutils.psf）

拟合重叠源

from photutils.psf import BasicPSFFitter, IntegratedGaussianPRF
from photutils.background import MMMBackground

# 1. Define the PSF model (Gaussian or custom)
psf_model = IntegratedGaussianPRF(sigma=1.5)

# 2. Initialize the fitter
fitter = BasicPSFFitter(psf_model, fitshape=(11, 11), 
                        finder=daofind, 
                        bkg_estimator=MMMBackground())

# 3. Perform fitting
# result_table contains fitted x, y, and flux
result_table = fitter(image)

# 4. Get the residual image (original - fitted stars)
residual_image = fitter.get_residual_image()

质心计算（photutils.centroids）

精确定位恒星中心

from photutils.centroids import centroid_com, centroid_2dg

# Center of mass (simple)
x1, y1 = centroid_com(image[y-5:y+6, x-5:x+6])

# 2D Gaussian fit (most precise for stars)
x2, y2 = centroid_2dg(image[y-5:y+6, x-5:x+6])

实用工作流

1. 完整图像处理流程（从检测到星等）

def extract_magnitudes(image, zero_point=25.0):
    # 1. Subtract 2D Background
    bkg = Background2D(image, box_size=(64, 64))
    data_sub = image - bkg.background
    
    # 2. Detect sources
    finder = DAOStarFinder(fwhm=2.5, threshold=5*bkg.background_rms_median)
    sources = finder(data_sub)
    
    # 3. Aperture Photometry
    pos = np.transpose((sources['xcentroid'], sources['ycentroid']))
    ap = CircularAperture(pos, r=4.0)
    phot = aperture_photometry(data_sub, ap)
    
    # 4. Convert to instrumental magnitudes
    phot['mag'] = -2.5 * np.log10(phot['aperture_sum']) + zero_point
    return phot

2. 使用残差检查测光质量

def check_psf_quality(image, sources):
    """Fits PSF and returns residual to check for artifacts."""
    # ... setup PSF fitter ...
    result = fitter(image, init_guesses=sources)
    residual = fitter.get_residual_image()
    
    # If residual has 'holes' or 'donuts', PSF model is bad.
    return residual

性能优化

在拟合中使用 n_jobs

当前版本的 Photutils 并非所有拟合器都原生支持 n_jobs，但可以将图像拆分成多个区块来实现并行处理。

孔径掩膜

对于超大图像，可使用 aperture.to_mask() 提取相关子数组，而不是对整幅图像进行计算。

常见问题与解决方案

log10 中的负通量

如果背景扣除过于激进，一些暗星的 aperture_sum 可能为负值。

# ✅ Solution: Filter results before log
valid = phot[phot['aperture_sum'] > 0]
valid['mag'] = -2.5 * np.log10(valid['aperture_sum'])

拥挤区域

当恒星相互重叠时，孔径测光会失效。

# ❌ Problem: One aperture contains 3 stars
# ✅ Solution: Switch to PSF Photometry (photutils.psf)
# This fits the shapes and accounts for overlapping light.

坐标顺序

photutils 使用 (x, y)，对应 (列, 行)。

# ❌ Error: image[x, y] # This is wrong for numpy!
# ✅ Solution: Use image[y, x] or image[row, col]

Photutils 是天文图像的“测量尺”。通过结合稳健的统计检测与精确的数学拟合，它使天文学家能够将数字图像转化为带有已知不确定度的科学星表。