103.5. Image display with matplotlib

103.5. Image display with matplotlib¶

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

Learning objective: How to use matplotlib as a backend for afwDisplay.

LSST data products: deep_coadd

Packages: lsst.daf.butler, lsst.afw.display, astropy.coordinates

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¶

This is a tutorial for displaying images with lsst.afw.display and matplotlib.

Firefly is the recommended backend for displaying LSST images, and is covered by another tutorial in this series.

Caveats. This tutorial demonstrates a pitfall of using matplotlib: the use of extent with imshow is necessary, because the pixel origin is not 1,1 for a deep_coadd. This is one of the reasons why matplotlib is not the recommended backend for image display; see the related Firefly tutorial for the recommended approach.

For Early Data Preview 2, this notebook covers deep_coadd display only, and visit_image support will be added once visit-level images become available in DP2.

Related tutorials: The 100-level series on the Butler demonstrates how to find and retrieve images, and the image display with Firefly tutorial.

1.1. Import packages¶

Import the Butler module from the lsst.daf.butler package, the display module from the lsst.afw package (for image display), the SkyCoord class from astropy.coordinates and the astropy.units module (for converting sky coordinates to pixel positions), the matplotlib.pyplot sublibrary for plotting, the general package numpy for analysis, and the gc (garbage collector) package to help clear the memory of large plots.

In [1]:

from lsst.daf.butler import Butler
import lsst.afw.display as afw_display
from lsst.images import get_legacy_deep_coadd_mask_planes
from astropy.coordinates import SkyCoord
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import gc

1.2. Define parameters and functions¶

1.2.1. Functions¶

matplotlib stores the data array associated with an image that is plotted. Since the LSST Charge-Coupled Device (CCD) detector images are large (~4k x 4k pixels), this can eventually lead to a memory overflow, which will cause the notebook kernel to die. To mitigate this issue, define a function to clean up after plotting them.

In [2]:

def remove_figure(fig):
    """
    Remove a figure to reduce memory footprint.

    Parameters
    ----------
    fig: matplotlib.figure.Figure
        Figure to be removed.

    Returns
    -------
    None
    """

    for ax in fig.get_axes():
        for im in ax.get_images():
            im.remove()

    fig.clf()
    plt.close(fig)
    gc.collect()

Define a function to draw a compass, showing the directions of North (N) and East (E), centered on the search location. The $cos(\rm{dec})$ term is included in the length of the compass arm for "East".

In [3]:

def draw_compass(ax, wcs, ra, dec, clen, x0=0, y0=0):
    """
    Draw a compass, showing the directions of North (N) and East (E),
      centered on the search location.

    Parameters
    ----------
    ax: matplotlib.axes._axes.Axes
        Axes for showing a plot.
    wcs: astropy.wcs.WCS
        WCS of an image.
    ra: float
        RA of a celestial location in degree.
    dec: float
        DEC of a celestial location in degree.
    clen: float
        Compass arm length in arcseconds.
    x0: float
        Pixel offset added to the WCS pixel coordinates, since
        world_to_pixel() returns coordinates local to the WCS rather
        than the bbox-based extent used to display the image.
    y0: float
        Pixel offset added to the WCS pixel coordinates, same as x0.

    Returns
    -------
    None
    """

    cosdec = np.cos(np.deg2rad(dec))
    x, y = wcs.world_to_pixel(SkyCoord(ra * u.deg, dec * u.deg))
    xN, yN = wcs.world_to_pixel(SkyCoord(ra * u.deg, (dec + clen / 3600.) * u.deg))
    xE, yE = wcs.world_to_pixel(SkyCoord((ra + (clen / 3600.) / cosdec) * u.deg, dec * u.deg))
    x, y = x + x0, y + y0
    xN, yN = xN + x0, yN + y0
    xE, yE = xE + x0, yE + y0
    ax.annotate('', xy=(xN, yN), xytext=(x, y),
                arrowprops=dict(color='yellow', linewidth=2, arrowstyle='->'))
    ax.text(xN, yN, 'N', color='yellow', fontsize=14)
    ax.annotate('', xy=(xE, yE), xytext=(x, y),
                arrowprops=dict(color='cyan', linewidth=2, arrowstyle='->'))
    ax.text(xE, yE, 'E', color='cyan', fontsize=14)

Define a function to zoom in on the search location.

In [4]:

def zoom_in(ax, wcs, ra, dec, side, x0=0, y0=0):
    """
    Zoom in on the search location.

    Parameters
    ----------
    ax: matplotlib.axes._axes.Axes
        Axes for showing a plot.
    wcs: astropy.wcs.WCS
        WCS of an image.
    ra: float
        RA of a celestial location in degree.
    dec: float
        DEC of a celestial location in degree.
    side: float
        Size of the display region in arcseconds.
    x0: float
        Pixel offset added to the WCS pixel coordinates, since
        world_to_pixel() returns coordinates local to the WCS rather
        than the bbox-based extent used to display the image.
    y0: float
        Pixel offset added to the WCS pixel coordinates, same as x0.

    Returns
    -------
    None
    """

    x, y = wcs.world_to_pixel(SkyCoord(ra * u.deg, dec * u.deg))
    x, y = x + x0, y + y0
    pixel_scale = 0.2
    side_pix = side / pixel_scale
    ax.set_xlim([x - side_pix / 2, x + side_pix / 2])
    ax.set_ylim([y - side_pix / 2, y + side_pix / 2])

1.2.2. Parameters¶

Set afwDisplay to use matplotlib.

In [5]:

afw_display.setDefaultBackend("matplotlib")

Instantiate the butler.

In [6]:

butler = Butler.from_config('dp2', collections='dp2')

2. Display an image with lsst.afw.display¶

Define an RA, Dec, and band (filter). These coordinates are near the center of the Extended Chandra Deep Field South (ECDFS).

In [7]:

ra = 53.076
dec = -28.110
band = 'r'

Define a query string using the coordinates and band as search constraints. Query the butler for matching deep_coadd images and retrieve the first on the list. Apply the 'pretty' background, which is tuned for visualization, since the raw pixel values otherwise still include an unsubtracted sky background.

In [8]:

query = "band.name = :band AND patch.region OVERLAPS POINT(:ra, :dec)"
bind = {'band': band, 'ra': ra, 'dec': dec}
deep_coadd_refs = butler.query_datasets("deep_coadd", where=query,
                                        bind=bind, order_by='patch')
deep_coadd = butler.get(deep_coadd_refs[0])
deep_coadd.apply_background('pretty')

Display the retrieved image in matplotlib. Define afw_display to show images. Recommend to set the scale to linear stretch and use the automatic algorithm zscale to select the white and black thresholds. Remove the underlying data from memory after creating and displaying the image.

In [9]:

fig, ax = plt.subplots()
display = afw_display.Display(frame=fig)
display.scale('linear', 'zscale')
display.image(deep_coadd.image)
plt.show()
remove_figure(fig)

Figure 1: A deep_coadd with the linear and zscale scaling and without the mask overlay, displayed in grayscale with a scale bar at right.

Option to read the API documentation about the above functions using the Jupyter notebook help() function.

In [10]:

# help(display.scale)

In [11]:

# help(display.image)

2.1. Scaling options¶

There are other options for scaling. Show the same image using the asinh and zscale scaling, and linear with explicit minimum and maximum values. Using plt.tight_layout() with multi-axis figures helps to avoid axis overlap or excessive white spaces and results in a nicer-looking plot.

In [12]:

fig, ax = plt.subplots(1, 2, figsize=(14, 7))
plt.sca(ax[0])
display1 = afw_display.Display(frame=fig)
display1.scale('asinh', 'zscale')
display1.image(deep_coadd.image)
plt.sca(ax[1])
display2 = afw_display.Display(frame=fig)
display2.scale('linear', min=-10., max=20)
display2.image(deep_coadd.image)
plt.tight_layout()
plt.show()
remove_figure(fig)

Figure 2: Left: The same deep_coadd with the asinh and zscale scaling. Right: The linear scaling with explicit minimum and maximum values.

2.2. Manipulate the mask display¶

Each image returned by the butler contains more than just the image pixel values. One other component is the mask associated with the image. A mask is composed of a set of "mask planes", 2D binary bit maps corresponding to pixels that are masked for various reasons.

The display framework renders each plane (each bit) of the mask in a different color. View the mapping between the mask plane and the color using the following code.

Notice: One of the changes with DP2's new image types is more informative mask names. However, the old mask names and colors are still being used by afw.display. The following is a temporary workaround to print new names with the colors defined by display. Once this is fixed, mask plane definitions will be able to be printed with print("Mask plane bit definitions:\n", display1.getMaskPlaneColor()).

In [13]:

mask = deep_coadd.mask
new_to_old_convert = {v.name: k for k, v in get_legacy_deep_coadd_mask_planes().items()}
for new_name in mask.schema.names:
    print('{} {}'.format(new_name,
                         display1.getMaskPlaneColor(new_to_old_convert[new_name])))

NO_DATA orange
INTERPOLATED green
COSMIC_RAY magenta
SATURATED green
DETECTION_EDGE yellow
CLIPPED None
REJECTED None
DETECTED blue
INEXACT_PSF None

Plot the image and mask plane side-by-side using matplotlib subplots.

Use plt.sca(ax[0]) to set the first axis as current, and then plt.sca(ax[1]) to switch to the second axis.

Plot the deep_coadd.

In [14]:

fig, ax = plt.subplots(1, 2, figsize=(14, 7))
plt.sca(ax[0])
display1 = afw_display.Display(frame=fig)
display1.scale('linear', 'zscale')
display1.image(deep_coadd.image)
plt.sca(ax[1])
display2 = afw_display.Display(frame=fig)
display2.image(deep_coadd.mask)
plt.tight_layout()
plt.show()
remove_figure(fig)

Figure 3: At left, the same deep_coadd. At right, the mask plane for the deep_coadd illustrates which pixels have a mask value.

The afwDisplay package also provides a nice interface for plotting the mask on top of the image using the maskedImage. This gives the same result as passing the image itself to display.image().

In [15]:

fig, ax = plt.subplots()
display = afw_display.Display(frame=fig)
display.scale('linear', 'zscale')
display.image(deep_coadd)
plt.show()
remove_figure(fig)

Figure 4: The deep_coadd with the mask plane overlaid.

To investigate the mask in a bit more detail, follow the same steps as above to display the image, but add a few modifications. Set the transparency of the overplotted mask to 10% (0 = opaque, 100 = transparent). Change the color of the 'DETECTED' mask plane from 'blue' (as above) to 'green' (i.e., all pixels associated with detected objects). Pass the full image object to display.image() instead of just the image plane.

In [16]:

fig, ax = plt.subplots()
display = afw_display.Display(frame=fig)
display.scale('asinh', 'zscale')
display.setMaskTransparency(10)
display.setMaskPlaneColor('DETECTED', 'green')
display.image(deep_coadd)
plt.show()
remove_figure(fig)

Figure 5: Similar to the previous figure of deep_coadd, but with the mask transparency reduced to 10% (mostly opaque), and the color representing 'DETECTED' pixels changed from blue to green.

2.3. Plot markers¶

In general the Table Access Protocol（TAP）service is the recommended access mechanism for catalog data.

However, in the case where all the objects detected in a given deep_coadd are desired, the Butler offers a quick way to retrieve them.

Use the dataId for the deep_coadd image to retrieve the object table for the tract.

In [17]:

dataId = deep_coadd_refs[0].dataId

Retrieve the patch, the objects' r-band sky coordinates, and the r-band Point Spread Function (PSF) flux (with uncertainty). The object table stores position as r_ra/r_dec rather than pixel coordinates, so these are converted to pixel positions using the deep_coadd WCS before plotting.

In [18]:

use_columns = ['objectId', 'patch',
               'r_ra', 'r_dec',
               'r_psfFlux', 'r_psfFluxErr']
objects = butler.get('object', tract=dataId.get('tract'),
                     parameters={'columns': use_columns})

Identify objects in the patch of the displayed deep_coadd image.

In [19]:

tx = np.where(objects['patch'] == dataId.get('patch'))[0]
print("Number of objects: ", len(tx))

Number of objects:  9058

Since there are many objects, filter objects by their Signal-to-Noise ratio (S/N) to make the following scatter plot more clear.

In [20]:

sel = objects['r_psfFlux'] / objects['r_psfFluxErr'] > 300
objects = objects[sel]
tx = np.where(objects['patch'] == dataId.get('patch'))[0]
print("Number of objects: ", len(tx))

Number of objects:  138

Redisplay the deep_coadd with no mask. Use buffering to display orange circles at the location of every object (this avoids re-drawing the image after each object is plotted.).

In [21]:

fig, ax = plt.subplots()
afw_display = afw_display.Display(frame=fig)
afw_display.scale('asinh', 'zscale')
afw_display.setMaskTransparency(100)
afw_display.image(deep_coadd.image)
with afw_display.Buffering():
    for i in tx:
        sky_point = SkyCoord(objects[i]['r_ra'] * u.deg, objects[i]['r_dec'] * u.deg)
        pix = deep_coadd.sky_projection.sky_to_pixel(sky_point)
        afw_display.dot('o', pix.x, pix.y, size=30, ctype='orange')
plt.show()
remove_figure(fig)

Figure 6: The same deep_coadd, but with orange circles marking the location of objects that have high S/N.

Why are no objects near the edges marked?

The deep_coadd images are per-patch, and patches and deep_coadd images overlap at their edges. The object table is by tract, and has no duplicates -- there is only one row per detected object. For every object the patch column is the patch for which they are closest to the center. The stars and galaxies near the edges of the displayed deep_coadd image are listed as belonging to the adjacent patch.

3. Display an image with imshow¶

In order to display the image axes in sky coordinates, use matplotlib.pyplot's subplot, imshow, and grid functions, along with astropy's WCS package.

3.1. Display the deep_coadd¶

The extent parameter in imshow defines the bounding box for the image to show. Define the extent using the boundary of BBox.

In [22]:

deep_coadd_extent = (deep_coadd.bbox.x.start, deep_coadd.bbox.x.stop,
                     deep_coadd.bbox.y.start, deep_coadd.bbox.y.stop)

Check the extent. Note the use of extent with imshow is necessary, because the pixel origin is not 1,1 for the deep_coadd.

In [23]:

print("Extent (xmin, xmax, ymin, ymax): ", deep_coadd_extent)
print("Number of rows and columns of the image array: ", np.shape(deep_coadd.image.array))

Extent (xmin, xmax, ymin, ymax):  (14850, 18150, 2850, 6150)
Number of rows and columns of the image array:  (3300, 3300)

Set the figure's projection to be the WCS of the image. Define the extent in pixel coordinates using the bounding box. Display the image data array using the gray colormap (cmap). Add solid white grid lines. Label the axes, and show the plot. Remove the underlying data from memory.

In [24]:

fig = plt.figure()
plt.subplot(projection=deep_coadd.fits_wcs)
im = plt.imshow(deep_coadd.image.array, cmap='gray',
                vmin=-75, vmax=125,
                extent=deep_coadd_extent,
                origin='lower')
plt.grid(color='white', ls='solid')
plt.xlabel('Right Ascension')
plt.ylabel('Declination')
plt.show()
remove_figure(fig)

Figure 7: The same deep_coadd displayed with the gray colormap at minimum -75 and maximum 125. Axes label and grid are in world coordinates. Notice that RA is in h:m:s, not degrees as in previous figures, when the fits_wcs is used for the sky projection.

3.2. Compass and zoom-in¶

Plot the image again. Show a compass (the arm length is 50 arcseconds) in one panel, and a zoom-in case (the box side is 140 arcseconds) in another panel.

In [25]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 7),
                               subplot_kw={'projection': deep_coadd.astropy_wcs})
for ax in [ax1, ax2]:
    im = ax.imshow(deep_coadd.image.array, cmap='gray',
                   vmin=-75, vmax=125,
                   extent=deep_coadd_extent,
                   origin='lower')
    ax.coords[0].set_major_formatter('hh:mm:ss.s')
    ax.coords[1].set_major_formatter('dd:mm:ss')
    ax.grid(color='white', ls='solid')
    ax.set_xlabel('Right Ascension')
    ax.set_ylabel('Declination')
x0 = deep_coadd.bbox.x.start
y0 = deep_coadd.bbox.y.start
clen = 50
draw_compass(ax1, deep_coadd.astropy_wcs, ra, dec, clen, x0=x0, y0=y0)
side = 140
zoom_in(ax2, deep_coadd.astropy_wcs, ra, dec, side, x0=x0, y0=y0)
plt.show()
remove_figure(fig)

Figure 8: The same deep_coadd displayed with the gray colormap at minimum -75 and maximum 125. Axes label and grid are in world coordinates. At left, the full deep_coadd with a compass. At right, a zoom-in of the deep_coadd on the search location.

3.3. Plot markers¶

Plot markers onto an image displayed with imshow. Repeat the previous to display the deep_coadd with matplotlib, and use a scatter plot to show objects.

In [26]:

fig = plt.figure()
ax = fig.add_subplot(111, projection=deep_coadd.astropy_wcs)
im = ax.imshow(deep_coadd.image.array, cmap='gray',
               vmin=-75, vmax=125,
               extent=deep_coadd_extent,
               origin='lower')
sky_points = SkyCoord(objects[tx]['r_ra'] * u.deg, objects[tx]['r_dec'] * u.deg)
pix_x, pix_y = deep_coadd.astropy_wcs.world_to_pixel(sky_points)
pix_x = pix_x + deep_coadd.bbox.x.start
pix_y = pix_y + deep_coadd.bbox.y.start
ax.scatter(pix_x, pix_y, facecolors='none', edgecolors='orange')
ax.coords[0].set_major_formatter('hh:mm:ss.s')
ax.coords[1].set_major_formatter('dd:mm:ss')
ax.set_xlabel('Right Ascension')
ax.set_ylabel('Declination')
ax.grid(color='white', ls='solid')
plt.show()
remove_figure(fig)

Figure 9: The deep_coadd, with orange circles marking the location of objects that have high S/N.

In [ ]: