306.2. Cross match queries#

306_2_Cross_match_queries

306.2. Cross match queries¶

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-10-16
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: Cross match to a list of extragalactic transients and retrieve light curves.

LSST data products: ForcedSourceOnDiaObject, DiaObject, Visit

Packages: lsst.rsp

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 tutorial demonstrates how to make efficient queries for the common question of "does LSST have difference-image light curves at these sky coordinates?"

For light curves of extragalactic transients, it is recommended to use forced flux measured on the difference image. These measurements are stored in the ForcedSourceOnDiaObject table.

Do not query source tables by coordinate. Queries by coordinate on the ForcedSourceOnDiaObject table are very inefficient, because this table contains forced measurements on every image for every row of the DiaObject table. It is more efficient to do the spatial cross-match with DiaObject once and then retrieve forced fluxes for successful matches based on their IDs. To enforce efficient queries, the ForcedSourceOnDiaObject table does not have coordinates - but the same advice applies to the ForcedSource and Object tables.

Avoid putting asynchronous queries into for loops. Although it is possible loop over a set of RA and Dec, and make a spatial cross-match query for each coordinate, this should be avoided. TAP queries have a rate limit to prevent the shared system from being overloaded by a single user with many queries executing simultaneously.

The correct method is to upload a table of coordinates and cross-match using a TAP table join. This notebook demonstrates how to upload a table of coordinates and use the TAP service to cross-match to the DiaObject table and then retrieve light curves for successful cross-matches.

Related tutorials: The 300-level tutorial on transient light curves. The 200-level tutorials on the DiaObject and ForcedSourceOnDiaObject catalogs. The 100-level tutorials on catalog access with TAP and ADQL query statements.

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 matplotlib.colors as mcolors
import skyproj
import numpy as np
from astropy.table import Table
from lsst.rsp import RSPDiscovery
from lsst.utils.plotting import (get_multiband_plot_colors,
                                 get_multiband_plot_symbols,
                                 set_rubin_plotstyle)

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

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

Define a function to truncate a matplotlib colormap. This can be used to cut off, e.g., the lightest or darkest colors of a sequential colormap, to help with data visibility on a light or dark background. For example, setting minval to 0.2 will cut off the lower 20% of the colormap; setting maxval to 0.7 will cut off the upper 30% of the colormap. See Section 3.2 for an example application to the "Blues" sequential colormap.

In [4]:

def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=256):
    """
    Customize a matplotlib colormap.

    Parameters
    ----------
    cmap: the matplotlib colormap
    minval: desired lower bound of colormap
    maxval: desired upper bound of colormap
    n: number of intervals (colormap resolution)

    Returns
    -------
    new_cmap: the truncated colormap
    """

    new_cmap = mcolors.LinearSegmentedColormap.from_list(
        'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval),
        cmap(np.linspace(minval, maxval, n)))
    return new_cmap

2. Load the coordinate list¶

A list of classified supernovae that were reported to the Transient Name Server (TNS) during the time when observations for Data Preview 2 were being obtained (April 2025 to January 2026) has been generated for this demonstration.

The TNS search query was executed by:

removing the default "reported within the last" constraint
checking the box "Limit to classified supernovae"
setting the Discovery Date Range to Apr 1 2025 to Jan 31 2026
setting the Discovery Mag Range to 18 to 23 mag
under "Advanced search", setting the Declination range to -90 to +30 deg

From the resulting TNS results table, four columns were stored in a file for this tutorial: the coordinates, name, and classification type. Read the file into an Astropy table.

In [5]:

path = '/rubin/cst_repos/tutorial-notebooks-data/data/'
filename = 'dp2_306_2_tns_search_results.ecsv'
user_table = Table.read(path + filename, format='ascii.ecsv')

Display the first five rows of the table.

In [6]:

user_table[:5]

Out[6]:

Table length=5

ra	dec	name	type
float64	float64	str10	str16
220.31468749999996	11.141855555555555	SN2026dej	SNIa
170.92526666666663	23.08451388888889	SN2026cnu	SNII
128.99505416666668	24.6421	SN2026cni	SLSN-II
222.71064583333333	12.544169444444444	SN2026cmq	SLSN-II
185.95345416666663	10.346022222222222	SN2026cfv	SNIa

3. Cross-match to the DiaObject table¶

Create the query that joins the TNS search results table with the DiaObject table using a radius of 1". Display the query.

In [7]:

search_radius = 1.0/3600.0

query = """
        SELECT diao.diaObjectId, diao.ra, diao.dec, diao.nDiaSources,
        ut.ra AS ut_ra, ut.dec AS ut_dec,
        ut.name AS ut_name, ut.type AS ut_type,
        DISTANCE(POINT('ICRS', diao.ra, diao.dec),
        POINT('ICRS', ut.ra, ut.dec)) AS distance
        FROM dp2.DiaObject AS diao
        JOIN TAP_UPLOAD.ut AS ut
        ON DISTANCE(POINT('ICRS', diao.ra, diao.dec),
        POINT('ICRS', ut.ra, ut.dec)) < {}
        """.format(search_radius)

print(query)

        SELECT diao.diaObjectId, diao.ra, diao.dec, diao.nDiaSources,
        ut.ra AS ut_ra, ut.dec AS ut_dec,
        ut.name AS ut_name, ut.type AS ut_type,
        DISTANCE(POINT('ICRS', diao.ra, diao.dec),
        POINT('ICRS', ut.ra, ut.dec)) AS distance
        FROM dp2.DiaObject AS diao
        JOIN TAP_UPLOAD.ut AS ut
        ON DISTANCE(POINT('ICRS', diao.ra, diao.dec),
        POINT('ICRS', ut.ra, ut.dec)) < 0.0002777777777777778

Execute the query, uploading the table to the TAP service.

In [8]:

job = tap_service.submit_job(query, uploads={"ut": user_table})
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

Retrieve the query results as an Astropy table, and print the number of rows.

In [9]:

assert job.phase == 'COMPLETED'
results = job.fetch_result().to_table()
print(len(results))

3.1. Check for duplicate matches¶

Get the number of unique matches between the SN in the table and rows of the DiaObject table (i.e., the number of unique elements of the ut_name (the SN name) column).

For every SN that had more than one cross-match result to the DiaObject table, print its name and the number of matches. Then for each match, print the diaObjectId, its number of detections, and its distance from the SN in arcseconds.

In [10]:

unique_matches, count_matches = np.unique(results['ut_name'], return_counts=True)
print(f"Of the {len(results)} cross-matches returned, only {len(unique_matches)} are unique.")

tx = np.where(count_matches > 1)[0]
for x in tx:
    print(' ')
    print(unique_matches[x], count_matches[x])
    ix = np.where(results['ut_name'] == unique_matches[x])[0]
    for i in ix:
        print(results['diaObjectId'][i],
              results['nDiaSources'][i],
              np.round(results['distance'][i], 6))
    del ix
del tx

Of the 140 cross-matches returned, only 136 are unique.

SN2025hlf 3
794819432252703796 9 0.0002
794819432252702927 17 0.000174
794819432252702737 147 5.4e-05

SN2025kdz 2
775517643066048695 2 0.000225
775517643066048567 27 8.4e-05

SN2025saa 2
778979592865055076 1 0.000208
778979592865054779 36 0.000161

In this case, for all three SNe with multiple matches to the DiaObject table, one of the matches looks to be the correct one: the distance is significantly lower and it has many more detections. The others are likely to be artifacts. For a real scientific analysis it is recommended to follow up by exploring the matches further and confirming they're extraneous (e.g., check for detection flags, look at the image cutouts). However for this demo, it is acceptable to remove the duplicates and move on.

Create a list of the row indexes for each of the duplicates and remove them from the table.

In [11]:

tx = np.where(count_matches > 1)[0]
temp = []
for x in tx:
    ix = np.where(results['ut_name'] == unique_matches[x])[0]
    mx = np.argmin(results['distance'][ix])
    for m, i in enumerate(ix):
        if m != mx:
            temp.append(i)
    del ix, mx

del_rows = np.asarray(temp, dtype='int')
print('Deleting the rows for:')
print(results['diaObjectId'][del_rows])

results.remove_rows([del_rows])

del temp, tx, del_rows

Deleting the rows for:
   diaObjectId
------------------
794819432252703796
794819432252702927
775517643066048695
778979592865055076

In [12]:

del unique_matches, count_matches

Also check how many diaObject rows might have been multiply matched to different SNe (it turns out that none were).

In [13]:

unique_objects, count_objects = np.unique(results['diaObjectId'], return_counts=True)
print(f"Of the {len(results)} cross-matches returned, only {len(unique_objects)} are unique.")
del unique_objects, count_objects

Of the 136 cross-matches returned, only 136 are unique.

3.2. Visualize the sample¶

Create a sky map of the locations of the cross-matched SNe from TNS, plotted over a 2-dimensional histogram of the number of observations (visits) in any filter.

First, retrieve the center coordinates of all visits from the Visit table.

In [14]:

query = """SELECT visit, ra, dec FROM dp2.Visit"""
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'
visit_table = job.fetch_result().to_table()
job.delete()
del query

Job phase is COMPLETED

Next, create the 2D histogram with a McBryde sky projection, using hex bins about the same size as the LSSTCam field of view (nside=19). Use a truncated version of the Blues colormap so that the low- and high-density regions are not too light or too dark. Overplot the RA, Dec of the TNS SNe as black crosses.

In [15]:

orig_cmap = plt.get_cmap('Blues')
new_cmap = truncate_colormap(orig_cmap, 0.2, 0.7)
fig, ax = plt.subplots(figsize=(10, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
vras = visit_table['ra'].value.data
vdecs = visit_table['dec'].value.data
sp.draw_hpxbin(vras, vdecs, nside=19, alpha=1, cmap=new_cmap)
sp.draw_colorbar(label='Number of DP2 visits (any filter)', shrink=0.5, pad=0.01)
sp.ax.plot(results['ra'], results['dec'], 'x', ms=5, color='black', label='TNS SN')
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
sp.ax.legend(loc="lower left")
plt.tight_layout()
plt.show()
del sp, vras, vdecs, new_cmap, orig_cmap

No description has been provided for this image

Figure 1: The sky location of all TNS SNe (black crosses) plotted over the 2D histogram of the DP2 visits (blue hexbins).

Delete the retrieved visit table, as it is no longer needed.

In [16]:

del visit_table

Plot a histogram of the number of detections (nDiaSources) per object.

In [17]:

fig = plt.figure(figsize=(6, 4))
plt.hist(results['nDiaSources'], bins=30)
plt.xlabel('Number of detections per object', fontsize=12)
plt.ylabel('Number of objects', fontsize=12)
plt.title('Distribution of the number of detections in any filter')
plt.show()

Figure 2: A histogram of the number of 5$\sigma$ difference-image detections, per cross-matched diaObject, showing that only a few were very well-observed with $>100$ detections.

4. Retrieve all light curves¶

Create a query to retrieve the light curves for all cross-matched objects.

The recommended data to use for transient light curves are the forced difference-image fluxes from the ForcedSourceOnDiaObject table. Observation dates are obtained by joining on the Visit table. The photometry for all identified diaObjectId values is retrieved by using a WHERE ... IN () statement, which requires a string-formatted tuple list of the diaObjectId.

Since there are 136 objects in the list, printing the query is left as optional.

In [18]:

id_array = np.asarray(results['diaObjectId'], dtype='long')
id_list = "(" + ", ".join(str(value) for value in id_array) + ")"
query = """SELECT fsodo.diaObjectId, 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 IN {}""".format(id_list)
# print(query)
del id_array, id_list

Execute the query.

In [19]:

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
7996

Display a few lines of the results to show that the photometry for multiple objects is stored in the same table.

In [20]:

flux_table.pprint(max_lines=10)

   diaObjectId     psfDiffFlux psfDiffFluxErr band    expMidptMJD
                       nJy          nJy                    d
------------------ ----------- -------------- ---- ------------------
761074595762536458    -13.9593        629.509    u 60878.207801085344
761074595762536458    -97.1393        1853.14    y  60864.15129082145
               ...         ...            ...  ...                ...
770565030377685241     16583.4        662.598    i  60939.16879158222
770565030377685241     16581.0        515.116    i  60940.10265472991
770565030377685241     15296.1        566.372    i  60940.11326337544
Length = 7996 rows

4.1. Plot one light curve¶

Identify the index (mx) of the SN with the highest number of detections, and then identify the indices (tx) for all rows of the flux_table that are associated with that object. Plot the light curve for that object only.

In [21]:

mx = np.argmax(results['nDiaSources'])
tx = np.where(flux_table['diaObjectId'] == results['diaObjectId'][mx])[0]
plot_title = """{} ({})    LSST-DP2-DO-{}
              """.format(results['ut_name'][mx], results['ut_type'][mx],
                         results['diaObjectId'][mx])

fig = plt.figure(figsize=(8, 4))
for f, filt in enumerate(filter_names):
    fx = np.where(flux_table['band'][tx] == filt)[0]
    plt.errorbar(flux_table['expMidptMJD'][tx[fx]], flux_table['psfDiffFlux'][tx[fx]],
                 yerr=flux_table['psfDiffFluxErr'][tx[fx]],
                 fmt=filter_symbols[filt], ms=10, mew=2, mec=filter_colors[filt],
                 ecolor=filter_colors[filt], alpha=0.5, color='none', label=filt)
    del fx
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.title(plot_title)
plt.show()

del mx, tx, plot_title

Figure 3: Forced fluxes measured on the difference image, as a function of MJD, for the diaObject cross-matched to SN2025mbp which was classified as a Type Ia supernova.

In [22]:

del results, flux_table

5. Single query with four-table join¶

The two individual queries in Section 3 and 4 can be combined into a single four-table join as shown below, but it is not recommended for most cases. The main drawback is that the resulting table is larger, includes photometry for the duplicates, and contains many duplicated elements because the per-object data is repeated for all rows of the photometry.

In [23]:

search_radius = 1.0/3600.0

query = """
        SELECT diao.diaObjectId, diao.ra, diao.dec, diao.nDiaSources,
        ut.ra AS ut_ra, ut.dec AS ut_dec, ut.name AS ut_name, ut.type AS ut_type,
        DISTANCE(POINT('ICRS', diao.ra, diao.dec), POINT('ICRS', ut.ra, ut.dec)) AS distance,
        fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, fsodo.band, v.expMidptMJD
        FROM dp2.DiaObject AS diao
        JOIN TAP_UPLOAD.ut AS ut
        ON DISTANCE(POINT('ICRS', diao.ra, diao.dec), POINT('ICRS', ut.ra, ut.dec)) < {}
        JOIN dp2.ForcedSourceOnDiaObject AS fsodo ON fsodo.diaObjectId = diao.diaObjectId
        JOIN dp2.Visit AS v ON v.visit = fsodo.visit
        """.format(search_radius)

print(query)

        SELECT diao.diaObjectId, diao.ra, diao.dec, diao.nDiaSources,
        ut.ra AS ut_ra, ut.dec AS ut_dec, ut.name AS ut_name, ut.type AS ut_type,
        DISTANCE(POINT('ICRS', diao.ra, diao.dec), POINT('ICRS', ut.ra, ut.dec)) AS distance,
        fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, fsodo.band, v.expMidptMJD
        FROM dp2.DiaObject AS diao
        JOIN TAP_UPLOAD.ut AS ut
        ON DISTANCE(POINT('ICRS', diao.ra, diao.dec), POINT('ICRS', ut.ra, ut.dec)) < 0.0002777777777777778
        JOIN dp2.ForcedSourceOnDiaObject AS fsodo ON fsodo.diaObjectId = diao.diaObjectId
        JOIN dp2.Visit AS v ON v.visit = fsodo.visit

Optionally, retrieve the query and show that the resulting table has 8651 rows, and that many of the elements are duplicated (an inefficient use of memory).

In [24]:

# job = tap_service.submit_job(query, uploads={"ut": user_table})
# 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_results = job.fetch_result().to_table()
# print(len(flux_table_results))

In [25]:

# flux_table_results

In [ ]: