306.1. Transient light curves#

306_1_Transient_light_curves

306.1. Transient light curves¶

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

Learning objective: An overview of plotting extragalactic transient light curves.

LSST data products: ForcedSourceOnDiaObject, DiaObject, Visit

Packages: matplotlib, numpy, lsst.rsp, lsst.utils.plotting

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 obtain and plot a difference-image light curve of an extragalactic transient.

Recommendations for extragalactic transient light curves.

For extragalactic transients, it is recommended to create light curves using forced difference-image photometry. Using the difference-image photometry ensures that any underlying host galaxy light does not contaminate the flux measurement, as can happen if science-image (or direct- or visit-image) photometry is used. Using the forced photometry ensures that a data point exists for every observation that overlapped the transient's coordinates, not just those in which it was detected at $>5\sigma$ in the difference image.

The forced difference-image photometry is stored in the ForcedSourceOnDiaObject table, and the relevant columns are:

diaObjectId: Unique DiaObject identifier
psfDiffFlux: Forced PSF flux measured on difference images at the diaObject position (nJy)
band: Filter associated with flux measurement
visit: Identifier of the visit for which the forced photometry was measured

Observation dates (MJDs) are stored in the Visit table, in column expMidptMJD (exposure midpoint MJD), and are obtained with a table join.

Difference-image fluxes can be negative and thus should not be converted to magnitudes, as without care this can lead to missing data points (especially observations with non-detections, which can be scientifically very useful for transients).

Related tutorials: There are 200-level tutorials on difference_images as well as the DiaSource, DiaObject and ForcedSourceOnDiaObject catalogs.

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 and image display functions from the LSST Science Pipelines (pipelines.lsst.io).

In [1]:

import matplotlib.pyplot as plt
import numpy as np
from lsst.rsp import RSPDiscovery
from lsst.utils.plotting import (get_multiband_plot_colors,
                                 get_multiband_plot_symbols,)

1.2. Define parameters and functions¶

Create an instance of the TAP service.

In [2]:

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

Define filter names, plot markers, and colors for plotting

In [3]:

filter_colors = get_multiband_plot_colors()
filter_names = filter_colors.keys()
filter_symbols = get_multiband_plot_symbols()

2. Retrieve the light curve data¶

Light curve data is retrieved from the ForcedSourceOnDiaObject by the associated diaObjectId, not by coordinates.

If the diaObjectId of the time-domain object for which the light curve is desired is already known, skip to Section 2.2.

2.1. Spatial query for the diaObject¶

Spatial queries require a table join: Spatial queries for time-domain objects are done on the DiaObject table, not on the tables containing the measured fluxes at the same location in every image. Spatial queries on ForcedSourceOnDiaObject (or DiaSource) are inefficient and for this reason, coordinates are not included in the ForcedSourceOnDiaObject table.

Cross-match radius: Although difference-image detections are associated within a 1" radius, and typically that is a fine search radius to use for point-source cross-match, a 2" radius is used for this example. It mimics a situation in which, e.g., the transient was detected in a lower-resolution survey, and demonstrates how to deal with multiple potential cross-matches. For scenarios in which cross-match for a set of coordinates is desired, instead of a single coordinate, see the tutorials on cross-match to user-uploaded tables.

Define the coordinates of a known transient and the search radius. Query the DiaObject table and retrieve the diaObjectId, the number of detections (nDiaSources), and the offset distance between all matches and the search coordinates. Retrieve the results as an Astropy table.

In [4]:

ra_targ = 151.592579
dec_targ = 1.158333
search_radius = 2.0/3600.0
query = """SELECT diaObjectId, nDiaSources,
        DISTANCE(POINT('ICRS', ra, dec), POINT('ICRS', {}, {}))
        FROM dp2.DiaObject
        WHERE CONTAINS(POINT('ICRS', ra, dec), CIRCLE('ICRS', {}, {}, {})) = 1
        """.format(ra_targ, dec_targ, ra_targ, dec_targ, search_radius)
job = tap_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'
diao_table = job.fetch_result().to_table()

Job phase is COMPLETED

Show the table.

In [5]:

diao_table

Out[5]:

Table length=2

diaObjectId	nDiaSources	DISTANCE3
int64	int64	float64
786352299365629972	201	2.0458469231481843e-07
786352299365631545	1	0.0004390324054815023

Although there are two diaObjects within the search table, the one with many more detections and a much smaller offset (distance) is clearly the object of interest.

2.2. Query for forced photometry¶

Create a query to retrieve the flux, band, and MJD of the diaObject using a table join on ForcedSourceOnDiaObject and Visit. Retrieve the results as an Astropy table.

In [6]:

use_diaObjectId = 786352299365629972
query = """SELECT fsodo.psfDiffFlux, fsodo.psfDiffFluxErr,
        fsodo.band, v.expMidptMJD
        FROM dp2.ForcedSourceOnDiaObject AS fsodo
        JOIN dp2.Visit AS v ON v.visit = fsodo.visit
        WHERE fsodo.diaObjectId = {}""".format(use_diaObjectId)
job = tap_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'
flux_table = job.fetch_result().to_table()
print(len(flux_table))

Job phase is COMPLETED
295

Option to view the table.

In [7]:

# flux_table

3. Plot the light curve¶

Plot the difference-image forced fluxes as open symbols with error bars.

In [8]:

fig = plt.figure(figsize=(8, 4))
for f, filt in enumerate(filter_names):
    fx = np.where(flux_table['band'] == filt)[0]
    plt.errorbar(flux_table['expMidptMJD'][fx], flux_table['psfDiffFlux'][fx],
                 yerr=flux_table['psfDiffFluxErr'][fx],
                 fmt=filter_symbols[filt], ms=10, mew=2, mec=filter_colors[filt],
                 ecolor=filter_colors[filt], alpha=0.5, color='none', label=filt)
plt.xlabel('Modified Julian Date', fontsize=12)
plt.ylabel('Difference-Image Forced Flux (nJy)', fontsize=12)
plt.legend(loc='upper right', ncol=2, fontsize=12)
plt.show()

No description has been provided for this image

Figure 1: Forced PSF difference-image photometry of an extragalactic transient from the ForcedSourceOnDiaObject table. Note that some difference-image fluxes are negative. This indicates that the template images that were subtracted from the science images likely contained some transient light. This is not surprising given that the templates were built out of the same images that were used in difference image analysis.

Zoom in to focus on the earlier days while the event was rising.

In [9]:

fig = plt.figure(figsize=(8, 4))
for f, filt in enumerate(filter_names):
    fx = np.where(flux_table['band'] == filt)[0]
    plt.errorbar(flux_table['expMidptMJD'][fx], flux_table['psfDiffFlux'][fx],
                 yerr=flux_table['psfDiffFluxErr'][fx],
                 fmt=filter_symbols[filt], ms=10, mew=2, mec=filter_colors[filt],
                 ecolor=filter_colors[filt], alpha=0.5, color='none', label=filt)
plt.xlabel('Modified Julian Date', fontsize=12)
plt.ylabel('Difference-Image Forced Flux (nJy)', fontsize=12)
plt.xlim([60809, 60829])
plt.legend(loc='upper left', ncol=2, fontsize=12)
plt.show()

Figure 2: Same as Figure 1, but focusing on the rising epochs.

3.1. Take care with magnitudes¶

Do not convert magnitudes to difference-image fluxes, as any negative fluxes will be lost. Negative difference-image fluxes are not erroneous -- they simply mean the time-varying object was brighter in the template than in the science image.

Re-plot the light curve in Figure 2, but with magnitudes, to show the error message and that some data points go missing. There will be a pink "RuntimeWarning" to represent negative values passed to np.log10.

In [10]:

flux_table['psfDiffMag'] = -2.5 * np.log10(flux_table['psfDiffFlux']) + 31.4
fig = plt.figure(figsize=(8, 4))
for f, filt in enumerate(filter_names):
    fx = np.where(flux_table['band'] == filt)[0]
    plt.plot(flux_table['expMidptMJD'][fx], flux_table['psfDiffMag'][fx],
             filter_symbols[filt], ms=10, mew=2, mec=filter_colors[filt],
             alpha=0.5, color='none', label=filt)
plt.xlabel('Modified Julian Date', fontsize=12)
plt.ylabel('Difference-Image Forced Magnitude', fontsize=12)
plt.xlim([60809, 60829])
plt.ylim([30, 18.0])
plt.legend(loc='upper left', ncol=2, fontsize=12)
plt.show()

/tmp/ipykernel_2913/2820271453.py:1: RuntimeWarning: invalid value encountered in log10
  flux_table['psfDiffMag'] = -2.5 * np.log10(flux_table['psfDiffFlux']) + 31.4

Figure 3: Similar to Figure 2, but only showing positive difference-image forced fluxes that have been converted to apparent magnitude. Data points for which the time-varying object is negative in the difference-image are missing from this plot.

3.2. Why forced flux is recommended¶

To demonstrate why forced difference-image photometry is recommended for transient light curves, retrieve the detection photometry from the DiaSource table for the same object and plot the light curve, for comparison.

In [11]:

use_diaObjectId = 786352299365629972
query = """SELECT psfFlux, psfFluxErr, band, midpointMjdTai
        FROM dp2.DiaSource AS dias
        WHERE diaObjectId = {}""".format(use_diaObjectId)
job = tap_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'
dias_table = job.fetch_result().to_table()
print(len(dias_table))

Job phase is COMPLETED
201

Plot the forced and detection difference-image flux light curves side by side.

In [12]:

fig = plt.figure(figsize=(8, 3))
for f, filt in enumerate(filter_names):
    fx = np.where(flux_table['band'] == filt)[0]
    plt.errorbar(flux_table['expMidptMJD'][fx], flux_table['psfDiffFlux'][fx],
                 yerr=flux_table['psfDiffFluxErr'][fx],
                 fmt=filter_symbols[filt], ms=10, mew=2, mec=filter_colors[filt],
                 ecolor=filter_colors[filt], alpha=0.5, color='none', label=filt)
plt.xlabel('Modified Julian Date', fontsize=12)
plt.ylabel('Forced Flux (nJy)', fontsize=12)
plt.xlim([60809, 60829])
plt.legend(loc='upper left', ncol=2, fontsize=12)
plt.show()

fig = plt.figure(figsize=(8, 3))
for f, filt in enumerate(filter_names):
    fx = np.where(dias_table['band'] == filt)[0]
    plt.errorbar(dias_table['midpointMjdTai'][fx], dias_table['psfFlux'][fx],
                 yerr=dias_table['psfFluxErr'][fx],
                 fmt=filter_symbols[filt], ms=10, mew=2, mec=filter_colors[filt],
                 ecolor=filter_colors[filt], alpha=0.5, color='none', label=filt)
plt.xlabel('Modified Julian Date', fontsize=12)
plt.ylabel('Detection Flux (nJy)', fontsize=12)
plt.xlim([60809, 60829])
plt.legend(loc='upper left', ncol=2, fontsize=12)
plt.show()

Figure 4: Notice that data points for some epochs, especially those in which the difference-image flux was close to 0, are missing from the bottom plot. For forced, difference-image photometry, a flux of zero is still an interesting and useful data point. But these points are missed if only the detections in the DiaSource table are used instead of all forced measurements in ForcedSourceOnDiaObject.

In [ ]: