205.2. PSF for visits#

205_2_PSF_for_visits

205.2. PSF for visits¶

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-22
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: To understand and characterize the Point Spread Function (PSF) model of a single visit in DP2.

LSST data products: visit_summary

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

Credit: Originally developed by the Rubin Community Science team. Section 3.2 on PSF size and shape uses the GalSim HSM interface. 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 obtain the PSF model of a single visit in DP2, and how to characterize its size, shape, and profile.

In DP2 the per-detector PSF model for a visit is stored in the visit_summary table (a legacy lsst.afw.detection.Psf, modeled with the "PSF in the Full Field of View" (PIFF) software). This notebook retrieves that legacy PSF and wraps it into the new lsst.images PointSpreadFunction with the from_legacy method, so that it can be used through the same interface as the deep_coadd PSF (see the DP2 tutorial on the PSF in coadded images, and the astrometric calibration tutorial, which wraps the visit WCS the same way).

The new PSF object provides compute_kernel_image and compute_stellar_image methods. It does not provide the shape, second-moment, aperture-flux, or peak methods of the DP1 afw PSF; the PSF size and shape are measured here instead with the GalSim HSM adaptive moments (Section 3.2).

Note: Early DP2 (EDP2) does not include the visit_image pixel data, so this notebook works with the PSF model only. A comparison of the model against observed stars (PSF residuals) requires the visit-image pixels and is deferred to a future update.

Related tutorials: See the DP2 tutorials on the PSF in deep_coadd images and on astrometric calibration.

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

From scipy, import tools for numerical optimization and statistics (scipy.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 display utilities, the Butler, and the new lsst.images PSF interface (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
import lsst.afw.display as afwDisplay
from lsst.daf.butler import Butler
from lsst.images import Box
from lsst.images.psfs import PointSpreadFunction

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 visit 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'

Find a visit and detector whose footprint overlaps the target position, using the visit_detector_region dimension records (registry metadata only; no image is read).

In [7]:

detector_records = butler.query_dimension_records(
    "visit_detector_region",
    where="band.name = :band AND visit_detector_region.region OVERLAPS POINT(:ra, :dec)",
    bind={"band": my_band, "ra": ra_cen, "dec": dec_cen})
visit_id = detector_records[0].dataId["visit"]
detector_id = detector_records[0].dataId["detector"]
print("Using visit:", visit_id, "| detector:", detector_id)

Using visit: 2025072200466 | detector: 94

Get the visit_summary for this visit, and retrieve the legacy PSF and bounding box for the chosen detector.

In [8]:

visit_summary = butler.get("visit_summary", instrument="LSSTCam", visit=visit_id)
record = visit_summary.find(detector_id)
legacy_psf = record.getPsf()
detector_bbox = Box.from_legacy(record.getBBox())

Wrap the legacy afw PSF into an lsst.images PointSpreadFunction with from_legacy, passing the detector bounding box as the region where the model is valid. The factory automatically selects the appropriate wrapper for the Piff PSF.

In [9]:

psf = PointSpreadFunction.from_legacy(legacy_psf, bounds=detector_bbox)

Note: The PSF is obtained here from visit_summary. Once the full DP2 release is available, the visit PSF will also be accessible directly as visit_image.psf (the same component-based access used for deep_coadd.psf), which will be the preferred method; the visit_summary approach shown here is used in the interim.

Define a pixel position within the detector at which to evaluate the PSF model.

In [10]:

x_pix, y_pix = 2000.0, 2000.0

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

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.

In [12]:

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.

In [13]:

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

Out[13]:

Box(y=Interval(start=1987, stop=2014), x=Interval(start=1987, stop=2014))

3. PSF size and shape¶

The new visit PSF (a PointSpreadFunction) does not provide the shape, second-moment, aperture-flux, or peak methods of the DP1 afw PSF. 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 [14]:

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: 27 x 27 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).

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

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

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):  4.480 pix
HSM FWHM:  10.550 pix
HSM ellipticity: e1 =  0.0248, e2 =  0.0604, |e| =  0.0653

4. PSF radial profile¶

Compute the center of the PSF image from its dimensions.

In [17]:

center_pix = np.array([xlen / 2, ylen / 2])

Create a coordinate grid and calculate the distance from each pixel to the center.

In [18]:

yy, xx = np.indices((ylen, xlen))
coords = np.stack((xx, yy), axis=-1)
distances = np.linalg.norm(coords - center_pix, axis=-1)

Apply a radial mask.

In [19]:

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

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

Fit a Gaussian profile.

In [21]:

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

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 a visit")
ax.set_aspect(1.0 / ax.get_data_ratio(), adjustable='box')
ax.legend()
plt.tight_layout()
plt.show()

No description has been provided for this image

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

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

Compute the normalized cumulative flux.

In [24]:

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

Plot the curve of growth.

In [25]:

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 a visit")
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 [26]:

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

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

Plot the contours.

In [28]:

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 a visit")
ax.set_xlim([-8, 8])
ax.set_ylim([-8, 8])
plt.tight_layout()
plt.show()

Figure 5: PSF contour plot.

In [ ]: