103.9. Image stamps

103.9. Image stamps¶

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: 2026-03-30
Repository: github.com/lsst/tutorial-notebooks

Learning objective: Generate many (bulk) image stamps using the Rubin image cutout service.

LSST data products: DiaObject, Visit, visit_image, template_coadd, difference_image

Packages: matplotlib, numpy, lsst.rsp,lsst.afw, lsst.geom, lsst.resources, lsst.utils.plotting, pyvo.dal.adhoc, reproject

Credit: Originally developed by the Rubin Community Science team with feedback from Gregory Dubois-Feldman. This notebook also utilizes a transient detection from LSSTComCam first identified by Dan Taranu. 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 produce bulk image stamps using the Rubin cutout service.

The cutout service performs image cutouts remotely on the server using a protocol for remote data processing operations provided by the International Virtual Observatory Alliance (IVOA). The International Virtual Observatory Alliance (IVOA) co-ordinates the community efforts of astronomical missions and archives to develop and maintain the Virtual Observatory (VO) standards. The VO standards enable interoperability between astronomical archives.

IVOA provides the Server-side Operations for Data Access (SODA) protocol to provide these remote data processing operations. This protocol allows users to perform computations (pixel operations, image transformations, etc) on the remote server, which avoids unnecessary data movement. The LSST architecture has a "VO-first" approach, meaning that VO standards are implemented in all applicable services, enabling the use of VO tools such as the image cutout service to access LSST data.

The procedure is to identify the remote web location of the image of interest (called a datalink), and use the web service to create a cutout from the linked data remotely, before transferring the cutout to the user on the Rubin Science Platform.

This notebook demonstrates the default cutout service that should be used to generate bulk image stamps: cutout-sync, which returns a small data package including basic header information and the science image pixel values. This service is ideal for bulk image stamps because it minimizes data transfer by just returning the image extension and header of minimal metadata, and if not needing to use LSST Pipelines for their analysis. See tutorial notebook 103.4 for a demonstration of how to generate cutout exposures that are suitable for analysis using LSST science pipeline packages.

Further details and information can be found at the IVOA data link documentation, where it says Access Data Services. Rubin-specific documentation for these can also be found in this document describing the RSP DataLink service implementation strategy.

Related tutorials: This is the second of two 100-level tutorials demonstrating use of the image cutout service. Tutorial notebook 103.4 demonstrates the return of exposure cutouts, meaning as Rubin LSST image type ExposureF and MaskedImageF.

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 the lsst package, import modules for accessing the Table Access Protocol (TAP) service, image display functions, and the Simple Image Access (SIA) service from the LSST Science Pipelines (pipelines.lsst.io).

Functions from pyvo are imported to perform and parse results from Datalink and SODA queries. This is used to obtain images from the SIA service.

The reproject package is an astropy-adjacent package that provides access to coordinate transformations between images with world coordinate systems (WCS), and will be used to align visit images with different WCS.

In [1]:

import math
import numpy as np
import io
import getpass
import glob
import os
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

from lsst.rsp import get_tap_service
from lsst.rsp.service import get_siav2_service
from lsst.rsp.utils import get_pyvo_auth

import lsst.afw.display as afwDisplay
from lsst.afw.image import ImageF
from lsst.afw.fits import MemFileManager

from pyvo.dal.adhoc import DatalinkResults, SodaQuery

from astropy import units as u
from astropy.wcs import WCS
from astropy.io import fits

from reproject import reproject_interp

from IPython.display import Image as dimg
from PIL import Image

1.2. Define parameters and functions¶

Create an instance of the TAP service, and assert that it exists.

In [2]:

service = get_tap_service("tap")
assert service is not None

Instantiate the Simple Image Access service (SIA).

In [3]:

sia_service = get_siav2_service("dp1")
assert sia_service is not None

Set the afwDisplay package backend to matplotlib.

In [4]:

afwDisplay.setDefaultBackend('matplotlib')

1.2.1. Set up a temporary directory¶

The Rubin cutout service allows the user to save many image stamps as fits files locally on the Rubin Science Platform (RSP).

RSP users should save temporary files in the shared scratch directory, /deleted-sundays/, in sub-directories named with their RSP username.

Get the username and ensure an appropriate /deleted-sundays/ sub-directory exists.

In [5]:

username = getpass.getuser()
print('username: ', username)

userdir = '/deleted-sundays/' + username
if not os.path.exists(userdir):
    os.makedirs(userdir)
    print('Created ', userdir)
else:
    print('Directory already existed: ', userdir)

username:  christinawilliams
Directory already existed:  /deleted-sundays/christinawilliams

Within that folder, create a sub-folder named dp1_103_9_temp in which to save files created by this tutorial. At the end of the notebook, the last cell will clear the contents and remove this temporary folder.

In [6]:

tempdir = userdir + '/dp1_103_9_temp'
if not os.path.exists(tempdir):
    os.makedirs(tempdir)
    print('Created ', tempdir)
else:
    print('Directory already existed: ', tempdir)

Directory already existed:  /deleted-sundays/christinawilliams/dp1_103_9_temp

Delete the username and userdir, but keep tempdir to be used.

In [7]:

del username, userdir

1.2.2. Generate a gif from pngs¶

The function below will generate a gif to visualize all of the image stamp files in png format produced in the notebook.

In [8]:

def make_gif(folder):
    """
    Generate a GIF for all *.png images contained in a folder

    Parameters
    ----------
    folder : 'string'
        string containing the path to the folder
        default filename is animation.gif

    Returns
    -------
    """
    frames = [Image.open(img) for img in sorted(glob.glob(f"{folder}/*.png"))]
    frame_one = frames[0]
    frame_one.save(folder+"/animation.gif", format="GIF",
                   append_images=frames, save_all=True, duration=500, loop=1)

1.2.3. Generate cutouts¶

Define a function to get an image cutout from the cutout service.

In [9]:

def get_cutout(dl_result, ra, dec, session, fov):
    """
    Wrapper function to generate a cutout using the cutout tool.

    Parameters
    ----------
    dl_result : 'pyvo.dal.DatalinkResults' The Datalink result object
        containing the 'cutout-sync' ad-hoc service descriptor.
    ra, dec : 'float'
        the ra and dec of the cutout center
    session : 'requests.Session'
        An authenticated session object required to
        access proprietary or restricted Rubin data services.
    fov : 'float'
        the edge size of the cutout in decimal degrees

    Returns
    -------
    ImageF(mem) : 'lsst.afw.image.ImageF' 
        The cutout image data instantiated directly in memory via MemFileManager.
    hdul : 'astropy.io.fits.HDUList'
        The identical cutout data parsed into 
        an Astropy FITS object for standard header and extension manipulation.
    """

    sq = SodaQuery.from_resource(dl_result,
                                 dl_result.get_adhocservice_by_id("cutout-sync"),
                                 session=session)

    sq.circle = (ra * u.deg, dec * u.deg, fov / 2. * u.deg)

    cutout_bytes = sq.execute_stream().read()
    sq.raise_if_error()
    mem = MemFileManager(len(cutout_bytes))
    mem.setData(cutout_bytes, len(cutout_bytes))
    hdul = fits.open(io.BytesIO(cutout_bytes))
    return ImageF(mem), hdul

1.2.4. Set up subplots¶

Define a function to organize the display with matplotlib of bulk image stamps.

In [10]:

def make_subplot_grid(n_subplots, figsize_per_plot=(4, 3)):
    """
    Create an optimal grid of matplotlib subplots.

    Parameters
    ----------
    n_subplots : int
        The total number of subplots required.
    figsize_per_plot : tuple of float or int, optional
        The (width, height) in inches for each individual subplot. 
        The total figure size is scaled automatically based on this 
        and the calculated grid size. Default is (4, 3).

    Returns
    -------
    fig : matplotlib.figure.Figure
        The generated matplotlib figure object.
    axes : numpy.ndarray of matplotlib.axes.Axes
        A flattened, 1D array of exactly `n_subplots` axes objects, 
        ready for plotting.
    """
    n_cols = math.ceil(math.sqrt(n_subplots))
    n_rows = math.ceil(n_subplots / n_cols)

    figsize = (figsize_per_plot[0] * n_cols,
               figsize_per_plot[1] * n_rows)
    fig, axes = plt.subplots(
        n_rows,
        n_cols,
        figsize=figsize,
        gridspec_kw={
            "wspace": 0,
            "hspace": 0})

    axes = axes.flatten()

    for ax in axes:
        ax.axis("off")

    for ax in axes[n_subplots:]:
        ax.remove()

    return fig, axes[:n_subplots]

2. Create single-epoch triplet stamps¶

2.1. Define target¶

Type Ia Supernovae (SNIa) have homogenous light curves with similar rise and decay rates. A candidate SNIa search can be conducted in DP1 by looking for fading transients that exhibit dimming a rate consistent with SNIa: ~0.05 mag / day = ~5% flux fading per day (See Papadogiannakis et al. 2019).

Note that the 7-week duration of the commissioning survey with ComCam in DP1 (24 Oct 2024 - 11 Dec 2024) should sufficiently sample that fading emission from a SNIa.

Define the ra and dec of the SNIa candidate that was identified as a DIA source in DP1 imaging (see DP1 notebook tutorial 306.2).

In [11]:

ra = 53.124767650110215
dec = -27.739814591611168

A triplet for a DIA source refers to the set of three stamps made from the science, template, and difference images. They are typically displayed side-by-side and serve as a visualization of the detection.

This section uses the SIA service to find images and the image cutout service to create a "cutout triplet". Both services are demonstrated in 100-level tutorials.

Define circle, the position (RA and Dec) and radius (in degrees) to search for science (visit_image), reference (template_coadd), and difference (difference_image) images with SIA.

In [12]:

circle = (ra, dec, 0.001)
print(ra, dec)

53.124767650110215 -27.739814591611168

The cell below will perform a query to make sure the images retrieved span the time during which the DIAObject was changing and thus detectable as a DIASource.

In [13]:

dia_query = "SELECT MIN(midpointMjdTai) as first_detect, DiaObjectId, " \
            "MAX(midpointMjdTai) as last_detect " \
            "FROM dp1.DiaSource " \
            "WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), " \
            "CIRCLE('ICRS',"+str(ra)+","+str(dec)+",0.001)) = 1 "

dia_job = service.submit_job(dia_query)
dia_job.run()
dia_job.wait(phases=['COMPLETED', 'ERROR'])
dia_results = dia_job.fetch_result().to_table()

DiaObjectId = dia_results['DiaObjectId'][0]
first_detect = dia_results['first_detect'][0]
last_detect = dia_results['last_detect'][0]

2.2. Query the SIA service for images¶

Use the SIA service to search for calibration level 2 images (visit_images) that overlap the search circle.

In [14]:

results_sci = sia_service.search(pos=circle, calib_level=2)

lvl2_table = results_sci.to_table()

Filter the results for r-band visit_images, and to only return the visit images during which there was a DIA detection. The if statement ensures that a DIA detection was identified for that object. Lastly, sort the query results by observation time.

In [15]:

if (not np.ma.is_masked(first_detect)
        and first_detect is not None and not np.ma.is_masked(last_detect)
        and last_detect is not None):

    tx = np.where(
        (lvl2_table['dataproduct_subtype'] == 'lsst.visit_image')
        & (lvl2_table['lsst_band'] == 'r')
        & (lvl2_table['t_min'] >= first_detect)
        & (lvl2_table['t_min'] <= last_detect)
    )[0]

    print(f"Found {len(tx)} visits taken between MJD {first_detect: .2f} and {last_detect: .2f}.")

else:
    print("No DIA sources were ever detected in this region. Returning empty table.")
    tx = []

scitab = lvl2_table[tx]
scitab.sort('t_max')

Found 114 visits taken between MJD  60623.26 and  60655.24.

Use the SIA service to search for calibration level 3 images (template_coadd images and difference_images) that overlap the search circle. Filter the results for r-band template_coadd images and difference_images and sort the latter by observation time. No sorting is performed for the reference template images because there is only one image.

In [16]:

results = sia_service.search(pos=circle, calib_level=3)

lvl3_table = results.to_table()

tx = np.where((lvl3_table['dataproduct_subtype'] == 'lsst.template_coadd')
              & (lvl3_table['lsst_band'] == 'r'))[0]
reftab = lvl3_table[tx]

tx = np.where((lvl3_table['dataproduct_subtype'] == 'lsst.difference_image')
              & (lvl3_table['lsst_band'] == 'r'))[0]
difftab = lvl3_table[tx]
difftab.sort('t_max')

Print the date of the earliest r-band image, for which bulk image stamps will be retrieved.

In [17]:

print('begin date:', scitab['t_max'][0], difftab['t_max'][0])
print('end date:', scitab['t_max'][len(scitab['t_max'])-1])
print('number of visits:', len(scitab['t_max']))

begin date: 60623.25907695602 60623.25907695602
end date: 60655.242405046294
number of visits: 114

Retrieve the datalink VOTable documents using the access URLs from the first entries in each of the table results. The image products being retrieved for the science and difference images correspond to the earliest r-band image products taken.

In [18]:

dl_result_sci = DatalinkResults.from_result_url(scitab['access_url'][0], session=get_pyvo_auth())

dl_result_ref = DatalinkResults.from_result_url(reftab['access_url'][0], session=get_pyvo_auth())

dl_result_diff = DatalinkResults.from_result_url(difftab['access_url'][0], session=get_pyvo_auth())

2.3. Retrieve stamps¶

Retrieve bulk image stamps with an edge size of fov (in degrees). The get_cutout function returns both the image array (first parameter) and the corresponding Header Data Unit List (HDUL), a Python object used by astropy to contain one or more Header Data Units (HDUs) that contain astronomical image data, metadata, or tables.

Define the image cutout to be 0.004 degrees (~14 arcseconds) on a side.

In [19]:

fov = 0.004
sci, scihdul = get_cutout(dl_result_sci, ra, dec, get_pyvo_auth(), fov)
ref, refhdul = get_cutout(dl_result_ref, ra, dec, get_pyvo_auth(), fov)
diff, diffhdul = get_cutout(dl_result_diff, ra, dec, get_pyvo_auth(), fov)

First, visualize the cutouts of the science, template, and difference images by displaying them side by side.

In [20]:

fig, ax = plt.subplots(1, 3, figsize=(8, 3))
plt.sca(ax[0])
display1 = afwDisplay.Display(frame=fig)
display1.scale('linear', 'zscale')
display1.mtv(sci)
plt.title('science')
plt.sca(ax[1])
display2 = afwDisplay.Display(frame=fig)
display2.scale('linear', 'zscale')
display2.mtv(ref)
plt.title('template')
plt.sca(ax[2])
display3 = afwDisplay.Display(frame=fig)
display3.scale('linear', 'zscale')
display3.mtv(diff)
ax[0].set_axis_off()
ax[1].set_axis_off()
ax[2].set_axis_off()
plt.title('difference')
plt.tight_layout()
fig.suptitle('DiaObject '+str(DiaObjectId)+', MJD='+str(np.round(scitab['t_max'][0], 2)))

plt.show()

Figure 1: The cutout triplet for the first detection of this DiaObject, a candidate SN Ia.

2.4. Rotate visit image¶

As Figure 1 shows, visit images are not aligned with the template image (North up, East left) since they are taken at many different telescope orientations. This section demonstrates how to align the visit image to match the orientation of the template and difference cutouts using the reproject python package. reproject requires the world coordinate system (WCS) of the existing and desired coordinate systems, and the image pixels of the image to reproject, in order to perform the calculation.

The cell below extracts these from the HUDL. The structure of the ImageF object returned by the cutout service has two HDUL extensions. The extension 0 contains a primary header of telescope metadata and no image pixels, while the extension 1 contains the science header and the science image pixels.

In [21]:

sci_primary_header = scihdul[0].header
sci_hdu = scihdul[1]
sci_header = scihdul[1].header
scidata = scihdul[1].data
sciwcs = WCS(scihdul[1].header)

ref_primary_header = refhdul[0].header
ref_hdu = refhdul[1]
ref_header = refhdul[1].header
refdata = refhdul[1].data
refwcs = WCS(refhdul[1].header)

Use the reproject_interp function to put the science visit image onto the WCS of the template image.

In [22]:

reproj_scidata, footprint = reproject_interp((scidata, sciwcs), refwcs)

Plot the reprojected science visit image together with the template image to compare.

In [23]:

ax1 = plt.subplot(1, 2, 1, projection=WCS(ref_hdu.header))
ax1.imshow(reproj_scidata, origin='lower', vmin=np.nanpercentile(reproj_scidata, 1),
           vmax=np.nanpercentile(reproj_scidata, 99))
ax1.coords['ra'].set_axislabel('Right Ascension')
ax1.coords['dec'].set_axislabel('Declination')
ax1.set_title('Reprojected Visit Image')

ax2 = plt.subplot(1, 2, 2, projection=WCS(ref_hdu.header))
ax2.imshow(refdata, origin='lower', vmin=np.nanpercentile(refdata, 1),
           vmax=np.nanpercentile(refdata, 99))
ax2.coords['ra'].set_axislabel('Right Ascension')
ax2.coords['dec'].set_axislabel('Declination')
ax2.coords['dec'].set_axislabel_position('r')
ax2.coords['dec'].set_ticklabel_position('r')
ax2.set_title('Template Image')
plt.show()

Figure 2: The science visit image (left) after being reprojected onto the same WCS as the template image cutout (right).

3. Bulk image stamps¶

Bulk image stamps will be a key visualization tool for DIA objects. Below, demonstrate how to retrieve many image stamps (and save as png files). Since many visits span the time during which the target was a DIA source, limit the number to visualize to the first 20 (to save time). The full set of stamps can be visualized by replacing the first line with num_visits = len(scitab['t_max'])

In [24]:

num_visits = 20
print('number of visits:', num_visits)
for i in range(num_visits):
    plt.title(f"MJD {np.round(scitab['t_max'][i], 4)}")

    dl_result_sci = DatalinkResults.from_result_url(scitab['access_url'][i],
                                                    session=get_pyvo_auth())

    sci, scihdul = get_cutout(dl_result_sci, ra, dec, get_pyvo_auth(), fov)
    sci_hdu = scihdul[1]
    sci_header = scihdul[1].header
    scidata = scihdul[1].data
    sciwcs = WCS(scihdul[1].header)
    reproj_scidata, footprint = reproject_interp((scidata, sciwcs), refwcs)

    plt.imshow(reproj_scidata, origin='lower', vmin=np.nanpercentile(reproj_scidata, 1),
               vmax=np.nanpercentile(reproj_scidata, 99))

    figname = os.path.join(tempdir, 'stamp_' + str(i) + '.png')
    if os.path.isfile(figname):
        os.remove(figname)
    plt.savefig(figname)

number of visits: 20

Below, visualize the bulk image stamps that were saved to disk.

In [25]:

fig, axes = make_subplot_grid(num_visits)
for i, ax in enumerate(axes):

    j = i
    figname = os.path.join(tempdir, 'stamp_' + str(j) + '.png')
    img = mpimg.imread(figname)

    ax.imshow(img)

plt.tight_layout()
plt.show()

Figure 3: Image stamps from the science visit images of the DIA object, in chronological order, after being reprojected onto the same WCS as the template image (right).

4. Visualize cutouts with gif¶

Below, the cell will use the make_gif function defined in Section 1.2.3 to make a single gif movie out of all the pngs in the folder from the visit images in chronological order.

In [26]:

print(tempdir)
make_gif(tempdir)
dimg(data=open(tempdir+'/animation.gif', 'rb').read())

/deleted-sundays/christinawilliams/dp1_103_9_temp

Out[26]:

<IPython.core.display.Image object>

Figure 4: A gif animation generated from the reprojected visit image stamps in Figure 3.