306.2. Candidate transient identification

306.2. Candidate transient identification¶

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

Learning objective: An overview of candidate transient (Ia Supernovae) identification in DP1.

LSST data products: DiaObject, ForcedSourceOnDiaObject, Visit, visit_image, template_coadd, difference_image

Packages: matplotlib, numpy, lsst.rsp,lsst.afw, lsst.geom, lsst.resources, lsst.utils.plotting, pyvo.dal.adhoc

Credit: Originally developed by the Rubin Community Science team with feedback from Eric Bellm. 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 how to use the light curve statistics in the DiaObject table to identify candidate transients. The specific case in this notebook is for Ia Supernova candidates.

Related tutorials: There are 200-level tutorials on difference_images as well as the DiaSource, DiaObject and ForcedSourceOnDiaObject catalogs. Use of the image cutout service is covered in the 100-level tutorials. There is another notebook in this 306 series (306.1) that presents how to make a light curve of an extragalactic transient.

1.1. Import packages¶

Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), and matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery).

From the lsst package, import modules for accessing the Table Access Protocol (TAP) service, image display functions, and the Simple Image Access (SIA) service from the LSST Science Pipelines (pipelines.lsst.io).

Functions from pyvo are imported to perform and parse results from Datalink and SODA queries. This is used to obtain image products cutouts from the SIA service.

In [1]:

import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u

from lsst.rsp import get_tap_service
from lsst.rsp.service import get_siav2_service
from lsst.rsp.utils import get_pyvo_auth

from lsst.utils.plotting import (get_multiband_plot_colors,
                                 get_multiband_plot_symbols)

import lsst.afw.display as afwDisplay
from lsst.afw.image import ExposureF
from lsst.afw.math import Warper, WarperConfig
from lsst.afw.fits import MemFileManager
import lsst.geom as geom
import lsst.resources

from pyvo.dal.adhoc import DatalinkResults, SodaQuery

1.2. Define parameters and functions¶

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

In [2]:

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

Instantiate the Simple Image Access service (SIA).

In [3]:

sia_service = get_siav2_service("dp1")
assert sia_service is not None

Define filter names, plot markers, linestyles, and colors for plotting

In [4]:

filter_names = ['u', 'g', 'r', 'i', 'z', 'y']
filter_colors = get_multiband_plot_colors()
filter_symbols = get_multiband_plot_symbols()

Set the afwDisplay package backend to matplotlib.

In [5]:

afwDisplay.setDefaultBackend('matplotlib')

2. Find candiate SN Ia¶

Type Ia Supernovae (SNIa) have homogenous light curves with similar rise and decay rates. A candidate SNIa search can be conducted in DP1 by looking for fading transients that exhibit dimming a rate consistent with SNIa: ~0.05 mag / day = ~5% flux fading per day (See Papadogiannakis et al. 2018).

Note that the 7-week duration of the commissioning survey with ComCam in DP1 (24 Oct 2024 - 11 Dec 2024) should sufficently sample that fading emission from a SNIa.

2.1. Get a sample of DiaObjects¶

Define the coordinates (RA and Dec) of the center of the Extended Chandra Deep Field-South field (ECDFS), which is the best-observed field in DP1.

In [6]:

ra = 53.13
dec = -28.10

Define the TAP query for a 2-degree cone search of the ECDFS for SNIa candidates in the DiaObject table, which contains useful light curve statistics like the slope of a linear model fit to the measured difference-image fluxes in each band (e.g. r_psfFluxLinearSlope).

Of the six bands, the ECDFS field was covered most by the r-band (See Table 3 of RTN-011). The SNIa candidate search will therefore utilize r-band light curve properties.

In [7]:

query = "SELECT ra, dec, diaObjectId, "\
        "r_psfFluxMean, nDiaSources, r_scienceFluxMean, "\
        "r_psfFluxNdata, r_psfFluxMax,r_psfFluxMin, r_psfFluxLinearSlope "\
        "FROM dp1.DiaObject "\
        "WHERE CONTAINS (POINT('ICRS', ra, dec), "\
        "CIRCLE('ICRS'," + str(ra) + ", "\
        + str(dec) + ", 2)) = 1 "
print(query)

SELECT ra, dec, diaObjectId, r_psfFluxMean, nDiaSources, r_scienceFluxMean, r_psfFluxNdata, r_psfFluxMax,r_psfFluxMin, r_psfFluxLinearSlope FROM dp1.DiaObject WHERE CONTAINS (POINT('ICRS', ra, dec), CIRCLE('ICRS',53.13, -28.1, 2)) = 1

Run the TAP search and fetch the results in table form.

In [8]:

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'
DiaObjsFull = job.fetch_result().to_table()

Job phase is COMPLETED

2.2. Review DiaObject characteristics¶

Plot the histogram distribution of DiaSources per DiaObject (i.e., the total number of $\geq5\sigma$ detections in difference images per DiaObject; at left), and the distribution of the number of r-band detections per DiaObject (at right).

Notice how the distribution is peaked at small numbers of DiaSources (and r-band detections) -- these are either artifacts from single difference images or time-variable sources that were only detected in a few difference images, and are most likely to have a faint time-variable flux.

In [9]:

fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)

ax[0].hist(DiaObjsFull['nDiaSources'], bins=50, log=True, color='gray')
ax[0].set_xlabel('nDiaSources')
ax[0].set_ylabel('Counts of DiaObjects')

ax[1].hist(DiaObjsFull['r_psfFluxNdata'], bins=50, log=True, color=filter_colors['r'])
ax[1].set_xlabel('r_psfFluxNdata')

plt.tight_layout()
plt.show()

Figure 1: At left, a histogram (grey) showing the distribution of the number of diaObjects in bins of number of diaSources (nDiaSources) across all filters. At right, a histogram (dark red) of the number of diaObjects in bins of number of r-band PSF flux detections (r_psfFluxNdata).

Plot the distribution of minimum and maximum r-band PSF flux from the difference-image detections and the distribution of the linear slope fit to r-band PSF flux measurements of the DiaObjects.

Notice how the PSF fluxes can be negative because they are measured on the difference images. Any source that is fainter in the direct image than in the template image will be detected as a negative-flux source in the difference image. Also note how these distributions are peaked at low time-variable flux components -- likely the same types of faint sources which were only detectable for a short duration (as shown above).

Warning: Take care when converting difference-image fluxes into magnitudes, because they can be negative.

In [10]:

fig, ax = plt.subplots(1, 3, figsize=(12, 3), sharey=False, sharex=False)

ax[0].hist(DiaObjsFull['r_psfFluxMin'], bins=100, log=True, color=filter_colors['r'])
ax[0].set_xlabel('r_psfFluxMin')
ax[0].set_ylabel('log(Number of DiaObjects)')

ax[1].hist(DiaObjsFull['r_psfFluxMax'], bins=100, log=True, color=filter_colors['r'])
ax[1].set_xlabel('r_psfFluxMax')

ax[2].hist(DiaObjsFull['r_psfFluxLinearSlope'], bins=100, log=True, color=filter_colors['r'])
ax[2].set_xlabel('r_psfFluxLinearSlope')

plt.tight_layout()
plt.show()

Figure 2: The DiaObject histogram distribution of the (left) minimum and (center) maximum r-band PSF flux measurements for diaObjects. (right) Histogram distribution of the linear slope of r-band PSF flux measurements.

2.3. Filter for SN Ia candidates¶

Filter the full DiaObject results (defined as DiaObjsFull) with the following conditions:

• nDiaSources > 5: Candidate has more than 5 detections (in all bands)
• r_scienceFluxMean < $10^{5.5}$ nJy: Candidate exhibits fluxes sufficiently below saturation limits
• r_psfFluxNdata > 10: Candidate has more than 10 r-band detections
• -0.07 < r_psfFluxLinearSlope/r_psfFluxMax < -0.03: r-band flux is fading at a rate between 3% - 7% per day from its max r-band flux.
• r_psfFluxMax - rpsfFluxMin > 0.1 * r_scienceFluxMean: The amplitude of the r-band flux is greater than 10% of the mean flux.

The filtered results (defined as DiaObjs) will be sorted by the r-band flux amplitude, which is defined as r_fluxamp = r_psfFluxMax - rpsfFluxMin.

In [11]:

print(str(len(DiaObjsFull))+' DiaObjects from initial query.')
DiaObjsFull['r_fluxamp'] = DiaObjsFull['r_psfFluxMax'] - DiaObjsFull['r_psfFluxMin']
DiaObjs = DiaObjsFull[DiaObjsFull['nDiaSources'] > 5]
DiaObjs = DiaObjs[DiaObjs['r_scienceFluxMean'] < 10**5.5]
DiaObjs = DiaObjs[DiaObjs['r_psfFluxNdata'] > 10]
DiaObjs = DiaObjs[(DiaObjs['r_psfFluxLinearSlope']/DiaObjs['r_psfFluxMax'] < -0.03)
                  & (DiaObjs['r_psfFluxLinearSlope']/DiaObjs['r_psfFluxMax'] > -0.07)]
DiaObjs = DiaObjs[DiaObjs['r_fluxamp'] > 0.1*DiaObjs['r_scienceFluxMean']]
DiaObjs.sort('r_fluxamp', reverse=True)
print(str(len(DiaObjs))+' SNIa candidates.')

210535 DiaObjects from initial query.
13 SNIa candidates.

Uncomment to view table.

In [12]:

# DiaObjs

3. Plot light curves for the SN Ia candidates¶

Obtain and plot the multi-band light curves of the top 10 SNIa candidates with the highest r-band flux amplitudes using their DiaObjectIds and querying the ForcedSourceOnDiaObject table with a JOIN on the Visit table.

The ForcedSourceOnDiaObject table to obtain forced PSF photometry for all difference images instead of the photometry from the DiaSource table which only measure difference-image sources with a signal-to-noise ratio (SNR) $\geq$ 5 detections. See tutorial 306.1. for more details on retrieving photometry from the ForcedSourceOnDiaObject table.

Note that the conditions used to filter the SNIa candidates in the previous section (Sec. 2.1.2.) were using statistics computed from the DiaSource table (i.e. on SNR $\geq5$ detections) and not the ForcedSourceOnDiaObject table.

In [13]:

cand_list = "(" + ", ".join(str(value) for value in DiaObjs['diaObjectId']) + ")"

query = "SELECT fsodo.coord_ra, fsodo.coord_dec, "\
        "fsodo.diaObjectId, fsodo.visit, fsodo.band, "\
        "fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, "\
        "vis.expMidptMJD "\
        "FROM dp1.ForcedSourceOnDiaObject as fsodo "\
        "JOIN dp1.Visit as vis ON vis.visit = fsodo.visit "\
        "WHERE diaObjectId IN {}" \
        "ORDER BY diaObjectId ASC".format(cand_list)
print(query)

SELECT fsodo.coord_ra, fsodo.coord_dec, fsodo.diaObjectId, fsodo.visit, fsodo.band, fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, vis.expMidptMJD FROM dp1.ForcedSourceOnDiaObject as fsodo JOIN dp1.Visit as vis ON vis.visit = fsodo.visit WHERE diaObjectId IN (611255141361779352, 611255759837069401, 611256447031836758, 609788805167185927, 611255072642302059, 611255141361779062, 609788942606139423, 611255278800732178, 609789423642476567, 611255759837069429, 611255072642301958, 611255003922825270, 611255003922825298)ORDER BY diaObjectId ASC

In [14]:

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

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

Create the multi-panel plot of light curves.

Warning: the following cell products a pink warning message that masks might be ignored. Columns with null forced source photometry are masked, but if they are masked they are still not plotted, so this is ok.

In [16]:

fig, ax = plt.subplots(5, 2, figsize=(10, 12), sharey=False, sharex=False)
x = 0
for i in range(5):
    for j in range(2):
        for f, filt in enumerate(filter_names):
            dx = np.where(FrcdSrc['diaObjectId'] == DiaObjs['diaObjectId'][x])
            lc = FrcdSrc[dx]
            fx = np.where(lc['band'] == filt)[0]
            ax[i, j].plot(lc['expMidptMJD'][fx], lc['psfDiffFlux'][fx],
                          filter_symbols[filt], ms=10, mew=0, alpha=0.5,
                          color=filter_colors[filt])
            del fx
        ax[i, j].set_title(DiaObjs['diaObjectId'][x])
        ax[i, j].set_ylim(np.median(lc['psfDiffFlux'])-3*np.std(lc['psfDiffFlux']),
                          np.median(lc['psfDiffFlux'])+3*np.std(lc['psfDiffFlux']))
        if i == 4:
            ax[i, j].xaxis.set_label_text('MJD (days)')
        if j == 0:
            ax[i, j].yaxis.set_label_text('DiffIm Flux (nJy)')
        x += 1
plt.tight_layout()
plt.show()

/opt/lsst/software/stack/conda/envs/lsst-scipipe-10.0.0/lib/python3.12/site-packages/numpy/_core/fromnumeric.py:868: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedColumn.
  a.partition(kth, axis=axis, kind=kind, order=order)

Figure 3: Multi-band forced-PSF difference-image photometry of 10 SNIa candidates from the ForcedSourceOnDiaObject table. Legend: $u$, light blue circle; $g$, green triangle; $r$, dark red inverted triangle; $i$, dark blue square; $z$, pink star; $y$, brown pentagon.

4. Use an alert broker to investigate a candidate SN Ia¶

The candidate in the top right plot (611255759837069401) exhibits a fairly compelling Ia-like light curve. The Zwicky Transient Facility has coverage of this same field and thus a search for ZTF alerts via the Alert Brokers could yield valuable insight.

A simple coordinate search in the brokers (e.g. in ANTARES) reveals that a ZTF alert was generated for this source in Sept 2024 (MJD 60565). Here are results from several brokers for this transient: ALERCE, ANTARES, Fink, Lasair. Note that brokers will also be streaming Rubin/LSST alerts.

If the transient peaked on MJD 60565 when it was detected by ZTF, the observations from DP1 would cover 60 - 90 days after the peak and would be consistent with the fading timescale of a SNIa.

5. Create a cutout triplet¶

A "cutout triplet" for a DIA source refers to the set of three stamps made from the science, template, and difference images. They are typically displayed side-by-side and serve as a visualization of the detection.

5.1. Define a function to get cutouts¶

This section uses the SIA service to find images and the image cutout service to create a "cutout triplet". Both services are demonstrated in 100-level tutorials.

Define a function to get an image cutout from the cutout service.

In [17]:

def get_cutout(dl_result, spherePoint, session, fov):
    sq = SodaQuery.from_resource(dl_result,
                                 dl_result.get_adhocservice_by_id("cutout-sync"),
                                 session=session)
    sphereRadius = fov * u.deg
    sq.circle = (spherePoint.getRa().asDegrees() * u.deg,
                 spherePoint.getDec().asDegrees() * u.deg,
                 sphereRadius)
    cutout_bytes = sq.execute_stream().read()
    sq.raise_if_error()
    mem = MemFileManager(len(cutout_bytes))
    mem.setData(cutout_bytes, len(cutout_bytes))
    return ExposureF(mem)

Use the DiaObjectID of the SNIa candidate from above, obtain its coordinates, and use the lsst.geom package to convert the coordinates to a SpherePoint format.

In [18]:

SNcandID = 611255759837069401
ra = DiaObjs[DiaObjs['diaObjectId'] == SNcandID]['ra'][0]
dec = DiaObjs[DiaObjs['diaObjectId'] == SNcandID]['dec'][0]
spherePoint = lsst.geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)

5.2. Query the SIA service for images¶

Define circle, the position (RA and Dec) and radius (in degrees) to search for science (visit_image), reference (template_coadd), and difference (difference_image) images with SIA.

In [19]:

circle = (ra, dec, 0.001)

Use the SIA service to search for calibration level 2 images (visit_images) that overlap the search circle. Filter the results for r-band visit_images and sort them by observation time.

In [20]:

results_sci = sia_service.search(pos=circle, calib_level=2)

lvl2_table = results_sci.to_table()

tx = np.where((lvl2_table['dataproduct_subtype'] == 'lsst.visit_image')
              & (lvl2_table['lsst_band'] == 'r'))[0]
scitab = lvl2_table[tx]
scitab.sort('t_max')

Use the SIA service to search for calibration level 3 images (template_coadd images and difference_images) that overlap the search circle. Filter the results for r-band template_coadd images and difference_images and sort the latter by observation time. No sorting is performed for the reference template images because there is only one image.

In [21]:

results = sia_service.search(pos=circle, calib_level=3)

lvl3_table = results.to_table()

tx = np.where((lvl3_table['dataproduct_subtype'] == 'lsst.template_coadd')
              & (lvl3_table['lsst_band'] == 'r'))[0]
reftab = lvl3_table[tx]

tx = np.where((lvl3_table['dataproduct_subtype'] == 'lsst.difference_image')
              & (lvl3_table['lsst_band'] == 'r'))[0]
difftab = lvl3_table[tx]
difftab.sort('t_max')

Print the date of the earliest r-band image, for which image cutouts will be retrieved.

In [22]:

print(scitab['t_max'][0], difftab['t_max'][0])

60623.25907695602 60623.25907695602

Retrieve the datalink VOTable documents using the access URLs from the first entries in each of the table results. The image products being retrieved for the science and difference images correspond to the earliest r-band image products taken.

In [23]:

dl_result_sci = DatalinkResults.from_result_url(scitab['access_url'][0], session=get_pyvo_auth())

dl_result_ref = DatalinkResults.from_result_url(reftab['access_url'][0], session=get_pyvo_auth())

dl_result_diff = DatalinkResults.from_result_url(difftab['access_url'][0], session=get_pyvo_auth())

5.3. Get cutouts from the images¶

Retrive the cutouts with a radius of fov (in degrees).

In [24]:

fov = 0.003
sci = get_cutout(dl_result_sci, spherePoint, get_pyvo_auth(), fov)
ref = get_cutout(dl_result_ref, spherePoint, get_pyvo_auth(), fov)
diff = get_cutout(dl_result_diff, spherePoint, get_pyvo_auth(), fov)

5.4. Warp the template¶

The template_coadd needs to be warped to match the orientation of the science and difference cutouts. Perform the image warping of the template_coadd cutout using the Warper from lsst.afw.math. The warped reference image is defined as warped_ref.

In [25]:

warper_config = WarperConfig()
warper = Warper.fromConfig(warper_config)

sci_wcs = sci.getWcs()
sci_bbox = sci.getBBox()

warped_ref = warper.warpExposure(sci_wcs, ref, destBBox=sci_bbox)

5.5. Display the cutout triplet¶

Display the science, reference, and difference cutouts for the SNIa candidate associated with the DiaObject 611255759837069401.

In [26]:

fig, ax = plt.subplots(1, 3, figsize=(8, 3))
plt.sca(ax[0])
display1 = afwDisplay.Display(frame=fig)
display1.scale('linear', 'zscale')
display1.mtv(sci.image)
plt.title('science')
plt.sca(ax[1])
display2 = afwDisplay.Display(frame=fig)
display2.scale('linear', 'zscale')
display2.mtv(warped_ref.image)
plt.title('tempate')
plt.sca(ax[2])
display3 = afwDisplay.Display(frame=fig)
display3.scale('linear', 'zscale')
display3.mtv(diff.image)
ax[0].set_axis_off()
ax[1].set_axis_off()
ax[2].set_axis_off()
plt.title('difference')
plt.tight_layout()
fig.suptitle('DiaObject 611255759837069401, MJD='+str(np.round(scitab['t_max'][0], 2)))
plt.show()

Figure 4: The cutout triplet for the first detection of this DiaObject, a candidate SN Ia.

In [ ]: