104.5. Image subsets with the Butler

104.5. Image subsets with the Butler¶

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-21
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: How to make image subsets with the Butler.

LSST data products: deep_coadd

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

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 make image subsets (a smaller piece of a parent LSST image) using the Butler. This method for generating small image subsets is the right tool when small images are desired from a single visit image or deep coadd patch, or, when the analysis workflow is already using the Butler for image analysis. Instead, the Rubin cutout service is a better tool when wanting to generate image cutouts at a range of coordinates or observations times. The Rubin cutout service is demonstrated in notebooks in the DP1 or DP2 103 notebook series.

Butler-related documentation:

• pipelines middleware Frequently Asked Questions
• Butler python module documentation
• Butler query expressions and operators

This tutorial demonstrates how to generate image subsets using the Butler using a variety of methods.

Related tutorials: The earlier 100-level Butler tutorials in this series show how to retrieve data with the Butler and how to explore and discover the dataset types and their properties. The 103 series demonstrate the Rubin cutout service to produce image cutouts.

1.1. Import packages¶

Import the butler module from the lsst.daf package, and the display module from the lsst.afw package (for image display). Import the lsst.sphgeom packages to enable spatial and subset functionality. The lsst.images package imports handling for LSST image types for releases DP2 and later.

In [1]:

import matplotlib.pyplot as plt
import math
from lsst.daf.butler import Butler
import lsst.sphgeom as sphgeom
import lsst.afw.display as afw_display
from lsst.images import Box, Interval
from astropy.coordinates import SkyCoord
from astropy.visualization import ZScaleInterval

1.2. Define parameters¶

Create an instance of the Butler with the repository and collection for DP2.

In [2]:

butler = Butler('dp2', collections="dp2")

Set afwDisplay to use matplotlib.

In [3]:

afw_display.setDefaultBackend('matplotlib')

Coordinates: Use coordinates RA, Dec = $53.076, -28.110$ deg, which are near the center of the Extended Chandra Deep Field South (ECDFS).

Region: Define a circle with a radius of 0.1 deg, centered on the coordinates.

In [4]:

ra = 53.076
dec = -28.110
band = 'r'

radius = 0.1
region = sphgeom.Region.from_ivoa_pos(f"CIRCLE {ra} {dec} {radius}")

2. The Bounding Box¶

A bounding box is useful for defining boundaries of images (either in absolute pixel coordinates of the parent tract, the local pixel coordinates of an individual deep_coadd image, or to make a subset from an image). This section will retrieve an image from the Butler and demonstrate how to define a bounding box (bbox) for a subset that is centered on pixel coordinate (y, x) and is 1000 pixels wide and 1000 pixels high.

Note that the bounding box is indexed as (y, x) with the y direction first.

2.1. Identify deep_coadds¶

First, query the Butler for r-band images that overlap the sky region defined above.

In [5]:

dataset_refs = butler.query_datasets("deep_coadd",
                                     where="band.name='r' AND\
                                     patch.region OVERLAPS POINT(ra, dec)",
                                     bind={"ra": ra, "dec": dec},
                                     with_dimension_records=True,
                                     order_by=["patch.tract"])

Get the dataId for the first returned dataset reference.

In [6]:

coadd = butler.get(dataset_refs[0])
dataId = dataset_refs[0].dataId
print(dataId)

{band: 'r', skymap: 'lsst_cells_v2', tract: 5063, patch: 15}

2.2. Access image boundaries¶

To retrieve the bounding box that defines the full deep_coadd, use the bbox method. By default, bbox returns pixel coordinates in the absolute coordinate system of the entire image tract. The first and last pixels in the x and y direction of the image can be accessed using the .start and .stop methods.

In [7]:

bbox = coadd.bbox
print('bbox defined using .bbox of existing image', bbox)

print(f"X-pixel range: [{bbox.x.start}, {bbox.x.stop}]")
print(f"Y-pixel range: [{bbox.y.start}, {bbox.y.stop}]")

bbox defined using .bbox of existing image [y=2850:6150, x=14850:18150]
X-pixel range: [14850, 18150]
Y-pixel range: [2850, 6150]

Alternatively one can also use the .min and .max which are inclusive bounds. The .start and .stop shown above are half inclusive (start and min identify the same pixel but stop is one less than max).

In [8]:

print(f"X-pixel range: [{bbox.x.min}, {bbox.x.max}]")
print(f"Y-pixel range: [{bbox.y.min}, {bbox.y.max}]")

X-pixel range: [14850, 18149]
Y-pixel range: [2850, 6149]

2.3. Define subset boundaries¶

Below, construct a new bounding box using the Box function to define a subset region of the full deep_coadd.

The Box function is part of the lsst.images package, and can also take explicit x and y keyword arguments, so that there is no confusion about the axes. In this case, an interval of x and y coordinates must be constructed using the Interval function in order for bbox to be in the correct format.

In [9]:

x_interval = Interval(bbox.x.start+100, bbox.x.stop-100)
y_interval = Interval(bbox.y.start+100, bbox.y.stop-100)

print("correct pixel range format: ", x_interval)

bbox = Box(x=x_interval, y=y_interval)
print(bbox)

correct pixel range format:  14950:18050
[y=2950:6050, x=14950:18050]

Equivalently, this can be achieved using the .padded method that shrinks the bounds of the existing bounding box by 100 pixels, but achieves the same format without needing the Interval or Box functions.

In [10]:

smaller_box = coadd.bbox.padded(-100)
x_interval = smaller_box.x
y_interval = smaller_box.y
print(smaller_box)
print(x_interval, y_interval)

[y=2950:6050, x=14950:18050]
14950:18050 2950:6050

2.4 Local vs absolute coordinates¶

By default, all pixel coordinates are in the absolute coordinate system (i.e. that of the parent tract to which the deep_coadd belongs).

If one wants to specify using the local array pixel coordinates (rather than absolute coordinates) then this can be done using the local method and pass the range of x and y pixels in local coordinates directly. Note however that the bounding box itself will still be in the absolute (parent patch/tract) coordinate system, even if defined initially in local coordinates (which is demonstrated below using the start and stop methods). Find documentation on the CellCoadd.local property at images.lsst.io.

In [11]:

coadd_local = coadd.local[10:1000,10:1000]
print(coadd_local.bbox)
print(f"X-pixel range: [{coadd_local.bbox.x.start}, {coadd_local.bbox.x.stop}]")
print(f"Y-pixel range: [{coadd_local.bbox.y.start}, {coadd_local.bbox.y.stop}]")

[y=2860:3850, x=14860:15850]
X-pixel range: [14860, 15850]
Y-pixel range: [2860, 3850]

2.5. Sky coordinates¶

The bounding box can also be defined starting with sky coordinates (converted to pixels) in the following way. Retrieving the world coordinate system (WCS) of the deep_coadd with the sky_projection method enables converting an ra, dec to pixel coordinates of the desired image.

In [12]:

point = SkyCoord(ra, dec, unit='deg')
xy = coadd.sky_projection.sky_to_pixel(point)
print(point, xy)

<SkyCoord (ICRS): (ra, dec) in deg
    (53.076, -28.11)> XY(x=15182.507513158846, y=4390.517978014746)

Then, define a spatial extent for the image subset (e.g. 60 arcseconds) and convert those to pixel coordinates to define the bounding box. Box.factory is a special property that can indexed to create a new box, is more concise, and harder to get wrong than passing x and y arguments to Box() as was done in Section 2.3.

In [13]:

pixel_scale_arcsec = 0.2
subset_side_arcsec = 60.0
subset_side_pixels = int(subset_side_arcsec / pixel_scale_arcsec)
half_edge = subset_side_pixels // 2
print(half_edge)
print(xy.y)

bbox_sky = Box.factory[round(xy.y):round(xy.y)+1, round(xy.x):round(xy.x)+1].padded(half_edge)
print(bbox_sky)

150
4390.517978014746
[y=4241:4542, x=15033:15334]

Alternatively, if the pixel range is already known it can be entered into Box.factory directly as below.

In [14]:

bbox = Box.factory[4241:4542, 15033:15334]

3. Generate image subsets¶

The last section explored how to access and construct bounding boxes. This section demonstrates how to use the bounding boxes to generate the image subsets.

3.1 Pass bounding box to the butler¶

Below, pass the bounding box to the butler in the parameters dictionary directly. This saves memory and time by only loading the requested pixels, rather than retrieving the full image first and generating the subset after loading it into memory.

Below, pass the bounding box to the Butler with the parameters keyword to request a subset of the deep_coadd image defined using the bbox. The output is an image with dimensions that match the bounding box in CellCoadd format (the new equivalent of ExposureF in the lsst.images image format).

In [15]:

subset = butler.get('deep_coadd',
                    dataId=dataId,
                    parameters={'bbox': bbox}
                   )

print(subset)

CellCoadd([y=4241:4542, x=15033:15334], tract=5063)

Display the image subset returned by the butler.

In [16]:

display = afw_display.Display(frame=1)
display.scale('asinh', 'zscale')
display.image(subset.image.to_legacy())
plt.show()

Figure 1: A subset of a sky region in the ECDFS returned after passing a bounding box to the Butler

3.2. Subsets from a full exposure¶

It is also possible to make subsets after retrieving the full deep_coadd from the butler. This may be preferable if many subsets are needed from the same patch and tract image.

First, retrieve the full deep_coadd of interest.

In [17]:

full_exposure = butler.get('deep_coadd', dataId=dataId)

Extract the region of interest using direct bounding box indexing on the full deep_coadd exposure. A simple example would look like:

subset_exposure = full_exposure[bbox_sky]

Note that subset is not creating a "deep copy" (equivalent of deep_copy from python's built-in copy module) of the image by default. The subset points to the same block of memory as the parent full_exposure. No new pixel data is duplicated into memory. While this is efficient in that it is quick and uses less memory, because the subset and the parent share the same memory, modifying one modifies the other. Any change to the pixel values in the subset will apply also to the full_exposure stored in memory. Any operations that are performed on the pixel array generated by the subset method, subset_subset, will also happen to full_exposure. Note that this means that even if the full image is deleted to save memory, the whole original memory block will continue to remain allocated until the image subset is also deleted.

Alternatively, use the .copy() method to make a copy of the subset from the full coadd that does not overwrite the memory location of the original:

subset_copy = full_exposure[bbox_sky].copy()

Using the copy method is instead is like deep_copy in that it allocates new memory and physically duplicates the subset into a new image object. This allows the user to apply mathematical operations on the subset without altering the the full_exposure image also stored in memory. However, it takes longer to execute and uses more memory (though for small or few subsets this is negligible).

The use case for copy: When the user already has the full image from the butler in memory, and wants to make one or more subsets from it. Generate such a subset to independent memory below, and leave full_exposure stored in memory in its original form.

In [18]:

subset_exposure = full_exposure[bbox_sky].copy()

Plot the image extension of the subset.

In [19]:

display = afw_display.Display(frame=1)
display.scale('asinh', 'zscale')
display.image(subset_exposure.image.to_legacy())
plt.show()

Figure 2: A subset of the ECDFS but generated using the bounding box indexing method

3.3. Convenience methods¶

Some properties of the subset can be accessed with methods such as .shape, and image extensions using .image or .variance

Values of interest to defining the subset, including including the dimensions of the bounding box can be accessed through the shape method.

In [20]:

print(f"Successfully created subset of shape: {subset.image.array.shape}")

Successfully created subset of shape: (301, 301)

In [21]:

print("The pixel dimensions of the coadd are ", bbox.shape.x, bbox.shape.y)

The pixel dimensions of the coadd are  301 301

It can be useful to first check whether an image subset defined by a bounding box is inside the full deep_coadd image before attempting to generate the subset, to avoid errors. Do this using the contains method.

In [22]:

print(full_exposure.bbox.contains(bbox_sky))

True

The intersection method generates a new bounding box for an overlapping region between two bounding boxes. In the case below, bbox_sky is fully contained within the bounding box of full_exposure and so bbox_sky is returned.

In [23]:

print(bbox_sky)
print(full_exposure.bbox)
bbox_overlap = full_exposure.bbox.intersection(bbox_sky)
print(bbox_overlap)

[y=4241:4542, x=15033:15334]
[y=2850:6150, x=14850:18150]
[y=4241:4542, x=15033:15334]

4. Retrieve multiple subsets from the butler¶

This section demonstrates an efficient way to retrieve multiple subsets from the butler. Instead of first retrieving the deep_coadd image into memory (which is many gigabytes is size), an alternative is to use the butler.getDeferred() method. This identifies the images on the remote server and returns a lightweight "pointer" or "handle" (specifically, a DeferredDatasetHandle) instead of the actual data. getDeferred is fast and uses nearly no memory, while queuing up a bunch of subsets to retrieve at a later time all at once without transferring the deep_coadd(s) themselves. The call handle.get(parameters=params) moves the subsets to the local disk. Passing the bbox parameter enables the butler to transfer only the pixels inside the bounding box.

First, define a list of targets in a dictionary containing bright (V band magnitude < 22) variability selected galaxies hosting active galactic nuclei (AGN) in the ECDFS field from Boutsia et al. 2009. Use the sky_projection.sky_to_pixel methods to convert the ra and dec to pixel coordinates.

In [24]:

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

First set up some useful parameters to generate nice looking and organized plots.

In [25]:

zscale = ZScaleInterval()
subset_size = 20
ncols = 4
nrows = math.ceil(len(targets) / ncols)

The cell below loops through the targets, first identifying the patch and tract containing the galaxy (since the galaxies are not necessarily all in the same patch and tract) and then constructing a 20x20 pixel bounding box (subset_size) centered at the target. As the loop iterates to the next galaxy, the previous subset object is overwritten so that the memory footprint never exceeds the size of a single subset.

In [26]:

fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(16, 4 * nrows))
axes = axes.flatten()

idx = 0
for target in targets:

    ra = target['ra']
    dec = target['dec']

    dataset_refs = butler.query_datasets("deep_coadd",
                                         where="band.name='r' AND\
                                         patch.region OVERLAPS POINT(ra, dec)",
                                         bind={"ra": ra, "dec": dec},
                                         with_dimension_records=True,
                                         order_by=["patch.tract"])

    dataId = dataset_refs[0].dataId
    point = SkyCoord(ra, dec, unit='deg')

    wcs = butler.get('deep_coadd.sky_projection', dataId=dataId)
    xy = wcs.sky_to_pixel(point)

    selection_bbox = Box.factory[round(xy.y):round(xy.y)+1, round(xy.x):round(xy.x)+1].padded(subset_size)
    params = {'bbox': selection_bbox}
    subset = butler.get(dataset_refs[0], parameters=params)

    print(f"Processed ID {target['id']}")

    img_array = subset.image.array
    vmin, vmax = zscale.get_limits(img_array)
    ax = axes[idx]

    ax.imshow(img_array, origin='lower')
    ax.set_title(f"{target['id']}")
    idx = idx+1

plt.tight_layout()
plt.show()

Processed ID 3

Processed ID 8

Processed ID 9

Processed ID 11

Processed ID 12

Processed ID 14

Processed ID 15

Processed ID 18

Processed ID 20

Processed ID 22

Processed ID 24

Processed ID 25

Processed ID 26

Processed ID 28

Processed ID 29

Processed ID 30

Processed ID 31

Processed ID 32

Processed ID 33

Processed ID 34

Processed ID 35

Processed ID 36

Processed ID 37

Processed ID 39

Processed ID 41

Processed ID 43

Processed ID 44

Figure 5: A grid of subsets for 27 AGNs in the ECDFS field generated by the Butler using deferred handles.