202.1. Deep coadds

202.1. Deep coadds¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 1
Container Size: large
LSST Science Pipelines version: r29.1.1
Last verified to run: 2025-06-21
Repository: github.com/lsst/tutorial-notebooks

Learning objective: To understand the format and metadata of a deep_coadd.

LSST data products: deep_coadd

Packages: lsst.daf.butler, lsst.rsp, 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¶

A deep_coadd is a combination of multiple processed, calibrated, and background-subtracted images, for a patch of sky, for one of the six LSST filters.

• butler dataset type: deep_coadd
• number of images: 2644

Related tutorials: The 100-level tutorials on how to use the Butler, TAP, SIA, Firefly, and the cutout service demonstrate how to query for and retrieve images, how to create small image cutouts ("stamps"), and how to manipulate the Firefly display interface.

1.1. Import packages¶

From the LSST Science Pipelines import the packages for the Butler, 2-dimensional sky geometry, and for image display. Also import the numpy package.

In [1]:

from lsst.daf.butler import Butler
import lsst.afw.display as afwDisplay
import lsst.geom
import numpy as np

1.2. Define parameters and functions¶

Instantiate the Butler.

In [2]:

butler = Butler("dp1", collections="LSSTComCam/DP1")
assert butler is not None

Set afw_display to use Firefly and open the Firefly display tab.

In [3]:

afwDisplay.setDefaultBackend('firefly')
afw_display = afwDisplay.Display(frame=1)

2. Data access¶

It is recommended to use the Rubin data Butler to retrieve image data within Jupyter Notebooks, but they are also available via the TAP and SIA services (Table Access Protocol and Simple Image Access).

Recommended access method: Butler

2.1. Butler¶

Show that the dimensions for the deep_coadd dataset type are the band, and the skymap's tract and patch, and that the returned image will be of type ExposureF. The LSST Science Piplines tools for image manipulation and visualization work best with the ExposureF format. This format includes pixel data and metadata.

In [4]:

butler.get_dataset_type('deep_coadd')

Out[4]:

DatasetType('deep_coadd', {band, skymap, tract, patch}, ExposureF)

In [5]:

butler.get_dataset_type('deep_coadd').dimensions.required

Out[5]:

{band, skymap, tract, patch}

2.1.1. Demo query¶

Query for deep_coadd imges that overlap coordinates RA, Dec near the center of the ECDFS field and were obtained with the r-band filter.

In [6]:

ra = 53.076
dec = -28.110
band = 'r'
query = f"band.name = '{band}' AND patch.region OVERLAPS POINT({ra}, {dec})"
dataset_refs = butler.query_datasets("deep_coadd", where=query)
assert len(dataset_refs) == 2
print(len(dataset_refs))
del ra, dec, band, query

2

The deep_coadd images do overlap at their edges, and for the above query the fact that 2 deep_coadd images match the query constraints mean that the RA, Dec coordinates are in the overlap region.

2.1.2. Get an image¶

The deep_coadd can be retrieved from the Butler using the (a) the dataset_refs; (b) the dataId; or (c) the tract and patch numbers.

(a) Use the dataset reference¶

For the first dataset reference returned by the query, get the corresponding deep_coadd.

In [7]:

ref = dataset_refs[0]
deep_coadd = butler.get(ref)

See that it is of type ExposureF.

In [8]:

deep_coadd

Out[8]:

<lsst.afw.image._exposure.ExposureF at 0x7a8b0517b9f0>

Display the image in the Firefly window. Set the mask transparency to 100 (fully transparent, i.e., do not show the mask).

In [9]:

afw_display.mtv(deep_coadd)
afw_display.setMaskTransparency(100)

Clean up.

In [10]:

del ref

(b) Use the dataId¶

For the first dataset reference returned by the query, get the dataId.

In [11]:

ref = dataset_refs[0]
use_dataId = ref.dataId
use_dataId

Out[11]:

{band: 'r', skymap: 'lsst_cells_v1', tract: 5063, patch: 14}

Option to use the defined dataId to retrieve the corresponding deep_coadd.

In [12]:

# deep_coadd = butler.get('deep_coadd', dataId = use_dataId)

Option to display the deep_coadd.

In [13]:

# afw_display.mtv(deep_coadd)

Clean up.

In [14]:

del ref, use_dataId

(c) Use the tract, patch, and band¶

For the first dataset reference returned by the query, get the dataId and extract the band and the tract and patch numbers.

In [15]:

ref = dataset_refs[0]
dataId = ref.dataId
use_band = dataId.get('band')
use_tract = dataId.get('tract')
use_patch = dataId.get('patch')
print(use_band, use_tract, use_patch)

r 5063 14

Even if the Butler is not used to query for the images (e.g., if the SIA or TAP service is used to find images matching a set of constraints), images can be retrieved from the Butler by passing the band, tract, and patch.

Option to use the band, tract, and patch to get the corresponding deep_coadd.

In [16]:

# deep_coadd = butler.get('deep_coadd', band=use_band,
#                         tract=use_tract, patch=use_patch)

Option to display the deep_coadd.

In [17]:

# afw_display.mtv(deep_coadd)

Clean up (but keep the deep_coadd and the dataset_refs to use below).

In [18]:

del ref, dataId, use_band, use_tract, use_patch

2.2. SIA (Simple Image Access)¶

The Simple Image Access (SIA) service is provides a standardized model for image metadata, and the capability to query and retrieve image datasets. Use of SIA, and data retrieval via the access_url in the SIA results table, is demonstrated in the 100-level tutorials.

To return deep_coadd images from a SIA query, set:

• calib_level = 3
• dpsubtype = 'lsst.deep_coadd'

See Section 2.1.2. for how to retrieve images via the Butler using the lsst_band, lsst_tract, and lsst_patch columns in the SIA results table.

2.3. TAP (Table Access Protocol)¶

The Table Access Protocol (TAP) service provides a standardized model for accessing image metadata and catalog data with the Astronomical Data Query Language (ADQL). Use of TAP, and data retrieval via the access_url in the TAP results table, is demonstrated in the 100-level tutorials.

Image metadata is stored in the ObsCore table.

To return deep_coadd images from a TAP query, set:

• Calibration level 3.
• Dataproduct subtype: lsst.deep_coadd

See Section 2.1.2. for how to retrieve images via the Butler using the lsst_band, lsst_tract, and lsst_patch columns in the SIA results table.

3. Pixel data¶

The deep_coadd images have three planes: image, variance, and mask.

Use the deep_coadd retrieved in Section 2.1.2.

Display the deep_coadd in frame 1, and set the mask transparency to 0 (not transparent; to show the mask).

In [19]:

afw_display.mtv(deep_coadd)
afw_display.setMaskTransparency(0)

With the mask plane fully opaque, it should look like this:

Figure 1: The mask plane of the r-band deep_coadd image for tract 5063, patch 14. The colors correspond to an informative mask flag value (Section 3.3). Almost every pixel has a color (has a nonzero mask value).

3.1. Image plane¶

Sky pixel data in units of nJy (nanoJanskys).

The background has been subtracted so values can be negative.

Extract the image data from the visit_image, and convert it to an array.

In [20]:

image = deep_coadd.getImage()
image_array = image.array

Option to show the image and the image_array.

In [21]:

# image

In [22]:

# image_array

Print statistics on the image_array.

In [23]:

print(image_array.min())
print(image_array.mean())
print(image_array.max())

-6.4893517
29.395851
1527183.1

Display the image in Firefly frame 2.

In [24]:

afw_display = afwDisplay.Display(frame=2)
afw_display.mtv(image)

Mouse-over display frame 2. As the image contains only pixel data, the image is displayed without a WCS: the coordinates are in pixels, not RA, Dec (as in frame 1).

3.2. Variance plane¶

Uncertainty (noise) in the sky pixel data in units of nJy^2 (nanoJanskys squared).

Values are always positive.

Extract the variance and variance_array.

In [25]:

variance = deep_coadd.getVariance()
variance_array = variance.array

Print statistics on the variance_array.

In [26]:

print(variance_array.min())
print(variance_array.mean())
print(variance_array.max())

1.3933268
171.81703
20352022.0

Display the variance.

In [27]:

afw_display = afwDisplay.Display(frame=3)
afw_display.mtv(variance)

3.3. Mask plane¶

An integer bitmask of representative flag values that indicate processing status or issues, similar to the SDSS bitmasks.

Extract the mask.

In [28]:

mask = deep_coadd.getMask()
mask_array = mask.array

Mask keys (names) are defined by the mask and interpreted as colors by afw_display.

Print the list of mask keys and color.

In [29]:

for mask_key, bit in mask.getMaskPlaneDict().items():
    print('{} ({}): {}'.format(mask_key, bit,
                               afw_display.getMaskPlaneColor(mask_key)))

BAD (0): red
CLIPPED (11): blue
CR (3): magenta
CROSSTALK (12): cyan
DETECTED (5): blue
DETECTED_NEGATIVE (6): cyan
EDGE (4): yellow
INEXACT_PSF (13): magenta
INTRP (2): green
ITL_DIP (14): yellow
NOT_DEBLENDED (15): orange
NO_DATA (8): orange
REJECTED (16): red
SAT (1): green
SENSOR_EDGE (17): green
STREAK (10): green
SUSPECT (7): yellow
UNMASKEDNAN (18): blue
VIGNETTED (9): red

Display the mask in Firefly frame 4.

In [30]:

afw_display = afwDisplay.Display(frame=4)
afw_display.mtv(mask)

The colors are different because the Firefly display interprets the pixel values as flux values (mouse-over and see).

Set the mask transparency to opaque to see the same colors as in frame 1.

In [31]:

afw_display.setMaskTransparency(0)

Set all mask key names to transparent except the "DETECTED" mask bit.

In [32]:

afw_display.setMaskTransparency(100)
afw_display.setMaskTransparency(0, 'DETECTED')

Option to print a list of the unique mask array values, the number of pixels with that value, and the overlapping masks for pixels with that value.

In [33]:

# values, counts = np.unique(mask_array, return_counts=True)
# for value, count in zip(values, counts):
#     print(value, '('+str(count)+')', mask.interpret(value))
# del values, counts

Clean up.

In [34]:

del image, image_array, variance, variance_array, mask, mask_array

4. Metadata and attributes¶

4.1. Header¶

The informational equivalent of a FITS image header is the metadata.

Get the deep_coadd metadata.

In [35]:

metadata = deep_coadd.getMetadata()

Option to display the metadata.

In [36]:

# metadata

Convert the metadata to a python dictionary.

In [37]:

md_dict = metadata.toDict()

Show any metadata key that contains the string "UNIT".

In [38]:

for key in md_dict.keys():
    if key.find('UNIT') >= 0:
        print(key)

BUNIT

Print the "BUNIT" from the dictionary (the flux units).

In [39]:

md_dict['BUNIT']

Out[39]:

'nJy'

Clean up.

In [40]:

del metadata, md_dict

4.2. Information¶

Information about the deep_coadd.

Get the info.

In [41]:

info = deep_coadd.getInfo()

Get the list of input images that were combined to generate the deep_coadd.

In [42]:

inputs = info.getCoaddInputs()

Extract the lists of visits and detectors (ccds) as astropy tables, with the option to display the tables.

In [43]:

visits = inputs.visits.asAstropy()
print(len(visits))
# visits

236

In [44]:

ccds = inputs.ccds.asAstropy()
print(len(ccds))
# ccds

726

In [45]:

del info, inputs, visits, ccds

4.3. World Coordinate System (WCS)¶

The reference system that maps pixel coordinate (x, y) to sky coordinate (RA, Dec).

Get the WCS and show it.

In [46]:

wcs = deep_coadd.getWcs()
wcs

Out[46]:

FITS standard SkyWcs:
Sky Origin: (53.0875576037, -27.5206611570)
Pixel Origin: (14999, 14999)
Pixel Scale: 0.2 arcsec/pixel

Convert pixel coordinates to sky coordinates, and vice versa.

In [47]:

coord = lsst.geom.SpherePoint(53.13*lsst.geom.degrees, -28.10*lsst.geom.degrees)
xy = wcs.skyToPixel(coord)

In [48]:

print(xy)

(14325, 4570.4)

In [49]:

xy2 = lsst.geom.Point2D(13000, 4000)
coord2 = wcs.pixelToSky(xy2)

In [50]:

print(coord2)

(53.2134826738, -28.1316360087)

In [51]:

del xy, xy2, coord, coord2

4.4. Bounding box (corners)¶

The bounding box defines the extent (corners) of the image.

Get the bounding box.

In [52]:

bbox = deep_coadd.getBBox()

In [53]:

bbox

Out[53]:

Box2I(corner=Point2I(11800, 2800), dimensions=Extent2I(3400, 3400))

Print the deep_coadd corners in pixels and sky coordinates.

In [54]:

corners = (lsst.geom.Point2D(bbox.beginX, bbox.beginY),
           lsst.geom.Point2D(bbox.beginX, bbox.endY),
           lsst.geom.Point2D(bbox.endX, bbox.endY),
           lsst.geom.Point2D(bbox.endX, bbox.beginY))
for corner in corners:
    print(corner, wcs.pixelToSky(corner))

(11800, 2800) (53.2891979304, -28.1982040136)
(11800, 6200) (53.2888500329, -28.0093360268)
(15200, 6200) (53.0749099188, -28.0094820512)
(15200, 2800) (53.0748880593, -28.1983511912)

Print the area of the bounding box.

In [55]:

print(bbox.area, 'square pixels')
scale = wcs.getPixelScale().asArcminutes()
area = bbox.area * scale**2
print(np.round(area, 2), 'square arcminutes')

11560000 square pixels
128.44 square arcminutes

Check if pixel values are within the bounding box.

In [56]:

assert bbox.contains(100, 100) is False
print(bbox.contains(100, 100))
assert bbox.contains(12000, 4000) is True
print(bbox.contains(12000, 4000))

False
True

Check if sky coordinates are within the bounding box.

In [57]:

coord = lsst.geom.SpherePoint(63.0*lsst.geom.degrees, -38.0*lsst.geom.degrees)
xy = wcs.skyToPixel(coord)
assert bbox.contains(xy.getX(), xy.getY()) is False
print(bbox.contains(xy.getX(), xy.getY()))

False

Clean up.

In [58]:

del wcs, bbox, corners, scale, area, coord, xy

4.5. Point spread function (PSF)¶

The PSF is the empirical estimate of the shape of a point source as a function of position in the image. The PSF encapsulates all of the atmospheric and optical effects that blur point sources into their measured shapes.

Get the PSF.

In [59]:

psf = deep_coadd.getPsf()

At pixel coordinates x, y = 12000, 4000, compute the PSF and display it in frame 5.

In [60]:

xy = lsst.geom.Point2D(12000, 4000)
psf_image = psf.computeImage(xy)
afw_display = afwDisplay.Display(frame=5)
afw_display.mtv(psf_image.convertF())

Compute the shape of the PSF. Print the size in $\sigma$ and in FWHM (full-width at half-max), and the second moments (all units are pixels).

In [61]:

psf_shape = psf.computeShape(xy)
psf_sigma = psf_shape.getDeterminantRadius()
psf_fwhm = psf_sigma * 2.0*np.sqrt(2.0*np.log(2.0))
print(f"PSF sigma: {psf_sigma}, PSF FWHM: {psf_fwhm}")
print("PSF shape (ixx, iyy, ixy): ",psf_shape.getIxx(), psf_shape.getIyy(), psf_shape.getIxy())

PSF sigma: 2.2006036420127026, PSF FWHM: 5.182025567379624
PSF shape (ixx, iyy, ixy):  4.866299254106444 4.820035315344007 0.0664330910196337

Clean up.

In [62]:

del psf, xy, psf_image, psf_shape, psf_sigma, psf_fwhm

5. Detected objects¶

In general the TAP service is the recommended access mechanism for catalog data.

However, in the case where all the objects (or sources, for visit_images) detected in a given image are desired, the Butler offers a quick way to retrieve them.

Get the tract and patch for the deep_coadd.

In [63]:

ref = dataset_refs[0]
dataId = ref.dataId
use_tract = dataId.get('tract')
use_patch = dataId.get('patch')
print(use_tract, use_patch)

5063 14

Show that the only dimension for the object dataset type is tract.

In [64]:

butler.registry.queryDatasetTypes('object')

Out[64]:

[DatasetType('object', {skymap, tract}, ArrowAstropy)]

Since a butler.get call on the object table will return all of the objects in a tract, not the single patch of the deep_coadd, return only a select few of the >1000 columns and include column patch.

In [65]:

use_columns = ['objectId', 'patch',
               'r_centroid_x', 'r_centroid_y']

Retrive the object table as objects.

In [66]:

objects = butler.get('object', tract=use_tract,
                     parameters={'columns': use_columns})

Print the number of sources in the detector for this deep_coadd (for its patch).

In [67]:

tx = np.where(objects['patch'] == use_patch)[0]
print(len(tx))

11123

Option to show the table.

In [68]:

# objects[tx]

Redisplay the deep_coadd in frame 1 with no mask.

In [69]:

afw_display = afwDisplay.Display(frame=1)
afw_display.mtv(deep_coadd)
afw_display.setMaskTransparency(100)

Overplot sources in the deep_coadd with orange circles.

In [70]:

with afw_display.Buffering():
    for i in tx:
        afw_display.dot('o', objects[i]['r_centroid_x'],
                        objects[i]['r_centroid_y'],
                        size=20, ctype='orange')

Notice that no objects near the edges are marked with an orange circle; they are listed as belonging to the adjacent patch.

Clean up.

In [71]:

del ref, dataId, use_tract, use_patch, use_columns, objects, tx
del dataset_refs, deep_coadd

In [ ]: