Photutils – Astronomical Photometry

Photutils is the modern standard for extracting quantitative data from astronomical images. It replaces legacy tools like SExtractor or DAOPHOT with a modular, Pythonic interface integrated with the Astropy ecosystem.

When to Use

• Detecting stars or galaxies in an astronomical image (DAOStarFinder, IRAFStarFinder).
• Measuring the brightness of objects using circular, elliptical, or rectangular apertures.
• Performing PSF-fitting photometry (essential for crowded stellar fields).
• Estimating and subtracting complex 2D backgrounds (sky background).
• Calculating the centroids (precise centers) of astronomical sources.
• Performing “isophotal analysis” (measuring shapes of galaxies).
• Segmentation of images based on thresholding.

Reference Documentation

Official docs: https://photutils.readthedocs.io/
GitHub: https://github.com/astropy/photutils
Search patterns: photutils.aperture, photutils.detection, photutils.background, photutils.psf

Core Principles

Background First

Accurate photometry is impossible without correct background subtraction. Photutils provides tools to estimate global and local sky levels.

Aperture vs. PSF

• Aperture Photometry: Summing pixels within a fixed shape. Best for isolated sources or galaxies.
• PSF Photometry: Fitting a mathematical model of a star to the data. Best for overlapping stars or high-precision work.

Units and WCS

Photutils is fully units-aware. If your image has WCS metadata, apertures can be defined in sky coordinates (degrees) rather than just pixels.

Quick Reference

Installation

pip install photutils

Standard Imports

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

Basic Pattern – Detection and Aperture Photometry

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)

Critical Rules

✅ DO

• Subtract background first – Detection and photometry always work better on a background-subtracted image.
• Use Background2D – For images with varying sky levels, a 2D mesh estimation is more accurate than a single median value.
• Specify the Gain – For accurate error (uncertainty) estimation, provide the detector gain to aperture_photometry.
• Match FWHM – Ensure the fwhm parameter in finders matches the actual PSF of your images.
• Vectorize apertures – One CircularAperture object can hold thousands of positions; don’t create them in a loop.
• Use SkyAperture – If you have a WCS, define apertures in RA/Dec to ensure they stay on target if the image is rotated or reprojected.

❌ DON’T

• Ignore mask pixels – If your image has bad pixels (cosmic rays, saturated stars), pass a mask to the photometry functions.
• Use small apertures – If an aperture is too small, you lose light (aperture correction required). If too large, you add noise.
• Blindly trust defaults – Always visualize your apertures on top of the image to verify they align with the stars.
• Mix pixel and sky units – Be careful whether you are passing (x, y) or (RA, Dec).

Anti-Patterns (NEVER)

# ❌ 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)

Background Estimation (photutils.background)

Robust 2D Background Mapping

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

Source Detection (photutils.detection)

Star Finding Algorithms

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']))

Aperture Photometry (photutils.aperture)

Handling Multiple Shapes

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 Photometry (photutils.psf)

Fitting Overlapping Sources

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()

Centroiding (photutils.centroids)

Precise Star Centering

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])

Practical Workflows

1. Full Image Pipeline (Detection to Magnitude)

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. Checking Photometry Quality with Residuals

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

Performance Optimization

Using n_jobs in fitting

While current Photutils doesn’t natively support n_jobs in all fitters, you can parallelize by splitting the image into tiles.

Aperture Masks

For very large images, use aperture.to_mask() to extract only relevant sub-arrays instead of performing calculations on the full image.

Common Pitfalls and Solutions

Negative Flux in log10

If the background subtraction is too aggressive, some faint stars might have negative aperture_sum.

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

Crowded Fields

Aperture photometry fails when stars overlap.

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

Coordinate Order

photutils uses (x, y) which corresponds to (column, row).

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

Photutils is the “measuring stick” for astronomical imagery. By combining robust statistical detection with precise mathematical fitting, it allows astronomers to turn digital images into scientific catalogs with known uncertainties.