205.1. PSF for visit images

205.1. The Point Spread Function for Visit Images¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 1
Container Size: large
LSST Science Pipelines version: r29.1.1
Last verified to run: 2025-06-20
Repository: github.com/lsst/tutorial-notebooks

Learning objective: To understand the format and metadata of the Point Spread Function in a visit image.

LSST data products: source, visit_image

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

Credit: Originally developed by Andrés A. Plazas Malagón and the Rubin Community Science team. Section 7 on "PSF model residuals" is based on a notebook by Pierre-François Léget. 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 visit_image at a particular location, as well as how to explore the PSF model methods to calculate and visualize its properties. The final section shows an example of PSF model residual visualization.

A visit_image is a fully processed, science-ready astronomical image from a single detector, sometimes also referred to as a Processed Visit Image (pvi). A visit_image object includes attributes such as a World Coordinate System (WCS) astrometric solution and an associated Point Spread Function (PSF) model.

For LSSTComCam, the PSF modeling was performed using the "PSF in the Full Field of View" (PIFF) software, with a PixelGrid model and BasisPolynomial interpolation on a per-CCD (Charge-Coupled Device) basis, in either pixel or sky coordinates. Currently, the model does not include features such as tree ring distortions or chromatic PSF corrections, which are planned for future implementation.

Related tutorials: DP1 tutorial on visit_images and PSF in coadded images.

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; scipy docs).

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

From lsst.rsp, import functions to interact with the Rubin Science Platform’s TAP (Table Access Protocol) service (rsp.lsst.io; dm.lsst.org).

In [1]:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy.stats import skew

import lsst.afw.display as afwDisplay
import lsst.geom
from lsst.daf.butler import Butler
from lsst.rsp import get_tap_service

1.2. Define parameters and functions¶

Constant: SIGMA_TO_FWHM

Defines the convertion factor between the standard deviation sigma of a Gaussian function a the full-withd 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

Defines 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 function evaluated at `x`.
    """
    return a * np.exp(-(x - x0)**2 / (2 * sigma**2))

Set afwDisplay to use matplotlib for image display.

In [4]:

afwDisplay.setDefaultBackend('matplotlib')

Create an instance of the Table Access Protocol (TAP) service, and assert that it exists.

In [5]:

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

Instantiate the Butler with the appropiate DP1 repository and collection.

In [6]:

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

2. Visualize the PSF model¶

For a given visit ID and a given ComCam detector, get the visit_image via the Butler.

First, retrieve a visit reference object around a chosen coordinate on the sky. Use the center location (International Celestial Reference System coordinates in degrees) of one of the fields observed by LSSTComCam during commissioning, the Extended Chandra Deep South Field.

In [7]:

ra, dec = 53.13, -28.10
my_band = 'i'

In [8]:

visitimage_refs = butler.query_datasets('visit_image',
                                        where="visit.region OVERLAPS POINT(ra, dec) "
                                        f"and BAND='{my_band}'",
                                        bind={'ra': ra, 'dec': dec},
                                        order_by=['visit', 'detector'],
                                        with_dimension_records=True,
                                        limit=1)

In [9]:

visitimage_refs

Out[9]:

[DatasetRef(DatasetType('visit_image', {band, instrument, day_obs, detector, physical_filter, visit}, ExposureF), {instrument: 'LSSTComCam', detector: 0, visit: 2024110800245, band: 'i', day_obs: 20241108, physical_filter: 'i_06'}, run='LSSTComCam/runs/DRP/DP1/DM-51335', id=1e1cd322-4ddd-4918-b2ac-c38e94fc8af5)]

Retrieve the visit_image object from the reference.

In [10]:

visit_image = butler.get(visitimage_refs[0])

Get the PSF model associated to the visit_image.

In [11]:

psf = visit_image.getPsf()

Define a central point where the PSF model should be evaluated.

In [12]:

xy = lsst.geom.Point2D(300, 300)

Use the computeImage() method to display the PSF model centered around the chosen point.

In [13]:

psf_image = psf.computeImage(xy)

In [14]:

afw_display = afwDisplay.Display(frame=1)
afw_display.mtv(psf_image.convertF())

Figure 1: PSF image using the computeImage() method.

Use the computeKernelImage() method. Note that with this method the postage stamp's coordinates are centered at the origin of the image. The coordinates of this origin point are (0, 0), resulting in negative coordinates for the lower left point.

In [15]:

psf_image_kernel = psf.computeKernelImage(xy)

In [16]:

afw_display = afwDisplay.Display(frame=1)
afw_display.mtv(psf_image_kernel.convertF())

Figure 2: PSF image using the computeKernelImage() method.

3. PSF model methods and properties¶

Option: explore the methods available for psf

In [17]:

# psf.

3.1. PSF size¶

Obtain the size of the PSF, as a $\sigma$ and full-width at half maximum (FWHM).

In [18]:

psf_shape = psf.computeShape(xy)

In [19]:

psf_sigma = psf_shape.getDeterminantRadius()
psf_fwhm = psf_sigma * SIGMA_TO_FWHM
print(f"PSF size: {psf_sigma: .3f} pix.\nPSF FWHM {psf_fwhm: .3f} pix.")

PSF size:  1.939 pix.
PSF FWHM  4.566 pix.

3.2. Second moments¶

Display the PSF second moments.

In [20]:

print(f"PSF Ixx: {psf_shape.getIxx(): .3f} pix^2")
print(f"PSF Iyy: {psf_shape.getIyy(): .3f} pix^2")
print(f"PSF Ixy: {psf_shape.getIxy(): .3f} pix^2")

PSF Ixx:  3.836 pix^2
PSF Iyy:  3.753 pix^2
PSF Ixy:  0.505 pix^2

3.3. Aperture flux¶

Display the PSF flux from aperture photometry weighted by a sinc function within 1, 2, and 3$\sigma$. Recall that the total flux of the PSF kernel has been normalized to an integrated flux of 1.

In [21]:

psf_apflux_1s = psf.computeApertureFlux(radius=1.0*psf_sigma, position=xy)
psf_apflux_2s = psf.computeApertureFlux(radius=2.0*psf_sigma, position=xy)
psf_apflux_3s = psf.computeApertureFlux(radius=3.0*psf_sigma, position=xy)

In [22]:

print("PSF aperture fluxes: ")
print(f"  within 1σ: {psf_apflux_1s: .3f}")
print(f"  within 2σ: {psf_apflux_2s: .3f}")
print(f"  within 3σ: {psf_apflux_3s: .3f}")

PSF aperture fluxes:
  within 1σ:  0.348
  within 2σ:  0.738
  within 3σ:  0.894

3.4. Peak flux¶

Print the peak flux of the PSF kernel.

In [23]:

psf_peak = psf.computePeak(position=xy)

In [24]:

print(f"PSF peak value at center: {psf_peak: .3f}")

PSF peak value at center:  0.040

3.5. Postage stamp dimensions¶

Calculate the PSF postage stamp dimensions, in pixels.

In [25]:

psf_dims = psf.computeImage(xy).getDimensions()

In [26]:

print(f"PSF postage stamp dimensions: {psf_dims} pix ")

PSF postage stamp dimensions: (25, 25) pix

4. PSF radial profile¶

Extract the data array from the PSF image, and compute its center from the PSF dimensions.

In [27]:

psf_image_array = psf_image.array
xlen, ylen = psf_dims
center = np.array([xlen / 2, ylen / 2])

In [28]:

print(f"PSF center position: ({center[0]: .1f}, {center[1]: .1f}) pix")

PSF center position: ( 12.5,  12.5) pix

Create a coordinate grid and calculate the distances from all the pixels to the center.

In [29]:

y, x = np.indices((xlen, ylen))
coords = np.stack((x, y), axis=-1)
distances = np.linalg.norm(coords - center, axis=-1)

Apply a radial mask.

In [30]:

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

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

In [31]:

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

Fit a Gaussian profile.

In [32]:

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

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)
    y_fit = gauss(x_fit, *pars)
    fitline = gauss(dists, *pars)
    skewness_val = skew(y_fit)

    ax.plot(dists, fitline, 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 (ADU)")
ax.set_title("Azimuthally-averaged radial profile")
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¶

Sort the distances calculated above.

In [34]:

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

Compute the normalized cumulative flux.

In [35]:

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

Plot the PSF curve of growth.

In [36]:

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

In [37]:

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

xx, yy = np.meshgrid(np.linspace(-xlen / 2, xlen / 2, xlen),
                     np.linspace(-ylen / 2, ylen / 2, ylen))

Plot the contours.

In [39]:

fig, ax = plt.subplots(figsize=(6, 5))
ax.contour(xx, yy, 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")
ax.set_xlim([-8, 8])
ax.set_ylim([-8, 8])
plt.tight_layout()
plt.show()

Figure 5: PSF countour plot.

7. PSF model residuals¶

Find a list of stars used for PSF modeling in the source catalog by setting the calib_psf_used flag to True, for the visit_id and detector_id defined before. For one star, display the residuals between the Piff PSF model at the star location and the observed star, using the visit_image object.

Get a table of stars used for PSF modeling by querying the source table with the TAP service.

In [40]:

visit_id = visitimage_refs[0].dataId['visit']
detector_id = visitimage_refs[0].dataId['detector']

In [41]:

query = (
    "SELECT sourceId, x, y, "
    "coord_ra, coord_dec, "
    "psfFlux, "
    "calib_psf_used "
    "FROM dp1.Source "
    f"WHERE visit = {visit_id} "
    f"AND detector = {detector_id} "
    "AND calib_psf_used = 1 "
    "ORDER BY sourceId "
)

In [42]:

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'
psf_used_table = job.fetch_result().to_table()

Job phase is COMPLETED

In [43]:

print(len(psf_used_table))

76

Choose one star.

In [44]:

x_star, y_star = psf_used_table['x'][1], psf_used_table['y'][1]

In [45]:

print(f"Star position: ({x_star: .1f}, {y_star: .1f}) pix")

Star position: ( 1947.8,  208.3) pix

Use the visit_image for the visit_id and detector_id obtained previously.

Get a postage stamp image of the star around its location, using as size the PSF model dimensions calculated previously.

In [46]:

position_star = lsst.geom.Point2D(x_star, y_star)
star_cutout = visit_image.getCutout(position_star, lsst.geom.Extent2I(xlen, ylen))

Get the observed star image.

In [47]:

star_masked_image = star_cutout.getMaskedImage()
star_image_array = star_masked_image.image.array / np.sum(star_masked_image.image.array)

Get the PSF model image.

In [48]:

psf_model = psf.computeImage(position_star)
psf_array = psf_model.array

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

In [49]:

max_star_image_array = np.max(np.abs(star_image_array))
PSF_model_name = "Piff"

fig, axes = plt.subplots(1, 3, figsize=(14, 5))
fig.subplots_adjust(wspace=0.3, left=0.07, right=0.95, bottom=0.15, top=0.8)
fig.suptitle(f"(x, y): {position_star} pix", fontsize=12)

images = [
    (star_image_array, 'Observed star', max_star_image_array),
    (
        psf_array,
        f"measured PSF\n(PSF model = {PSF_model_name})",
        max_star_image_array
    ),
    (
        star_image_array - psf_array,
        f"Residuals\n(Star - PSF model ({PSF_model_name}))",
        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')
    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 from Piff (center), and residuals (right).