103.9. Custom color images

103.9. Custom color images¶

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

Learning objective: To make a custom color composite image.

LSST data products: deep_coadd

Packages: lsst.daf.butler, lsst.geom, lsst.sphgeom, lsst.pipe.tasks.prettyPictureMaker, lsst.ip.diffim, lsst.pipe.tasks.makeDiscreteSkyMap, lsst.pipe.base, matplotlib.pyplot.

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.

Related tutorials: See the 103-series tutorial on multi-coadd mosaic cutouts, which introduces the multi-patch cutout technique used in Section 4.

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 tutorial demonstrates how to use the lsst.pipe.tasks.prettyPictureMaker code from the LSST Science Pipelines to produce high-quality RGB composite images from multi-band data.

The prettyPictureMaker task in the LSST Science Pipelines is specifically optimized for LSST images, incorporating algorithms that are hue preserving so that astrophysical colors remain consistent across brightness levels. Unlike many external packages that emphasize contrast at the expense of color fidelity, Rubin's approach ensures that faint objects retain distinguishable colors and bright objects do not lose astrophysical color information. By keeping hue constant across all brightness levels, the algorithm enables consistent comparison of astrophysical objects regardless of their distance or intrinsic brightness. In addition, it corrects for the non-uniform way the human eye perceives color, producing perceptually uniform hues that more accurately convey the information contained in astrophysical colors.

Section 2 builds and explains the PrettyPictureTask configuration object that controls image stretching, tone mapping, color balance, and other visual aspects. In Section 3, a cutout of an Euclid Quick Release 1 strong lens candidate in the Euclid Deep Field South (a Data Preview 2 [DP2] field) is made and rendered as a color composite. Finally, in Section 4, a color composite of the galaxy cluster PSZ2 G309.43-72.86 is displayed.

1.1. Import packages¶

The PrettyPictureTask and ChannelRGBConfig classes from lsst.pipe.tasks.prettyPictureMaker are used for creating color composites, while GetTemplateTask from lsst.ip.diffim, MakeDiscreteSkyMapTask from lsst.pipe.tasks.makeDiscreteSkyMap, and InMemoryDatasetHandle from lsst.pipe.base support working with image templates and sky maps.

In [1]:

import numpy as np
import matplotlib.pyplot as plt

from astropy.coordinates import SkyCoord
import astropy.units as u

from lsst.daf.butler import Butler

from lsst.geom import Box2I, Point2I, Extent2I, SpherePoint, degrees
from lsst.sphgeom import Box as SphBox

from lsst.pipe.tasks.prettyPictureMaker import PrettyPictureTask
from lsst.pipe.tasks.prettyPictureMaker._task import ChannelRGBConfig
from lsst.ip.diffim import GetTemplateTask
from lsst.pipe.tasks.makeDiscreteSkyMap import MakeDiscreteSkyMapTask
from lsst.pipe.base import InMemoryDatasetHandle

1.2. Define parameters and functions¶

Function: display_image_array

Displays an RGB color image. PrettyPictureTask returns its result as an lsst.images.ColorImage, whose .array attribute is the underlying HxWx3 numpy array; the function accepts either a ColorImage or a plain numpy array.

In [2]:

def display_image_array(image_array, ax=None, figsize=(8, 8), name=None):
    """Display an RGB color image.

    Parameters
    ----------
    image_array : `lsst.images.ColorImage` or `numpy.ndarray`
        An RGB color image, either as an ``lsst.images.ColorImage`` (as
        returned by ``PrettyPictureTask``) or as an HxWx3 numpy array.
    ax : matplotlib.axes.Axes, optional
        Axis to draw on; if None, a new Figure and Axes are created.
    figsize : tuple, optional
        Size of the created `matplotlib.Figure` (ignored if `ax` is provided).
    name : str, optional
        Label to annotate the image.
    """
    # PrettyPictureTask returns an lsst.images.ColorImage; use its numpy array.
    if hasattr(image_array, "array"):
        image_array = image_array.array

    if ax is None:
        fig = plt.figure(figsize=figsize)
        ax = fig.add_subplot(1, 1, 1)

    ax.imshow(image_array, interpolation="nearest", origin="lower")
    ax.axis("off")

    if name:
        ax.text(0, 0.95, name, color="white", fontsize=12, ha="left", va="top",
                transform=ax.transAxes)

    return

Create an instance of the butler, and assert that it exists.

In [3]:

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

Get the skymap. DP2 coadds use the lsst_cells_v2 skymap.

In [4]:

skymap = butler.get("skyMap", skymap="lsst_cells_v2")

2. PrettyPictureTask configuration¶

Define the configuration file for PrettyPictureTask, and set some parameters.

In [5]:

pretty_pic_config = PrettyPictureTask.ConfigClass()

2.1. imageRemappingConfig.absMax¶

This parameter sets the input data value (for Rubin data, in nJy) that will be treated as "white" in the output image. In this context, "white" means the maximum display intensity in the output RGB image. When an input pixel has a flux equal to imageRemappingConfig.absMax, it will be mapped to pure white (full intensity in all channels) after stretching and scaling. Brighter fluxes will also be clipped to white.

If left unset, the brightest pixel in the image might be used as the white point --- but this can be unstable due to rare, very bright pixels in the long tail of the brightness histogram.

In [6]:

pretty_pic_config.imageRemappingConfig.absMax = 11000

2.2. luminanceConfig.stretch¶

This parameter controls the strength of the arcsinh stretch applied to the luminance channel, as defined in Lupton et al. 2004. The arcsinh function compresses the dynamic range of the image while still preserving faint structures, making it well-suited for astronomical images where bright sources and faint details coexist. Smaller values of stretch produce a more aggressive compression that enhances faint features but can cause bright regions to saturate, while larger values approach a nearly linear scaling.

In [7]:

pretty_pic_config.luminanceConfig.stretch = 750

2.3. luminanceConfig.highlight, luminanceConfig.shadow, luminanceConfig.midtone¶

These parameters control a tone-mapping function applied after the arcsinh stretch, similar to the highlight, shadow, and midtone sliders in many photo-editing applications. They determine how the luminance values, now scaled to the 0–1 range, are mapped to the final display brightness.

The highlight value sets the luminance level that will be mapped to white, while shadow sets the level that will be mapped to black. The midtone parameter defines the balance point between white and black, affecting the overall contrast and brightness perception of the image. Adjusting these values allows fine-tuning of the tonal balance to emphasize particular brightness ranges without altering the color mapping.

In [8]:

pretty_pic_config.luminanceConfig.highlight = 0.905882
pretty_pic_config.luminanceConfig.shadow = 0.5
pretty_pic_config.luminanceConfig.midtone = 0.5

2.4. doPsfDeconvolve¶

This parameter controls whether the image will undergo Point Spread Function (PSF) deconvolution as part of the processing. Setting it to False disables PSF deconvolution, which can simplify and speed up the processing, especially for tutorial or exploratory work. When enabled, PSF deconvolution can sharpen features in the image by compensating for the blurring effects of the point-spread function, but it also increases computational complexity and may amplify noise in low signal-to-noise regions.

In [9]:

pretty_pic_config.doPsfDeconvolve = False

2.5. exposureBracketerConfig.exposureBrackets¶

This parameter enables the use of exposure bracketing, an advanced technique for further compressing the dynamic range of an image. Exposure brackets combine differently scaled versions of the same image to retain detail in both the brightest and faintest regions, similar to high dynamic range (HDR) processing in photography. While powerful, this feature can require careful tuning to achieve the desired visual effect. For simplicity and consistency in this tutorial, it is set to None, which disables bracketing and uses a single stretch for the entire image.

In [10]:

pretty_pic_config.exposureBracketerConfig.exposureBrackets = None

2.6. colorConfig.maxChroma and colorConfig.saturation¶

These parameters control how color intensity changes with brightness. Increasing maxChroma (> 0) and saturation produces richer, more vibrant colors compared to the default settings, enhancing visual appeal while preserving the relative color balance of objects in the image. The saturation factor ranges from 0 to 1, but values near 1 are not recommended because they can make bright pixels appear overly saturated.

In [11]:

pretty_pic_config.colorConfig.maxChroma = 80
pretty_pic_config.colorConfig.saturation = 0.6

2.7. channelConfig¶

This setting controls how each input band contributes to the red, green, and blue channels of the final RGB image. In the default configuration, the g band maps entirely to blue, the r band to green, and the i band to red, producing a straightforward g–r–i color composite (with the rgb values in the [0,1] range).

It's possible to assign any floating-point weights, allowing a band to contribute to multiple color channels. This can be used to fine-tune color balance or to include additional bands beyond three by mixing them into the RGB channels (e.g., adding u to blue and red to simulate violet, or mixing z into red). The mapping is applied uniformly across all brightness levels, so object colors remain consistent regardless of brightness. Keeping each band mapped to a single channel is simplest to interpret, but mixing can provide more flexibility and information.

In [12]:

pretty_pic_config.channelConfig = {}
pretty_pic_config.channelConfig["g"] = ChannelRGBConfig(r=0, g=0, b=1)
pretty_pic_config.channelConfig["r"] = ChannelRGBConfig(r=0, g=1, b=0)
pretty_pic_config.channelConfig["i"] = ChannelRGBConfig(r=1, g=0, b=0)

Define an instance of PrettyPictureTask with the configuration file.

In [13]:

task = PrettyPictureTask(config=pretty_pic_config)

3. Strong-lens candidate cutout¶

Make a cutout of a strong lens candidate from Euclid data in the Euclid Deep Field South (EDFS) DP2 field.

Define the approximate coordinates of the lens, in degrees. The coordinates are from "Euclid Quick Data Release (Q1): The Strong Lensing Discovery Engine". This example uses one target, but the workflow is structured to accept a list of positions.

In [14]:

ra_coordinates = [59.626134]
dec_coordinates = [-49.06175]

Define the stamp size, in pixels.

In [15]:

stamp_size = 32

Define the photometric bands of the input images used to create the composite color image.

In [16]:

bands = ["g", "r", "i"]

Format the bands into a query string.

In [17]:

bands_str = ",".join(f"'{b}'" for b in bands)

In [18]:

bands_str

Out[18]:

"'g','r','i'"

Find the tract, patch for each of the RA and Dec coordinate pairs (in this case, just one), using the coord_list_to_tract_patch function.

In [19]:

def coord_list_to_tract_patch(coord_list, skymap):
    """Create a list of tract, patch pairs for each coordinate in the input list.

    Parameters
    ----------
    coord_list : `list` of `lsst.geom.SpherePoint`
        List of coordinates for which the closest tract, patch region is to be
        found.
    skymap : `lsst.skymap.BaseSkymap`
        The skymap from where the tract, patch pairs are to be determined.

    Results
    -------
    results : `list` of `tuple` of `int`, `int`
        A list the same length as the input which each element is the closest
        tract patch to the corresponding coordinate. A tuple of -1, -1 will
        be returned if a coordinate can't be attributed to a tract, patch.

    Notes
    ----
    The patch number returned is the sequential index used to identify patches
    in the butler.
    """
    results = []
    for coord in coord_list:
        tractInfo = skymap.findTract(coord)
        for patchInfo in tractInfo:
            if patchInfo.getInnerSkyPolygon().contains(coord.getVector()):
                results.append((tractInfo.getId(), patchInfo.getSequentialIndex()))
                break
        else:
            results.append((-1, -1))
    return results

In [20]:

sph_points = [SpherePoint(r, d, degrees) for r, d in zip(ra_coordinates, dec_coordinates)]
skymap_tract_patch = coord_list_to_tract_patch(sph_points, skymap)

Sort into tract–patch keys so that all points can be processed in a single tract–patch combination for efficiency. Because of patch overlap regions, patch is a list of one or more patches, and the first one is used.

In [21]:

sorted_tract_patch = {}
for i, (tract, patch) in enumerate(skymap_tract_patch):
    container = sorted_tract_patch.setdefault((tract, patch), list())
    container.append(sph_points[i])

Loop over tract–patch pairs to create RGB images from coadd cutouts. For each patch, the relevant coadd datasets for the requested bands are fetched through the Butler, with duplicates removed using set. Each coadd is retrieved with storageClass="Exposure" so that getCutout returns an ExposureF (in DP2, the default deep_coadd storage class is a cell-based coadd). Since all bands in a tract–patch share the same World Coordinate System (WCS), one WCS is used to transform sky coordinates into pixel coordinates. For every target point, a bounding box is defined to extract a cutout from each band image. These cutouts are passed to the PrettyPictureTask, which produces an RGB image for display.

Note: The PrettyPictureTask run method natively takes a dictionary mapping of bands as keys to exposures. However, the task also supports arrays via task.makeInputsFromArrays, which can be used if the band data have already been preprocessed into numpy arrays.

Warning: Fetching cutouts individually can be inefficient because the Butler re-reads the images for every request; for multiple cutouts, it is usually faster to fetch the full image once and then slice it locally.

In [22]:

for (tract, patch), points in sorted_tract_patch.items():
    lens_refs = list(
        set(
            butler.query_datasets(
                "deep_coadd",
                patch=patch,
                tract=tract,
                where=f"band in ({bands_str})",
            )
        )
    )

    wcs = skymap[tract].wcs
    coadd_dict = {}
    for ref in lens_refs:
        coadd_dict[ref.dataId["band"]] = butler.get(ref, storageClass="Exposure")

    for point in points:
        xy = Point2I(wcs.skyToPixel(point))
        cutoutSize = Extent2I(stamp_size, stamp_size)
        bbox = Box2I(xy - cutoutSize // 2, cutoutSize)
        cutouts = {k: exp.getCutout(bbox) for k, exp in coadd_dict.items()}
        rgbImage = task.run(cutouts).outputRGB
        display_image_array(rgbImage, figsize=(8, 4), name="EUCLID QR1 Lens Candidate")

Figure 1: gri color-composite cutout of a Euclid QR1 strong lens candidate in DP2 data.

4. Galaxy cluster big cutout¶

PSZ2 G309.43-72.86 is a Sunyaev-Zel'dovich-selected galaxy cluster at redshift 0.35 that is included in the DP2 data release, in the ELAIS-S1 field. Create a color composite of the cluster to illustrate the use of PrettyPictureTask on a dense field of galaxies, combined with the multi-patch cutout technique introduced in the 103-series tutorial on multi-coadd mosaic cutouts.

Define the coordinates of the cluster center (the Sunyaev-Zel'dovich coordinates), in degrees.

In [23]:

ra_bcg = 10.196
dec_bcg = -44.134

Define a half-width of 0.05 degrees from the cluster center in declination, so the cutout spans 0.1 degrees on the sky (twice the half-width).

In [24]:

box_half_size_deg = 0.05

To obtain a cutout that is (approximately) square on the sky, widen the extent in right ascension (RA) by a factor of 1/cos(Dec), because lines of constant RA converge toward the poles. Without this correction, equal ranges in RA and Dec degrees map to different numbers of pixels, and the resulting image is rectangular rather than square.

In [25]:

cos_dec = np.cos(np.radians(dec_bcg))
ra_half_deg = box_half_size_deg / cos_dec
dec_half_deg = box_half_size_deg

In [26]:

box_region = SphBox.fromDegrees(ra_bcg - ra_half_deg, dec_bcg - dec_half_deg,
                                ra_bcg + ra_half_deg, dec_bcg + dec_half_deg)

dataset_refs = set(
    butler.query_datasets(
        "deep_coadd",
        where=f"patch.region OVERLAPS(region) and band in ({bands_str})",
        bind={"region": box_region}
    )
)

Prepare the inputs for building a custom skymap. For each dataset reference, the WCS and bounding box of its patch are stored, so the collection of positions can later define a skymap that includes only the data in use. At the same time, the dataset references are organized into a nested dictionary: first grouped by band, then within each band grouped by tract. Each tract key contains all the dataset refs for that tract band combination.

In [27]:

position_container = []
sorted_data = {}
for ds_ref in dataset_refs:
    did = ds_ref.dataId

    patchInfo = skymap[did["tract"]][did["patch"]]
    position_container.append((patchInfo.wcs, patchInfo.outer_bbox))

    band_container = sorted_data.setdefault(did["band"], dict())
    tract_container = band_container.setdefault(did["tract"], list())
    tract_container.append(ds_ref)

The goal is to create an image of an arbitrary sky region that is not aligned with the native LSST sky tessellation. This is similar to the process in alert production, where past LSST observations are combined to generate static-sky templates for detecting transients in new exposures. The same template-generation code can be used to transform input data aligned on the LSST tessellation into the required geometry.

To use the template code, a single WCS that encompasses all input data is needed, instead of having a separate tangent-plane WCS for each tract. A discrete skymap is therefore created for the input data, producing one tract with a single WCS projection. A bounding box is then defined in these new coordinates, centered on the cluster and spanning the specified width.

With this setup, images are generated per band by combining the inputs that fall within the bounding box. The GetTemplateTask creates coadds for each band, which are then combined into a final RGB mosaic.

Note: In DP2, the default storage class of a deep_coadd is a cell-based coadd (lsst.images.cells.CellCoadd), while GetTemplateTask expects afw Exposure objects. Retrieve each coadd with storageClass="Exposure" and wrap it in an InMemoryDatasetHandle, then pass the handles to the task through the coaddExposureHandles argument. This is a validated workaround for a temporary issue with passing storageClass="Exposure" to deferred (lazy) Butler handles; see the Rubin Community Forum discussion.

First, make a skymap which contains all input data using the MakeDiscreteSkyMapTask() task.

Override the default MakeDiscreteSkyMapTask() configuration values and replace them with those appropriate for LSSTCam, in particular the pixel scale and the World Coordinate System (WCS) projection.

Use the cluster location to obtain the pixel scale. First, get the WCS of the patch and tract where the cluster is located.

In [28]:

bcg_point = SpherePoint(ra_bcg * degrees,
                        dec_bcg * degrees)
bcg_tract = skymap.findTract(bcg_point)
bcg_patch = bcg_tract.findPatch(bcg_point)

bcg_tract_number = bcg_tract.tract_id
bcg_patch_number = bcg_patch.getSequentialIndex()
print(bcg_tract_number, bcg_patch_number)

2877 34

Get the data reference for the coadd in the g band that contains the cluster center.

In [29]:

bcg_data_ref = next(
    ref
    for ref in sorted_data[bands[0]][bcg_tract_number]
    if ref.dataId["patch"] == bcg_patch_number
)

Get the pixel scale at the cluster location from the coadd's sky_projection component. In DP2, the astrometric solution of a cell-based deep_coadd is stored as an lsst.images SkyProjection, retrieved through the deep_coadd.sky_projection component (see the DP2 astrometric calibration tutorial). The pixel scale is measured as the on-sky separation between two pixels that are one pixel apart.

In [30]:

coadd_wcs = butler.get("deep_coadd.sky_projection",
                       band=bands[0], skymap="lsst_cells_v2",
                       tract=bcg_tract_number, patch=bcg_patch_number)

target = SkyCoord(ra=ra_bcg * u.deg, dec=dec_bcg * u.deg)
xy = coadd_wcs.sky_to_pixel(target)
sky_00 = coadd_wcs.pixel_to_sky(x=xy.x, y=xy.y)
sky_10 = coadd_wcs.pixel_to_sky(x=xy.x + 1, y=xy.y)
plate_scale_bcg = sky_00.separation(sky_10).to(u.arcsec).value
print(f"Pixel scale: {plate_scale_bcg: .4f} arcsec/pixel")

Pixel scale:  0.2000 arcsec/pixel

Define an instance of the configuration class and override the plate scale.

In [31]:

discrete_skymap_config = MakeDiscreteSkyMapTask().ConfigClass()

In [32]:

discrete_skymap_config.skyMap.pixelScale = plate_scale_bcg

Set the WCS projection. The WCS of deep_coadd images are TAN projections (see the DP2 documentation).

In [33]:

discrete_skymap_config.skyMap.projection = 'TAN'

In [34]:

new_skymap = MakeDiscreteSkyMapTask(config=discrete_skymap_config).run(position_container).skyMap

lsst.makeDiscreteSkyMap INFO: Extracting bounding boxes of 6 images

lsst.makeDiscreteSkyMap INFO: Computing spherical convex hull

lsst.makeDiscreteSkyMap INFO: tract 0 has corners (10.480, -44.356), (9.859, -44.356), (9.861, -43.912), (10.478, -43.912) (RA, Dec deg) and 2 x 2 patches

Get the new WCS from the single tract discrete skymap.

In [35]:

new_wcs = new_skymap[0].wcs

Create a bounding box in the new WCS coordinates that is centered on the cluster and has the width defined above.

In [36]:

small_box = Box2I(
    Point2I(new_wcs.skyToPixel(SpherePoint(ra_bcg - ra_half_deg,
                                           dec_bcg - dec_half_deg, degrees))),
    Point2I(new_wcs.skyToPixel(SpherePoint(ra_bcg + ra_half_deg,
                                           dec_bcg + dec_half_deg, degrees))),
)

Loop over bands to create images that are combinations of the inputs contained within the defined box. Each coadd is retrieved with storageClass="Exposure" and wrapped in an InMemoryDatasetHandle, then passed to the task through coaddExposureHandles. Note that the run method modifies the bounding box in place, so a copy is made.

In [37]:

getTemplateTask = GetTemplateTask()
coadds = {}
for band, mapping in sorted_data.items():
    small_box_copy = Box2I(small_box)
    new_image = getTemplateTask.run(
        coaddExposureHandles={
            tract: [InMemoryDatasetHandle(butler.get(ref, storageClass="Exposure"))
                    for ref in refs]
            for tract, refs in mapping.items()
        },
        bbox=small_box_copy,
        wcs=new_wcs,
        dataIds={tract: [ref.dataId for ref in refs] for tract, refs in mapping.items()},
        physical_filter=band,
    ).template
    coadds[band] = new_image

Run the task to make the image.

In [38]:

mosaic_image = task.run(coadds).outputRGB

Display the image.

In [39]:

display_image_array(mosaic_image, name="PSZ2 G309.43-72.86")

Figure 2: gri color-composite big cutout of galaxy cluster PSZ2 G309.43-72.86.