302.1. PSF Photometry for stars

302.1. PSF Photometry for stars¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: Small
LSST Science Pipelines version: r30.0.9
Last verified to run: 2026-07-19
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: To understand the Point Spread Function (PSF) photometry for stars.

LSST data products: object, deep_coadd

Packages: lsst.daf.butler, lsst.rsp.RSPDiscovery, lsst.images.

Credit: Originally developed by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

This notebook demonstrates PSF photometry for stars in the ECDFS field. Among many other available photometric measurements from the LSST science pipelines, PSF photometry is recommended for point-like sources. In DP2, each visit image is modeled with a spatially varying PSF using Piff ("PSF in the Full Field of View"). The coadd PSF stored with each deep_coadd is the position-dependent weighted sum of the contributing PSF models from the input visit_images with seeing FWHM < 1".7. Building on this, the Object table provides forced PSF measurements on the deep_coadds at the locations of all deblended objects detected with signal-to-noise ratio > 5 in a deep_coadd of any filter.

This notebook uses the Object table to examine the spatial and magnitude distributions of bright, isolated stars selected for PSF modeling, with ~10% reserved for validation. This tutorial visualizes the coadd PSF at the location of a selected star and display residuals between the model and the observed star. Finally, it compares PSF fluxes to aperture fluxes as a function of aperture size.

Related tutorials: See the 200-level DP1 tutorials on PSFs in visit_image and deep_coadd.

1.1. Import packages¶

Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery), re, a module to enable regular expression operations for pattern matching and text manipulation, and astropy.coordinates as well as astropy.units modules for coordinate manipulation.

From the LSST pacakge (pipelines.lsst.io), import modules for accessing the Table Access Protocol (TAP) service, for retrieving datasets from the Butler, and for defining image cutout.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import re
from astropy.coordinates import SkyCoord
import astropy.units as u

from lsst.rsp import RSPDiscovery
from lsst.daf.butler import Butler
from lsst.images import Box

1.2. Define parameters and functions¶

Instantiate RSPDiscovery with the DP2 release, create an instance of the TAP service, and assert that it exists.

In [2]:

discovery = RSPDiscovery("dp2")
rsp_tap = discovery.get_tap_client()

Instantiate the Butler with the appropiate DP2 repository and collection.

In [3]:

butler = Butler('dp2', collections='dp2')

Define a circular region with a 1-degree radius centered on the ECDFS field (coordinates in degrees).

In [4]:

ra_cen = 53.13
dec_cen = -28.10
radius = 1.0

Define my_band as the $r$ band.

In [5]:

my_band = 'r'

Set font sizes globally for all subsequent matplotlib figures in this notebook.

In [6]:

plt.rcParams["font.size"] = 15
plt.rcParams["figure.titlesize"] = 15
plt.rcParams["legend.fontsize"] = 12
plt.rcParams["axes.titlesize"] = 12

2. PSF photometry in the ECDFS field¶

The Object table contains forced flux measurements on the deep coadded images at the locations of all objects detected with signal-to-noise ratio > 5 in a deep_coadd of any filter. Examine PSF photometry for stars in the ECDFS field.

2.1. Query for point-like objects¶

Query the Object table in the field for coordinates, PSF photometry, multiple aperture fluxes, and flags indicating whether objects were used for PSF modeling or reserved for PSF evaluation. Restrict to point sources (extendedness = 0 in my_band) with signal-to-noise ratio in my_band > 5, and require my_band_pixelFlags_inexact_psfCenter = 0 to exclude objects with a discontinuous PSF at the centroid.

In [7]:

