206.2. Deblender footprints

Deblender footprints¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: Large
LSST Science Pipelines version: v30.0.11
Last verified to run: 2026-10-23
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: This notebook demonstrates how to use the deblended models generated by scarlet lite.

LSST data products: Object table, deep_coadd, object_scarlet_models

Packages: lsst.rsp, lsst.daf.butler, lsst.scarlet.lite, lsst.meas.extensions.scarlet

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. The deblender packages are based on Melchior et al. 2018, and the Scarlet software written by Peter Melchior and Fred Moolekamp. Further documentation on Scarlet is available at that site.

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 reconstruct per-pixel models of objects that are deblended by the LSST Science Pipeline when building the object table. The focus is on how to recover the model of an object from its information in the object catalog, and how to investigate the blend as a whole.

To measure the photometry of objects in deep_coadd images, the LSST Science Pipeline uses the multi-band deblending scarlet algorithm Melchior et al. 2018 implemented with optimizations for LSST data in Scarlet Lite. The algorithm models objects in a blend as a collection of components that have the same spectrum in each pixel of their morphology, assuming that flux monotonically decreases from the peak (maximum pixel value) of the component.

This tutorial will begin with the coordinates of an object of interest and show how to obtain its model and the model of the other sources that are extracted from the same blend.

Note: There have been updates to the deblending implementation since DP1 so expect the API to be slightly different (improved) in future data releases, including having more information available to understand how blended an object is with its neighbors and useful cuts to remove potentially more difficult blends from analysis.

References:

• A general discussion of blending impacts to LSST can be found in this article titled "The challenge of blending in large sky surveys" by Melchior et al. 2021.
• Some more in-depth presentations on the deblending implementation in the LSST Science Pipeline are available from recorded talks during the Rubin Project and Community Workshop 2022's session on "Deblending Plans and Challenges".
• Further information about the available deblending flags, and the processes that produce the deblender data products discussed in this tutorial can be found in this overview of the deblending flags.

Related tutorials: A notebook introducing the deblender data products stored in the Object table is available as 206.1 Introduction to deblender outputs. A notebook introducing image subsets and defining bounding boxes is available as 104.5 titled Image subsets with the butler.

1.1. Import packages¶

The matplotlib, numpy, and astropy libraries are widely used Python libraries for plotting, scientific computing, and astronomical data analysis.

lsst.afw.image provides access to some of the extended products created by the deblender pipeline. The lsst.afw.display and lsst.images libraries provides access to image visualization routines and the lsst.daf.butler library is used to access data products via the butler. The RSPDiscovery class in the lsst.rsp package is for accessing the TAP service. Finally, lsst.scarlet and lsst.meas.extensions.scarlet are two packages containing the scarlet software infrastructure, to enable displaying and analyzing deblended models and object footprints.

In [1]:

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

from astropy.coordinates import SkyCoord
from astropy.visualization import ZScaleInterval

from lsst.daf.butler import Butler
from lsst.afw.image import MultibandExposure
import lsst.afw.display as afw_display
from lsst.images import Box, Interval
from lsst.rsp import RSPDiscovery

import lsst.scarlet.lite as scl
import lsst.meas.extensions.scarlet as mes

1.2. Define parameters¶

Set matplotlib parameters to define the lower left corner of images as the origin.

In [2]:

matplotlib.rcParams['image.origin'] = 'lower'

Set the backend for afwDisplay to matplotlib.

In [3]:

afw_display.setDefaultBackend('matplotlib')

Get an instance of the TAP service.

In [4]:

discovery = RSPDiscovery("dp2")
service = discovery.get_tap_client()

Instantiate the butler.

In [5]:

butler = Butler('dp2', collections=["dp2"])
assert butler is not None

2. Identify a blend¶

First, identify a blended system in the DP2 data. This example will load data for the deblended object that is located at (151.78888574607532, 1.55284113175051) along with the relevant fields required to access the blend models associated with it.

• tract: the LSST track that contains the object
• patch: the patch in the tract that contains the object
• parentObjectId: the unique ID of the object's parent blend
• deblend_blendId: the index of the sub-blend inside the parent blend
• deblend_blendNChild: Describes the number of children that were in the parent blend
• deblend_peak_center_x,y: Not required but will be used later to mark the location of the object

2.1. Query for deblended object¶

Search for all objects within 1 arcsec of the object of interest. There should only be a single result.

In [6]:

ra = 151.78888574607532
dec = 1.55284113175051

query = f"""SELECT objectId, coord_ra, coord_dec, parentObjectId, 
            deblend_blendId, tract, patch, 
            deblend_blendNChild, deblend_peak_center_x, deblend_peak_center_y 
            FROM dp2.Object 
            WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), 
            CIRCLE('ICRS', {ra}, {dec}, 0.000278)) = 1"""

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
if job.phase == 'ERROR':
    job.raise_if_error()

results = job.fetch_result().to_table()
print(len(results))

1

Identify the parent object from which the target was deblended (using parentObjectId from the object catalog) and identify all the deblended children. This parent object was deblended into 11 separate objects.

In [7]:

parentId = results['parentObjectId'][0]

query = f"""SELECT objectId, coord_ra, coord_dec, parentObjectId, footprintArea,
            deblend_peak_center_x, deblend_peak_center_y
            FROM dp2.Object
            WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec ),
            CIRCLE('ICRS',  {ra}, {dec},  .1)) = 1
            AND parentObjectId = {parentId}"""


job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

assert job.phase == 'COMPLETED'
children = job.fetch_result().to_table()
children

Job phase is COMPLETED

Out[7]:

Table length=11

objectId	coord_ra	coord_dec	parentObjectId	footprintArea	deblend_peak_center_x	deblend_peak_center_y
	deg	deg		pix	pix	pix
int64	float64	float64	int64	int32	int32	int32
788017509725927499	151.7876581354054	1.5520202792623203	788017509725899960	1103	14061	2769
788017509725927503	151.78799725447507	1.5541580723790012	788017509725899960	648	14055	2808
788017509725927502	151.78940061329126	1.5513401386815453	788017509725899960	879	14030	2757
788017509725927501	151.7873381242345	1.5527339536780616	788017509725899960	1025	14067	2782
788017509725927500	151.78783048296668	1.5526027394298354	788017509725899960	1228	14058	2780
788017509725927498	151.78830479725625	1.5537066759440619	788017509725899960	1081	14049	2800
788017509725927497	151.78888574607532	1.55284113175051	788017509725899960	1588	14039	2784
788017509725927496	151.7888307832984	1.5518285256916162	788017509725899960	1647	14040	2766
788017509725927495	151.78923138010427	1.5490635184869188	788017509725899960	804	14033	2716
788017509725927494	151.78887003102082	1.5501733885410403	788017509725899960	1096	14039	2736
788017509725927493	151.78821708813376	1.549245434908877	788017509725899960	1156	14051	2720

2.2. Visualize the blend¶

To visualize the blend, plot a deep_coadd image subset that contains the blend with children identified. First, retrieve the patch and tract for the coadd containing the blend, and set the band to be i-band.

Next, define the dimensions of the image subset to be 5x the square root of the total footprintArea of the blend, in arcseconds. Since footprintArea is in pixels, multiply by the pixel scale 0.2 arcsec / pixel to define an appropriate area for the cutout in arcsec.

In [8]:

total_footprint_area = np.sum(children['footprintArea'])
cutoutSideLength = 5.0 * np.sqrt(total_footprint_area) * 0.2

coadd = butler.get("deep_coadd", tract=results['tract'][0], patch=results['patch'][0], band='i')

Define the center to be at the deblended target, and define the bounding box to construct the image subset.

In [9]:

center = SkyCoord(ra, dec, unit='deg')
xy = coadd.sky_projection.sky_to_pixel(center)
bbox = Box.factory[round(xy.y):round(xy.y)+1, round(xy.x):round(xy.x)+1].padded(cutoutSideLength)
subset = coadd[bbox]

Finally, visualize the blend and its child objects. First display the cutout, and then over plot circles to identify the child objects that were deblended.

In [10]:

display = afw_display.Display(frame=1)
display.scale('asinh', 'zscale')

display.mtv(subset.image)
with display.Buffering():
    for ci in range(len(children)):
        display.dot('o', children['deblend_peak_center_x'][ci],
                    children['deblend_peak_center_y'][ci], size=10)
plt.gca().axis('off')
display.show_colorbar(False)

plt.show()

Figure 1: The image subset of the i-band deep_coadd image displayed in greyscale, with green circles marking the deblended children.

3. Blend and object models¶

Scarlet Lite uses the multi-wavelength images to identify peaks. The scarlet models produced during deblending are available as a data product via the object_scarlet_model catalog. This contains all of the information needed to reconstruct the multi-band model for each source.

3.1. Load parent blend data¶

Scarlet lite models are stored in a zip archive with keys for metadata and each blend, allowing efficient retrieval of data for a single parent blend.

Below, find the object_scarlet_model data product in the butler by searching for all dataset types with the word "scarlet" in the name.

In [11]:

for dt in sorted(butler.registry.queryDatasetTypes('*scarlet*')):
    print(dt)

DatasetType('object_scarlet_models', {skymap, tract, patch}, LsstScarletModelData)

Use the butler to retrieve the object_scarlet_models data product that corresponds to the image that contains this blend.

In [12]:

obj = results[0]

model_data = butler.get(
    "object_scarlet_models",
    tract=obj["tract"],
    patch=obj["patch"],
    parameters={"blend_id": [obj["parentObjectId"]]}
)

3.2. Extracting the blend data¶

In DP2, all blends use the same model PSF (a Gaussian with sigma=0.8) and the same observed PSF (the PSF at the center of the patch), and the same bands (usually all bands, ugrizy, except on the edges of the survey or crowded fields where not all bands successfully processed). Their models can thus be extracted directly from the metadata (new in DP2). LsstScarletModelData is a container for all of the blends in a patch. Omitting the parameters={"blend_id": [obj["parentObjectId"]]} in the command above would load all of the models.

model_data (above) only loaded a single blend, but it is still returned from the Butler in the full container, which is needed because it contains the metadata required to reconstruct the scarlet lite Blend object.

Below, extract the model_psf, psf, and the bands in the blend. The LSST science pipelines use persistable dataclass objects as intermediaries between the storage containers and python classes. In the case of Scarlet Lite blends, this is the ScarletBlendData class. Some attributes of the blend model can be retrieved from model_data.

In order to reconstruct the blend, the model PSF from the deconvolved model space is also required. The model PSF is the same for all blends, so extract that here.

Also extract the set of blends associated with this parentObjectId.

In [13]:

model_psf = model_data.metadata["model_psf"]
psf = model_data.metadata["psf"]
bands = model_data.metadata["bands"]

Print the bands metadata associated with the blend model.

In [14]:

print(bands)

['g', 'i', 'r', 'u', 'y', 'z']

Below, extract the parent data from model_data. parent_data holds the full blend from the object_detection catalog. Then, extract the blend data from parent_data. blend_data holds the deconvolved blend of the object. A model can be extracted for each deblended child from the parent blend.

Multiple parent IDs can also be passed in the list in order to load multiple blends.

In [15]:

parent_data = model_data.blends[obj["parentObjectId"]]

blend_data = parent_data.children[obj["deblend_blendId"]]

Finally, convert the blend_data storage class into a scarlet lite Blend instance called blend. The minimal_data_to_blend method is used to access the blend model.

In [16]:

blend = blend_data.minimal_data_to_blend(model_psf[None], psf, bands)

Print the object IDs for all of the objects in the blend by looking at Blend.sources. Notice that the target object is the second one in the list.

In [17]:

for source in blend.sources:
    print(source.metadata["id"])

788017509725927496
788017509725927497
788017509725927498
788017509725927499
788017509725927500
788017509725927501
788017509725927502
788017509725927503

3.3. Plot the blend model¶

3.3.1. Extracting a model¶

To compare the extracted blend model to observations, first convolve it with the difference kernel so that it matches the observed seeing. This happens automatically when passing convolve=True to Blend.get_model.

In [18]:

model = blend.get_model(convolve=True)

3.3.2. Index slicing¶

It is useful to be able to change the order of the filters stored in the blend model.

Print the observed bands used in the deblending process. Notice that the observed bands are not in wavelength order, and they are still in the same order from model_data.bands above. This is because the LSST science pipelines process each band in alphabetical order, not wavelength order.

In [19]:

print(blend.observation.bands)

('g', 'i', 'r', 'u', 'y', 'z')

All of the band-dependent scarlet lite data products can be re-ordered using indexing by passing slices or tuples of bands, similar to numpy arrays, with the exception that the last band in the slice is included in the resulting sliced model.

In [20]:

blend = blend[tuple("ugrizy")]
print(blend.observation.bands)
print(blend["g":"i"].observation.bands)
print(blend["z"].observation.bands)

('u', 'g', 'r', 'i', 'z', 'y')
('g', 'r', 'i')
('z',)

To change to wavelength order, specify model[tuple("ugrizy")] or blend[tuple("ugrizy")] when retrieving the blend model to display the image with the expected color assignment. If fewer bands of coverage exist, simply replace this by the combination of observed bands, e.g. tuple("griz").

Demonstrate the new filter order below.

In [21]:

print(blend.observation.bands)

('u', 'g', 'r', 'i', 'z', 'y')

3.3.3. Visualize¶

To display multi-band image data as a color (red green blue; RGB) image, use Scarlet lite helper functions. Below, use the Lupton RGB color mapping to preserve colors and display with an asinh normalization. The asinh stretch preserves colors independent of brightness (for more information see Lupton et al. 2024).

In [22]:

stretch = 50
Q = 10
norm = scl.display.AsinhMapping(minimum=0, stretch=stretch, Q=Q)

By default, Scarlet Lite assumes the input filters are in wavelength order, and thus mapping into an RGB scaling where the first filter map to blue and last filter map to red. This relative weighting between filters to RGB can be altered with the channel_map keyword to img_to_rgb (see documentation here).

Below, convert the blend model into an RGB image and display it, and mark the object of interest with an x. (Scarlet lite also has helper functions to display a full blend, mark the location of peaks, show the residuals automatically. This will be demonstrated in Section 4.3).

In [23]:

extent = scl.display.get_extent(model.bbox)
rgb = scl.display.img_to_rgb(model, norm=norm)
plt.imshow(rgb, extent=extent)
plt.title("blend model")
plt.plot(obj["deblend_peak_center_x"], obj["deblend_peak_center_y"], "wx")
plt.show()

Figure 2: A ugrizy RGB color image of the blend model. The original target object from Section 2 is marked with a white x. By default, Scarlet Lite maps the filters into an RGB scaling where the first bands contribute to the blue color and last bands contribute to red.

3.4. Displaying a single source model¶

To display the model for the object, first extract the object from the blend.

In [24]:

source = next((src for src in blend.sources if
    src.metadata["id"] == obj["objectId"]), None)
print(source)

Source<1>

Recall that the object model is in the deconvolved space, so displaying it will show the internal scarlet lite model (intrinsic model) for the object.

In [25]:

rgb = scl.display.img_to_rgb(source.get_model(), norm=norm)
plt.imshow(rgb, extent=extent)
plt.title("Intrinsic (deconvolved) model")
plt.show()

Figure 3: An RGB (PSF-deconvolved) image of the initial deblended source of interest.

Now, to view the model as directly comparable to the data, first convolve it with the observation using the convolve method.

In [26]:

model = blend.observation.convolve(source.get_model())
rgb = scl.display.img_to_rgb(model, norm=norm)
plt.imshow(rgb, extent=extent)
plt.title("Intrinsic (convolved) model")
plt.show()

Figure 4: An RGB image of the model of the same initial deblended source in Figure 3, but now convolved with the PSF.

4. Compare to observation¶

This section will make visual comparisons between the blend model and the observed data, and demonstrate some of the Scarlet Lite functionality.

4.1. Loading a MultibandExposure¶

First, load a multi-band exposure for the region that covers the blend. To load the observation, the following information is required: 1) the tract and patch that defines the region of the sky imaged in the deep_coadd; 2) the bounding box of the parent blend (bbox), which allows extraction a smaller region from the full deep_coadd; and 3) the bands desired for display.

In [27]:

bbox = mes.utils.scarletBoxToBBox(blend.observation.bbox)

mCoadd = MultibandExposure.fromButler(
    butler,
    bands,
    "deep_coadd",
    tract=obj["tract"],
    patch=obj["patch"],
    storageClass="ExposureF",
    parameters={'bbox': bbox}
)

4.2. Build a full Blend with an Observation¶

To compare the Scarlet Lite deblended models with the observed data, first create an observation object of the Observation class. This bundles the actual images in each band (mCoadd) with the telescope's PSF (modelPsf). Scarlet Lite uses this to convolve the intrinsic model of the blend for comparison to the observed images.

In [28]:

observation = mes.utils.buildObservation(
    modelPsf=model_psf,
    psfCenter=bbox.getCenter(),
    mExposure=mCoadd,
)

Pass observation into the blend_data, to build full_blend, which is an instance of the scarlet.lite.Blend class that contains all this information together. This enables both the model and real observations to be compared directly using Scarlet Lite's built-in show_scene function.

In [29]:

observed_bands = tuple("ugrizy")
full_blend = blend_data.to_blend(observation)[observed_bands]

4.3. Display full blend residuals¶

Now, use the Scarlet Lite built in display function called show_scene to show the image, blend model, convolved blend model, and residuals between the model and observation.

In [30]:

scl.display.show_scene(
    full_blend,
    norm=norm,
    show_model=True,
    show_rendered=True,
    show_observed=True,
    show_residual=True,
    linear=False,
    add_boxes=True,
    add_labels=True,
    figsize=(10, 10),
)

plt.show()

Figure 5: Visualization of the model generated by Scarlet Lite, with x and y axes indicating pixel number from the patch/tract. The panels show the blend model (not PSF-convolved; top left), the blend model rendered (PSF-convolved; top right), the observation (bottom left) and the residual image (observation - model rendered; bottom right). All deblended sources are identified with white boxes and a number indicating its ID within the set of deblended objects.

4.4. Display source models¶

Now, use the Scarlet Lite built-in display function called show_sources to display the models of the individual sources compared to the observed image and its spectral energy distribution.

In [31]:

scl.display.show_sources(
    full_blend,
    norm=norm,
    show_observed=True,
    add_boxes=True
)
plt.show()

Figure 6: Visualization of the individual deblended sources (top to bottom show ID 0 to 7). For each source, left panel shows the PSF-deconvolved model footprint (where the color scale corresponds to the weighted RGB color), middle panel shows the observation of the blend with the source identified by a white box, and right panel shows the flux in each band. The bands are shown all 6 LSST bands, and are in order of increasing wavelength, defined in observed_bands. The units of intensity is nJy, as in the image (nJy/pixel).

4.5. Display the measurement models¶

In the LSST science pipeline, the final properties saved in the Object catalog are not measured on the idealized Scarlet models directly. Instead, the pipeline uses the Scarlet models as templates to determine what fraction of light in each pixel belongs to which object. The algorithm uses these fractions to redistribute the raw image flux, creating an isolated, flux-conserved image for each source. Finally, the pipeline runs its standard measurement tools on these isolated images to calculate the object's centroid position, total flux, and shape.

The first step is to conserve the flux by applying conserve_flux() to modify the full_blend object. Doing so will create a Source.flux_weighted_image for each Source in the full_blend, and save those internally.

In [32]:

full_blend.conserve_flux()

Now, pass use_flux=True to show_sources to see the flux weighted models (left panel) that were used in all of the measurement tasks. Notice that the spectrum in Figure 7 is now also different from in Figure 6. This is because Figure 7 is based on the summing the flux from each band the real images rather than based on the idealized models.

In [33]:

scl.display.show_sources(
    full_blend,
    norm=norm,
    show_observed=True,
    add_boxes=True,
    use_flux=True,
)
plt.show()

Figure 7: Same as Figure 6, except after applying the conserve_flux function, to show the flux-redistributed model that the science pipelines use for measurements. For each source, left panel shows the PSF-deconvolved model footprint (where the color scale corresponds to the weighted RGB color), middle panel shows the observation of the blend with the source identified by a white box, and right panel shows the flux in each band. The bands are shown all 6 LSST bands, and are in order of increasing wavelength, defined in observed_bands.

5. Analyzing a single object¶

Individual objects can be inspected. Here, select the second brightest object in the blend (listed as Source 1 in the output in Section 4.1) and access its objectId.

In [34]:

object1 = full_blend.sources[1]
print(object1.metadata["id"])

788017509725927497

As in Section 4.5, conserve the flux of the overall scene.

In [35]:

full_blend.conserve_flux()

Three different models for the object are available for inspection. 1) The deconvolved scarlet model; 2) the convolution of the scarlet model with the difference kernel; and 3) the flux-redistributed model that the science pipelines use for measurement. The section demonstrates how to plot these three different models for source 1 in the blend.

First, extract the deconvolved model from the object. The output is a Scarlet Lite image.

In [36]:

source_model = object1.get_model()
display_bands = tuple("griz")

Next, project the model into a box large enough to fit the convolved model.

In [37]:

psf_radius = int((np.max(blend.observation.psfs.shape)+1)//2)
enlarged_box = source_model.bbox.grow(psf_radius)
convolved_model = source_model.copy().project(bbox=enlarged_box)

Finally, plot the three models for the object for inspection.

In [38]:

fig, ax = plt.subplots(1, 3, figsize=(15, 5))

extent = scl.display.get_extent(object1.bbox)
rgb = scl.display.img_to_rgb(source_model[display_bands], norm=norm)
ax[0].imshow(rgb, origin='lower', extent=extent)
ax[0].set_title("Deconvolved Model")

convolved_model = blend.observation.convolve(convolved_model)
extent = scl.display.get_extent(enlarged_box)
rgb = scl.display.img_to_rgb(convolved_model[display_bands], norm=norm)
ax[1].imshow(rgb, origin='lower', extent=extent)
ax[1].set_title("Convolved Model")

flux_model = object1.flux_weighted_image
extent = scl.display.get_extent(flux_model.bbox)
rgb = scl.display.img_to_rgb(flux_model[display_bands], norm=norm)
ax[2].imshow(rgb, origin='lower', extent=extent)
ax[2].set_title("Flux Re-distributed Model")

plt.show()

Figure 8: Visualization of deblended source 1. Left panel shows the PSF-deconvolved model footprint (where the color scale corresponds to the weighted RGB color), middle panel shows the convolution of the scarlet model with the difference kernel and right panel shows the flux-redistributed model that the science pipelines use for measurements.