103.6. Cutout Exposures

103.6. Cutout Exposures¶

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: 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 use the image cutout service with DP2 images. 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 two types of cutout services: cutout-sync-exposure, which returns the full set of metadata and image extensions that are contained in the deep_coadd data type as a lsst.images.cells.CellCoadd storage class (hereafter CellCoadd; Section 3), and cutout-sync-maskedimage which returns the mask extension with the science and variance images, but no extended metadata (lsst.images.MaskedImage storage class; hereafter MaskedImage; Section 4). These two cutout services are ideal for generating a small number of cutouts if needing to use LSST pipelines to run analysis on the cutouts. Another cutout service cutout-sync exists to generate bulk cutouts, which minimizes data transfer by just returning only the image extension and header of minimal metadata (lsst.images.Image storage class). This service is best for users who do not need to use LSST Pipelines for their analysis. See tutorial notebook 103.9 for a demonstration of how to generate bulk cutouts with the cutout service.

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.

To identify the remote location of the image, use Simple Image Access (SIA; see 100-level tutorial notebook on the SIA service). SIA is a protocol of the International Virtual Observatory Alliance (IVOA). It provides a standardized model for image metadata, and the capability to query and retrieve image datasets. Learn more in the IVOA SIA documentation.

Related tutorials: See also the 103-series tutorial on image stamps, which demonstrates how to perform bulk cutouts using the cutout service that returns only the image extension to minimize the data transferred.

1.1. Import packages¶

Import common scientific analysis packages numpy and astropy.

Import LSST Science Pipelines packages for image display lsst.afw.display, and utilities for remote data access from lsst.rsp.

Import pyvo packages for working with the virtual observatory cutout service.

In [1]:

import io
import matplotlib.pyplot as plt
import numpy as np

import lsst.afw.display as afwDisplay
import lsst.geom as geom
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.coordinates import Angle
from astropy.wcs import WCS

1.2. Define parameters¶

Set the backend for afwDisplay to matplotlib.

In [2]:

afwDisplay.setDefaultBackend('matplotlib')

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.

In [3]:

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.

In [4]:

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.

In [5]:

circle = (target_ra, target_dec, 0.01)

This query will return 1 deep_coadd (by design).

In [6]:

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.

In [7]:

results.to_table()

Out[7]:

Table length=1

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.

In [8]:

dl_result = discovery.get_datalink_results(results[0])

f"Datalink status: {dl_result.status}."

Out[8]:

"Datalink status: ('OK', 'QUERY_STATUS not specified')."

Lastly, call the Rubin Image Cutout Service. This section will demonstrate cutout-sync-exposure, which returns the full set of metadata and image extensions that are available for deep_coadd data type, as contained in the CellCoadd storage class. The return has the same format and data structure as the parent CellCoadd.

To use the cutout service in this example, the IVOA procedure cutout-sync-exposure 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.

In [9]:

sq = SodaQuery.from_resource(dl_result,
                             dl_result.get_adhocservice_by_id("cutout-sync-exposure"),
                             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. In the case of a polygon, either a square or a rectangular cutout will result, depending on whether the length and width edge dimensions are different values. Only cutouts with 4 corners and 90 degree angles are supported.

.circle defaults to assuming the units are degrees; this notebook demonstrates its use when specifying the units with astropy.

In [10]:

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 in the CellCoadd format using read_archive from the lsst.images.serialization package.

In [11]:

cutout = read_archive(io.BytesIO(cutout_bytes))
print(cutout)

CellCoadd([y=10862:11223, x=14229:14590], tract=5063)

Display the image extension of the CellCoadd cutout by appending .image to the exposure. Similarly, the variance can be displayed appending .variance. To display the image and mask together, pass the cutout variable directly to display.image.

In [12]:

display = afwDisplay.Display()
display.scale('asinh', 'zscale')
display.image(cutout.image)
plt.show()

Figure 1: The cutout image, displayed in pixel coordinates using LSST pipeline tools in grayscale with a scale bar at right.

3.2.1. Option to save cutout to disk¶

The cutout can be saved to disk as a FITS file with the f.write() function. The following commands demonstrate how to save it in a temporary folder in the user's home directory.

In [13]:

# tempdir = os.path.join(os.getenv('HOME'), 'cutouts_temp/')
# if not os.path.exists(tempdir):
#     os.makedirs(tempdir)
#     print('Created ', tempdir)
# else:
#     print('Directory already existed: ', tempdir)

# sodaCutout = os.path.join(tempdir, 'cutout-circle.fits')
# with open(sodaCutout, 'bw') as f:
#     f.write(sq.execute_stream().read())

Note: Each CellCoadd contains multiple Header-Data Units (HDUs) with extensive exposure metadata. To reduce file size, one option is to save only the image plane, along with the WCS recorded in the Primary HDU. This can be done by writing cutout.image into a FITS file using lsst.images.fits.write.

In [14]:

# sodaCutout_small = os.path.join(tempdir, 'cutout-circle_small.fits')

# if os.path.exists(sodaCutout_small):
#     os.remove(sodaCutout_small)

# fits_write = lsst.images.fits.write(cutout.image, sodaCutout_small)

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.

3.3. Use polygon to define shape¶

It is also possible to define the cutout geometry using a polygon, which enables the cutout to be rectangular, but not necessarily be square. For this, use polygon, which takes as input the four corners in celestial coordinates. A minimum of 3 vertices are required (the line from the last vertex back to the first is implicit). Vertices must be ordered in the counter-clockwise direction. For example: a polygon is defined as a set of 4 (x,y) coordinates from (12,34) to (14,34) to (14,36) to (12,36) and (implicitly) back to (12,34) as:

POLYGON=12 34 14 34 14 36 12 36

Since the center of the cutout is already defined in unit degrees in the cells above, this example will define each x,y set as RA+/-edge and Dec+/-edge.

In [15]:

sqp = SodaQuery.from_resource(dl_result,
                              dl_result.get_adhocservice_by_id("cutout-sync-exposure"),
                              session=get_pyvo_auth())

ra_edge = 0.02 * u.deg
de_edge = 0.005 * u.deg

sqp.polygon = (cutout_ra - ra_edge,
               cutout_dec - de_edge,
               cutout_ra - ra_edge,
               cutout_dec + de_edge,
               cutout_ra + ra_edge,
               cutout_dec + de_edge,
               cutout_ra + ra_edge,
               cutout_dec - de_edge)

cutout_bytes = sqp.execute_stream().read()
sqp.raise_if_error()
polygon = read_archive(io.BytesIO(cutout_bytes))

In [16]:

display = afwDisplay.Display()
display.scale('asinh', 'zscale')
display.mtv(polygon.image)
plt.show()

Figure 2: A rectangular polygon cutout from a deep_coadd.

3.4. Correcting for cos(d)¶

There is an important difference to note between the circle and polygon shape definitions. The angular distance on the sky that defines the circular cutout size already accounts for the difference in angular distance in the RA direction is smaller by a factor of cos(declination), where declination is in units radians. The difference increases with higher declination. However, the polygon definition does not automatically account for this cosine factor. Thus, circle and polygon cutout definitions using the same cutout edge length will not match size in the RA direction. The 2 cells below demonstrate how to make this correction to the polygon cutout definition to create symmetric cutouts with polygon. Here, reset the edge sizes to be the same as Radius from the circle definition above.

First, generate a polygon cutout without factoring in cos(dec).

In [17]:

sq2 = SodaQuery.from_resource(dl_result,
                              dl_result.get_adhocservice_by_id("cutout-sync-exposure"),
                              session=get_pyvo_auth())

sq2.polygon = (cutout_ra - Radius,
               cutout_dec - Radius,
               cutout_ra - Radius,
               cutout_dec + Radius,
               cutout_ra + Radius,
               cutout_dec + Radius)

cutout_bytes = sq2.execute_stream().read()
sq2.raise_if_error()
polygon2 = read_archive(io.BytesIO(cutout_bytes))

Second, generate a polygon cutout and include the factor of cos(dec), to match the area generated by circle.

In [18]:

spherePoint = geom.SpherePoint(target_ra*geom.degrees, target_dec*geom.degrees)
a = Angle(spherePoint.getDec().asDegrees(), u.deg)
cosd = np.cos(a.radian)

sq3 = SodaQuery.from_resource(dl_result,
                              dl_result.get_adhocservice_by_id("cutout-sync-exposure"),
                              session=get_pyvo_auth())

sq3.polygon = (cutout_ra - Radius/cosd,
               cutout_dec - Radius,
               cutout_ra - Radius/cosd,
               cutout_dec + Radius,
               cutout_ra + Radius/cosd,
               cutout_dec + Radius,
               cutout_ra + Radius/cosd,
               cutout_dec - Radius)

cutout_bytes = sq3.execute_stream().read()
sq3.raise_if_error()
polygon3 = read_archive(io.BytesIO(cutout_bytes))

Plot the three cutouts as a comparison below. Setting width_ratios makes sure the y-axes (declination direction) span the same extent in the figure to emphasize that without the cos(dec) factor, the R.A. direction (x-axis) is truncated using polygon relative to circle.

In [19]:

fig, ax = plt.subplots(1, 3, width_ratios=[0.35, 0.3, 0.35], figsize=(10, 14))

plt.sca(ax[0])
display1 = afwDisplay.Display(frame=fig, useWcsAxes=False)
display1.scale('linear', 'zscale')
display1.mtv(cutout.image, title='Cutout defined with circle')

plt.sca(ax[1])
display2 = afwDisplay.Display(frame=fig, useWcsAxes=False)
display2.scale('linear', 'zscale')
display2.mtv(polygon2.image, title='Cutout defined with polygon')

plt.sca(ax[2])
display3 = afwDisplay.Display(frame=fig, useWcsAxes=False)
display3.scale('linear', 'zscale')
display3.mtv(polygon3.image, title='Polygon including cos(dec)')
plt.tight_layout()
plt.show()

Figure 3: A comparison of a cutout generated with the circle function (same as Figure 1; left panel) with a cutout defined using the polygon functionality defined using the same edge size (middle panel). The right panel uses polygon but accounts for the $cos({\rm dec})$ term, to replicate the same cutout that was made using circle.

For citizen science applications, the zooniverse package should be used instead of this procedure.

4. Masked image cutouts¶

As mentioned earlier, there are three types of cutout services. Demonstrated above in Section 3 is the cutout-sync-exposure service that returns an image as a CellCoadd storage class with all LSST image extensions. The cutout-sync-maskedimage service returns an image as a MaskedImage storage class that contains the LSST image, variance, and bitmask. But, unlike CellCoadd, the MaskedImage is missing the full range of metadata that may be needed to perform science on the images, including the PSF and aperture correction map, and information about the background that was subtracted.

Reuse the data link dl_result that was defined in Section 3 to define a new cutout, this time using cutout-sync-maskedimage service. Call this soda query sqM to differentiate from the earlier ones. Define the cutout as a square subtended by a circle of radius 0.01 degrees, as done in Section 3.

In [20]:

sqM = SodaQuery.from_resource(dl_result,
                              dl_result.get_adhocservice_by_id("cutout-sync-maskedimage"),
                              session=get_pyvo_auth())

sqM.circle = (cutout_ra, cutout_dec, Radius)

cutout_bytes = sqM.execute_stream().read()
sqM.raise_if_error()

In [21]:

masked_im = read_archive(io.BytesIO(cutout_bytes))

Extract the attributes associated to the masked image and the cutout into sets, then use the .difference function to see what is present in the CellCoadd that is not in the MaskedImage.

In [22]:

masked_im_attributes = set(dir(masked_im))
cell_coadd_attributes = set(dir(cutout))
print('Present in CellCoadd but not MaskedImage:')
cell_coadd_attributes.difference(masked_im_attributes)

Present in CellCoadd but not MaskedImage:

Out[22]:

{'_aperture_corrections',
 '_backgrounds',
 '_band',
 '_mask_fractions',
 '_noise_realizations',
 '_patch',
 '_provenance',
 '_psf',
 'aperture_corrections',
 'apply_background',
 'backgrounds',
 'band',
 'bounds',
 'from_legacy_cell_coadd',
 'grid',
 'mask_fractions',
 'noise_realizations',
 'patch',
 'provenance',
 'psf',
 'skymap',
 'to_legacy_cell_coadd',
 'tract'}

masked_im is an LSST MaskedImage object. In addition to the mask, it also contains the image and variance extensions, which can then be accessed using the .image or .variance attributes.

In [23]:

display = afwDisplay.Display()
display.scale('asinh', 'zscale')
display.image(masked_im)
plt.show()

Figure 4: The cutout image retrieved as a MaskedImage object, displayed in pixel coordinates using LSST pipeline tools in grayscale with a scale bar at right. The mask is plotted on top of the image pixels in colors.

5. Exercise for the learner¶

Reproduce the cutout below, whose center is (ra, dec) = 59.1, -48.8 with 0.06 degrees on a side.

In [ ]: