202.3. Difference images

202.3. Difference images¶

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

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

LSST data products: difference_image

Packages: lsst.daf.butler, lsst.rsp, lsst.afw.display

Credit: Originally developed by the Rubin Community Science team with feedback from Eric Bellm. This notebook also highlights a transient from LSSTComCam first identified by Dan Taranu. 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 difference_image is the pixel-by-pixel difference of a visit_image and a deep, coadded template image of the same position in the sky from a single detector (a single CCD) in a single band.

Before subtracting the visit and coadd template images, they are warped to the same pixel grid, scaled to the same photometric response, matched to the same point spread function (PSF) shape, and corrected for differential chromatic refraction.

• butler dataset type: difference_image
• number of images: 15972

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. The 200-level tutorials on the photometric and astrometric calibration, and the point spread function, describe these aspects of the images in more detail.

1.1. Import packages¶

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

In [1]:

from lsst.daf.butler import Butler, Timespan
from lsst.rsp import get_tap_service
import lsst.geom
import lsst.afw.display as afwDisplay
from astropy.time import Time
import numpy as np

1.2. Define parameters and functions¶

Instantiate the Butler.

In [2]:

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

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

In [3]:

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

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

In [4]:

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 difference_image dataset type are the filter, instrument, observation date, visit, and detector, and that the returned image will be of type ExposureF. The LSST Science Pipelines tools for image manipulation and visualization work best with the ExposureF format. This format includes pixel data and metadata.

In [5]:

butler.get_dataset_type('difference_image')

Out[5]:

DatasetType('difference_image', {band, instrument, day_obs, detector, physical_filter, visit}, ExposureF)

Show that only instrument, detector, and visit are required to uniquely identify a difference_image.

In [6]:

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

Out[6]:

{instrument, detector, visit}

2.1.1. Demo query¶

Query for difference_images that overlap coordinates RA, Dec in the ECDFS field with a known transient, were obtained with the r-band filter, and within a timespan near the start of the LSSTComCam campaign.

In [7]:

ra = 53.125
dec = -27.740
band = "r"
time1 = Time("2024-11-09T00:00:00.0", format="isot", scale="tai")
time2 = Time(60624.0, format="mjd", scale="tai")
timespan = Timespan(time1, time2)
del time1, time2

In [8]:

dataset_refs = butler.query_datasets("difference_image",
                                     where="band.name = :band AND \
                                     visit.timespan OVERLAPS :timespan AND \
                                     visit_detector_region.region OVERLAPS POINT(:ra, :dec)",
                                     bind={"band": band, "timespan": timespan,
                                           "ra": ra, "dec": dec},
                                     order_by=["visit.timespan.begin"])
print(len(dataset_refs))

16

In [9]:

del band, timespan

2.1.2. Get a difference image¶

The difference_image can be retrieved from the Butler using the (a) the dataset_refs; (b) the dataId; or (c) the visit and detector numbers.

(a) Use the dataset reference¶

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

In [10]:

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

See that it is of type ExposureF.

In [11]:

difference_image

Out[11]:

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

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

In [12]:

afw_display.mtv(difference_image)
afw_display.setMaskTransparency(100)

Clean up.

In [13]:

del ref, difference_image

(b) Use the dataId¶

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

In [14]:

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

Out[14]:

{instrument: 'LSSTComCam', detector: 5, visit: 2024110800247, band: 'r', day_obs: 20241108, physical_filter: 'r_03'}

Use the defined dataId to retrieve the corresponding difference_image.

In [15]:

difference_image = butler.get('difference_image', dataId=use_dataId)

Option to display the difference_image.

In [16]:

# afw_display.mtv(difference_image)

Clean up.

In [17]:

del ref, use_dataId, difference_image

(c) Use the visitId and detector¶

For the third dataset reference returned by the query, get the dataId and extract the visit and detector numbers.

In [18]:

ref = dataset_refs[2]
dataId = ref.dataId
use_visit = dataId.get('visit')
use_detector = dataId.get('detector')
print(use_visit, use_detector)

2024110800251 5

Use the visit and detector to get the corresponding difference_image.

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 visit and detector.

In [19]:

difference_image = butler.get('difference_image', visit=use_visit,
                              detector=use_detector)

Option to display the difference_image.

In [20]:

# afw_display.mtv(difference_image)

Clean up.

In [21]:

del ref, dataId, use_visit, use_detector, difference_image

2.2. SIA (Simple Image Access)¶

The Simple Image Access (SIA) service 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 difference_images from a SIA query, set:

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

See Section 2.1.2. for how to retrieve images via the Butler using the lsst_visit and lsst_detector 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.

2.3.1. ObsCore table¶

Image metadata is stored in the ObsCore table.

To return only visit_images from a TAP query, set:

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

See Section 2.1.2. for how to retrieve images via the butler using the lsst_visit and lsst_detector columns in the TAP results table.

3. Pixel data¶

The difference_images have three planes of pixel data: image, variance, and mask.

Retrieve the first difference_image from the Butler query in Section 2.1.1.

In [22]:

ref = dataset_refs[0]
print(ref.dataId)

{instrument: 'LSSTComCam', detector: 2, visit: 2024110800246, band: 'r', day_obs: 20241108, physical_filter: 'r_03'}

In [23]:

difference_image = butler.get(ref)

Display the difference_image in frame 1, and set the mask tranparency to 50 (semi-transparent; to overlay the mask).

In [24]:

afw_display.mtv(difference_image)
afw_display.setMaskTransparency(50)

3.1. Image plane¶

Sky pixel data in units of nJy (nanoJanskys).

Since a template image has been subtracted, sources that are fainter in the visit image than they were in the template will have negative flux values.

Extract the image data from the difference_image, and extract the pixel values as an array.

In [25]:

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

Option to show the image and the image_array.

In [26]:

# image

In [27]:

# image_array

Print statistics on the image_array.

In [28]:

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

-74727.85
-0.025696142
13028.367

Display the image in Firefly frame 2.

In [29]:

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

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

3.2. Variance plane¶

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

Values are always positive.

Extract the variance and variance_array.

In [30]:

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

Print statistics on the variance_array.

In [31]:

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

164.21577
309.26575
78008.25

Display the variance.

In [32]:

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.

Extract the mask.

In [33]:

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

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

Print the list of mask keys and color.

In [34]:

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

BAD: red
CLIPPED: blue
CR: magenta
CROSSTALK: cyan
DETECTED: blue
DETECTED_NEGATIVE: cyan
EDGE: yellow
INEXACT_PSF: magenta
INJECTED: yellow
INJECTED_TEMPLATE: orange
INTRP: green
ITL_DIP: red
NOT_DEBLENDED: green
NO_DATA: orange
REJECTED: blue
SAT: green
SAT_TEMPLATE: cyan
SENSOR_EDGE: magenta
STREAK: green
SUSPECT: yellow
UNMASKEDNAN: yellow
VIGNETTED: red

Display the mask in Firefly frame 4.

In [35]:

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

afw_display.setMaskTransparency(0)

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

In [37]:

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

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

In [38]:

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

Clean up.

In [39]:

del image, image_array, variance, variance_array, mask
del ref, dataset_refs

4. Metadata and attributes¶

4.1. Header¶

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

Get the difference_image metadata.

In [40]:

metadata = difference_image.getMetadata()

Option to display the long list of metadata.

In [41]:

# metadata

Convert the metadata to a python dictionary.

In [42]:

md_dict = metadata.toDict()

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

In [43]:

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

LSST ISR UNITS
LSST ISR OVERSCAN SERIAL UNITS
LSST ISR READNOISE UNITS
LSST ISR OVERSCAN PARALLEL UNITS
BUNIT

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

In [44]:

md_dict['BUNIT']

Out[44]:

'nJy'

Clean up.

In [45]:

del metadata, md_dict

4.2. Visit information¶

Information about the observation associated with this difference_image, such as pointing, airmass, and weather.

Get the visit info.

In [46]:

info = difference_image.getInfo().getVisitInfo()

In [47]:

print(info.getBoresightAirmass())
print(info.getBoresightParAngle())
print(info.getDate())
print(info.getExposureTime())

1.03814024464709
1.78496 rad
DateTime("2024-11-09T06:13:26.028993809", TAI)
30.0

Get weather information.

In [48]:

weather = info.getWeather()

In [49]:

print(weather.getAirPressure())
print(weather.getAirTemperature())
print(weather.getHumidity())

74135.0
14.375
17.0499992370605

In [50]:

del info, weather

4.3. World Coordinate System (WCS)¶

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

Get the WCS and show it.

In [51]:

wcs = difference_image.getWcs()
wcs

Out[51]:

Non-standard SkyWcs (Frames: PIXELS, NORMEDPIXELS, INTERMEDIATE0, IWC, SKY):
Sky Origin: (53.3274246637, -28.0723466054)
Pixel Origin: (-3193.23, 6504.86)
Pixel Scale: 0.200609 arcsec/pixel

Convert pixel coordinates to sky coordinates, and vice versa.

In [52]:

coord = lsst.geom.SpherePoint(53.125*lsst.geom.degrees, -27.740*lsst.geom.degrees)
xy = wcs.skyToPixel(coord)

In [53]:

print(xy)

(1905.3, 2040.7)

In [54]:

xy2 = lsst.geom.Point2D(1905.3, 2040.7)
coord2 = wcs.pixelToSky(xy2)

In [55]:

print(coord2)

(53.1249985150, -27.7400011325)

In [56]:

del xy, xy2, coord, coord2

4.4. Bounding box (corners)¶

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

Get the bounding box.

In [57]:

bbox = difference_image.getBBox()

In [58]:

bbox

Out[58]:

Box2I(corner=Point2I(0, 0), dimensions=Extent2I(4072, 4000))

Print the difference_image corners in pixels and sky coordinates.

In [59]:

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

(0, 0) (52.9730157867, -27.8178541310)
(0, 4000) (53.2184235349, -27.8677139888)
(4072, 4000) (53.2753713343, -27.6467395750)
(4072, 0) (53.0303569372, -27.5968876224)

Print the area of the bounding box, first in pixels squared and then in arcminutes squared. The pixel scale is calculated at the center of the bounding box, as it is recommended not to call the wcs.getPixelScale() method without an argument.

In [60]:

print(bbox.area, 'square pixels')
scale = wcs.getPixelScale(lsst.geom.Point2D(
    bbox.centerX, bbox.centerY)).asArcminutes()

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

16288000 square pixels
181.68 square arcminutes

Check if pixel values are within the bounding box.

In [61]:

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

True
False

Check if sky coordinates are within the bounding box.

In [62]:

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

False

Clean up.

In [63]:

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

4.5. Point spread function (PSF)¶

The PSF refers to the "shape" of a point source in the image due due to diffraction in the optical system and sky conditions. Note that the PSF of the difference_image should be the exact same as that of the visit_image.

Get the PSF.

In [64]:

psf = difference_image.getPsf()

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

In [65]:

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

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

In [66]:

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

PSF sigma: 1.89762088083317, PSF FWHM: 4.468555688055235
PSF shape (ixx, iyy, ixy):  3.363654908117207 3.859563794899629 0.12365968451024756

Clean up.

In [67]:

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

5. Detected DIA sources¶

Use the TAP service to retrieve the DiaSources detected in the difference_image.

Get the visit and detector numbers for the difference_image.

In [68]:

visit = difference_image.getInfo().getVisitInfo().id
detector = difference_image.getInfo().getDetector().getId()
print(visit, detector)

2024110800246 2

Define the TAP query from the DiaSource table to retrieve the x and y pixel coordinate positions of the sources in the visit and detector of the difference image.

In [69]:

query = "SELECT diaSourceId, x, y " \
        "FROM dp1.DiaSource " \
        "WHERE visit = " + str(visit) + " " \
        "AND detector = " + str(detector)

Perform the TAP query and define the result in table form.

In [70]:

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()
assert job.phase == 'COMPLETED'

diasource = job.fetch_result().to_table()

Job phase is COMPLETED

Option to show the pixel positions of the DiaSources.

In [71]:

# diasource

Redisplay the difference_image in frame 1 with no mask.

In [72]:

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

Overplot diasources in the difference_image with orange circles.

In [73]:

with afw_display.Buffering():
    for i in np.arange(len(diasource)):
        afw_display.dot('o', diasource[i]['x'], diasource[i]['y'], size=20, ctype='orange')

Clean up.

In [74]:

del difference_image
del visit, detector, diasource

In [ ]: