304.1. Cosmological Applications of Point Spread Function Analysis#

304_1_PSF_cosmology

304.1. Cosmological Applications of Point Spread Function Analysis¶

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

Learning objective: To use coadded PSFs in cosmological analyses.

LSST data products: object, deep_coadd.

Packages: lsst.daf.butler, lsst.geom, treecorr, astropy, numpy, matplotlib, scipy.

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¶

Accurate modeling of the Point Spread Function (PSF) is a fundamental component of modern cosmological analyses. In particular, weak gravitational lensing (cosmic shear, galaxy-galaxy lensing, and galaxy cluster mass calibration) relies on precise knowledge of the PSF to infer subtle distortions in the shapes of background galaxies. Imperfections in PSF modeling can bias shear measurements and, consequently, cosmological parameter estimation.

This notebook explores diagnostic tools and validation techniques for coadd PSF modeling using DP2 (LSSTCam) data, focusing on the galaxy cluster PSZ2 G309.43-72.86. This Sunyaev-Zel'dovich-selected cluster, at a redshift of $z = 0.35$ and with a Planck SZ mass of $M_{500c} = 5.82 \times 10^{14}\,M_\odot$, lies in the ELAIS-S1 deep drilling field.

All diagnostics here use the coadd-level object catalog, which provides both the measured star shapes and the PSF model shapes at the star positions.

Related tutorials: DP2 tutorials on the PSF in deep_coadd images and in visit_images.

1.1. Import packages¶

Import numpy (numpy.org) and matplotlib (matplotlib.org).

Import treecorr, for two-point correlation functions (TreeCorr docs).

From scipy, import binned_statistic; from astropy, import SkyCoord and units; and from the LSST Science Pipelines, import the Butler and geometry utilities (pipelines.lsst.io).

In [1]:

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

from scipy.stats import binned_statistic
from matplotlib.patches import Circle

from astropy.coordinates import SkyCoord
import astropy.units as u

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

1.2. Define parameters and functions¶

Function: get_ellipticity

Compute the star and PSF model ellipticity components (e1, e2) from second moments.

In [2]:

def get_ellipticity(catalog, band, obj_type='star'):
    """Compute the ellipticity components (e1, e2) from second moments
    for stars or PSF models.

    Parameters
    ----------
    catalog : `astropy.table.Table`
        Catalog containing second moment columns for the
        specified band and object type.
    band : `str`
        Photometric band (e.g., 'i', 'r') used to construct the column names.
    obj_type : `str`, optional
        Type of object: 'star' (default) or 'PSF'.

    Returns
    -------
    e1 : `np.array`
        First ellipticity component, (Ixx - Iyy) / (Ixx + Iyy).
    e2 : `np.array`
        Second ellipticity component, 2 * Ixy / (Ixx + Iyy).
    """
    if obj_type == 'star':
        suffix = ''
    elif obj_type == 'PSF':
        suffix = 'PSF'
    else:
        raise ValueError("obj_type must be either 'star' or 'PSF'.")

    ixx = catalog[f'{band}_ixx{suffix}']
    iyy = catalog[f'{band}_iyy{suffix}']
    ixy = catalog[f'{band}_ixy{suffix}']
    denom = ixx + iyy
    e1 = (ixx - iyy) / denom
    e2 = 2. * ixy / denom
    return e1, e2

Function: compute_et

Compute the tangential ellipticity residual in radial bins around a point.

In [3]:

def compute_et(cat, e1_star, e2_star, e1_psf, e2_psf,
               center_ra, center_dec, r_bins):
    """Compute the tangential ellipticity residual in radial bins.

    Parameters
    ----------
    cat : `astropy.table.Table`
        Catalog with 'coord_ra' and 'coord_dec' columns (in degrees).
    e1_star, e2_star : `np.ndarray`
        Measured star ellipticity components.
    e1_psf, e2_psf : `np.ndarray`
        PSF model ellipticity components at the star positions.
    center_ra, center_dec : `float`
        Center coordinates, in degrees.
    r_bins : `array_like`
        Bin edges (in arcminutes).

    Returns
    -------
    mid : `np.ndarray`
        Midpoints of the radial bins (arcminutes).
    e_t_mean : `np.ndarray`
        Mean tangential ellipticity residual per bin.
    e_t_err : `np.ndarray`
        Standard error per bin.
    count : `np.ndarray`
        Number of objects per bin.
    """
    coord0 = SkyCoord(center_ra, center_dec, frame='icrs', unit='deg')
    coord1 = SkyCoord(cat['coord_ra'], cat['coord_dec'],
                      frame='icrs', unit='deg')

    pos_angle = coord0.position_angle(coord1).rad + np.pi / 2.
    r = coord0.separation(coord1).arcmin

    e1 = e1_star - e1_psf
    e2 = e2_star - e2_psf

    e_t = -e1 * np.cos(2. * pos_angle) - e2 * np.sin(2. * pos_angle)
    mid = 0.5 * (r_bins[1:] + r_bins[:-1])
    count, _, _ = binned_statistic(r, e_t, bins=r_bins, statistic='count')
    e_t_mean, _, _ = binned_statistic(r, e_t, bins=r_bins)
    e_t_std, _, _ = binned_statistic(r, e_t, bins=r_bins, statistic='std')
    e_t_err = e_t_std / np.sqrt(count)

    return mid, e_t_mean, e_t_err, count

Instantiate the Butler with the DP2 repository and collection, and get the skymap.

In [4]:

butler = Butler('dp2', collections=['dp2'])
assert butler is not None

In [5]:

skymap = butler.get('skyMap', skymap='lsst_cells_v2')

Define the cluster center (the Sunyaev-Zel'dovich coordinates of PSZ2 G309.43-72.86, in degrees) and the band.

In [6]:

ra_cl, dec_cl = 10.196, -44.134
my_band = 'i'

Find the tract containing the cluster center.

In [7]:

center_point = lsst.geom.SpherePoint(
    ra_cl*lsst.geom.degrees, dec_cl*lsst.geom.degrees)
my_tract = skymap.findTract(center_point).getId()
print('tract:', my_tract)

tract: 2877

2. PSF whisker plots from the `object` catalog¶

Whisker plots represent the shape and orientation of PSFs or galaxy ellipticities as small line segments ("whiskers") at their sky positions. Each whisker's direction encodes the ellipticity angle, and its length the ellipticity magnitude. They help visualize spatial patterns and systematics in the PSF across the field, which is essential for weak-lensing shear measurements.

Retrieve the object catalog for the tract, selecting the columns needed for the star and PSF model ellipticities and the flags for stars used in, and reserved from, PSF modeling.

In [8]:

columns = ["patch", "refExtendedness",
           f"{my_band}_calib_psf_used", f"{my_band}_calib_psf_reserved",
           f"{my_band}_pixelFlags_inexact_psfCenter",
           "coord_ra", "coord_dec",
           f"{my_band}_ixx", f"{my_band}_ixxPSF",
           f"{my_band}_ixy", f"{my_band}_ixyPSF",
           f"{my_band}_iyy", f"{my_band}_iyyPSF"]
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.

Define the angular separation of each object from the cluster center, and a base selection of point-like sources with a valid PSF centroid within 0.75 degrees of the center.

In [9]:

obj_coords = SkyCoord(np.asarray(objects["coord_ra"])*u.deg,
                      np.asarray(objects["coord_dec"])*u.deg)
separation = obj_coords.separation(
    SkyCoord(ra_cl*u.deg, dec_cl*u.deg)).deg

base = ((np.asarray(objects["refExtendedness"]) == 0.0)
        & ~np.asarray(objects[f"{my_band}_pixelFlags_inexact_psfCenter"])
        & (separation < 0.75))

Split into stars used in, and reserved from, PSF modeling. The reserved stars are about 10 percent of the total.

In [10]:

psf_used_table = objects[base & np.asarray(objects[f"{my_band}_calib_psf_used"])]
psf_reserved_table = objects[base & np.asarray(objects[f"{my_band}_calib_psf_reserved"])]
print(len(psf_used_table), len(psf_reserved_table))

2709 304

Calculate ellipticities and their Cartesian components for the used PSF models and observed stars, and the residuals.

For the PSF model, at the location of the used stars.

In [11]:

e1_psf_used, e2_psf_used = get_ellipticity(psf_used_table, my_band, obj_type='PSF')
e_psf_used = np.sqrt(e1_psf_used**2 + e2_psf_used**2)
theta_psf_used = 0.5 * np.arctan2(e2_psf_used, e1_psf_used)
cx_psf_used = np.asarray(e_psf_used * np.cos(theta_psf_used))
cy_psf_used = np.asarray(e_psf_used * np.sin(theta_psf_used))

For the observed stars used in PSF modeling.

In [12]:

e1_star_used, e2_star_used = get_ellipticity(psf_used_table, my_band, obj_type='star')
e_star_used = np.sqrt(e1_star_used**2 + e2_star_used**2)
theta_star_used = 0.5 * np.arctan2(e2_star_used, e1_star_used)
cx_star_used = np.asarray(e_star_used * np.cos(theta_star_used))
cy_star_used = np.asarray(e_star_used * np.sin(theta_star_used))

For the reserved stars and PSF models (used for validation below).

In [13]:

e1_psf_res, e2_psf_res = get_ellipticity(psf_reserved_table, my_band, obj_type='PSF')
e1_star_res, e2_star_res = get_ellipticity(psf_reserved_table, my_band, obj_type='star')

The ellipticity residuals for the used stars.

In [14]:

e1_residual_used = e1_star_used - e1_psf_used
e2_residual_used = e2_star_used - e2_psf_used
theta_residual_used = 0.5 * np.arctan2(e2_residual_used, e1_residual_used)
e_residual_used = np.sqrt(e1_residual_used**2 + e2_residual_used**2)
cx_residual_used = np.asarray(e_residual_used * np.cos(theta_residual_used))
cy_residual_used = np.asarray(e_residual_used * np.sin(theta_residual_used))

Define the plotting limits and the reference whisker, relative to the cluster center.

In [15]:

scale = 0.6
ref_norm = scale * 0.1
ref_x, ref_y = ra_cl + 0.5, dec_cl - 0.6
ra_limits = [ra_cl + 0.55, ra_cl - 0.55]
dec_limits = [dec_cl - 0.65, dec_cl + 0.65]
x_centroid = np.asarray(psf_used_table['coord_ra'])
y_centroid = np.asarray(psf_used_table['coord_dec'])

Prepare the whisker data. Switch the sign of the abscissa component so RA increases to the left.

In [16]:

whiskers = [
    (scale * -cx_star_used, scale * cy_star_used,
     f'Star ellipticity (from {my_band}_ixx, {my_band}_ixy, {my_band}_iyy)'),
    (scale * -cx_psf_used, scale * cy_psf_used,
     f'PSF model ellipticity (from {my_band}_ixxPSF, {my_band}_ixyPSF, {my_band}_iyyPSF)'),
    (scale * -cx_residual_used, scale * cy_residual_used,
     'Star - PSF ellipticity residuals')
]

Plot the three whisker panels side by side: measured star ellipticities, PSF model ellipticities, and star–PSF residuals, each with a reference whisker of e = 0.1.

In [17]:

fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 5))

for i, ax in enumerate(axes):
    cx, cy, title = whiskers[i]
    ax.quiver(x_centroid, y_centroid, cx, cy, angles='xy', color='black',
              scale_units='xy', scale=1, headlength=0, headwidth=0,
              headaxislength=0)
    ax.scatter([ra_cl], [dec_cl], marker='+', s=100, c='darkviolet')
    ax.add_patch(Circle((ra_cl, dec_cl), 0.5, color='darkviolet',
                        fill=False, linewidth=1))
    ax.quiver(ref_x, ref_y, ref_norm, 0, angles='xy', color='red',
              scale_units='xy', scale=1, headlength=0, headwidth=0,
              headaxislength=0, label=f'Ref: {ref_norm}')
    ax.text(ref_x + ref_norm / 2, ref_y + 0.02, 'e=0.1', color='red', ha='center')
    ax.invert_xaxis()
    ax.set_title(title)
    ax.set_xlabel('RA (deg)')
    ax.set_xlim(ra_limits)
    ax.set_ylim(dec_limits)
    if i == 0:
        ax.set_ylabel('Dec (deg)')

fig.tight_layout()

No description has been provided for this image

Figure 1: Whisker plot for measured stars used in PSF modeling (left), PSF models from Piff (Jarvis et al. 2021, center), and residuals (right) in a region with a half-degree radius around the center of PSZ2 G309.43-72.86.

Plot the histograms of the ellipticity residuals for used and reserved stars.

In [18]:

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 5))

ax[0].hist(np.asarray(e1_star_used - e1_psf_used), bins=30, range=[-0.04, 0.04],
           histtype='step', density=True, label='used')
ax[0].hist(np.asarray(e1_star_res - e1_psf_res), bins=30, range=[-0.04, 0.04],
           linestyle='--', histtype='step', density=True, label='reserved')
ax[0].set_xlabel(r'$\delta e_1 = $ Star e1 - PSF e1')
ax[0].legend()

ax[1].hist(np.asarray(e2_star_used - e2_psf_used), bins=30, range=[-0.04, 0.04],
           histtype='step', density=True, label='used')
ax[1].hist(np.asarray(e2_star_res - e2_psf_res), bins=30, range=[-0.04, 0.04],
           linestyle='--', histtype='step', density=True, label='reserved')
ax[1].set_xlabel(r'$\delta e_2 = $ Star e2 - PSF e2')
ax[1].legend()

plt.show()

Figure 2: Ellipticity residuals (e1 on the left; e2 on the right) for used and reserved stars.

3. PSF size correlation function from the `object` catalog¶

Calculate the two-point correlation function of the PSF size in the chosen band, using treecorr and the psf_used_table.

In [19]:

psf_size = np.sqrt((np.asarray(psf_used_table[f'{my_band}_ixxPSF'])
                    + np.asarray(psf_used_table[f'{my_band}_iyyPSF'])) / 2)
obj_ra = np.asarray(psf_used_table['coord_ra'])
obj_dec = np.asarray(psf_used_table['coord_dec'])

Build a treecorr.Catalog using the normalized PSF size (mean subtracted, in pixels) as the scalar field k, with angular positions in degrees; then run a KKCorrelation with separations between min_sep = 0.0001 deg and max_sep = 0.06 deg over 12 log-spaced bins to obtain xi (ξ), the two-point autocorrelation of PSF size as a function of angular separation θ.

In [20]:

cat = treecorr.Catalog(ra=obj_ra, dec=obj_dec,
                       k=psf_size - np.mean(psf_size),
                       ra_units='deg', dec_units='deg')
kk_config = {'max_sep': .06, 'min_sep': .0001, 'nbins': 12}
kk = treecorr.KKCorrelation(kk_config)
kk.process(cat)
xi = kk.xi
bins = kk.rnom

Plot the spatial map of PSF size across the field (left) and the corresponding two-point size correlation function ξ(θ) as a function of angular separation (right).

In [21]:

_, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.5),
                             gridspec_kw={'wspace': .3})

ax1.set_title(f'PSF size. Band: {my_band}')
scatter_plot = ax1.scatter(obj_ra, obj_dec, c=psf_size, s=1, cmap='cividis')
ax1.set_xlabel('RA [deg]')
ax1.set_ylabel('DEC [deg]')
plt.colorbar(scatter_plot, ax=ax1, label='PSF size [pixels]')

ax2.set_title(f'PSF size correlation. Band: {my_band}')
ax2.plot(np.degrees(bins), xi*1e4, 'o-', color='darkblue')
ax2.axhline(0, linestyle='--', color='lightgrey')
ax2.set_xscale('log')
ax2.set_ylabel(r'$\xi \times 10^{4}$', labelpad=10)
ax2.set_xlabel(r'$\theta$ [degree]')

plt.show()

Figure 3: Left: spatial distribution of PSF size (pixels, normalized by subtracting the mean) for a half-degree-radius region around the center of PSZ2 G309.43-72.86, using the object catalog. Right: PSF size correlation function in the same field.

4. Tangential PSF ellipticity residuals around the cluster center¶

An important systematic check for weak lensing is whether there are significant tangential ellipticity residuals around the cluster center, using both used and reserved stars. Residual tangential patterns in the stellar ellipticities—after subtracting the PSF model—can mimic or contaminate the tangential shear signal expected around massive clusters.

Calculate the angular separation of the used stars from the cluster center, and define radial bins.

In [22]:

r_used = SkyCoord(psf_used_table['coord_ra'], psf_used_table['coord_dec'],
                  unit='deg').separation(
    SkyCoord(ra_cl, dec_cl, unit='deg')).arcmin
r_max = np.max(r_used)
bins = np.linspace(0, r_max, 8)

Calculate the tangential residual for used stars.

In [23]:

mid_u, e_t_u, err_u, count_u = compute_et(
    psf_used_table, e1_star_used, e2_star_used, e1_psf_used, e2_psf_used,
    ra_cl, dec_cl, bins)

Calculate the tangential residual for reserved stars.

In [24]:

mid_r, e_t_r, err_r, count_r = compute_et(
    psf_reserved_table, e1_star_res, e2_star_res, e1_psf_res, e2_psf_res,
    ra_cl, dec_cl, bins)

Plot the profiles.

In [25]:

plt.figure(figsize=(6, 4))
plt.errorbar(mid_u, e_t_u, err_u, fmt='o', label='Used stars')
plt.errorbar(mid_r, e_t_r, err_r, fmt='s', label='Reserved stars')
plt.axhline(0, color='k', linestyle=':')
plt.xlabel('Separation [arcmin]', fontsize=14)
plt.ylabel(r'Residuals: $\langle e_t \rangle$', fontsize=14)
plt.ylim([-0.006, 0.006])
plt.legend(fontsize=12)
plt.tight_layout()
plt.show()

Figure 4: Tangential PSF ellipticity residuals around the center of PSZ2 G309.43-72.86, for used and reserved stars. Residuals consistent with zero indicate that the PSF modeling does not introduce a spurious tangential signal.

5. Exercises for the learner¶

Change my_band from 'i' to 'r' and rerun the notebook. Notice that the number of PSF stars and the amplitude of the PSF size correlation function change with the band.
Reduce the selection radius in Section 2 from 0.75 to 0.5 degrees and rerun. The whisker plots will cover a smaller area, and the number of stars will decrease.
In Section 4, change the number of radial bin edges from 8 to 5 (i.e., from 7 to 4 bins) and rerun. The tangential residual profile will have larger bins with smaller error bars per bin.

In [ ]: