103.8. Image display with matplotlib

103.8. Image display with matplotlib¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 1
Container Size: large
LSST Science Pipelines version: r29.2.0
Last verified to run: 2025-11-20
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, visit_image

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

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 afwDisplay and matplotlib. In fact, Firefly is the recommended backend, but matplotlib can also be used as the backend.

This tutorial demonstrates the pitfalls of using matplotlib, how use of extent with imshow is necessary, and how Right Ascension/Declination (RA/DEC) gridlines can appear to fail -- those are the reasons why matplotlib is not the recommended backend, but Firefly is.

The tutorial includes both visit_images, where the rotation can be tricky to interpret and RA/DEC gridlines can fail, and deep_coadds, where the pixel origin is not 1,1 and the extent is required.

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), and the lsst.geom package (for dealing with sky coordinates).

Import the matplotlib.pyplot sublibrary for plotting, and the astropy.wcs package for dealing with World Coordinate Systems (WCS). Import astropy.unit and astropy.coordinates for making celestial coordinates. Import the general package numpy for analysis.

Import 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 afwDisplay
import lsst.geom as geom
import matplotlib.pyplot as plt
from astropy.wcs import WCS
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):
    """
    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: lsst.afw.geom.SkyWcs
        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.

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

    cosdec = np.cos(np.deg2rad(dec))
    x, y = wcs.skyToPixel(geom.SpherePoint(ra, dec, geom.degrees))
    xN, yN = wcs.skyToPixel(geom.SpherePoint(ra, dec + clen/3600., geom.degrees))
    xE, yE = wcs.skyToPixel(geom.SpherePoint(ra + (clen/3600.)/cosdec, dec, geom.degrees))
    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):
    """
    Zoom in on the search location.

    Parameters
    ----------
    ax: matplotlib.axes._axes.Axes
        Axes for showing a plot.
    wcs: lsst.afw.geom.SkyWcs
        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.

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

    x, y = wcs.skyToPixel(geom.SpherePoint(ra, dec, geom.degrees))
    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])

Set the x/y-axis label properly based on the rotation angle of a visit_image.

In [5]:

def set_label(ax, rotation_angle):
    """
    Set the x/y-axis label properly based on the rotation angle of a visit image.

    Parameters
    ----------
    ax: matplotlib.axes._axes.Axes
        Axes for showing a plot.
    rotation_angle: float
        Rotation angle of the visit image in degree.

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

    case1 = (rotation_angle > 45) & (rotation_angle < 135)
    case2 = (rotation_angle > 225) & (rotation_angle < 315)
    if (case1 | case2):
        ax.set_ylabel('Right Ascension')
        ax.set_xlabel('Declination')
    else:
        ax.set_xlabel('Right Ascension')
        ax.set_ylabel('Declination')

1.2.2. Parameters¶

Set afwDisplay to use matplotlib.

In [6]:

afwDisplay.setDefaultBackend("matplotlib")

Instantiate the butler.

In [7]:

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

2. Display an image with AFWDisplay¶

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

In [8]:

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.

In [9]:

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])

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

fig, ax = plt.subplots()
afw_display = afwDisplay.Display(frame=fig)
afw_display.scale('linear', 'zscale')
afw_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. Axes labels are in pixels.

Read the API documentation about the above functions using the Jupyter notebook help() function.

In [11]:

help(afw_display.scale)

Help on method scale in module lsst.afw.display.interface:

scale(algorithm, min, max=None, unit=None, **kwargs) method of lsst.afw.display.interface.Display instance
    Set the range of the scaling from DN in the image to the image
    display.

    Parameters
    ----------
    algorithm : `str`
        Desired scaling (e.g. "linear" or "asinh").
    min
        Minimum value, or "minmax" or "zscale".
    max
        Maximum value (must be `None` for minmax|zscale).
    unit
        Units for min and max (e.g. Percent, Absolute, Sigma; `None` if
        min==minmax|zscale).
    **kwargs
        Optional keyword arguments to the backend.

In [12]:

help(afw_display.image)

Help on method image in module lsst.afw.display.interface:

image(data, title='', wcs=None, metadata=None) method of lsst.afw.display.interface.Display instance
    Display an image on a display, with semi-transparent masks
    overlaid, if available.

    Parameters
    ----------
    data : `lsst.afw.image.Exposure` or `lsst.afw.image.MaskedImage` or `lsst.afw.image.Image`
        Image to display; Exposure and MaskedImage will show transparent
        mask planes.
    title : `str`, optional
        Title for the display window.
    wcs : `lsst.afw.geom.SkyWcs`, optional
        World Coordinate System to align an `~lsst.afw.image.MaskedImage`
        or `~lsst.afw.image.Image` to; raise an exception if ``data``
        is an `~lsst.afw.image.Exposure`.
    metadata : `lsst.daf.base.PropertySet`, optional
        Additional FITS metadata to be sent to display.

    Raises
    ------
    RuntimeError
        Raised if an Exposure is passed with a non-None wcs when the
        ``wcs`` kwarg is also non-None.
    TypeError
        Raised if data is an incompatible type.

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

fig, ax = plt.subplots(1, 2, figsize=(14, 7))

plt.sca(ax[0])
display1 = afwDisplay.Display(frame=fig)
display1.scale('asinh', 'zscale')
display1.image(deep_coadd.image)

plt.sca(ax[1])
display2 = afwDisplay.Display(frame=fig)
display2.scale('linear',
               min=np.nanmin(deep_coadd.image.array),
               max=0.00001*np.nanmax(deep_coadd.image.array))
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.

In [14]:

print("Mask plane bit definitions:\n", display1.getMaskPlaneColor())

Mask plane bit definitions:
 {'BAD': 'red', 'CR': 'magenta', 'EDGE': 'yellow', 'INTERPOLATED': 'green', 'SATURATED': 'green', 'DETECTED': 'blue', 'DETECTED_NEGATIVE': 'cyan', 'SUSPECT': 'yellow', 'NO_DATA': 'orange', 'INTRP': 'green', 'SAT': 'green'}

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

fig, ax = plt.subplots(1, 2, figsize=(14, 7))
plt.sca(ax[0])
display1 = afwDisplay.Display(frame=fig)
display1.scale('linear', 'zscale')
display1.image(deep_coadd.image)
plt.sca(ax[1])
display2 = afwDisplay.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.

Repeat this for a visit_image. Query the butler, order the exposures by time, and retrieve the first exposure.

In [16]:

query = "band.name = :band AND visit_detector_region.region OVERLAPS POINT(:ra, :dec)"
bind = {'band': band, 'ra': ra, 'dec': dec}
visit_image_refs = butler.query_datasets("visit_image", where=query,
                                         bind=bind, order_by=["visit.timespan.begin"])
visit_image = butler.get(visit_image_refs[0])

Plot the visit_image.

In [17]:

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

Figure 4: At left, a visit_image. At right, the mask plane for the visit_image 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(). Show the deep_coadd first.

In [18]:

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

Figure 5: The deep_coadd with the mask plane overlaid.

Repeat this for the visit_image.

In [19]:

fig, ax = plt.subplots()
display = afwDisplay.Display(frame=fig)
display.scale('linear', 'zscale')
display.image(visit_image.maskedImage)
plt.show()
remove_figure(fig)

Figure 6: The visit_image 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. First show the deep_coadd.

In [20]:

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

Figure 7: 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.

Repeat it on the visit_image.

In [21]:

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

Figure 8: Similar to the previous figure of visit_image, 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 (or sources for visit_image), 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 [22]:

dataId = deep_coadd_refs[0].dataId

Retrieve the patch, the objects' x and y coordinates, and the r-band Point Spread Function (PSF) flux (with uncertainty).

In [23]:

use_columns = ['objectId', 'patch',
               'r_centroid_x', 'r_centroid_y',
               '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 [24]:

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

Number of objects:  11123

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

In [25]:

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:  410

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

fig, ax = plt.subplots()
afw_display = afwDisplay.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:
        afw_display.dot('o', objects[i]['r_centroid_x'],
                        objects[i]['r_centroid_y'],
                        size=30, ctype='orange')
plt.show()
remove_figure(fig)

Figure 9: 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.

Repeat the step for plotting sources on the visit_image. First, get the sources.

In [27]:

dataId = visit_image_refs[0].dataId

use_columns = ['sourceId', 'visit',
               'x', 'y',
               'psfFlux', 'psfFluxErr']

sources = butler.get('source', dataId=dataId,
                     parameters={'columns': use_columns})

print("Number of sources: ", len(sources))

Number of sources:  27708

Filter the sources by S/N.

In [28]:

sel = sources['psfFlux'] / sources['psfFluxErr'] > 300
sources = sources[sel]
print("Number of sources: ", len(sources))

Number of sources:  592

Display the image and the sources.

In [29]:

fig, ax = plt.subplots()
afw_display = afwDisplay.Display(frame=fig)
afw_display.scale('asinh', 'zscale')
afw_display.setMaskTransparency(100)
afw_display.image(visit_image.image)
with afw_display.Buffering():
    for i in sources:
        afw_display.dot('o', i['x'],
                        i['y'],
                        size=30, ctype='orange')
plt.show()
remove_figure(fig)

Figure 10: The same visit_image, but with orange circles marking the location of sources that have high S/N.

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

deep_coadd_extent = (deep_coadd.getBBox().beginX,
                     deep_coadd.getBBox().endX,
                     deep_coadd.getBBox().beginY,
                     deep_coadd.getBBox().endY)

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

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):  (11800, 15200, 2800, 6200)
Number of rows and columns of the image array:  (3400, 3400)

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

fig = plt.figure()
plt.subplot(projection=WCS(deep_coadd.getWcs().getFitsMetadata()))
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 11: The same deep_coadd displayed with the gray colormap at minimum -75 and maximum 125. Axes label and grid are in world coordinates.

3.2. Display the visit_image¶

Repeat the steps for the visit_image. Define and check the extent first.

In [33]:

visit_image_extent = (visit_image.getBBox().beginX,
                      visit_image.getBBox().endX,
                      visit_image.getBBox().beginY,
                      visit_image.getBBox().endY)

Print the extent, and note the use of extent with imshow is not necessary here, since the origin is 1,1 for the visit_image.

In [34]:

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

Extent (xmin, xmax, ymin, ymax):  (0, 4072, 0, 4000)
Number of rows and columns of the image array:  (4000, 4072)

In astronomy, the convention is to display images that are north-up. However, the camera obtains visits at a range of rotation angles. When the rotation angle of the visit_image is between 45 deg to 135 deg, or between 225 deg to 315 deg, the east/west direction in RA is along the y-axis, and the north/south direction in DEC is along the x-axis, which is different from the convention, and the x-axis label and y-axis label need to be switch to DEC and RA; otherwise, use the normal RA and DEC for the x- and y-axis labels.

Get the rotation angle and print the rotation angle. It is between 45 deg and 135 deg, which causes the difference of the RA and DEC directions from the convention.

In [35]:

rotation_angle = visit_image.getInfo().getVisitInfo().boresightRotAngle.asDegrees()
print("Rotation angle (deg): ", rotation_angle)

Rotation angle (deg):  102.79934100531

Plot the image. Set the x/y-axis label properly based on the rotation angle. Show a compass (the arm length is 50 arcseconds) in one panel. Show a zoom-in case (the box side is 140 arcseconds) in another panel.

In [36]:

fig = plt.figure(figsize=(14, 7))

ax1 = fig.add_subplot(121, projection=WCS(visit_image.getWcs().getFitsMetadata()))
ax2 = fig.add_subplot(122, projection=WCS(visit_image.getWcs().getFitsMetadata()))

for ax in [ax1, ax2]:
    im = ax.imshow(visit_image.image.array, cmap='gray',
                   vmin=-75, vmax=125,
                   extent=visit_image_extent,
                   origin='lower')
    ax.grid(color='white', ls='solid')
    set_label(ax, rotation_angle)

clen = 50
draw_compass(ax1, visit_image.getWcs(), ra, dec, clen)

side = 140
zoom_in(ax2, visit_image.getWcs(), ra, dec, side)

plt.tight_layout()
plt.show()
remove_figure(fig)

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

In addition, when the rotation angle is close to 45 deg, 135 deg, 225 deg, or 315 deg, the grid lines of RA/DEC can look similar to each other, which can be confusing. Therefore, it is recommended to use Firefly (presented in the related tutorials mentioned in Section 1).

3.3. Plot markers¶

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

In [37]:

fig = plt.figure(figsize=(14, 7))

ax1 = fig.add_subplot(121, projection=WCS(deep_coadd.getWcs().getFitsMetadata()))
ax2 = fig.add_subplot(122, projection=WCS(visit_image.getWcs().getFitsMetadata()))


im = ax1.imshow(deep_coadd.image.array, cmap='gray',
                vmin=-75, vmax=125,
                extent=deep_coadd_extent,
                origin='lower')
ax1.scatter(objects[tx]['r_centroid_x'], objects[tx]['r_centroid_y'],
            facecolors='none', edgecolors='orange')
ax1.set_xlabel('Right Ascension')
ax1.set_ylabel('Declination')

im = ax2.imshow(visit_image.image.array, cmap='gray',
                vmin=-75, vmax=125,
                extent=visit_image_extent,
                origin='lower')
ax2.scatter(sources['x'], sources['y'], edgecolors='orange', facecolors='none')
set_label(ax2, rotation_angle)

for ax in [ax1, ax2]:
    ax.grid(color='white', ls='solid')

plt.tight_layout()
plt.show()
remove_figure(fig)

Figure 13: The deep_coadd (left) and visit_image (right), with orange circles marking the location of objects (left) and sources (right) that have high S/N.

In [ ]: