205.1. PSF for deep coadd imags

205.1. PSF for deep coadd imags¶

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

Learning objective: To understand the format and metadata of the Point Spread Function (PSF) in a coadded image, and how to access and characterize it in the DP2 deep_coadd data model.

LSST data products: object, deep_coadd

Packages: lsst.daf.butler, lsst.afw.display, lsst.geom, astropy, numpy, matplotlib, scipy, galsim.

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 how to visualize the PSF model of a deep_coadd at a particular location, and how to explore the PSF model methods to calculate and visualize its properties. The final section shows an example of PSF model residual visualization.

In DP2, a deep_coadd is stored in the new lsst.images data model as a CellCoadd object. Its PSF is accessed through the deep_coadd.psf attribute, and its astrometric mapping through deep_coadd.sky_projection (see the DP2 tutorial on deep_coadd images and astrometric calibration). The PSF modeling was performed with the "PSF in the Full Field of View" (PIFF) software; the coadd PSF is a position-dependent weighted sum of the contributing single-visit PSF models.

The new PSF object provides compute_kernel_image and compute_stellar_image methods.

Related tutorials: See the DP2 tutorials on deep_coadd images, on the PSF in visit_images, and on cosmological applications of PSF analysis.

1.1. Import packages¶

Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), and matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery).

From scipy, import tools for numerical optimization and statistics, including curve fitting and computation of distribution skewness (scipy.org).

From astropy, import celestial coordinate and unit handling (astropy.org).

Import galsim, the modular galaxy image simulation toolkit, for its HSM adaptive-moments interface (github.com/GalSim-developers/GalSim).

From the LSST Science Pipelines, import modules for display utilities, geometry handling, and the Butler data access system (pipelines.lsst.io).

In [1]:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy.stats import skew
import galsim
from astropy.coordinates import SkyCoord
import astropy.units as u
import lsst.afw.display as afwDisplay
from lsst.images import Box
import lsst.geom
from lsst.daf.butler import Butler

1.2. Define parameters and functions¶

Constant: SIGMA_TO_FWHM

The conversion factor between the standard deviation sigma of a Gaussian function and the full-width at half-maximum (FWHM): $2\sqrt{2\ln(2)}$.

In [2]:

SIGMA_TO_FWHM = 2.0*np.sqrt(2.0*np.log(2.0))

Function: gauss

A one-dimensional Gaussian profile.

In [3]:

def gauss(x, a, x0, sigma):
    """Evaluate a 1D Gaussian function.

    Parameters
    ----------
    x : `array_like`
        Input values where the Gaussian is evaluated.
    a : `float`
        Amplitude of the Gaussian peak.
    x0 : `float`
        Mean (center) of the Gaussian.
    sigma : `float`
        Standard deviation (width) of the Gaussian.

    Returns
    -------
    y : `array_like`
        Values of the Gaussian evaluated at `x`.
    """
    return a * np.exp(-(x - x0)**2 / (2 * sigma**2))

Set afwDisplay to use firefly for image display.

In [4]:

afwDisplay.setDefaultBackend('firefly')
afw_display = afwDisplay.Display(frame=1)

Instantiate the Butler with the DP2 repository and collection.

In [5]:

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

2. Retrieve the coadd and visualize the PSF model¶

Use a location near the center of the Extended Chandra Deep Field South (ECDFS), in degrees, and the r band.

In [6]:

ra_cen, dec_cen = 53.076, -28.110
my_band = 'r'

Query the Butler for a deep_coadd that overlaps this location, and retrieve it. The returned object is a CellCoadd.

In [7]:

query = f"band.name=:band AND patch.region OVERLAPS POINT(:ra, :dec)"
bind_params = {'band':my_band, 'ra':ra_cen, 'dec':dec_cen}

dataset_refs = butler.query_datasets("deep_coadd", where=query,
                                     bind=bind_params)
deep_coadd = butler.get(dataset_refs[0])
deep_coadd

Out[7]:

CellCoadd([y=2850:6150, x=14850:18150], tract=5063)

Get the PSF model and the sky projection (WCS) from the coadd.

In [8]:

psf = deep_coadd.psf
wcs = deep_coadd.sky_projection

Note: Here psf and wcs are taken from the deep_coadd object loaded above. When the full image is not needed, each can instead be read directly from the Butler as a component—for example butler.get("deep_coadd.psf", ...) and butler.get("deep_coadd.sky_projection", ...)—without loading the whole coadd.

Convert the central sky coordinate to pixel coordinates (in the tract system used by the coadd).

In [9]:

center = SkyCoord(ra=ra_cen*u.deg, dec=dec_cen*u.deg, frame='icrs')
xy = wcs.sky_to_pixel(center)
x_pix, y_pix = xy.x, xy.y
print(f"pixel position: ({x_pix: .1f}, {y_pix: .1f})")

pixel position: ( 15182.5,  4390.5)

Use the compute_kernel_image method to evaluate the PSF at the chosen point. This returns an lsst.images.Image centered on the origin (0, 0), suitable for convolution.

In [10]:

psf_kernel_image = psf.compute_kernel_image(x=x_pix, y=y_pix)
afw_display = afwDisplay.Display(frame=1)
afw_display.scale('asinh', 'zscale')
afw_display.mtv(psf_kernel_image)

Figure 1: PSF model from the compute_kernel_image method, centered on the origin.

Use the compute_stellar_image method. This evaluates the PSF centered on the given coordinates, just like the postage stamp of a star would appear, which is convenient for comparison with observed stars.

In [11]:

psf_stellar_image = psf.compute_stellar_image(x=x_pix, y=y_pix)

afw_display = afwDisplay.Display(frame=2)
afw_display.scale('asinh', 'zscale')
afw_display.mtv(psf_stellar_image)

Figure 2: PSF model from the compute_stellar_image method, centered on the star position.

The compute_stellar_bbox method returns the bounding box of the image that compute_stellar_image would produce at a given position, without computing the image itself. This is useful to know the size and location of the stamp in advance.

In [12]:

stellar_bbox = psf.compute_stellar_bbox(x=x_pix, y=y_pix)
stellar_bbox

Out[12]:

Box(y=Interval(start=4374, stop=4409), x=Interval(start=15166, stop=15201))

3. PSF size and shape¶

The CellCoadd PSF (deep_coadd.psf, a CellPointSpreadFunction) does not provide the shape, second-moment, aperture-flux, or peak methods that the DP1 afw PSF had. This tutorial measures the PSF size and shape with the GalSim HSM adaptive moments (Section 3.2); the postage-stamp dimensions are read directly from the PSF image array (Section 3.1).

3.1. Postage stamp dimensions¶

Get the PSF postage stamp dimensions, in pixels, from the kernel image array.

In [13]:

psf_image_array = psf_kernel_image.array
ylen, xlen = psf_image_array.shape
print(f"PSF postage stamp dimensions: {xlen} x {ylen} pix")

PSF postage stamp dimensions: 35 x 35 pix

3.2. PSF size and shape with GalSim HSM¶

An independent way to characterize the PSF size and shape is with adaptive moments, using the HSM algorithm (Hirata & Seljak 2003; Mandelbaum et al. 2005) as implemented in GalSim. galsim.hsm.FindAdaptiveMom fits an elliptical Gaussian to the image and returns the adaptive-moments size moments_sigma (the determinant radius $|\det M|^{1/4}$, in pixels) and the observed_shape (a Shear with ellipticity components e1 and e2). The HSM algorithm is the same algorithm that is used to measure PSF moments that are tabulated in various DP2 catalogs.

Wrap the PSF kernel image array in a galsim.Image with a pixel scale of 1 (so the size is measured in pixels), and run the adaptive-moments measurement.

In [14]:

gs_image = galsim.Image(np.ascontiguousarray(psf_kernel_image.array), scale=1.0)
hsm_result = galsim.hsm.FindAdaptiveMom(gs_image, guess_sig=3.0)

Report the HSM size, full-width at half maximum, and ellipticity.

In [15]:

hsm_sigma = hsm_result.moments_sigma
hsm_fwhm = hsm_sigma * SIGMA_TO_FWHM
shape = hsm_result.observed_shape
print(f"HSM size (moments_sigma): {hsm_sigma: .3f} pix")
print(f"HSM FWHM: {hsm_fwhm: .3f} pix")
print(f"HSM ellipticity: e1 = {shape.e1: .4f}, e2 = {shape.e2: .4f}, |e| = {shape.e: .4f}")

HSM size (moments_sigma):  2.161 pix
HSM FWHM:  5.090 pix
HSM ellipticity: e1 =  0.0135, e2 = -0.0051, |e| =  0.0145

4. PSF radial profile¶

Create a coordinate grid and calculate the distance from each pixel to the center. Using the meshgrid method associated with the psf bbox retrieves coordinates with (0, 0) at the center of the kernel image.

In [16]:

xx, yy = psf_kernel_image.bbox.meshgrid()
coords = np.stack((xx, yy), axis=-1)
distances = np.linalg.norm(coords, axis=-1)

Apply a radial mask.

In [17]:

mask = distances <= xlen / 2
values = psf_image_array[mask]
dists = distances[mask]

Set initial parameters and bounds for a Gaussian fit of the radial profile.

In [18]:

p0 = [np.max(values), 0, 10]
bounds = ((0, 0, 0), (np.inf, np.inf, np.inf))

Fit a Gaussian profile.

In [19]:

try:
    pars, _ = curve_fit(gauss, dists, values, p0=p0, bounds=bounds)
    pars[0] = np.abs(pars[0])
    pars[2] = np.abs(pars[2])
except RuntimeError:
    pars = None

Plot the radial profile.

In [20]:

fig, ax = plt.subplots(figsize=(6, 5))
ax.plot(dists, values, 'x', label='Radial data')

if pars is not None:
    fitAmp, fitMean, fitSigma = pars
    fitFwhm = fitSigma * SIGMA_TO_FWHM
    x_fit = np.linspace(-np.max(dists), np.max(dists), 200)
    skewness_val = skew(gauss(x_fit, *pars))
    ax.plot(dists, gauss(dists, *pars), label=(
        f"Gaussian Fit\n"
        f"Amp: {fitAmp: .3f}\n"
        f"Position: {fitMean: .2f}\n"
        f"FWHM: {fitFwhm: .2f}\n"
        f"Skewness: {skewness_val: .2f}"))

ax.set_xlabel("Radius (pix)")
ax.set_ylabel("Flux (normalized)")
ax.set_title("Azimuthally-averaged radial profile - PSF in deep_coadd")
ax.set_aspect(1.0 / ax.get_data_ratio(), adjustable='box')
ax.legend()
plt.tight_layout()
plt.show()

Figure 3: Azimuthally-averaged PSF radial profile, fitted with a Gaussian function.

5. PSF curve of growth (encircled energy)¶

Sort the distances calculated above.

In [21]:

sorted_indices = np.argsort(dists)
dists_sorted = dists[sorted_indices]
values_sorted = values[sorted_indices]

Compute the normalized cumulative flux.

In [22]:

cum_fluxes = np.cumsum(values_sorted)
cum_fluxes_norm = cum_fluxes / np.max(cum_fluxes)

Plot the curve of growth.

In [23]:

fig, ax = plt.subplots(figsize=(6, 5))
ax.plot(dists_sorted, cum_fluxes_norm, markersize=10)
ax.set_ylabel('Encircled flux (normalized)')
ax.set_xlabel('Radius (pix)')
ax.set_title("Encircled Flux - PSF in deep_coadd")
ax.grid(True)
plt.tight_layout()
plt.show()

Figure 4: PSF curve of growth.

6. Make a 2D PSF contour plot¶

Define intensity percentiles for the contour levels.

In [24]:

vmin = np.percentile(psf_image_array, 0.1)
vmax = np.percentile(psf_image_array, 99.9)
nContours = 10
lvls = np.linspace(vmin, vmax, nContours)

Define a coordinate grid centered around (0, 0).

In [25]:

xg, yg = np.meshgrid(np.linspace(-xlen / 2, xlen / 2, xlen),
                     np.linspace(-ylen / 2, ylen / 2, ylen))

Plot the contours.

In [26]:

fig, ax = plt.subplots(figsize=(6, 5))
ax.contour(xg, yg, psf_image_array, levels=lvls)
ax.tick_params(which="both", direction="in", top=True, right=True, labelsize=8)
ax.set_aspect("equal")
ax.set_xlabel('x (pix)')
ax.set_ylabel('y (pix)')
ax.set_title("Contour plot - PSF in deep_coadd")
ax.set_xlim([-8, 8])
ax.set_ylim([-8, 8])
plt.tight_layout()
plt.show()

Figure 5: PSF contour plot.

7. PSF model residuals¶

Find stars used for PSF modeling in the object catalog by selecting the {band}_calib_psf_used flag, retrieved from the Butler for the coadd's tract. For one star, display the residuals between the PSF model at the star location and the observed star.

Get the tract and patch of the deep_coadd.

In [27]:

data_id = dataset_refs[0].dataId
my_tract = data_id['tract']
my_patch = data_id['patch']
print('tract:', my_tract, '| patch:', my_patch)

tract: 5063 | patch: 15

Retrieve the object catalog for the tract, selecting only the columns needed. The object catalog contains all objects in the tract, so restrict to the coadd's patch.

In [28]:

columns = ["patch", "refExtendedness",
           f"{my_band}_calib_psf_used",
           f"{my_band}_pixelFlags_inexact_psfCenter",
           f"{my_band}_psfFlux", f"{my_band}_psfFluxErr",
           "coord_ra", "coord_dec"]
objects = butler.get("object", tract=my_tract, skymap="lsst_cells_v2",
                     parameters={"columns": columns})

Note: Here refExtendedness = 0 selects point-like sources, and is used for illustration. LSST DP2 provides several derived quantities that can be used for star–galaxy separation, described in the DP2 schema:

• {band}_extendedness (binary): classifies an object as a point source (0) if its PSF flux exceeds 0.985 times its cModel flux, and as extended (1) otherwise.
• refExtendedness (binary): the {band}_extendedness measured in the object's reference band (chosen per object by detection significance, in the priority order i, r, z, y, g, u).
• {band}_sizeExtendedness (binary): a likelihood-based comparison of the object's HSM moments and PSF size (measured on the g band by default).
• {band}_model_extendedness (floating point): a Sérsic-plus-PSF flux and size estimate of whether an object is point-like or extended.
• griz_model_extendedness (floating point): a multi-band (griz) Sérsic-model flux and size estimate.

The most suitable classifier depends on the science case.

Select point-like sources used in PSF modeling, in this patch, with a high signal-to-noise ratio.

In [29]:

snr = (np.asarray(objects[f"{my_band}_psfFlux"])
       / np.asarray(objects[f"{my_band}_psfFluxErr"]))
mask = ((np.asarray(objects["patch"]) == my_patch)
        & (np.asarray(objects["refExtendedness"]) == 0.0)
        & np.asarray(objects[f"{my_band}_calib_psf_used"])
        & ~np.asarray(objects[f"{my_band}_pixelFlags_inexact_psfCenter"])
        & (snr > 50))
star_ra = np.asarray(objects["coord_ra"])[mask]
star_dec = np.asarray(objects["coord_dec"])[mask]
print(len(star_ra))

52

Get the bounding box (bbox) for the deep coadd image, and find the first PSF star that falls within it.

In [30]:

bbox = deep_coadd.bbox

position_star = None
for ra_star, dec_star in zip(star_ra, star_dec):
    coord = SkyCoord(ra=ra_star*u.deg, dec=dec_star*u.deg, frame='icrs')
    pixel_point = wcs.sky_to_pixel(coord)
    if bbox.contains(x=pixel_point.x, y=pixel_point.y):
        position_star = pixel_point
        break

position_star

Out[30]:

XY(x=16130.208090528844, y=3291.1001875072616)

Evaluate the PSF model at the star location with compute_stellar_image, and normalize it.

In [31]:

psf_model = psf.compute_stellar_image(x=position_star.x,
                                      y=position_star.y)
psf_array = psf_model.array
psf_array = psf_array / np.sum(psf_array)
stamp_ny, stamp_nx = psf_array.shape

Use the dimensions of the PSF model image as the cutout size, then make a cutout around the star from the deep coadd and normalize its flux values.

In [32]:

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

In [33]:

star_cutout = deep_coadd[cutout_bbox].copy()

In [34]:

star_image_array = star_cutout.image.array
star_image_array /= star_image_array.sum()

Plot the observed star, the PSF model at that location, and the residuals.

In [35]:

max_star = np.max(np.abs(star_image_array))
PSF_model_name = "Piff"
print(max_star)

0.028992428

In [36]:

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)): {np.round(position_star.x, 1)} {np.round(position_star.y, 1)}", fontsize=12)

images = [
    (star_image_array, 'Observed star', max_star),
    (psf_array, 'model PSF', max_star),
    (star_image_array - psf_array,
     'Star - PSF model', max_star/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 6: Measured star used for PSF modeling (left), PSF model (center), and residuals (right).

In [ ]:

In [ ]: