103.7. Bulk cutout stamps#
103.7. Bulk cutout stamps¶
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-26
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20
Learning objective: How to use the Rubin image cutout service to make cutout exposures with DP2.
LSST data products: deep_coadd
Packages: lsst.images, lsst.rsp.utils, pyvo, lsst.rsp.RSPDiscovery, 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 notebook demonstrates how to produce bulk image stamps using the Rubin cutout service, using DP2 deep_coadd images. Details on the Rubin cutout service can be found in the introduction of notebook tutorial 103.6 on image cutout exposures.
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. See tutorial notebook 103.6 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.6 demonstrates the return of exposure cutouts. See also 104.5 that demonstrates how to generate subsets of images using the butler.
1.1. Import packages¶
Import common scientific analysis packages numpy, matplotlib and astropy. Import os and getpass to set up a temporary directory to store cutouts.
Import LSST Science Pipelines packages for image display lsst.afw.display, lsst.images for working with images, and utilities for remote data access from lsst.rsp.
Import pyvo packages for working with the virtual observatory cutout service.
import io
import matplotlib.pyplot as plt
import numpy as np
import getpass
import matplotlib.image as mpimg
import lsst.afw.display as afwDisplay
import lsst.images
from lsst.images.serialization import read_archive
from lsst.rsp import RSPDiscovery
from lsst.rsp.utils import get_pyvo_auth
from pyvo.dal.adhoc import SodaQuery
from astropy import units as u
from astropy.table import Table
1.2. Define parameters¶
Set the backend for afwDisplay to matplotlib.
afwDisplay.setDefaultBackend('matplotlib')
plt.rcParams.update({'font.size': 16})
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.
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: plazas Directory already existed: /deleted-sundays/plazas
Within that folder, create a sub-folder named dp2_103_7_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.
tempdir = userdir + '/dp2_103_7_temp'
if not os.path.exists(tempdir):
os.makedirs(tempdir)
print('Created ', tempdir)
else:
print('Directory already existed: ', tempdir)
Created /deleted-sundays/plazas/dp2_103_7_temp
Delete the username and userdir, but keep tempdir to be used.
del username, userdir
1.2.2. Generate cutouts¶
Define a function to get an image cutout from the cutout service.
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
-------
cutout : 'lsst.image'
The cutout image data instantiated directly in memory via read_archive.
"""
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()
cutout = read_archive(io.BytesIO(cutout_bytes))
return cutout
1.3. Initiate the SIA service¶
SIAv2 is an IVOA standard for querying and retrieving image data from astronomical archives. This is used to retrieve a datalink that uniquely identifies DP2 images (in the format of a web URL identifying where the data is hosted).
Load the RSPDiscovery tool, and use it to instantiate the SIA service.
discovery = RSPDiscovery("dp2")
sia_client = discovery.get_sia_client()
2. Find the coadd image¶
The cutout service needs the access_url for the image from which a cutout is desired.
The SIA or TAP services can be used to find the desired image and retrieve its access_url.
For this example, make an $r$-band (effective wavelength 622.1 nm) cutout centered on a set of coordinates in the ECDFS field.
Define the coordinates right ascension (target_ra) and declination (target_dec) in degrees, and the band via its effective wavelength.
target_ra = 53.1246023
target_dec = -27.7404715
eff_wl = 622.1e-09
2.1. Query for images with SIA¶
It is recommended to tightly constrain image queries, so that they return only the image data products needed for a given scientific analysis.
Define the search position as a 0.01 degree circle, centered on the target.
circle = (target_ra, target_dec, 0.01)
This query will return 1 deep_coadd (by design).
results = sia_client.search(pos=circle, calib_level=3,
dpsubtype='lsst.deep_coadd',
band=eff_wl)
print(len(results))
1
Display the results as an Astropy table.
results.to_table()
| dataproduct_type | dataproduct_subtype | facility_name | calib_level | target_name | obs_id | obs_collection | obs_publisher_did | access_url | access_format | s_resolution | s_xel1 | s_xel2 | t_xel | t_min | t_max | t_exptime | t_resolution | em_xel | em_min | em_max | em_res_power | em_filter_name | o_ucd | pol_xel | instrument_name | lsst_visit | lsst_detector | lsst_tract | lsst_patch | lsst_band | lsst_filter | obs_title | s_ra | s_dec | s_fov | s_region |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| arcsec | d | d | s | s | m | m | deg | deg | deg | |||||||||||||||||||||||||||
| object | object | object | int32 | object | object | object | object | object | object | float64 | int64 | int64 | int64 | float64 | float64 | float64 | float64 | int64 | float64 | float64 | float64 | object | object | int64 | object | int64 | int64 | int64 | int64 | object | object | object | float64 | float64 | float64 | object |
| image | lsst.deep_coadd | Rubin:Simonyi | 3 | lsst_cells_v2-5063-34 | LSST.DP2 | ivo://org.rubinobs/usdac/lsst-dp2?repo=dp2&id=019ed8bf-9458-7317-baf3-e92bca961d1e | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp2%3Frepo%3Ddp2%26id%3D019ed8bf-9458-7317-baf3-e92bca961d1e | application/x-votable+xml;content=datalink | -- | 3300 | 3300 | -- | -- | -- | -- | -- | -- | 5.51e-07 | 6.891e-07 | -- | r | phot.flux.density | -- | LSSTCam | -- | -- | 5063 | 34 | r | deep_coadd - r - tract=5063 patch=34 | 53.18170603020091 | -27.7705990503479 | 0.2592691563673628 | POLYGON ICRS 53.078100 -27.862296 53.285469 -27.862155 53.285139 -27.678826 53.078116 -27.678966 |
In the table, the access_url contains the web URL datalink for the image. This datalink will be needed to generate the image cutout.
2.2. Query for images with TAP¶
It is also possible to get the access_url from the ObsCore table in the TAP service.
An example TAP query using this method that retrieves the access_url is:
SELECT dataproduct_type,dataproduct_subtype,calib_level,lsst_band,em_min,em_max,lsst_tract,lsst_patch,
lsst_filter,lsst_visit,lsst_detector,t_exptime,t_min,t_max,s_ra,s_dec,s_fov,obs_id,
obs_collection,o_ucd,facility_name,instrument_name,obs_title,s_region,access_url,
access_format
FROM ivoa.ObsCore
WHERE CONTAINS(POINT('ICRS', 53.1567053, -27.7815854), s_region)=1
AND obs_collection = 'LSST.DP2' AND calib_level = 3
AND dataproduct_type = 'image' AND instrument_name = 'LSSTCam'
AND dataproduct_subtype = 'lsst.deep_coadd'
AND ( 622e-9 BETWEEN em_min AND em_max )
3. Generating an image cutout¶
Use the discovery.get_datalink_results method to create a DatalinkResults object to be able to access the datalink URL, which will be stored as dl_result and available for approximately 15 minutes, in a format that can be used by the IVOA tools below. The datalink is a VOTable document, stored as dl_result.
dl_result = discovery.get_datalink_results(results[0])
f"Datalink status: {dl_result.status}."
"Datalink status: ('OK', 'QUERY_STATUS not specified')."
Lastly, call the Rubin Image Cutout Service. In this example, the IVOA procedure cutout-sync is called using get_adhocservice_by_id. It is done by feeding the data link created above (called dl_result) to from_resource. Since the Rubin DP2 imaging is proprietary it is necessary to again provide the authorization for the current RSP session. Do this using the get_pyvo_auth function.
sq = SodaQuery.from_resource(dl_result,
dl_result.get_adhocservice_by_id("cutout-sync"),
session=get_pyvo_auth())
The variable sq now holds the result of the SODA query using the data link (which currently still points to the full LSST deep_coadd, at its remote location in the database). The cell below will now demonstrate how to extract a cutout from sq.
3.1. Define cutout center and edge¶
Only two shape definitions are supported: a circle function, and a polygon function can be used to define the cutout dimensions. These shape definitions do not produce circle or polygon cutouts, but rather are methods for defining the edges of cutouts with 4 sides. In the case of circle, the resulting cutout is always a square, with edge size that is the same as the circle diameter. Only cutouts with 4 corners and 90 degree angles are supported. See notebook 103.6 for more information on how to use these different definitions.
.circle defaults to assuming the units are degrees; this notebook demonstrates its use when specifying the units with astropy.
cutout_ra = target_ra * u.deg
cutout_dec = target_dec * u.deg
Radius = 0.01 * u.deg
sq.circle = (cutout_ra, cutout_dec, Radius)
cutout_bytes = sq.execute_stream().read()
sq.raise_if_error()
3.2. Retrieve the cutout¶
Read the cutout into memory using read_archive from the lsst.images.serialization package.
cutout = read_archive(io.BytesIO(cutout_bytes))
print(cutout)
Image([y=10862:11223, x=14229:14590], float32)
Display the image cutout. cutout-sync returns only the image extension of an lsst.images object and thus it is not required to append .image as demonstrated in 103.6.
display = afwDisplay.Display()
display.scale('asinh', 'zscale')
display.image(cutout)
plt.show()
Figure 1: The cutout image, displayed using LSST pipeline tools in grayscale with a scale bar at right.
Future planned options for the Rubin cutout service, including the potential to retrieve other image formats such as jpeg, are listed at the Rubin Science Platform image cutout implementation strategy document.
4. Bulk image stamps¶
Bulk image stamps will be a key visualization tool. Below, it is demonstrated how to retrieve many image stamps (and save as png files). Use the target list already used in tutorial notebook 104.5 for bright (V band magnitude < 22) variability selected galaxies hosting active galactic nuclei (AGN) in the ECDFS field from Boutsia et al. 2009. Define the targets as a dictionary containing the ID, ra, and dec, store it as an astropy table called ut1.
targets = [{'id': '3', 'ra': 53.335875, 'dec': -27.819472},
{'id': '8', 'ra': 53.194833, 'dec': -28.146306},
{'id': '9', 'ra': 53.411917, 'dec': -27.671222},
{'id': '11', 'ra': 53.318792, 'dec': -27.844306},
{'id': '12', 'ra': 53.191458, 'dec': -27.962583},
{'id': '14', 'ra': 53.133333, 'dec': -28.052750},
{'id': '15', 'ra': 52.977708, 'dec': -28.176583},
{'id': '18', 'ra': 53.380708, 'dec': -27.942833},
{'id': '20', 'ra': 53.049333, 'dec': -28.153056},
{'id': '22', 'ra': 53.370542, 'dec': -27.944750},
{'id': '24', 'ra': 53.290458, 'dec': -27.937222},
{'id': '25', 'ra': 52.819542, 'dec': -27.724889},
{'id': '26', 'ra': 53.144042, 'dec': -28.053889},
{'id': '28', 'ra': 53.337875, 'dec': -27.653278},
{'id': '29', 'ra': 53.155375, 'dec': -28.146389},
{'id': '30', 'ra': 53.302625, 'dec': -27.931000},
{'id': '31', 'ra': 52.797417, 'dec': -27.692194},
{'id': '32', 'ra': 53.344958, 'dec': -27.923278},
{'id': '33', 'ra': 53.224583, 'dec': -27.898361},
{'id': '34', 'ra': 53.359333, 'dec': -27.974917},
{'id': '35', 'ra': 53.084583, 'dec': -28.037444},
{'id': '36', 'ra': 53.132417, 'dec': -28.119556},
{'id': '37', 'ra': 53.371750, 'dec': -27.990750},
{'id': '39', 'ra': 53.184083, 'dec': -28.174583},
{'id': '41', 'ra': 52.992208, 'dec': -28.044861},
{'id': '43', 'ra': 53.333375, 'dec': -27.986778},
{'id': '44', 'ra': 52.812667, 'dec': -27.921833}]
ut1 = Table(targets)
Loop over the targets, use the SIA client to identify the datalink of the relevant deep_coadd, and use the get_cutout function defined in section 1 to generate the cutout. Then use imshow to internally display the image so it can be written to disk as a png file.
for i in range(len(ut1['id'])):
plt.title("ID = " + ut1['id'][i])
cutout_ra = ut1['ra'][i]
cutout_dec = ut1['dec'][i]
circle = (cutout_ra, cutout_dec, 0.001)
eff_wl = 622.1e-09
results_sci = sia_client.search(pos=circle, calib_level=3,
dpsubtype='lsst.deep_coadd',band=eff_wl)
results_sci.to_table()
dl_result = discovery.get_datalink_results(results_sci[0])
f"Datalink status: {dl_result.status}."
cutout = get_cutout(dl_result, cutout_ra, cutout_dec, get_pyvo_auth(), 0.01)
plt.imshow(cutout.array, origin='lower', vmin=np.nanpercentile(cutout.array, 1),
vmax=np.nanpercentile(cutout.array, 99))
figname = os.path.join(tempdir, 'stamp_' + str(i) + '.png')
if os.path.isfile(figname):
os.remove(figname)
plt.savefig(figname)
Figure 2: The last image stamp generated from the loop that was written to disk as a png file.
Define a function to organize the display with matplotlib of bulk image stamps.
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 = np.ceil(np.sqrt(n_subplots)).astype(int)
n_rows = np.ceil(n_subplots / n_cols).astype(int)
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]
Below, visualize the bulk image stamps that were saved to disk.
fig, axes = make_subplot_grid(len(targets))
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 generated from the
deep_coaddimages of each target.