308.2. Solar System object observations#

308_02_Observations

308.2. Solar System object observations¶

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

Learning objective: An overview of the Data Preview 2 (DP2) Solar System object observations.

LSST data products: SSSource, SSObject

Packages: lsst.rsp, lsst.daf.butler, 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¶

This notebook examines the 8,136,150 observations of the 299,343 Solar System objects detected by Rubin in the DP2 data release.

Related tutorials: The 100-level tutorials demonstrate how to use the TAP service. The 200-level tutorials introduce the types of catalog data.

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 for retrieving datasets from the butler from the LSST Science Pipelines (pipelines.lsst.io). Additional modules support standardized multiband plotting (lsst.utils.plotting) for LSST data analysis and visualization.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker

from lsst.rsp import RSPDiscovery
from lsst.daf.butler import Butler
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")
service = discovery.get_tap_client()

Create an instance of the Rubin data butler.

In [3]:

butler = Butler("dp2", collections="dp2")

Define parameters to use colorblind-friendly colors with matplotlib.

In [4]:

plt.style.use('seaborn-v0_8-colorblind')
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']

Define the colors and symbols to represent the LSST filters in plots.

In [5]:

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

2. Sky coverage¶

Query the SSSource table to retrieve the on-sky coordinates for all Solar System sources in DP2. The query may take several minutes to complete.

In [6]:

query = "SELECT ssObjectId, "\
        "eclLambda, eclBeta, "\
        "ephRa, ephDec "\
        "FROM dp2.SSSource " \
        "ORDER BY ssObjectId ASC "
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

Fetch the results as an astropy table. There are 8,136,150 Solar System source measurements (detections) in DP2.

In [7]:

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

Option to print the results.

In [8]:

# results

Plot the sky coverage in ecliptic and equatorial coordinates for all DP2 Solar System sources.

In [9]:

fig, axes = plt.subplots(1, 2, figsize=(16, 5),
                         subplot_kw={'projection': 'mollweide'})

ax = axes[0]
lam = np.deg2rad(results['eclLambda'])
lam = (lam + np.pi) % (2 * np.pi) - np.pi
bet = np.deg2rad(results['eclBeta'])
ax.scatter(lam, bet, s=0.02, alpha=0.03,
           rasterized=True)
ax.set_title('Ecliptic coordinates', pad=15)
ax.grid(True, alpha=0.3)

ax = axes[1]
ra = np.deg2rad(results['ephRa'])
ra = (ra + np.pi) % (2 * np.pi) - np.pi
dec = np.deg2rad(results['ephDec'])
ax.scatter(ra, dec, s=0.02, alpha=0.03,
           rasterized=True)
ax.set_title('Equatorial coordinates', pad=15)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

No description has been provided for this image

Figure 1: Sky coverage in ecliptic coordinates (left) and equatorial coordinates (right) for the ~6% subset of DP2 Solar System sources.

3. Observations per object¶

3.1. Multiple detections¶

Find the SS objects with more than one SS source detected in the SSSource table. The vast majority (281,856) of the 299,343 DP2 Solar System objects have more than one detection. 17,487 DP2 Solar System objects have a single detection.

In [10]:

values, counts = np.unique(results['ssObjectId'], return_counts=True)
tx = np.where(counts > 1)[0]
print('Solar System objects with >1 detection: ', len(tx))

Solar System objects with >1 detection:  281856

Plot a histogram of the number of sources for each object with more than one detection.

In [11]:

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

ax.hist(counts[tx].clip(0, 100), bins=50,
        edgecolor='none', alpha=0.85)
plt.xlabel('Number of sources per object (clipped at 100)')
plt.ylabel('Count')
plt.title('Number of sources with more than one detection')
ax.axvline(np.median(counts[tx]), color='black', linestyle='--', lw=1.5)
ax.text(0.95, 0.92,
        f'median = {np.median(counts[tx]): .0f}',
        transform=ax.transAxes, ha='right', fontsize=11,
        fontstyle='italic')

plt.tight_layout()
plt.show()

Figure 2: Histogram of the number of source detections for every Solar System object with more than one detection.

Clean up.

In [12]:

job.delete()
del query, results, values, counts, tx

3.2. Bands¶

The number of observations for each ssObjectId can also be retrieved from the SSObject table.

TAP is the recommended way to access the DP2 tables, but the butler can be used to retrieve the same table information. Define the columns to retrieve: observation arc, number of observations, and number of observations in each band.

In [13]:

col_list = ['ssObjectId', 'arc', 'nObs',
            'u_nObs', 'g_nObs', 'r_nObs',
            'i_nObs', 'z_nObs', 'y_nObs']

Get the data from the butler.

In [14]:

results = butler.get('ss_object', parameters={'columns': col_list})
print(len(results))

Option to display the results.

In [15]:

# results

Plot histograms of the number of observations, observing arcs, and total number of per-band Solar System source detections for all Solar System object sources.

In [16]:

fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))

ax = axes[0]
nobs = results['nObs']
ax.hist(nobs.clip(0, 100), bins=100,
        edgecolor='none', alpha=0.85)
ax.set_xlabel('Total nObs (clipped at 100)')
ax.set_ylabel('Count')
ax.set_title('Observation count')

ax = axes[1]
arc = results['arc']
ax.hist(arc[arc > 0], bins=100,
        edgecolor='none', alpha=0.85)
ax.set_xlabel('Arc (days)')
ax.set_ylabel('Count')
ax.set_title('Observing arc')

ax = axes[2]
filter_nobs = [results[f'{b}_nObs'].sum() for b in filter_names]
ax.bar(filter_names, filter_nobs,
       color=[filter_colors[b] for b in filter_names])
ax.set_ylabel('Total observations')
ax.set_title('Observations by band')

plt.tight_layout()
plt.show()

Figure 3: Histograms of the number of observations, observing arcs, and total number of per-band Solar System source detections for all DP2 Solar System object sources.

Clean up.

In [17]:

del col_list, results

4. Distances¶

Query the SSSource table to retrieve the heliocentric and topocentric distances (helioRange and topoRange, respectively), phase angle, solar elongation, and on-sky motion rate (ephRate) for all Solar System sources in DP2. The query may take several minutes to complete.

In [18]:

query = "SELECT ssObjectId, "\
        "helioRange, topoRange, "\
        "phaseAngle, elongation, "\
        "ephRate "\
        "FROM dp2.SSSource " \
        "ORDER BY ssObjectId ASC "
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

Fetch the results as an astropy table.

In [19]:

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

Option to print the results.

In [20]:

# results

Plot histograms of the heliocentric and topocentric distances for all DP2 Solar System sources, comparing the inner, outer, and distant Solar System.

In [21]:

fig, axes = plt.subplots(1, 3, figsize=(12, 3))

ax = axes[0]
ax.hist(results['helioRange'].clip(0, 5), bins=100,
        edgecolor='none', alpha=0.85, color='blue',
        label='Heliocentric')
ax.hist(results['topoRange'].clip(0, 5), bins=100,
        edgecolor='none', alpha=0.6, color='red',
        label='Topocentric')
ax.set_xlabel('On-sky rate (deg/day)')
ax.set_ylabel('Count')
ax.set_title('Distance distributions (< 5 au)')
ax.legend()

ax = axes[1]
ax.hist(results['helioRange'].clip(5, 65), bins=50,
        edgecolor='none', alpha=0.85, color='blue',
        label='Heliocentric')
ax.hist(results['topoRange'].clip(5, 65), bins=100,
        edgecolor='none', alpha=0.6, color='red',
        label='Topocentric')
ax.set_xlim(5., 65.)
ax.set_ylim(0., 2600)
ax.set_xlabel('On-sky rate (deg/day)')
ax.set_ylabel('Count')
ax.set_title('Distance distributions (5-65 au)')
ax.legend()

ax = axes[2]
ax.hist(results['helioRange'].clip(65, 100), bins=50,
        edgecolor='none', alpha=0.85, color='blue',
        label='Heliocentric')
ax.hist(results['topoRange'].clip(65, 100), bins=50,
        edgecolor='none', alpha=0.6, color='red',
        label='Topocentric')
ax.set_xlim(64., 91.)
ax.set_ylim(0., 51.)
ax.set_xlabel('On-sky rate (deg/day)')
ax.set_ylabel('Count')
ax.set_title('Distance distributions (>65 au)')
ax.legend()

plt.tight_layout()
plt.subplots_adjust(wspace=0.3)
plt.show()

Figure 4: Histograms of the heliocentric and topocentric distances for all DP2 Solar System sources, comparing the inner (left), outer (middle), and distant (right) Solar System.

5. On-sky motion rate¶

Plot a histogram of the on-sky motion rate and a scatter plot of the on-sky motion rate as a function of phase angle.

In [22]:

fig, axes = plt.subplots(1, 2, figsize=(12, 3))

ax = axes[0]
ax.hist(results['ephRate'].clip(0, 0.5), bins=100,
        edgecolor='none', alpha=0.85)
ax.set_xlabel('On-sky rate (deg/day)')
ax.set_ylabel('Count')
ax.set_title('On-sky motion rate')
ax.text(0.95, 0.92,
        f'median = {np.median(results['ephRate']): .3f} deg/day',
        transform=ax.transAxes, ha='right', fontsize=11,
        fontstyle='italic')

ax = axes[1]
ax.scatter(results['phaseAngle'], results['ephRate'], s=0.3)
ax.set_xlabel('Phase angle (deg)')
ax.set_ylabel('On-sky rate (deg/day)')

plt.subplots_adjust(wspace=0.3)
plt.show()

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

Figure 5: Histogram of on-sky motion rates with the median rate labeled (left) and scatter plot of phase angle versus on-sky motion rates (right) for all DP2 Solar System source measurements.

6. Phase angle and solar elongation¶

Plot histograms of the phase angle and solar elongation distributions.

In [23]:

fig, axes = plt.subplots(1, 2, figsize=(12, 3))

ax = axes[0]
ax.hist(results['phaseAngle'], bins=200,
        edgecolor='none', alpha=0.85)
ax.set_xlabel('Phase angle (deg)')
ax.set_ylabel('Count')
ax.set_title('Phase angle')
ax.text(0.95, 0.92,
        f'{results["phaseAngle"].min(): .1f}\u00b0 \u2013 '
        f'{results["phaseAngle"].max(): .1f}\u00b0',
        transform=ax.transAxes, ha='right', fontsize=11,
        fontstyle='italic')

ax = axes[1]
ax.hist(results['elongation'], bins=200,
        edgecolor='none', alpha=0.85)
ax.set_xlabel('Solar elongation (deg)')
ax.set_ylabel('Count')
ax.set_title('Solar elongation')
ax.text(0.95, 0.92,
        f'min = {results["elongation"].min(): .1f}\u00b0',
        transform=ax.transAxes, ha='right', fontsize=11,
        fontstyle='italic')

plt.subplots_adjust(wspace=0.3)
plt.show()

Figure 6: Histograms of the phase angle (left) and solar elongation (right) distributions for all DP2 Solar System source measurements. The phase angle range and minimum solar elongation are labeled on each respective plot.

Plot solar elongation as a function of heliocentric distance for every Solar System source.

In [24]:

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

ax.scatter(results['helioRange'], results['elongation'], s=0.5)
plt.xlabel('Heliocentric Distance (au)')
plt.ylabel('Solar elongation (deg)')
plt.tight_layout()
ax.xaxis.set_minor_locator(ticker.MultipleLocator(5))
ax.grid(which='major', color='gray', linestyle='-')
ax.grid(which='minor', color='gray', linestyle='--')

plt.show()

Figure 7: Solar elongation as a function of heliocentric distance for every DP2 Solar System source. Notice a small cluster of difference image detections at low solar elongation (SE; roughly 38 - 48 deg). All but one of the low-SE detections are at a heliocentric distance < 5 au; the single low-SE outlier was observed at a heliocentric distance of nearly 30 au out near Neptune's semimajor axis. The highest concentrations in heliocentric distance (r) are at r < 5 au and roughly 30 au < r < 48 au (the Kuiper belt beyond Neptune); a handful of objects were detected at large heliocentric distances out to almost 90 au.

Clean up.

In [25]:

del query, results