photutils

# 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 ```bash pip install photutils ``` ### Standard Imports ```python 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 ```python 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) ```python # ❌ 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 ```python 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 ```python 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 ```python 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 ```python 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 ```python 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) ```python 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 ```python 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`. ```python # ✅ 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. ```python # ❌ 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). ```python # ❌ 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.