query = f"""
SELECT
    coord_ra, coord_dec,
    {my_band}_psfMag, {my_band}_psfMagErr, {my_band}_psfFlux, 
    {my_band}_ixxPSF, {my_band}_iyyPSF, {my_band}_ixyPSF, 
    {my_band}_ap03Flux, {my_band}_ap06Flux, {my_band}_ap09Flux, 
    {my_band}_ap12Flux, {my_band}_ap17Flux, {my_band}_ap25Flux, 
    {my_band}_ap35Flux, {my_band}_ap50Flux, {my_band}_ap70Flux, 
    {my_band}_calib_psf_used, {my_band}_calib_psf_reserved
FROM
    dp2.Object
WHERE
    CONTAINS(POINT('ICRS', coord_ra, coord_dec), 
    CIRCLE('ICRS', {ra_cen}, {dec_cen}, {radius})) = 1
    AND {my_band}_extendedness = 0
    AND {my_band}_psfFlux / {my_band}_psfFluxErr > 5
    AND {my_band}_pixelFlags_inexact_psfCenter = 0
"""

In [8]:

job = rsp_tap.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()

Job phase is COMPLETED

Fetch the results, store them as a table, then sort by my_band in ascending order, i.e., brighter stars first.

In [9]:

table = job.fetch_result().to_table()
table.sort(f'{my_band}_psfMag')
print(f'The query returned {len(table)} objects.')

The query returned 75879 objects.

Option to display the table of results.

In [10]:

# table

2.2. PSF modeling and evaluation stars¶

The LSST science pipelines use isolated, bright stars for PSF modeling, with roughly 10% reserved to evaluate the PSF fit. Examine the number of stars used for PSF modeling and reserved for evaluation, along with their spatial and magnitude distributions.

In [11]:

psf_used = table['r_calib_psf_used'] == 1
psf_reserved = table['r_calib_psf_reserved'] == 1

print(f"Number of PSF stars used: {psf_used.sum()}")
print(f"Number of PSF stars reserved: {psf_reserved.sum()}")

Number of PSF stars used: 5684
Number of PSF stars reserved: 657

Plot the spatial and magnitude distributions of PSF stars relative to all retrieved point sources (extendedness = 0 in my_band) in the field.

In [12]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

hb = ax1.hexbin(table['coord_ra'], table['coord_dec'],
                gridsize=100, bins='log', cmap='gray_r')
ax1.scatter(table['coord_ra'][psf_used], table['coord_dec'][psf_used],
            s=3, alpha=0.3, c='tab:orange', label='PSF used')
ax1.scatter(table['coord_ra'][psf_reserved], table['coord_dec'][psf_reserved],
            s=7, marker="x", c='tab:green', label='PSF reserved')
ax1.legend(loc=2)
ax1.set_xlabel('RA (deg)')
ax1.set_ylabel('Dec (deg)')
cb = fig.colorbar(hb, ax=ax1, label='Counts of all point-like objects')

ax2.hist(table[f'{my_band}_psfMag'], bins=100, range=(16, 27),
         density=True, histtype='step', lw=2, label='All stars')
ax2.hist(table[f'{my_band}_psfMag'][psf_used], bins=100, range=(16, 27),
         density=True, histtype='step', lw=2, label='PSF used')
ax2.hist(table[f'{my_band}_psfMag'][psf_reserved], bins=100, range=(16, 27),
         density=True, histtype='step', ls='--', lw=2, label='PSF reserved')
ax2.legend()
ax2.set_xlabel(f'{my_band}_psfMag (mag)')
ax2.set_ylabel('Normalized density')

plt.tight_layout()
plt.show()

Figure 1: Spatial (left) and magnitude (right) distributions of all stars, PSF-used stars, and PSF-reserved stars, with PSF-selected stars being brighter and more isolated.

2.3. Photometric uncertainty¶

Visualize photometric uncertainty as a function of magnitude for stars used in PSF modeling and for those reserved to evaluate the PSF modeling. Also plot the average photometric uncertainty in 0.5-magnitude bins for all stars in the field.

In [13]:

fig, ax = plt.subplots(figsize=(5, 4))

mag = table[f'{my_band}_psfMag']
magerr = table[f'{my_band}_psfMagErr']

bins = np.arange(mag.min(), mag.max(), 0.5)
bin_centers = 0.5 * (bins[:-1] + bins[1:])
bin_pairs = zip(bins[:-1], bins[1:])
avg_magerr = [magerr[(mag >= b1) & (mag < b2)].mean() for b1, b2 in bin_pairs]

ax.plot(bin_centers, avg_magerr, c='k', marker="+", label='All stars')
ax.scatter(table[f'{my_band}_psfMag'][psf_used],
           table[f'{my_band}_psfMagErr'][psf_used],
           s=5, alpha=0.3, label="PSF used")
ax.scatter(table[f'{my_band}_psfMag'][psf_reserved],
           table[f'{my_band}_psfMagErr'][psf_reserved],
           s=5, alpha=0.6, label="PSF reserved")
ax.set_xlabel(f"{my_band}_psfMag (mag)")
ax.set_ylabel(f"{my_band}_psfMagErr (mag)")
ax.set_yscale("log")
ax.legend()

plt.tight_layout()
plt.show()

Figure 2: Photometric uncertainty versus magnitude for PSF-used and PSF-reserved stars, and average uncertainties in 0.5-mag bins for all stars in the field. PSF stars are bright and exhibit <1% uncertainty, consistent with the precision achieved by all bright stars.

In [14]:

del mag, magerr, bins, bin_centers, avg_magerr

2.4. Star-PSF comparison¶

This section compares an observed star and PSF model for the observed star.

2.4.1. Retrieve a deep_coadd image¶

Select the brightest PSF-used star and display the residuals between the observed image and the PSF model at its location.

In [15]:

ra_star = table['coord_ra'][psf_used][0]
dec_star = table['coord_dec'][psf_used][0]

Query the butler for a deep_coadd image near the selected star.

In [16]:

coadd_datasetrefs = butler.query_datasets(
    "deep_coadd",
    where=f"band.name='{my_band}' AND\
    patch.region OVERLAPS POINT(ra, dec)",
    bind={"ra": ra_star, "dec": dec_star},
    with_dimension_records=True,
    order_by=["patch.tract"])

Option to display the dataset references.

In [17]:

# coadd_datasetrefs[0]

Get the first deep_coadd image in the list and its World Coordinate System (WCS) information.

In [18]:

coadd = butler.get(coadd_datasetrefs[0])
wcs_coadd = coadd.sky_projection

2.4.2. Make image cutouts¶

Find the pixel location of the selected star on the deep_coadd image using the WCS function.

In [19]:

sky_coord = SkyCoord(ra=ra_star*u.deg, dec=dec_star*u.deg, frame='icrs')
pixel_coord = wcs_coadd.sky_to_pixel(sky_coord)
print(pixel_coord)

XY(x=11329.410524412953, y=25422.19512515213)

Retrieve the PSF from the deep_coadd, compute its image at the target position using the compute_stellar_image method and extract the pixel array.

In [20]:

psf_coadd = coadd.psf
psf_image = psf_coadd.compute_stellar_image(x=pixel_coord.x, y=pixel_coord.y)
psf_array = psf_image.array

Extract a cutout centered on the target star with the same dimensions as the PSF image, obtain its masked image, and normalize the pixel array by its total flux.

In [21]:

psf_bbox = psf_image.bbox
bbox = Box.factory[psf_bbox.y.start:psf_bbox.y.stop, psf_bbox.x.start:psf_bbox.x.stop]

star_cutout = coadd[bbox].copy()
star_image_array = star_cutout.image.array
star_image_array /= star_image_array.sum()

2.4.3. Display PSF residuals¶

Plot the observed star, the PSF model at its location, and their residual.

In [22]:

max_star_image_array = np.max(np.abs(star_image_array))

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
fig.subplots_adjust(wspace=0.3, left=0.07, right=0.95, bottom=0.15, top=0.8)
fig.suptitle(f"{pixel_coord}", fontsize=12)

images = [
    (star_image_array, 'Observed star', max_star_image_array),
    (psf_array, 'model PSF', max_star_image_array),
    (star_image_array - psf_array,
     'Star - PSF model', max_star_image_array/10)
]

for ax, (img, title, v) in zip(axes, images):
    im = ax.imshow(img, vmin=-v, vmax=v, cmap='viridis',
                   origin='lower', aspect='auto')
    fig.colorbar(im, ax=ax)
    ax.set_xlabel('x (pixel)', fontsize=14)
    ax.set_ylabel('y (pixel)', fontsize=14)
    ax.set_title(title, fontsize=12)

Figure 3: Observed star (left), PSF model at the same location (middle), and residual image (right). The small residuals demonstrate the high quality of the PSF fit.

In [23]:

del max_star_image_array, images, img

2.5. PSF performance¶

The LSST Science Pipelines quantify the PSF performance via the trace of the second-moment matrix T and the ellipticity components e$_1$ and e$_2$ (see DP1 paper). Compute the object size estimator T (the trace of the second-moment matrix) and the ellipticity components e$_1$, and e$_2$ for the PSF-reserved stars.

In [24]:

ixx = table[f"{my_band}_ixxPSF"][psf_reserved]
iyy = table[f"{my_band}_iyyPSF"][psf_reserved]
ixy = table[f'{my_band}_ixyPSF'][psf_reserved]

T_psf = ixx + iyy
e1_psf = (ixx - iyy) / T_psf
e2_psf = 2 * ixy / T_psf

Display their distributions.

In [25]:

fig, axes = plt.subplots(1, 3, figsize=(14, 4))

mt = np.ma.median(T_psf)
h_t = axes[0].hist(T_psf, bins=50, density=True, histtype="step")
axes[0].axvline(mt, c="tab:orange", linestyle="--")
axes[0].set_xlabel(f"$T$ in the {my_band}-band [pixel$^2$]")
axes[0].set_ylabel("Normalized density")
axes[0].set_title(f"median = {mt: .4f}")

m1 = np.ma.median(e1_psf)
h_e1 = axes[1].hist(e1_psf, bins=50, density=True, histtype="step")
axes[1].axvline(m1, c="tab:orange", linestyle="--")
axes[1].set_xlabel(f"$e_1$ in the {my_band}-band")
axes[1].set_title(f"median = {m1: .4f}")

m2 = np.ma.median(e2_psf)
h_e2 = axes[2].hist(e2_psf, bins=50, density=True, histtype="step")
axes[2].axvline(m2, c="tab:orange", linestyle="--")
axes[2].set_xlabel(f"$e_2$ in the {my_band}-band")
axes[2].set_title(f"median = {m2: .4f}")

plt.show()

Figure 4: Histograms of the trace of the second-moment matrix T (left) and the ellipticity components e$_1$ (middle) and e$_2$ (right) for the PSF-reserved stars. The orange dashed lines mark the median values, which is printed as the figure title for each panel. The median T = 9.2282 pixel$^2$ suggests a full-width-half-maximum (FWHM) of $\sim$1.01" for the 0.2"/pixel pixel scale (FWHM = 2.355$\times\sqrt(T/2)\times$0.2"). The median e1 = 0.0323 and e2 = -0.0087 give |e| = $\sqrt(e_1^2 + e_2^2) \simeq$ 0.0335 and an axis ratio q = $\sqrt(1 - e)/(1 + e) \simeq$ 0.97 (q = 1 is circular).

In [26]:

del h_t, h_e1, h_e2

3. Compare PSF and Aperture fluxes¶

In principle, a properly normalized PSF model integrates to 100% of a point source’s flux, whereas a circular aperture measures only the encircled flux within its radius—requiring an aperture correction to recover the total. In practice, small apertures miss flux; large ones admit neighbor/sky contamination. This section compares circular-aperture fluxes (r = 3–70 pixels centered on stellar positions) to PSF fluxes for the PSF-used stars.

Define a function binned_stats that computes quantile values within each magnitude bin in the next cell.

In [27]:

def binned_stats(x, y, nbins=30, q=(0.1, 0.5, 0.9)):
    """
    Compute binned quantiles of `y` as a function of `x`.

    Parameters
    ----------
    x : array-like
        Independent variable used for binning. Non-finite values are ignored.
    y : array-like
        Dependent variable; quantiles are computed within bins of `x`.
    nbins : int, optional
        Number of equal-width bins spanning [min(x), max(x)]. Default is 30.
    q : float or sequence of float, optional
        Quantiles to compute in each bin, expressed in [0, 1]. Default is
        (0.1, 0.5, 0.9).

    Returns
    -------
    xc : ndarray of shape (nbins,)
        Bin centers (midpoints of the bin edges).
    qs : list[ndarray]
        List of arrays, one per requested quantile, each of shape (nbins,).
        Each element is the quantile of `y` among samples whose `x` falls in
        that bin. Bins with no samples are filled with NaN.

    Examples
    --------
    >>> xc, (q10, q50, q90) = binned_stats(x, y, nbins=20, q=(0.1, 0.5, 0.9))
    """
    x = np.asarray(x)
    y = np.asarray(y)
    m = np.isfinite(x) & np.isfinite(y)
    if not np.any(m):
        return np.array([]), [np.array([]) for _ in q]
    x, y = x[m], y[m]
    bins = np.linspace(np.nanmin(x), np.nanmax(x), nbins + 1)
    idx = np.digitize(x, bins) - 1
    xc = 0.5 * (bins[:-1] + bins[1:])
    qs = []
    for qq in q:
        arr = np.full(nbins, np.nan)
        for i in range(nbins):
            yy = y[idx == i]
            if yy.size > 0:
                arr[i] = np.nanpercentile(yy, 100 * qq)
        qs.append(arr)
    return xc, qs

Collect the table’s column names, select the aperture-flux columns, then sort them numerically by aperture radius and store the result in ap_cols.

In [28]:

colnames = table.colnames
ap_cols = [c for c in colnames if re.match(r"^r_ap\d+Flux$", c)]
ap_cols.sort(key=lambda s: int(re.search(r"ap(\d+)", s).group(1)))

Extract the PSF flux and PSF magnitude for the psf_used stars in band my_band, then build a boolean mask m that selects entries with finite, positive PSF flux.

In [29]:

psf = table[f"{my_band}_psfFlux"][psf_used]
psfmag = table[f"{my_band}_psfMag"][psf_used]
m = np.isfinite(psf) & (psf > 0)

Make a figure comparing PSF flux to aperture flux across different aperture sizes for stars in the field.

In [30]:

ncols = 3
nrows = int(np.ceil(len(ap_cols) / ncols)) if ap_cols else 1
fig, axs = plt.subplots(nrows=nrows, ncols=ncols, figsize=(4*ncols, 3.5*nrows))
axs = axs.ravel()

for ax, ap_col in zip(axs, ap_cols):
    ap = table[psf_used][ap_col]
    mm = m & np.isfinite(ap) & (ap > 0)
    x = psfmag[mm]
    ratio = psf[mm] / ap[mm]

    ax.scatter(x, ratio, s=3, alpha=0.25, c='gray')
    xc, (q10, q50, q90) = binned_stats(x, ratio, nbins=30, q=(0.1, 0.5, 0.9))
    ax.plot(xc, q50, lw=1.5, label="Median")
    ax.plot(xc, q10, lw=1, linestyle="--", label="10th quantile")
    ax.plot(xc, q90, lw=1, linestyle="--", label="90th quantile")

    ax.axhline(1.0, lw=1, linestyle=":")
    ax.set_ylabel(f"r_psfFlux / {ap_col}")
    ax.set_xlabel("PSF mag")
    ax.set_ylim(0, 2)
    ax.set_title(ap_col)

axs[0].legend()

for j in range(len(ap_cols), len(axs)):
    fig.delaxes(axs[j])

fig.tight_layout()
fig.suptitle("PSF vs. aperture flux ratios", y=1.02)
plt.show()

Figure 5: Each panel shows the ratio of PSF flux to a given aperture flux as a function of PSF magnitude for stars in the field. The aperture size ranges from 3 pixels to 70 pixels. Ratio values > 1 indicate apertures too small to capture the total stellar flux, while values < 1 suggest contamination from neighboring sources due to too large aperture sizes.

4. Exercise for the learner¶

Display the PSF residuals for the faintest PSF-used star and compare them with those of the brightest PSF-used star, as shown in Figure 3. Also compare these with the residuals of the brightest and faintest PSF-reserved stars.

In [ ]: