308.6. SSO lightcurves

308.6. SSO Lightcurves¶

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

Learning objective: How to retrieve Solar System Object photometry for lightcurve analysis in Data Preview 1.

LSST data products: SSObject, MPCORB, DiaSource, SSSource

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

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¶

The observed brightness of a Solar System Object (SSO) changes primarily due to the following factors:

• Heliocentric and observer (topocentric) distances
• Viewing geometry, i.e. the Sun-SSO-observer phase angle
• Rotational variation, most often shape-driven

For active SSOs (comets and active asteroids), the observed photometry is also affected by the photometric properties of the emitted dust and gas and is more complicated to model.

This notebook demonstrates how to retrieve SSO photometry from Data Preview 1 (DP1) and how to analysis the lightcurve, accounting for viewing distances and geometry.

Warning: DP1 presents only a fraction of the data that will be released with the full LSST. As such the lightcurve and phase angle coverage of SSOs in DP1 is limited* and this science analysis notebook is for demonstration purposes.

As part of Rubin First Look (RFL) and Solar System Processing (SSP) testing there has been a significant amount of LSSTCam SSO photometry submitted to the Minor Planet Center (MPC). Furthermore the pre-LSST alert stream contains suitable coverage of a number of SSOs. These data are more suitable for full lightcurve analysis, as demonstrated by recent publications by Greenstreet et al. 2026 and Kurlander et al. 2026, but require different methods to access and analyse compared to the upcoming LSST prompt products and data releases. See Prompt Products documentation for information on alerts and an MPC photometry tutorial.

Related tutorials: The 100-level tutorials demonstrate how to use the TAP service. The 200-level tutorials introduce the types of catalog data; see 201.5, 201.7, 201.8 and 201.9 for information on the DiaSource, SSObject, SSSource and MPCORB tables respectively. The 300-level tutorials demonstrate science analysis; see 308.4 for details on SSO photometry.

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). Import astropy.units to convert between fluxes and magnitudes.

From the lsst package, import modules for accessing the Table Access Protocol (TAP) service from the LSST Science Pipelines (pipelines.lsst.io).

In [1]:

import numpy as np
import matplotlib.pyplot as plt

import astropy.units as u

from lsst.rsp import get_tap_service
from lsst.utils.plotting import (get_multiband_plot_colors,
                                 get_multiband_plot_symbols,
                                 get_multiband_plot_linestyles)

1.2. Define parameters and functions¶

Define colors, symbols, and linestyles to represent the six LSST filters, $ugrizy$.

In [2]:

filter_names = ['u', 'g', 'r', 'i', 'z', 'y']
filter_colors_white = get_multiband_plot_colors()
filter_colors_black = get_multiband_plot_colors(dark_background=True)
filter_symbols = get_multiband_plot_symbols()
filter_linestyles = get_multiband_plot_linestyles()

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

In [3]:

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

The DP1 MPCORB table provides the orbital elements for each ssObjectId from the Minor Planet Center database. The perihelion q and eccentricity e are provided but semi-major axis must be calculated.

Define a function to calculate the semi-major axis of an asteroid's orbit.

In [4]:

def calc_semimajor_axis(q, e):
    """
    Given a perihelion distance and orbital eccentricity,
      calculate the semi-major axis of the orbit.

    Parameters
    ----------
    q: float
        Distance at perihelion, in au.
    e: float
        Orbital eccentricity.

    Returns
    -------
    a: float
        Semi-major axis of the orbit, in au.
        q = a(1-e), so a = q/(1-e)
    """

    return q / (1.0 - e)

2. Define target¶

The data specific to SSOs in DP1 are spread across four tables (see DP1 schema). Object level information is contained in SSObject and MPCORB tables, where each row corresponds to a unique SSO and is indexed by ssObjectId. Source level information is contained in DiaSource and SSSource, where each row is a unique observation of an SSO indexed by diaSourceId and associated with an ssObjectId to identify the object. As such querying data for SSOs will often involve joining tables on ssObjectId and/or diaSourceId.

Warning: DP1 does not contain all the fields as listed in the Data Products Definition Document. This is especially true for the SSObject table.

Warning: Future data releases will treat the association of diaSourceId and ssObjectID differently to DP1 - a given ssObjectId will be assigned to the diaSourceId nearest the predicted SSO position at the time of exposure (this is encapsulated by SSSource:diaDistanceRank; see Source Assocation prompt products documentation for details).

Warning: For static transients the ForcedSourceOnDiaObject table is recommended for retrieving photometry (see tutorial 205.1). As SSO DiaSources are not static it is recommended to use the DiaSource table directly.

2.1. Select an ssObjectId¶

Select an SSO with a large number of observations to look at its lightcurve. The following query retrieves ssObjectId, sorted by the total number of observations: numObs. Order the results by numObs in descending order.

In [5]:

query = "SELECT ssObjectId, numObs "\
        "FROM dp1.SSObject as sso "\
        "ORDER BY numObs DESC"

Submit the query as a job to the TAP service. Run the job and wait until the job phase is either "COMPLETED" or "ERROR", then print the job phase and the error message, if applicable.

In [6]:

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

Job URL is https://data.lsst.cloud/api/tap/async/x40a7alwncuv4bum
Job phase is COMPLETED

If the status returned was "COMPLETED", execute this cell to retrieve the query results as an astropy table using the to_table method, and assert that the table length is 173 as expected.

In [7]:

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

Print the first 10 results showing the number of observations in descending order.

In [8]:

results[:10]

Out[8]:

Table length=10

ssObjectId	numObs
int64	int32
21163620217073748	133
23133931615301680	82
21163637482928473	80
21164728252512342	66
23133931615301681	65
23133931615301682	63
21164741002015298	39
21163645768577093	33
23133931615301687	32
23133931615301942	30

Save space in memory by deleting the job and the retrieved results.

In [9]:

job.delete()
del results

2.2. Retrieve orbit information¶

Select the DP1 SSO with the most observations for lightcurve analysis (ssObjectId = 21163620217073748). Get information such as name and orbital parameters for the selected ssObjectId. This requires a JOIN between the MPCORB and SSObject tables.

In [10]:

ssoid = 21163620217073748

In [11]:

query = "SELECT mpc.mpcDesignation, "\
        "mpc.q, mpc.e, mpc.incl, "\
        "sso.numObs "\
        "FROM dp1.MPCORB as mpc "\
        "INNER JOIN dp1.SSObject as sso "\
        "ON mpc.ssObjectId = sso.ssObjectId "\
        "WHERE mpc.ssObjectId = 21163620217073748"

In [12]:

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

Job URL is https://data.lsst.cloud/api/tap/async/tq4u9tzpz59sgl7j
Job phase is COMPLETED

In [13]:

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

Calculate the orbital semi-major axis a and look at the object's orbital parameters:

In [14]:

results['a'] = calc_semimajor_axis(results['q'], results['e'])
results

Out[14]:

Table length=1

mpcDesignation	q	e	incl	numObs	a
	AU		deg
str13	float64	float64	float64	int32	float64
2003 QT30	1.70936064352526	0.09492430251121	26.6295069323684	133	1.88863831861581

The full Minor Planet Center (MPC) designation of an object is retrieved from the MPCORB table. This is useful for looking up the object outside the Rubin databases. In this case ssObjectId= 21163620217073748 corresponds to mpcDesignation = 2003 QT30, an inner main-belt asteroid.

In [15]:

mpcDes = results['mpcDesignation'][0]
mpcDes

Out[15]:

np.str_('2003 QT30')

Clean up.

In [16]:

job.delete()
del results

2.3. Retrieve observation data¶

To get all information for observations of the SSO, the query must JOIN the tables DiaSource and SSSource.

Flag columns are used to identify observations that the photometry algorithms may have had issues with, these are also retrieved. See tutorial 308.4 Photometry for more details.

In [17]:

diaFlags = np.array(['isDipole', 'dipoleFitAttempted', 'pixelFlags',
                     'pixelFlags_offimage', 'pixelFlags_edge',
                     'pixelFlags_interpolated', 'pixelFlags_saturated', 'pixelFlags_cr',
                     'pixelFlags_bad', 'pixelFlags_suspect',
                     'pixelFlags_interpolatedCenter', 'pixelFlags_saturatedCenter',
                     'pixelFlags_crCenter', 'pixelFlags_suspectCenter', 'centroid_flag',
                     'apFlux_flag', 'apFlux_flag_apertureTruncated', 'psfFlux_flag',
                     'psfFlux_flag_noGoodPixels', 'psfFlux_flag_edge',
                     'forced_PsfFlux_flag', 'forced_PsfFlux_flag_noGoodPixels',
                     'forced_PsfFlux_flag_edge', 'shape_flag', 'shape_flag_no_pixels',
                     'shape_flag_not_contained', 'shape_flag_parent_source',
                     'trail_flag_edge', 'pixelFlags_streak', 'pixelFlags_streakCenter',
                     'pixelFlags_injected', 'pixelFlags_injectedCenter',
                     'pixelFlags_injected_template',
                     'pixelFlags_injected_templateCenter', 'pixelFlags_nodata',
                     'pixelFlags_nodataCenter'])

In [18]:

query = """SELECT dia.diaSourceId, dia.midpointMjdTai, dia.band,
           dia.psfFlux, dia.psfFluxErr, dia.scienceFlux, dia.scienceFluxErr,
           sss.topocentricDist, sss.heliocentricDist, sss.phaseAngle, {}
           FROM dp1.DiaSource as dia
           INNER JOIN dp1.SSSource as sss
           ON dia.diaSourceId = sss.diaSourceId
           WHERE dia.ssObjectId = {}
        """.format(','.join(['dia.'+x for x in diaFlags]), ssoid)
print(query)

SELECT dia.diaSourceId, dia.midpointMjdTai, dia.band,
           dia.psfFlux, dia.psfFluxErr, dia.scienceFlux, dia.scienceFluxErr,
           sss.topocentricDist, sss.heliocentricDist, sss.phaseAngle, dia.isDipole,dia.dipoleFitAttempted,dia.pixelFlags,dia.pixelFlags_offimage,dia.pixelFlags_edge,dia.pixelFlags_interpolated,dia.pixelFlags_saturated,dia.pixelFlags_cr,dia.pixelFlags_bad,dia.pixelFlags_suspect,dia.pixelFlags_interpolatedCenter,dia.pixelFlags_saturatedCenter,dia.pixelFlags_crCenter,dia.pixelFlags_suspectCenter,dia.centroid_flag,dia.apFlux_flag,dia.apFlux_flag_apertureTruncated,dia.psfFlux_flag,dia.psfFlux_flag_noGoodPixels,dia.psfFlux_flag_edge,dia.forced_PsfFlux_flag,dia.forced_PsfFlux_flag_noGoodPixels,dia.forced_PsfFlux_flag_edge,dia.shape_flag,dia.shape_flag_no_pixels,dia.shape_flag_not_contained,dia.shape_flag_parent_source,dia.trail_flag_edge,dia.pixelFlags_streak,dia.pixelFlags_streakCenter,dia.pixelFlags_injected,dia.pixelFlags_injectedCenter,dia.pixelFlags_injected_template,dia.pixelFlags_injected_templateCenter,dia.pixelFlags_nodata,dia.pixelFlags_nodataCenter
           FROM dp1.DiaSource as dia
           INNER JOIN dp1.SSSource as sss
           ON dia.diaSourceId = sss.diaSourceId
           WHERE dia.ssObjectId = 21163620217073748

In [19]:

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

assert job.phase == 'COMPLETED'
results = job.fetch_result()

Store the results in a pandas dataframe as df_obs.

In [21]:

df_obs = results.to_table().to_pandas()

Option to print results.

In [22]:

# df_obs

Clean up.

In [23]:

job.delete()
del results

2.4. Calculate Magnitudes¶

Rubin reports fluxes, this is especially important for forced photometry and non-detections of static transients. However it is common to analyze SSO lightcurves with magnitude, therefore we use astropy.units to convert flux nJy to ABmag.

Warning DiaSource detections can have negative flux, in these cases the RuntimeWarning: invalid value encountered in log10 will be raised.

It is also possible to perform the flux to magnitude conversion in the ADQL query using scisql_nanojanskyToAbMag (see tutorial 103.3) but this function is not recommended when negative fluxes are present.

In [24]:

df_obs['psfABmag'] = (np.array(df_obs['psfFlux']) * u.nJy).to(u.ABmag).value
df_obs['scienceABmag'] = (np.array(df_obs['scienceFlux']) * u.nJy).to(u.ABmag).value

df_obs['psfABmagErr'] = ((2.5 / np.log(10)) * (df_obs['psfFluxErr'] / df_obs['psfFlux']))
df_obs['scienceABmagErr'] = ((2.5 / np.log(10)) * (df_obs['scienceFluxErr'] / df_obs['scienceFlux']))

2.5. Reduced Magnitude¶

The reduced magnitude of an SSO accounts for viewing distances at the time of observation, where apparent magnitude is corrected to a heliocentric and topocentric distance of 1 au. By removing the effects of observing distance, the reduced magnitude lightcurve is determined by phase angle effects and rotational brightness variations.

In [25]:

thdist = df_obs['topocentricDist']*df_obs['heliocentricDist']
df_obs['red_psfABmag'] = df_obs['psfABmag'] - 5.0*np.log10(thdist)
df_obs['red_scienceABmag'] = df_obs['scienceABmag'] - 5.0*np.log10(thdist)

3. SSO Lightcurve Analysis¶

Select the type of magnitude measurement to use for the SSO lightcurve. For DP1 one should generally use the PSF measurements of flux.

In this example psfAbmag which was measured on the difference image gives a more accurate observation on the first night compared to scienceABmag, due to background contamination of an extended source in the science image. Take a closer look at the cutouts for diaSourceId = 600438911280349188 and 600438911951437875 by using tutorial 308.5 Image stamps.

In [26]:

mag_col = 'psfABmag'

Use the DiaSource flag columns to identify observations that should probably be excluded from scientific analysis.

In [27]:

any_flag_mask = df_obs[diaFlags].any(axis=1)
df_obs[any_flag_mask][['diaSourceId'] + list(diaFlags[df_obs[diaFlags].any()])]

Out[27]:

	diaSourceId	isDipole	dipoleFitAttempted	centroid_flag	apFlux_flag	psfFlux_flag	forced_PsfFlux_flag	shape_flag
8	600447733812690946	True	True	False	False	False	False	False
13	600447736631263242	True	True	False	False	False	False	False
39	600447734351134753	True	True	False	False	False	False	False
77	600452146489458689	False	False	True	True	True	True	True

Plot the magnitude brightness of the SSO as a function of time denoting observations in different filters. Indicate observations with at least one DiaSource flag; these observations should generally be excluded from further analysis.

In [28]:

x_plot = 'midpointMjdTai'
y_plot = 'red_{}'.format(mag_col)
yErr_plot = '{}Err'.format(mag_col)

df_plot = df_obs
df_plot2 = df_obs[any_flag_mask]

fig, ax = plt.subplots(figsize=(6, 4))

for i, filt in enumerate(np.unique(df_plot['band'])):
    _df_plot = df_plot[df_plot['band'] == filt]
    ax.errorbar(_df_plot[x_plot], _df_plot[y_plot], _df_plot[yErr_plot],
                c=filter_colors_white[filt],
                fmt=filter_symbols[filt],
                label=filt)

ax.scatter(df_plot2[x_plot], df_plot2[y_plot],
           c='k', marker='x',
           zorder=5, label='flagged')

ax.legend()
ax.set_xlabel(x_plot)
ax.set_ylabel(y_plot)
ax.invert_yaxis()

plt.title(f"{mpcDes} (ssObjectId = {ssoid})")

plt.show()

Figure 1: Magnitude brightness of an SSO as a function of time. Marker colour and symbol denote observations in different filters. Observations with at least one DiaSource flag are indicated with black crosses; these observations should generally be excluded from further analysis.

3.1. Phase Angle Correction¶

As an asteroid and Earth move on their heliocentric orbits the Sun-asteroid-observer viewing geometry, or phase angle, changes and causes the observed reduced magnitude to get brighter towards peak brightness as phase angle approaches zero degrees. Due to the limited time range of DP1 observations the range of phase angles for SSOs in DP1 is very small. Typical phase curve analysis for main-belt asteroids generally requires observations spanning zero to tens of degrees in phase angle. Furthermore a typical phase curve changes most rapidly at phase angles < 5 degrees, this object has been observed at large phase angle where brightness changes are less significant.

Determine the range of phase angles using numpy's peak-to-peak (np.ptp) function. The range of phase angle coverage for this SSO is 0.5 degrees.

In [29]:

np.ptp(df_obs['phaseAngle'])

Out[29]:

np.float32(0.51647377)

The mean observed phase angle is ~25 degrees.

In [30]:

np.mean(df_obs['phaseAngle'])

Out[30]:

np.float32(24.867834)

Plot the change in reduced magnitude as a function of phase angle, indicating observations in each filter.

In [31]:

x_plot = 'phaseAngle'
y_plot = 'red_{}'.format(mag_col)
yErr_plot = '{}Err'.format(mag_col)
filt_ref = 'r'

df_plot = df_obs[~any_flag_mask]

fig, ax = plt.subplots(figsize=(6, 4))

for i, filt in enumerate(np.unique(df_plot['band'])):
    _df_plot = df_plot[df_plot['band'] == filt]

    ax.errorbar(_df_plot[x_plot], _df_plot[y_plot], _df_plot[yErr_plot],
                c=filter_colors_white[filt],
                fmt=filter_symbols[filt],
                label=filt)

ax.legend()
ax.set_xlabel(x_plot)
ax.set_ylabel('{}'.format(y_plot))
ax.invert_yaxis()

plt.title(f"{mpcDes} (ssObjectId = {ssoid})")

plt.show()

Figure: Scatter plot of the change in reduced magnitude as a function of phase angle, where observations in each filter are indicated. Normally the phase curve increases in brightness as phase angle decreases, but for this object the narrow phase angle coverage means that the probable rotational variation signal outweighs the expected phase curve trend (N.B. flagged outlying photometry has already been dropped). Analysis of the phase curve is reserved for when new observations extend the phase angle coverage.

4. Exercise for the learner¶

Repeat the lightcurve analysis with visit image magnitudes scienceABmag instead of difference image magnitudes psfABmag.