105.5. Background subtraction

105.5. Background subtraction¶

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-26
Repository: github.com/lsst/tutorial-notebooks

Learning objective: Re-estimate the background for a processed visit image using LSST Science Pipelines tasks.

LSST data products: visit_image, Object

Packages: lsst.daf.butler, lsst.rsp, lsst.meas.algorithms, lsst.afw.display, lsst.afw.image, lsst.pipe.tasks

Credit: Originally developed by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

This notebook demonstrates one approach to custom background subtraction at the invidiual visit image level using the LSST Science Pipelines. By default, visit images are background-subtracted; retrieving both the visit image and its corresponding background image is necessary to apply a custom subtraction. Custom background subtraction is valuable when the pipeline-estimated background is insufficient, particularly near low-surface-brightness features or extended galaxies in crowded or structured fields.

The notebook presents a basic workflow as an example; it does not evaluate alternative background subtraction methods. Users should adopt or develop strategies suited to their specific scientific goals.

Related tutorials: The 100-level butler tutorial to retrieve image data with the butler. The 100-level image reprocessing tutorial to perform forced photometry on sources using IDs and coordinates.

1.1. Import packages¶

Import matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery).

Also import lsst packages. These imports provide access to LSST data (Butler), TAP queries (get_tap_service), background subtraction (SubtractBackgroundTask), image handling (afwImage, afwDisplay), and forced photometry (ForcedMeasurementDriverTask).

In [1]:

import matplotlib.pyplot as plt

from lsst.daf.butler import Butler
from lsst.rsp import get_tap_service
from lsst.meas.algorithms import SubtractBackgroundTask
import lsst.afw.display as afwDisplay
import lsst.afw.image as afwImage
from lsst.pipe.tasks.measurementDriver import (
    ForcedMeasurementDriverConfig,
    ForcedMeasurementDriverTask,
)

1.2. Define parameters¶

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

In [2]:

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

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

In [3]:

service = get_tap_service("tap")
assert service is not None

objectId: Use objectId = 592913157106704436 for an example galaxy in this tutorial.

In [4]:

objectId = 592913157106704436

Set some parameters of matplotlib.pyplot to use a color-blind–friendly palette and improve overall plot aesthetics.

In [5]:

plt.style.use('tableau-colorblind10')

In [6]:

params = {'axes.labelsize': 20,
          'axes.titlesize': 12,
          'font.size': 15,
          'legend.fontsize': 15
          }
plt.rcParams.update(params)

2. Query images and coordinates¶

2.1. Retrieve images corvering the target galaxy¶

Retrieve the visit_image and visit_image_background using the visit identifier 2024111800093 and detector 2 for this example.

In [7]:

vi_id = {'visit': 2024111800093,
         'skymap': 'lsst_cells_v1',
         'detector': 2}

In [8]:

vi = butler.get("visit_image", vi_id)
vi_bg = butler.get("visit_image_background", vi_id)

2.2. Retrieve coordinates for the target galaxy¶

Obtain the coordinates of the target galaxy.

In [9]:

query = """
        SELECT objectId, coord_ra, coord_dec
        FROM dp1.Object
        WHERE objectId = {}
        """.format(objectId)

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

Job phase is COMPLETED

In [10]:

results = job.fetch_result().to_table()

Convert the galaxy’s sky coordinates to pixel coordinates on the image.

In [11]:

wcs = vi.getWcs()
xy = wcs.skyToPixelArray(results['coord_ra'], results['coord_dec'], degrees=True)

2.3. Display images¶

The afwDisplay package provides a nice interface for plotting the mask on top of the image using the maskedImage property.

Plot the zoomed-in processed visit image and its corresponding zoomed-in background image, and overplot the galaxy position on both.

In [12]:

x_slice = slice(2500, 3500)
y_slice = slice(2500, 3500)

Make a copy of the exposure to avoid altering the original. Extract its MaskedImage. Add the default background model back in to restore the image to its state with the original background subtraction applied.

In [13]:

vi_copy = vi.clone()
mi_default = vi_copy.maskedImage
mi_default += vi_bg.getImage()

In [14]:

fig, axs = plt.subplots(1, 2, figsize=(14, 7), sharey=True)

plt.sca(axs[0])
display1 = afwDisplay.Display(frame=fig, backend='matplotlib')
display1.scale('linear', 'zscale')
display1.mtv(vi.maskedImage[x_slice, y_slice])
axs[0].scatter(xy[0], xy[1], c='r', marker='o', label='galaxy')
axs[0].set_title('Background-subtracted visit image')
axs[0].set_xlabel('X')
axs[0].set_ylabel('Y')
axs[0].legend()

plt.sca(axs[1])
display2 = afwDisplay.Display(frame=fig, backend='matplotlib')
display2.scale('linear', 'zscale')
display2.mtv(vi_copy.maskedImage[x_slice, y_slice])
axs[1].scatter(xy[0], xy[1], c='r', marker='o')
axs[1].set_title('Visit image with restored background')
axs[1].set_xlabel('X')

plt.show()

Figure 1: Zoomed-in visit_image around the target galaxy (left) and its background-restored version using the default background (right). The red dot denotes the center of the target galaxy.

3. Re-estimate background¶

3.1. Define the configuration¶

The bkgConfig instance defined below specifies the background modeling and subtraction settings used in LSST image processing.

In [15]:

bkgConfig = SubtractBackgroundTask.ConfigClass()

Option to display the default configuration.

In [16]:

# bkgConfig

The galaxy spans approximately 100 pixels along the x-axis and 200 pixels along the y-axis. Increase the binSize to 256 for the new background estimation. Set the useApprox to False to do a direct 2D interpolation across the binned image.

In [17]:

bkgConfig.binSize = 256
bkgConfig.useApprox = False

3.2. Measure and subtract new backgound¶

Create the background subtraction task using the configuration.

In [18]:

bkgTask = SubtractBackgroundTask(config=bkgConfig)

Run the task on the copied image with the default background restored, since the background subtraction task modifies the data in-place. This task re-estimates the background and subtracts it from the vi_copy ExposureF.

In [19]:

bkgResult = bkgTask.run(vi_copy)
newbkg = bkgResult.background

Compute image and background differences as NumPy arrays, then wrap them as ImageF objects for displaying.

In [20]:

diff = vi.getImage().array - vi_copy.getImage().array
diff_image = afwImage.ImageF(diff)

diffbg = vi_bg.getImage().array - newbkg.getImage().array
diffbg_image = afwImage.ImageF(diffbg)

In [21]:

fig, axs = plt.subplots(1, 3, figsize=(15, 5), sharey=True)
axs = axs.ravel()

plt.sca(axs[0])
display1 = afwDisplay.Display(frame=fig, backend='matplotlib')
display1.scale('linear', min=757, max=762)
display1.mtv(newbkg.getImage()[x_slice, y_slice])
axs[0].scatter(xy[0], xy[1], c='r', marker='o', label='galaxy')
axs[0].set_title('New Background')
axs[0].set_xlabel('X')
axs[0].set_ylabel('Y')
axs[0].legend()

plt.sca(axs[1])
display2 = afwDisplay.Display(frame=fig, backend='matplotlib')
display2.scale('linear', 'zscale')
display2.mtv(diffbg_image[x_slice, y_slice])
axs[1].scatter(xy[0], xy[1], c='r', marker='o')
axs[1].set_title(r'$\Delta$Background (Original - New)')
axs[1].set_xlabel('X')

plt.sca(axs[2])
display3 = afwDisplay.Display(frame=fig, backend='matplotlib')
display3.scale('linear', 'zscale')
display3.mtv(vi_copy.maskedImage[x_slice, y_slice])
axs[2].scatter(xy[0], xy[1], c='r', marker='o')
axs[2].set_title(r'New Background-subtracted Image')
axs[2].set_xlabel('X')

plt.tight_layout()
plt.show()

Figure 2: Zoomed-in view of the new background estimate (left), the difference between the default and new estimates (middle), and the visit_image after subtraction of the new background (right). The red dot denotes the center of the target galaxy. The new background estimate is higher than the default around the galaxy.

4. Forced measurements with different backgrounds¶

4.1. Run the forced measurement driver¶

Configure the forced measurement driver to use base_TransformedCentroidFromCoord for centroids transformed from the reference catalog. Set other config parameters to disable shape measurement, avoid replacing other detected footprints with noise, and enable aperture correction.

In [22]:

measConfig = ForcedMeasurementDriverConfig()

In [23]:

measConfig.measurement.slots.centroid = "base_TransformedCentroidFromCoord"
measConfig.measurement.slots.shape = None
measConfig.measurement.doReplaceWithNoise = False
measConfig.doApCorr = True

Create the forced photometry driver task using the configuration.

In [24]:

driver = ForcedMeasurementDriverTask(config=measConfig)

Run the task using the input results table containing the galaxy's objectId and coordinates and the visit image with its default background subtracted.

In [25]:

meas_on_default = driver.runFromAstropy(
    results,
    vi,
    id_column_name="objectId",
    ra_column_name="coord_ra",
    dec_column_name="coord_dec",
)

lsst.forcedMeasurementDriver INFO: Measuring 1 sources in a single band using 'ForcedMeasurementTask'

lsst.forcedMeasurementDriver.measurement INFO: Performing forced measurement on 1 source

lsst.forcedMeasurementDriver INFO: Applying aperture corrections to a single band

lsst.forcedMeasurementDriver.applyApCorr INFO: Applying aperture corrections to 1 instFlux fields

lsst.forcedMeasurementDriver INFO: Finished processing for a single band; output catalog has 88 fields and 1 records

Option to display the measurement results.

In [26]:

# meas_on_default

Run the same task on the image with the new background subtracted.

In [27]:

meas_on_new = driver.runFromAstropy(
    results,
    vi_copy,
    id_column_name="objectId",
    ra_column_name="coord_ra",
    dec_column_name="coord_dec",
)

lsst.forcedMeasurementDriver INFO: Measuring 1 sources in a single band using 'ForcedMeasurementTask'

lsst.forcedMeasurementDriver.measurement INFO: Performing forced measurement on 1 source

lsst.forcedMeasurementDriver INFO: Applying aperture corrections to a single band

lsst.forcedMeasurementDriver.applyApCorr INFO: Applying aperture corrections to 1 instFlux fields

lsst.forcedMeasurementDriver INFO: Finished processing for a single band; output catalog has 88 fields and 1 records

Option to display the measurement results.

In [28]:

# meas_on_new

4.2. Compare the measurements¶

Compare the aperture photometry measurements (using the largest, 70-pixel radius) from the two images with different background subtractions applied, given the large size of the target galaxy.

Convert their instFlux to magnitudes.

In [29]:

photcalib = vi.getPhotoCalib()

flux_default = meas_on_default['base_CircularApertureFlux_70_0_instFlux'][0]
mag_default = photcalib.instFluxToMagnitude(flux_default)

flux_new = meas_on_new['base_CircularApertureFlux_70_0_instFlux'][0]
mag_new = photcalib.instFluxToMagnitude(flux_new)

In [30]:

print(f"Aperture photometry with the default background: {mag_default:.3f}")

Aperture photometry with the default background: 16.323

In [31]:

print(f"Aperture photometry with the new background: {mag_new:.3f}")

Aperture photometry with the new background: 16.331

Note: The resulting magnitudes differ slightly. Because the new background estimate is higher (i.e., more background is subtracted), the measured magnitude is correspondingly fainter.

5. Exercises for the learner¶

Experiment with other parameters in bkgConfig to see how they influence the background estimation.

In [ ]: