205.2. PSF for deep coadd images

205.2. The Point Spread Function for Coadded 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 coadded image.

LSST data products: object, deep_coadd

Packages: lsst.daf.butler, lsst.afw.display, lsst.rsp, lsst.geom, 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 Leget. 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 object 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 deep_coadd object has attributes such as a World Coordinate System (WCS) astrometric solution and a Point Spread Function (PSF) model associated with it.

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

Related tutorials: See the 200-level DP1 tutorials on coadd images and PSF in visit 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, collection, and skymap.

In [6]:

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

2. Visualize the PSF model¶

Use the center location (International Celestial Reference System coordinates in degrees) of one of the fields observed by ComCam during commissioning, the Extended Chandra Deep South Field (ECDS).

In [7]:

ra_cen, dec_cen = 53.13, -28.10

Query the butler for an i-band deep_coadd image that is near the center of the ECDS field.

In [8]:

my_band = 'i'

In [9]:

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

In [10]:

center_coadd_datasetrefs

Out[10]:

[DatasetRef(DatasetType('deep_coadd', {band, skymap, tract, patch}, ExposureF), {band: 'i', skymap: 'lsst_cells_v1', tract: 5063, patch: 14}, run='LSSTComCam/runs/DRP/DP1/DM-51335', id=63cf9a40-44ce-40d8-b1a9-827f236c3a8a)]

Get the coadd image.

In [11]:

coadd = butler.get(center_coadd_datasetrefs[0])

Get the PSF model from the coadd info.

In [12]:

info_coadd = coadd.getInfo()
psf_coadd = info_coadd.getPsf()

Find the pixel location of the image using the World Coordinate System (WCS) function, obatined from the coadd info.

In [13]:

my_spherePoint = lsst.geom.SpherePoint(ra_cen*lsst.geom.degrees,
                                       dec_cen*lsst.geom.degrees)

In [14]:

wcs_coadd = info_coadd.getWcs()
point_image = wcs_coadd.skyToPixel(my_spherePoint)
point_image

Out[14]:

Point2D(14325.053060337099, 4570.427844687465)

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

In [15]:

psf_image_coadd = psf_coadd.computeImage(point_image).convertF()

In [16]:

fig = plt.figure(figsize=(6, 4))
afw_display = afwDisplay.Display(frame=fig)
afw_display.scale('asinh', 'zscale')
afw_display.mtv(psf_image_coadd)
plt.xlabel('x (pix)')
plt.ylabel('y (pix)')
plt.tick_params(axis='both', labelsize=9)
plt.gca().axis('on')
plt.show()

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

psf_image_coadd_kernel = psf_coadd.computeKernelImage(point_image).convertF()

In [18]:

fig = plt.figure(figsize=(6, 4))
afw_display = afwDisplay.Display(frame=fig)
afw_display.scale('asinh', 'zscale')
afw_display.mtv(psf_image_coadd_kernel)
plt.xlabel('x (pix)')
plt.ylabel('y (pix)')
plt.tick_params(axis='both', labelsize=9)
plt.gca().axis('on')
plt.show()

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

3. PSF model methods and properties¶

Option: explore the methods available for psf_coadd

In [19]:

# psf_coadd.

3.1. PSF size¶

In [20]:

psf_shape = psf_coadd.computeShape(point_image)

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

In [21]:

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

PSF size:  2.039 pix.
PSF FWHM  4.801 pix.

3.2. Second moments¶

Display the PSF second moments.

In [22]:

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:  4.210360 pix^2
PSF Iyy:  4.102 pix^2
PSF Ixy: -0.032 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 [23]:

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

In [24]:

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.349
  within 2σ:  0.741
  within 3σ:  0.898

3.4. Peak flux¶

Print the peak flux of the PSF kernel.

In [25]:

psf_peak = psf_coadd.computePeak(position=point_image)

In [26]:

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

PSF peak value at center:  0.036

3.5. Postage stamp dimensions¶

Calculate the PSF postage stamp dimensions, in pixels

In [27]:

psf_dims = psf_coadd.computeImage(point_image).getDimensions()

In [28]:

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

PSF postage stamp dimensions: (35, 35) pix

4. PSF radial profile¶

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

In [29]:

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

In [30]:

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

PSF center position: ( 17.5,  17.5) pix

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

In [31]:

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

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

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

Fit a Gaussian profile.

In [34]:

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

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

Sort the distances calculated above.

In [36]:

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

Compute the normalized cumulative flux.

In [37]:

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

Plot the curve of growth

In [38]:

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

In [39]:

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

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

Plot the contours.

In [41]:

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 - PSF in deep_coadd")
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 object catalog by setting the calib_psf_used flag to True, using the TAP service. For one star, display the residuals between the Piff PSF model at the star location and the observed star.

Define a 2-degree radius around the central point.

In [42]:

radius = 2.0

Define the query.

In [43]:

query = (
    "SELECT objectId, x, y, "
    "refExtendedness, "
    f"{my_band}_calib_psf_used, "
    f"{my_band}_pixelFlags_inexact_psfCenter, "
    f"{my_band}_calibFlux "
    "FROM dp1.Object "
    "WHERE refExtendedness = 0.0 "
    f"AND {my_band}_calib_psf_used = 1 "
    f"AND {my_band}_pixelFlags_inexact_psfCenter = 0 "
    f"AND {my_band}_calibFlux > 360 "
    "AND CONTAINS(POINT('ICRS', coord_ra, coord_dec), "
    f"CIRCLE('ICRS', {ra_cen}, {dec_cen}, {radius})) = 1 "
    "ORDER BY objectId"
)

Use the TAP service to get a table with objects used for PSF modeling.

In [44]:

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

print(len(psf_used_table))

1281

Choose one star in the field of view of the coadd image. First, get the bounds of the image.

In [46]:

box = coadd.getBBox()
coadd_min_x, coadd_max_x = box.getMinX(), box.getMaxX()
coadd_min_y, coadd_max_y = box.getMinY(), box.getMaxY()

In [47]:

for x, y in zip(psf_used_table['x'], psf_used_table['y']):
    if x > coadd_min_x and x < coadd_max_x and y > coadd_min_y and y < coadd_max_y:
        print(x, y)
        break

11985.316271125881 2910.939376113397

In [48]:

position_star = lsst.geom.Point2D(x, y)

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

In [49]:

star_cutout = coadd.getCutout(position_star, lsst.geom.Extent2I(xlen, ylen))

Get the observed star image.

In [50]:

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

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

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

In [52]:

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 (pix), y (pix)): {position_star}", 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).