103.4. Cutout exposures

103.4. Cutout Exposures¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: DP1
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: How to use the Rubin image cutout service to make cutout exposures with DP1.

LSST data products: visit_image

Packages: lsst.rsp.utils, pyvo, lsst.rsp.service, 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 DP1 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 ExposureF data type (Section 3), and cutout-sync-maskedimage which returns the mask extension with the science and variance images (MaskedImageF data type; 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 the image extension and header of minimal metadata, and if not needing 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, and utilities for remote data access lsst.rsp.

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

In [1]:

import numpy as np
import matplotlib.pyplot as plt

import lsst.afw.display as afwDisplay
from lsst.afw.image import ExposureF, MaskedImageF
from lsst.afw.fits import MemFileManager

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

import lsst.geom as geom

from pyvo.dal.adhoc import DatalinkResults, SodaQuery

from astropy import units as u
from astropy.coordinates import Angle
from astropy.time import Time
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 DP1 images (in the format of a web URL identifying where the data is hosted).

Instantiate the SIA service.

In [3]:

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

2. Find the visit image¶

The cutout services 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, for a visit image that was obtained between MJD 60623.256 and 60623.259.

Define the coordinates right ascension (target_ra) and declination (target_dec) in degrees. Define the band and the start and end times.

In [4]:

target_ra = 53.1246023
target_dec = -27.7404715
eff_wl = 622.1e-09
time1 = Time(60623.256, format="mjd", scale="tai")
time2 = Time(60623.259, format="mjd", scale="tai")

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.05 degree circle, centered on the target.

In [5]:

circle = (target_ra, target_dec, 0.05)

This query will return 1 visit_image (by design).

In [6]:

results = service.search(pos=circle, calib_level=2,
                         dpsubtype='lsst.visit_image',
                         band=eff_wl, time=(time1, time2))
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.visit_image	Rubin:Simonyi	2	slew_icrs	CC_O_20241108_000246	LSST.DP1	ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=019b1962-3d55-76a7-b646-7714eb4c3c79	https://data-int.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1962-3d55-76a7-b646-7714eb4c3c79	application/x-votable+xml;content=datalink	--	4072	4000	--	60623.25872464106	60623.25907695602	30.0	--	--	5.51e-07	6.891e-07	--	r	phot.flux.density	--	LSSTComCam	2024110800246	2	--	--	r	r_03	visit_image - r - CC_O_20241108_000246-R22_S02 2024-11-09T06:12:33.808988Z	53.124250253081	-27.73236806603477	0.3179328995436719	POLYGON ICRS 53.030319 -27.596909 52.972978 -27.817875 53.218386 -27.867735 53.275334 -27.646761

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.DP1' AND calib_level = 2
      AND dataproduct_type = 'image' AND instrument_name = 'LSSTComCam'
      AND dataproduct_subtype = 'lsst.visit_image'
      AND ( t_min <= 60623.259 AND 60623.256 <= t_max )
      AND ( 622e-9 BETWEEN em_min AND em_max )

3. Generating an image cutout¶

First, extract the datalink (access URL) from the first row of the results as datalink_url, using the from_result_url method from the pyvo package.

Next, provide authorization for the current RSP session, since the Rubin DP1 imaging is proprietary. Do this using the get_pyvo_auth function.

Finally, create a DatalinkResults object to be able to access this 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 VOTable document, stored as dl_result.

In [8]:

datalink_url = results[0].access_url
dl_result = DatalinkResults.from_result_url(datalink_url,
                                            session=get_pyvo_auth())

f"Datalink status: {dl_result.status}. Datalink service url: {datalink_url}"

Out[8]:

"Datalink status: ('OK', 'QUERY_STATUS not specified'). Datalink service url: https://data-int.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1962-3d55-76a7-b646-7714eb4c3c79"

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 contained in the ExposureF data type.

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 DP1 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 the full LSST visit_image, 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¶

The cutout should be read into memory in the ExposureF format using LSST pipeline tools. To enable image display, first store the bytes from mem as an LSST ExposureF object. The image extension can then be accessed using the .image method. REPEATED BELOW FIX

In [11]:

mem = MemFileManager(len(cutout_bytes))
mem.setData(cutout_bytes, len(cutout_bytes))
exposure = ExposureF(mem)

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

In [12]:

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

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

The WCS information can be obtained using the getWcs method, and converted to an astropy WCS using the WCS function. Optionally display the WCS by uncommenting the last line. Values of specific wcs keywords can be retrieved as, e.g. exposureWCS['CRPIX2']

In [13]:

exposureWCS = exposure.getWcs().getFitsMetadata()
astropyWCS = WCS(exposureWCS)
# print(exposureWCS)

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

# 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: ExposureF contains multiple Header-Data Units (HDUs) with extensive exposure metadata. To reduce file size, one option is to save only the image, mask, and variance planes, along with the WCS recorded in the Primary HDU. This can be done by writing exposure.maskedImage-a three-plane object containing the image, mask, and variance-into a FITS file using writeFitsImage.

In [15]:

# fits.info(sodaCutout)

# sodaCutout_small = os.path.join(tempdir, 'cutout-circle_small.fits')
# afwDisplay.writeFitsImage(sodaCutout_small,
#                           exposure.maskedImage,
#                           wcs=exposure.wcs)

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.

Warning: Visit images are rotated, and although it is the "Declination edge" that is defined to be smaller, that corresponds to the x-axis when the cutout is displayed below, due to image rotation.

In [16]:

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()
mem = MemFileManager(len(cutout_bytes))
mem.setData(cutout_bytes, len(cutout_bytes))
polygon = ExposureF(mem)

In [17]:

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

Figure 2: A rectangular polygon cutout from a rotated visit_image.

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

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()
mem = MemFileManager(len(cutout_bytes))
mem.setData(cutout_bytes, len(cutout_bytes))
polygon2 = ExposureF(mem)

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

In [19]:

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()
mem = MemFileManager(len(cutout_bytes))
mem.setData(cutout_bytes, len(cutout_bytes))
polygon3 = ExposureF(mem)

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

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)
display1.scale('linear', 'zscale')
display1.mtv(exposure.image, title='Cutout defined with circle')

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

plt.sca(ax[2])
display3 = afwDisplay.Display(frame=fig)
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.

The zooniverse package should be used instead of this procedure for citizen science applications.

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 ExposureF type image with all LSST image extensions. The cutout-sync-maskedimage service that returns an MaskedImageF type containing the LSST image, variance, and bitmask. But, unlike ExposureF, the MaskedImageF is missing the full range of metadata that may be needed to perform science on the images, including the WCS model and bounding box, the photometric calibration and physical filter used, optical models including the PSF and aperture correction map, the visit info and detector used.

Reuse the data link dl_result that was defined in Section 3 to define a new cutout, this time using cutout-sync-maskedimage service. The returned 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 [21]:

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

mem = MemFileManager(len(cutout_bytes))
mem.setData(cutout_bytes, len(cutout_bytes))

To enable image display, first store the bytes from mem as an LSST MaskedImageF object. In addition to the mask, it also contains the image and variance extensions, which can then be accessed using the .image or .variance attribute.

In [23]:

mask = MaskedImageF(mem)
display = afwDisplay.Display()
display.scale('asinh', 'zscale')
display.image(mask)
plt.show()

Figure 4: The cutout image retrieved in MaskedImageF format, displayed in pixel coordinates using LSST pipeline tools in grayscale with a scale bar at right. 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.