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 1
Container Size: large
LSST Science Pipelines version: r29.2.0
Last verified to run: 2025-09-09
Repository: github.com/lsst/tutorial-notebooks

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

LSST data products: object, deep_coadd

Packages: lsst.daf.butler, lsst.rsp, lsst.geom.

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 Fornax dSph field, where the majority of stars are members of the Fornax dSph. Among many other available photometric measurements from the LSST science pipelines, PSF photometry is recommended for point-like sources. In DP1, 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 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), and re, a module to enable regular expression operations for pattern matching and text 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 geometry handling.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import re

from lsst.rsp import get_tap_service
from lsst.daf.butler import Butler
import lsst.geom

1.2. Define parameters and functions¶

Create an instance of the TAP service, and assert that it exists.

In [2]:

rsp_tap = get_tap_service("tap")
assert rsp_tap is not None

Instantiate the Butler with the appropiate DP1 repository and collection.

In [3]:

butler = Butler('dp1', collections='LSSTComCam/DP1')
assert butler is not None

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

In [4]:

ra_cen = 40.080
dec_cen = -34.450
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 Fornax dSph 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 Fornax dSph field.

2.1. Query for point-like objects¶

Query the Object table in the Fornax dSph 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) and require my_band_pixelFlags_inexact_psfCenter = 0 to exclude objects with a discontinuous PSF at the centroid.

In [7]:

query = f"""
SELECT objectId, x, y, 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 dp1.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}_pixelFlags_inexact_psfCenter = 0
ORDER BY objectId ASC
"""

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()
assert job.phase == 'COMPLETED'

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 90589 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 in DP1. 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: 3902
Number of PSF stars reserved: 475

Plot the spatial and magnitude distributions of PSF stars relative to all retrieved point sources (extendedness = 0 in my_band) in the Fornax dSph 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. The central hole in the spatial map arises from deblending being skipped in regions of high crowding.

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.

Note: This figure might raise the NumPy warning "converting a masked element to nan" and it can be safely ignored because subsequent steps handle NaNs explicitly.

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

/opt/lsst/software/stack/conda/envs/lsst-scipipe-10.1.0/lib/python3.12/site-packages/numpy/_core/shape_base.py:63: UserWarning: Warning: converting a masked element to nan.
  result = asanyarray(arys[0])

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 typically 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

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

In [18]:

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

2.4.2. Make image cutouts¶

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

In [19]:

my_sphere_point = lsst.geom.SpherePoint(ra_star*lsst.geom.degrees,
                                        dec_star*lsst.geom.degrees)
point_image = wcs_coadd.skyToPixel(my_sphere_point)

Retrieve the PSF from the coadd, compute its image at the target position using the computeImage() method, convert it to float, and extract the pixel array.

In [20]:

psf_coadd = coadd.getPsf()
psf_image_coadd = psf_coadd.computeImage(point_image).convertF()
psf_array = psf_image_coadd.getArray()

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

star_cutout = coadd.getCutout(point_image, psf_image_coadd.getDimensions())
star_image_array = star_cutout.getMaskedImage().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"(x (pix), y (pix)): {point_image}", 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 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 = 6.4456 pixel$^2$ suggests a full-width-half-maximum (FWHM) of $\sim$0.85" for the 0.2"/pixel pixel scale (FWHM = 2.355$\times$T/2$\times$0.2"). The median e1 = 0.0196 and e2 = 0.0251 give |e| = $\sqrt(e_1^2 + e_2^2)$ 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 magntidue 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 mangitude 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 [ ]: