202.1. Deep coadd images

202.1. Deep coadds¶

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-23
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: 925460

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, matplotlib plotting tools, and Astropy coordinate utilities.

In [1]:

from lsst.daf.butler import Butler
import lsst.afw.display as afwDisplay
from lsst.images._cell_grid import CellIJ

from astropy.coordinates import SkyCoord
import astropy.units as u
import matplotlib.pyplot as plt
from matplotlib.collections import PatchCollection
import numpy as np

1.2. Define parameters and functions¶

Instantiate the Butler.

In [2]:

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

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 CellCoadd. The LSST Science Pipelines tools for image manipulation and visualization work best with the CellCoadd 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}, CellCoadd)

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'
dataset_refs = butler.query_datasets("deep_coadd",
                                     where="band.name = :band AND \
                                     patch.region OVERLAPS POINT(:ra, :dec)",
                                     bind={"band": band, "ra": ra, "dec": dec})

print(len(dataset_refs))
del band

1

The query returns only a single $r$-band deep_coadd image overlapping the RA, Dec coordinates.

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 CellCoadd.

In [8]:

deep_coadd

Out[8]:

CellCoadd([y=2850:6150, x=14850:18150], tract=5063)

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_v2', tract: 5063, patch: 15}

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, band, and skymap¶

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

In [15]:

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

r 5063 15 lsst_cells_v2

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, patch, and skymap.

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

In [16]:

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

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 15. The colors correspond to an informative mask flag value (Section 3.3). A large fraction of the pixels have a color (have 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 deep_coadd, and convert it to an array.

In [20]:

image = deep_coadd.image
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())

-150.19159
3.5985272
48950.227

Display the image in Firefly frame 2.

In [24]:

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

Mouse-over display frame 2. The image plane contains less information than the full CellCoadd, but nonetheless displays coordinates in RA, Dec, as it does carry a WCS.

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.variance
variance_array = variance.array

Print statistics on the variance_array.

In [26]:

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

0.23627509
inf
inf

The mean and maximum values printed to the screen are inf (infinity). This is because masked pixels (e.g., the saturated pixels at the centers of bright stars) can have no data in the deep_coadd, which causes the variance for those pixels to be infinite.

Flag the pixels with inf values, count how many there are, and recalculate statistics with those pixels removed.

In [27]:

inf_variance = (variance_array == np.inf)
print(f"Number of pixels with infinite variance: {np.sum(inf_variance)}")

print("\nStatistics without inf pixels:")
print(variance_array[~inf_variance].min())
print(variance_array[~inf_variance].mean())
print(variance_array[~inf_variance].max())

Number of pixels with infinite variance: 9

Statistics without inf pixels:
0.23627509
25.55358
6961.53

Display the variance. Note that the default auto-scaling of the image stretch will cause the image to look "blank." Change the stretch (to, for example, "Z Scale Linear Stretch") to see detail in the variance image.

In [28]:

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

mask = deep_coadd.mask
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 their descriptions.

In [30]:

print(mask.schema)

NO_DATA [0@0x1]: No data was available for this pixel.
INTERPOLATED [0@0x2]: Pixel value is the result of interpolating nearby good pixels.
COSMIC_RAY [0@0x4]: A cosmic ray affected this pixel on at least one input image (and was interpolated).
SATURATED [0@0x8]: More than 10% of the potential input visits had a saturated pixel at this location ('potential' because saturated pixel values are not actually propagated to the coadd). SATURATED always implies REJECTED, and is often a reason for NO_DATA.
DETECTION_EDGE [0@0x10]: Pixel was too close to the edge of the patch to be considered for detection, due to the finite size of the detection kernel.
CLIPPED [0@0x20]: Region was identified as a probable artifact when comparing multiple single-visit warps. CLIPPED always implies REJECTED.
REJECTED [0@0x40]: At least one input visit was left out of the coadd for this pixel due to masking. REJECTED always implies INEXACT_PSF.
DETECTED [0@0x80]: Pixel was part of a detected source.
INEXACT_PSF [1@0x1]: The set of visits contributing to this pixel differs from the set of visits contributing to the PSF model for its cell.

Display the mask in Firefly frame 4.

In [31]:

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

afw_display.setMaskTransparency(0)

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

In [33]:

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

Option to print a list of the mask planes and the number of pixels with each mask plane set.

In [34]:

# for plane in mask.schema:
#     print(plane.name, np.sum(mask.get(plane.name)))

3.4. Background¶

Coadds have a final round of background subtraction applied (in addition to the background subtraction that was already applied to each of the input visits). Extract the backgrounds dict and examine the contents.

In [35]:

backgrounds = deep_coadd.backgrounds

for key in backgrounds.keys():
    print(key)

pretty
object

Print the description of each background type to the screen.

In [36]:

bg_obj = deep_coadd.backgrounds['object']
bg_pretty = deep_coadd.backgrounds['pretty']

print("object:\n", bg_obj.description)
print('\n')
print("pretty:\n", bg_pretty.description)

object:
 Background subtracted from the image when generating the Object catalog. This intentionally oversubtracts the background to reduce blending and ensure scattered light is subtracted. Restoring this background does not restore all original backgrounds, as the coadd was built from background-subtracted visit images; in most cases this background term is actually quite small.


pretty:
 An alternate background optimized for visually attractive RGB images. 'Subtracting' this background will generally *add* flux, correcting for most oversubtraction problems, but leaving in scattered light and some instrumental backgrounds in some cases. Because this background is fit to a difference of two wholly different coadds that may have different input images, it can also be pulled up or down by bright variable or transient objects.

Render the backgrounds as images and display them.

In [37]:

bg_obj_image = bg_obj.field.render()
bg_pretty_image = bg_pretty.field.render()

In [38]:

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

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

Make a copy of the deep_coadd, then use the apply_background method to restore some of the background that is often oversubtracted by the default behavior. Display this "restored" pretty_coadd in frame 1.

In [39]:

pretty_coadd = deep_coadd.copy()
pretty_coadd.apply_background('pretty')

In [40]:

afw_display = afwDisplay.Display(frame=1)
afw_display.mtv(pretty_coadd)

Use the Image alignment drop-down in Firefly to lock the images by pixel origin, then examine various regions of the images to compare the pretty_coadd (now in frame 1) to the deep_coadd image (in frame 2) and the backgrounds in frames 3 and 4.

Clean up.

In [41]:

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

4. Image attributes¶

4.1. Cell-based coadds¶

4.1.1. Input visits¶

Information about the input visits that contributed to the deep_coadd.

The deep_coadd images are made up of a grid of 22x22 150-pixel "cells", within which the PSF may be considered "edge-free" (meaning it was created in a way that avoids any discontinuities). The cells are created separately, then stitched together into a deep_coadd. Thus the set of inputs for each cell can be slightly different.

Get the provenance info about each cell and inspect it. Cells are identified by indices i and j.

In [42]:

contrib = deep_coadd.provenance.contributions.to_pandas()

In [43]:

contrib

Out[43]:

	cell_i	cell_j	visit	detector	overlaps_center	overlap_fraction	unmasked_fraction	weight	psf_shape_xx	psf_shape_yy	psf_shape_xy	psf_shape_flag	instrument
0	0	0	2025072800396	97	True	1.0	0.992044	0.004219	7.407177	5.427278	-0.020960	False	LSSTCam
1	0	0	2025072800397	97	True	1.0	0.986889	0.004169	8.689157	6.845690	-0.636649	False	LSSTCam
2	0	0	2025072800398	97	True	1.0	0.997067	0.004197	7.595631	6.659030	-0.131480	False	LSSTCam
3	0	0	2025072800399	97	True	1.0	0.978044	0.004235	6.482252	6.736140	-0.323755	False	LSSTCam
4	0	0	2025072800400	97	True	1.0	0.983867	0.004250	5.927198	5.833827	-0.329723	False	LSSTCam
...	...	...	...	...	...	...	...	...	...	...	...	...	...
10382	21	21	2026010500066	94	True	1.0	1.000000	0.001420	3.044833	3.736348	0.369175	False	LSSTCam
10383	21	21	2026010500067	94	True	1.0	1.000000	0.001420	2.789960	3.681143	0.449140	False	LSSTCam
10384	21	21	2026010500069	94	True	1.0	1.000000	0.001409	2.926084	3.519976	0.265914	False	LSSTCam
10385	21	21	2026010500070	94	True	1.0	1.000000	0.001411	3.288111	3.596816	0.277896	False	LSSTCam
10386	21	21	2026010500071	91	True	1.0	1.000000	0.001398	3.399342	4.159896	0.539613	False	LSSTCam

10387 rows × 13 columns

Make a histogram of the number of inputs contributing to each cell. The pandas "groupby" function makes it easy to group combinations of ij indices.

In [44]:

plt.hist(contrib.groupby(['cell_i', 'cell_j']).count()['visit'], bins=np.arange(11, 33, 1.0),
         histtype='step')
plt.minorticks_on()
plt.xlabel('number of input visits per cell')
plt.ylabel('number of cells')
plt.show()

Figure 2: Histogram showing the number of input visits contributing to each cell of the CellCoadd.

Get the list of input images that were combined to generate the deep_coadd. Note: this list contains all inputs that contributed anywhere in the patch.

Display the first 3 rows of the inputs table.

In [45]:

inputs = deep_coadd.provenance.inputs
print(f"Number of visits contributing to coadd: {len(inputs)}")
inputs[:3]

Number of visits contributing to coadd: 98

Out[45]:

Table length=3

visit	detector	day_obs	polygon
uint64	uint16	uint32	object
2025072800396	93	20250728	POLYGON ((18149.5 3789.35009765625, 18048.044921875 3895.27978515625, 18149.5 3992.398681640625, 18149.5 3789.35009765625))
2025072800396	94	20250728	POLYGON ((18149.5 4170.5751953125, 17960.4140625 3988.98828125, 15884.8779296875 6149.5, 18149.5 6149.5, 18149.5 4170.5751953125))
2025072800396	96	20250728	POLYGON ((16961.666015625 2849.5, 17888.99609375 3738.50341796875, 18149.5 3466.904052734375, 18149.5 2849.5, 16961.666015625 2849.5))

Make a spatial plot in tract xy coordinates, color-coding each cell based on the number of input visits that contributed to that cell's coadd.

In [46]:

unique_combinations = contrib[['cell_i', 'cell_j']].drop_duplicates()
ninputs = []
fig, ax = plt.subplots(figsize=(7, 7))

ninputs = contrib.groupby(['cell_i', 'cell_j']).size().to_numpy()

normal = plt.Normalize(np.min(ninputs), np.max(ninputs))
cmap = plt.cm.viridis
c = cmap(normal(ninputs))

patches = []

for i in range(len(unique_combinations)):
    cell_i = unique_combinations.iloc[i]['cell_i']
    cell_j = unique_combinations.iloc[i]['cell_j']
    pick_cell = (contrib['cell_i'] == cell_i) & (contrib['cell_j'] == cell_j)
    cell_bbox = deep_coadd.grid.bbox_of(CellIJ(cell_i, cell_j))
    patches.append(plt.Rectangle((cell_bbox.x.start, cell_bbox.y.start),
                                 150, 150, facecolor=c[i], fill=True))

col = PatchCollection(patches)
col.set(array=ninputs, cmap=cmap)

ax.add_collection(col)
ax.autoscale()
fig.colorbar(col, label='number of inputs')
ax.set_aspect('equal')

plt.xlabel('tract x (pix)')
plt.ylabel('tract y (pix)')
plt.show()

4.1.2. Info about the cells¶

Explore the cell grid associated with the deep_coadd to learn how to access cell information.

First, extract the xy coordinate of the patch corner, yx0.

In [47]:

xy0 = deep_coadd.yx0

Extract the grid indices of the cell corresponding to an xy coordinate, remembering to offset to tract xy.

In [48]:

xcoord = 1914 + xy0.x
ycoord = 1200 + xy0.y
print(f"x: {xcoord}, y: {ycoord}")

cij = deep_coadd.grid.index_of(x=xcoord, y=ycoord)
print(cij)

x: 16764, y: 4050
CellIJ(i=8, j=12)

Extract the bounding box of the cell indices retrieved above.

In [49]:

deep_coadd.grid.bbox_of(cij)

Out[49]:

Box(y=Interval(start=4050, stop=4200), x=Interval(start=16650, stop=16800))

Print the grid size, then the dict of overall parameters defining the grid.

In [50]:

deep_coadd.grid.grid_size

Out[50]:

CellIJ(i=22, j=22)

In [51]:

grid_dict = deep_coadd.grid.model_dump()

for key in grid_dict.keys():
    print(key, grid_dict[key], '\n')

bbox {'y': {'start': 2850, 'stop': 6150}, 'x': {'start': 14850, 'stop': 18150}}

cell_shape {'y': 150, 'x': 150}

In [52]:

del inputs, grid_dict, xy0, xcoord, ycoord

4.2. World Coordinate System (WCS)¶

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

4.2.1. sky_projection¶

There are multiple encodings of the mapping between sky and pixel coordinates. The most precise version is the sky_projection. Retrieve the sky_projection.

In [53]:

wcs = deep_coadd.sky_projection

Convert pixel coordinates to sky coordinates, and vice versa.

In [54]:

coord = SkyCoord(ra=ra*u.deg, dec=dec*u.deg, frame='icrs')
xy = wcs.sky_to_pixel(coord)
coord

Out[54]:

<SkyCoord (ICRS): (ra, dec) in deg
    (53.076, -28.11)>

In [55]:

print(xy)

XY(x=15182.507513158846, y=4390.517978014746)

Always pass xy values in tract coordinates when converting from pixel to sky with sky_projection.

In [56]:

xy2 = (15000, 4300)
print("x, y in tract coordinates: ", xy2[0], xy2[1])
print("Corner coordinates of patch: ", deep_coadd.yx0)
print("x, y in patch coordinates: ", xy2[0]-deep_coadd.yx0.x, xy2[1]-deep_coadd.yx0.y)
coord2 = wcs.pixel_to_sky(x=xy2[0], y=xy2[1])

x, y in tract coordinates:  15000 4300
Corner coordinates of patch:  YX(y=2850, x=14850)
x, y in patch coordinates:  150 1450

In [57]:

print(coord2)

<SkyCoord (ICRS): (ra, dec) in deg
    (53.08749462, -28.11502872)>

4.2.2. astropy_wcs¶

For users' convenience, an Astropy-compatible version of the WCS is also provided.

In [58]:

wcs_astropy = deep_coadd.astropy_wcs

Convert the SkyCoord from the previous section ("coord2") to xy using the Astropy WCS, and compare the results to see that Astropy returns xy in local (patch) pixel values.

In [59]:

print("Input xy: ", xy2)
xy2_astropy = wcs_astropy.world_to_pixel(coord2)
print("Output xy: ", xy2_astropy)
print("Output xy in tract coords: ",
      xy2_astropy.x+deep_coadd.yx0.x,
      xy2_astropy.y+deep_coadd.yx0.y)

Input xy:  (15000, 4300)
Output xy:  XY(x=149.99999999994907, y=1450.0000000000055)
Output xy in tract coords:  14999.999999999949 4300.0000000000055

Apply the Astropy pixel_to_world conversion, remembering to pass xy coordinates in local (patch) coordinates.

In [60]:

coord2_astropy = wcs_astropy.pixel_to_world(xy2[0]-deep_coadd.yx0.x,
                                            xy2[1]-deep_coadd.yx0.y)
print(coord2_astropy)

<SkyCoord (ICRS): (ra, dec) in deg
    (53.08749462, -28.11502872)>

The RA, Dec coordinate agrees with the result from the previous section.

4.2.3. fits_wcs¶

There is also a WCS in standard FITS format, provided as fits_wcs. For deep_coadd images this WCS is an exact representation of the LSST WCS for the image, but for other image types it may be an approximation to the WCS.

Read this WCS and print it to the screen.

In [61]:

wcs_fits = deep_coadd.fits_wcs

In [62]:

wcs_fits

Out[62]:

WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN' 'DEC--TAN'
CUNIT : 'deg' 'deg'
CRVAL : 53.08755760368661 -27.520661157024808
CRPIX : 150.0 12150.0
CD1_1 CD1_2  : -5.5555555555553004e-05 0.0
CD2_1 CD2_2  : 0.0 5.555555555555367e-05
NAXIS : 3300  3300

Print the object's associated methods to confirm that it also includes conversions between sky to pixel coordinates.

In [63]:

dir(wcs_fits)

Out[63]:

['__abstractmethods__',
 '__class__',
 '__copy__',
 '__deepcopy__',
 '__delattr__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__firstlineno__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__getitem__',
 '__getstate__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__iter__',
 '__le__',
 '__lt__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__static_attributes__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_abc_impl',
 '_all_pix2world',
 '_all_world2pix',
 '_array_converter',
 '_as_mpl_axes',
 '_denormalize_sky',
 '_det2im',
 '_fix_ctype',
 '_fix_pre2012_scamp_tpv',
 '_fix_scamp',
 '_get_components_and_classes',
 '_get_naxis',
 '_init_kwargs',
 '_naxis',
 '_normalize_sky',
 '_out_of_bounds_to_nan',
 '_p4_pix2foc',
 '_pix2foc',
 '_pixel_bounds',
 '_preserve_units',
 '_read_d2im_old_format',
 '_read_det2im_kw',
 '_read_distortion_kw',
 '_read_sip_kw',
 '_remove_sip_kw',
 '_write_det2im',
 '_write_distortion_kw',
 '_write_sip_kw',
 'all_pix2world',
 'all_world2pix',
 'array_index_to_world',
 'array_index_to_world_values',
 'array_shape',
 'axis_correlation_matrix',
 'axis_type_names',
 'calc_footprint',
 'celestial',
 'copy',
 'cpdis1',
 'cpdis2',
 'deepcopy',
 'det2im',
 'det2im1',
 'det2im2',
 'dropaxis',
 'fix',
 'footprint_contains',
 'footprint_to_file',
 'get_axis_types',
 'has_celestial',
 'has_distortion',
 'has_spectral',
 'has_temporal',
 'is_celestial',
 'is_spectral',
 'is_temporal',
 'low_level_wcs',
 'naxis',
 'p4_pix2foc',
 'pix2foc',
 'pixel_axis_names',
 'pixel_bounds',
 'pixel_n_dim',
 'pixel_scale_matrix',
 'pixel_shape',
 'pixel_to_world',
 'pixel_to_world_values',
 'preserve_units',
 'printwcs',
 'proj_plane_pixel_area',
 'proj_plane_pixel_scales',
 'reorient_celestial_first',
 'serialized_classes',
 'sip',
 'sip_foc2pix',
 'sip_pix2foc',
 'slice',
 'spectral',
 'sub',
 'swapaxes',
 'temporal',
 'to_fits',
 'to_header',
 'to_header_string',
 'wcs',
 'wcs_pix2world',
 'wcs_world2pix',
 'world_axis_names',
 'world_axis_object_classes',
 'world_axis_object_components',
 'world_axis_physical_types',
 'world_axis_units',
 'world_n_dim',
 'world_to_array_index',
 'world_to_array_index_values',
 'world_to_pixel',
 'world_to_pixel_values']

In [64]:

del xy, xy2, coord, coord2

4.3. Bounding box (corners)¶

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

Get the bounding box.

In [65]:

bbox = deep_coadd.bbox

Print the deep_coadd corners in pixels and sky coordinates.

In [66]:

bbox.min.x

Out[66]:

14850

In [67]:

corners = ((bbox.min.x, bbox.min.y),
           (bbox.min.x, bbox.max.y),
           (bbox.max.x, bbox.max.y),
           (bbox.max.x, bbox.min.y))
for corner in corners:
    print("corner (x, y): ", corner, '\n', wcs.pixel_to_sky(x=corner[0], y=corner[1]))

corner (x, y):  (14850, 2850)
 <SkyCoord (ICRS): (ra, dec) in deg
    (53.09694922, -28.19557406)>
corner (x, y):  (14850, 6149)
 <SkyCoord (ICRS): (ra, dec) in deg
    (53.09693349, -28.01231544)>
corner (x, y):  (18149, 6149)
 <SkyCoord (ICRS): (ra, dec) in deg
    (52.88934327, -28.01217359)>
corner (x, y):  (18149, 2850)
 <SkyCoord (ICRS): (ra, dec) in deg
    (52.88901088, -28.19543113)>

Check if pixel values are within the bounding box.

In [68]:

assert bbox.contains(x=100, y=100) is False
print(bbox.contains(x=100, y=100))
assert bbox.contains(x=16000, y=4000) is True
print(bbox.contains(x=16000, y=4000))

False
True

Check if sky coordinates are within the bounding box.

In [69]:

coord = SkyCoord(ra=63.0*u.deg, dec=-38.0*u.deg, frame='icrs')
coord2 = SkyCoord(ra=53.0*u.deg, dec=-28.1*u.deg, frame='icrs')
xy = wcs.sky_to_pixel(coord)
xy2 = wcs.sky_to_pixel(coord2)
print(bbox.contains(x=xy[0], y=xy[1]))
print(bbox.contains(x=xy2[0], y=xy2[1]))

False
True

Clean up.

In [70]:

del corners, coord, coord2, xy

4.4. 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 [71]:

psf = deep_coadd.psf

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

In [72]:

xy = (16000, 4000)
psf_image = psf.compute_kernel_image(x=xy[0], y=xy[1])
afw_display = afwDisplay.Display(frame=5)
afw_display.mtv(psf_image)

Clean up.

In [73]:

del psf, xy, psf_image

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

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

5063 15

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

In [75]:

butler.registry.queryDatasetTypes('object')

Out[75]:

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

use_columns = ['objectId', 'patch',
               'coord_ra', 'coord_dec']

Retrive the object table as objects.

In [77]:

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

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

In [78]:

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

9058

Option to show the table.

In [79]:

# objects[tx]

Redisplay the deep_coadd in frame 1 with no mask.

In [80]:

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

sc_patch = SkyCoord(ra=objects['coord_ra']*u.deg,
                    dec=objects['coord_dec']*u.deg,
                    frame='icrs')

with afw_display.Buffering():
    for i in tx:
        xy_tmp = wcs.sky_to_pixel(sc_patch[i])
        afw_display.dot('o',
                        xy_tmp.x,
                        xy_tmp.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 [82]:

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

In [ ]: