204.3. Astrometric Calibration#

204_3_Astrometric_calibration

204.3. Astrometric Calibration¶

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 astrometric calibration in DP2, including how to access the World Coordinate System (WCS) for coadd images and centroid uncertainties.

LSST data products: deep_coadd.sky_projection, object.

Packages: lsst.daf.butler, lsst.geom, lsst.utils.plotting, astropy, numpy, matplotlib.

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¶

The World Coordinate System (WCS), or astrometric solution, defines how pixel positions in an image map to celestial coordinates—right ascension and declination (RA, Dec).

The LSST Science Pipelines perform astrometric calibration in two main stages.

An initial WCS is estimated using a static camera model and refined per detector and exposure by matching sources to Gaia DR3, achieving an accuracy of approximately 10-20 milliarcseconds—sufficient for associating sources across visits.
A more precise solution is derived using overlapping visits with the GBDES (Bernstein, Armstrong, Plazas et al., 2017) algorithm. This model fits two components: a per-detector static map (capturing camera distortions and the physical layout of the camera) and a per-visit dynamic map (accounting for time-varying effects like atmospheric refraction). The solution can optionally include proper motion, parallax, and corrections for differential chromatic refraction (DCR).

For DP2, the astrometric calibration follows the DP1 method with two main additions (see the DP2 paper, RTN-115). A joint calibration is performed on all visits in a given band that overlap a region of the sky—an order-3 HEALPix pixel in DP2, rather than the tract used in DP1—by associating the isolated point sources shared among the overlapping visits and matching them to Gaia DR3 with the gbdes package. In addition, the wave-like residuals left by atmospheric turbulence are modeled with a Gaussian Process (Léget et al. 2021) and combined with the gbdes model, significantly improving the solution (see DMTN-324). The proper motions and parallaxes of the calibration stars are fit jointly (reference epoch 2024.9); differential chromatic refraction and lateral color are not yet included and will be added in future data releases.

This is a Data Preview 2 (DP2) tutorial, using data from the LSST Science Camera (LSSTCam). Two aspects differ from the DP1 (LSSTComCam) version:

The WCS is stored as an AST-based SkyWcs object, not a plain FITS WCS. A FITS-standard representation is also written to the image headers so that tools such as astropy can read the WCS. For a coadd this FITS WCS is exact, because the coadd is built on an exact tangent-plane (TAN) projection; for a visit image it is only an approximation to the full SkyWcs, slightly less accurate near the edges of the LSSTCam focal plane (see DMTN-339).
DP2 coadds use the new lsst.images data model, in which the WCS is an lsst.images SkyProjection (with pixel_to_sky, sky_to_pixel, and as_fits_wcs methods). The projection can be read on its own with a Butler component query (deep_coadd.sky_projection), without loading the full image.

References

Related tutorials: The Monster reference catalog, used as the astrometric reference, is covered in another DP2 200-level calibration tutorial.

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 astropy, import units and celestial-coordinate handling (astropy.org).

From the lsst package, import the Butler, geometry utilities, and standard plotting helpers from the LSST Science Pipelines (pipelines.lsst.io).

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.coordinates import SkyCoord

from lsst.daf.butler import Butler
from lsst.images import SkyProjection, DetectorFrame, Box
from lsst.utils.plotting import (get_multiband_plot_colors,
                                 get_multiband_plot_linestyles)

1.2. Define parameters and functions¶

Instantiate the Butler with the DP2 repository and collection.

In [2]:

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

Get the standard per-band colors and linestyles.

In [3]:

filter_colors = get_multiband_plot_colors()
filter_linestyles = get_multiband_plot_linestyles()

2. WCS for a coadd image¶

In DP2 a deep_coadd is stored in the new lsst.images data model (CellCoadd), and its WCS is an lsst.images SkyProjection. This section shows how to obtain a SkyProjection and use its methods; the visit-image case (Section 3) uses the same interface.

Define a particular location on the sky and a band. Use a position in the ELAIS-S1 deep drilling field, in degrees.

In [4]:

ra_cen, dec_cen = 10.196, -44.134
my_band = 'i'

Find the tract and patch that contain the target position with a spatial Butler query on the patch dimension records (registry metadata only).

In [5]:

patch_records = butler.query_dimension_records(
    "patch",
    where="skymap = :skymap AND patch.region OVERLAPS POINT(:ra, :dec)",
    bind={"skymap": "lsst_cells_v2", "ra": ra_cen, "dec": dec_cen})
my_tract = patch_records[0].dataId["tract"]
my_patch = patch_records[0].dataId["patch"]
print("tract:", my_tract, "| patch:", my_patch)

tract: 2877 | patch: 34

Read only the sky_projection and bbox components of the deep_coadd from the Butler; the full image is not loaded.

In [6]:

data_id = {"band": my_band, "skymap": "lsst_cells_v2",
           "tract": my_tract, "patch": my_patch}
coadd_wcs = butler.get("deep_coadd.sky_projection", **data_id)
coadd_bbox = butler.get("deep_coadd.bbox", **data_id)

2.1. Convert between pixel and sky coordinates¶

The SkyProjection maps between pixel coordinates (in the tract pixel grid) and sky positions, which are astropy SkyCoord objects. Convert the target sky position to pixels, and back.

In [7]:

target = SkyCoord(ra=ra_cen*u.deg, dec=dec_cen*u.deg)
xy = coadd_wcs.sky_to_pixel(target)
xy

Out[7]:

XY(x=14656.538058947777, y=10504.271749935344)

In [8]:

coadd_wcs.pixel_to_sky(x=xy.x, y=xy.y)

Out[8]:

<SkyCoord (ICRS): (ra, dec) in deg
    (10.196, -44.134)>

2.2. Pixel scale¶

Compute the pixel scale (in arcseconds per pixel) by transforming two pixels one pixel apart and measuring their angular separation.

In [9]:

sky_00 = coadd_wcs.pixel_to_sky(x=xy.x, y=xy.y)
sky_10 = coadd_wcs.pixel_to_sky(x=xy.x + 1, y=xy.y)
pixel_scale = sky_00.separation(sky_10).to(u.arcsec)
print(f"Pixel scale: {pixel_scale.value: .4f} arcsec/pixel")

Pixel scale:  0.2000 arcsec/pixel

2.3. FITS WCS¶

A FITS-standard WCS can be produced with as_fits_wcs, given the image bounding box. For a coadd this FITS WCS is exact, because the coadd is built on an exact tangent-plane (TAN) projection (see Section 1).

In [10]:

fits_wcs_coadd = coadd_wcs.as_fits_wcs(coadd_bbox)
fits_wcs_coadd

Out[10]:

WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN' 'DEC--TAN'
CUNIT : 'deg' 'deg'
CRVAL : 10.169491525423737 -43.88429752066116
CRPIX : 3150.0 6150.0
CD1_1 CD1_2  : -5.5555555555553024e-05 0.0
CD2_1 CD2_2  : 0.0 5.555555555555367e-05
NAXIS : 3300  3300

3. WCS for a visit image¶

The per-detector WCS for every detector in a visit is collected in the visit_summary table. In DP2 the WCS stored there is a legacy lsst.afw.geom.SkyWcs; wrap it into an lsst.images SkyProjection with from_legacy to use the same interface as the coadd (Section 2). Once a SkyProjection is in hand, pixel_to_sky, sky_to_pixel, and the pixel-scale calculation work exactly as in Section 2.

Important: Visit images are not available in "Early DP2" but they will be available with the "Full DP2" release.

Note: Once the Full DP2 release is available (see RTN-011), the visit WCS will be accessible directly as visit_image.sky_projection (the same component-based access used for the coadd in Section 2), which will be the preferred method; the visit_summary approach shown here is used in the interim.

Find a visit and detector whose footprint overlaps the target position, using the visit_detector_region dimension records.

In [11]:

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: 2026010100038 | detector: 82

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

In [12]:

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

Wrap the legacy SkyWcs into an lsst.images SkyProjection with from_legacy. Build a DetectorFrame from the visit ID, detector ID, and bounding box (all available from the same visit_summary row).

In [13]:

wcs_visit = SkyProjection.from_legacy(
    legacy_wcs,
    DetectorFrame(instrument="LSSTCam", visit=visit_id,
                  detector=detector_id, bbox=detector_bbox))

3.1. FITS approximation¶

The pixel_to_sky, sky_to_pixel, and pixel-scale methods behave exactly as for the coadd (Section 2). The one difference is the FITS WCS: for a visit image the true SkyWcs cannot be represented exactly in FITS, so as_fits_wcs returns a SIP approximation and must be called with allow_approximation=True (see Section 1).

In [14]:

fits_wcs_visit = wcs_visit.as_fits_wcs(detector_bbox, allow_approximation=True)
fits_wcs_visit

Out[14]:

WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN-SIP' 'DEC--TAN-SIP'
CUNIT : 'deg' 'deg'
CRVAL : 10.195640728405449 -44.15321830027622
CRPIX : 2048.5 2002.5
CD1_1 CD1_2  : -4.5664821768378106e-05 -3.165696114864475e-05
CD2_1 CD2_2  : -3.1649358566663904e-05 4.567304467816654e-05
NAXIS : 4096  4004

4. Sky coordinates, errors, and astrometric accuracy¶

In the source and object catalogs, centroids are measured using the SdssCentroid algorithm, which operates in pixel coordinates. Once the centroid and its uncertainty are determined in pixels, the WCS transforms these measurements into right ascension and declination, and the uncertainties are propagated through this transformation.

4.1. Centroid uncertainties¶

Retrieve the centroid uncertainties and magnitudes for PSF stars from the object dataset for the tract. In DP2 the Object catalog is read directly from the Butler, which returns all objects in the tract.

In [15]:

columns = ["refExtendedness",
           f"{my_band}_calib_psf_used",
           f"{my_band}_pixelFlags_inexact_psfCenter",
           f"{my_band}_raErr", f"{my_band}_decErr", f"{my_band}_psfFlux"]
object_data = butler.get("object", skymap="lsst_cells_v2",
                         tract=my_tract, parameters={"columns": columns})

Select point-like sources used in PSF modeling, keeping only those with finite centroid errors and positive flux.

In [16]:

mask = ((np.asarray(object_data["refExtendedness"]) == 0.0)
        & np.asarray(object_data[f"{my_band}_calib_psf_used"])
        & ~np.asarray(object_data[f"{my_band}_pixelFlags_inexact_psfCenter"])
        & np.isfinite(np.asarray(object_data[f"{my_band}_raErr"]))
        & np.isfinite(np.asarray(object_data[f"{my_band}_decErr"]))
        & (np.asarray(object_data[f"{my_band}_psfFlux"]) > 0))

Convert the centroid errors from degrees to milliarcseconds using astropy units, and convert the PSF fluxes to AB magnitudes.

Converting flux to AB magnitude. Fluxes in the Object table are in nanojansky (nJy). Convert a flux $f$ to an AB magnitude with

$$m_{\mathrm{AB}} = -2.5\,\log_{10}\left(\frac{f}{3631\,\mathrm{Jy}}\right) = -2.5\,\log_{10}(f_{\mathrm{nJy}}) + 31.4,$$

where the constant $31.4$ is the AB zero point for fluxes in nanojansky ($3631\,\mathrm{Jy} = 3.631\times10^{12}\,\mathrm{nJy}$). Only positive fluxes have a defined magnitude—forced-photometry fluxes can be negative for faint sources—so the selection above keeps only sources with positive {band}_psfFlux.

In [17]:

ra_err = (np.asarray(object_data[f"{my_band}_raErr"])[mask] * u.deg).to(u.mas).value
dec_err = (np.asarray(object_data[f"{my_band}_decErr"])[mask] * u.deg).to(u.mas).value
psf_flux = np.asarray(object_data[f"{my_band}_psfFlux"])[mask]
mag = -2.5 * np.log10(psf_flux) + 31.4

Plot the distribution of the propagated uncertainties on the centroids.

In [18]:

max_err = max(ra_err.max(), dec_err.max())
error_bins = np.linspace(0, max_err, 50)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3), sharey=True)

for ax, err, label in zip((ax1, ax2), (ra_err, dec_err), ('RA Error', 'Dec Error')):
    _, _, patches = ax.hist(err, bins=error_bins, histtype='step',
                            linewidth=2, color=filter_colors[my_band], label=label)
    for patch in patches:
        patch.set_linestyle(filter_linestyles[my_band])
    ax.set_xlabel(f'{label} [mas]')
    ax.set_title(f'{label} Distribution')
    ax.legend()

ax1.set_ylabel('Count')
plt.tight_layout()
plt.show()

No description has been provided for this image

Figure 1: Right ascension and declination SdssCentroid error distributions for a subset of stars used in PSF modeling from the object dataset.

4.2. Centroid error vs magnitude¶

Visualize the centroid error as a function of magnitude. The centroid error is proportional to the seeing and inversely proportional to the signal-to-noise ratio of each object.

In [19]:

centroid_err = np.sqrt(ra_err**2 + dec_err**2)

In [20]:

plt.figure(figsize=(8, 5))
plt.scatter(mag, centroid_err, s=10, alpha=0.7)
plt.xlabel(f"{my_band} PSF Magnitude", fontsize=12)
plt.ylabel("Centroid Error (mas)", fontsize=12)
plt.title(f"Centroid Error vs. Magnitude in {my_band}-band", fontsize=14)
plt.grid(True)
plt.tight_layout()
plt.show()

Figure 2: SdssCentroid errors as a function of i-band magnitude.

5. Notes on astrometric accuracy¶

The plots in Section 4 use the propagated uncertainties on the centroids from the SdssCentroid fit, which are not indicative of the overall astrometric accuracy. To assess the performance of the complete astrometric solution, internal and external consistency tests are performed.

For internal consistency, the repeatability of position measurements for the same object (the RMS of the fit positions per object) is examined. For external consistency, the separation between sources not included in the astrometric fit and their counterparts in Gaia DR3 is measured, with a design requirement of 50 milliarcseconds for the main survey. Remaining residuals are due to distortions not yet included in the astrometric model (planned for future inclusion), such as atmospheric, camera, and detector-level distortions.

In [ ]: