308.2. Orbits

308.2. Solar System orbit classifications¶

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

Learning objective: An overview of the Data Preview 1 (DP1) Solar System object orbits.

LSST data products: MPCORB, SSObject

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 notebook examines the orbital classifications of the DP1 Solar System object, including the DP1 discoveries. The DP1 data release contains 431 Solar System objects, 93 of those are discoveries.

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 pandas as pd
import matplotlib.pyplot as plt

from lsst.rsp import get_tap_service

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

Define parameters to use colorblind-friendly colors with matplotlib.

In [3]:

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

Define a function to calculate orbital classifications from object semimajor axes $a$, eccentricities $e$, and inclinations $incl$, where AMO = Amor NEO, APO = Apollo NEO, ATE = Aten NEO, IEO = inner-Earth object NEO, IMB = inner main-belt asteroid, MBA = main-belt asteroid, MCA = Mars-crosser asteroid, OMB = outer main-belt asteroid, TJN = Jupiter Trojan, CEN = Centaur, TNO = transneptunian object, PAA = parabolic asteroid, HYA = hyperbolic asteroid, closely following the orbit type definitions from the Minor Planet Center.

In [4]:

def calc_orbit_class(q, e):
    a = q / (1 - e)
    Q = a * (1 + e)

    orbit_class = np.full(len(a), "AST")

    orbit_class_dict = {
        "AMO": np.where((a > 1.0) & (q > 1.017) & (q < 1.3)),
        "APO": np.where((a > 1.0) & (q < 1.017)),
        "ATE": np.where((a < 1.0) & (Q > 0.983)),
        "CEN": np.where((a > 5.5) & (a < 30.1)),
        "IEO": np.where((Q < 0.983)),
        "IMB": np.where((a < 2.0) & (q > 1.666)),
        "MBA": np.where((a > 2.0) & (a < 3.2) & (q > 1.666)),
        "MCA": np.where((a < 3.2) & (q > 1.3) & (q < 1.666)),
        "OMB": np.where((a > 3.2) & (a < 4.6)),
        "TJN": np.where((a > 4.6) & (a < 5.5) & (e < 0.3)),
        "TNO": np.where((a > 30.1)),
        "PAA": np.where((e == 1)),
        "HYA": np.where((e > 1)),
    }
    for c, v in orbit_class_dict.items():
        orbit_class[v] = c

    return orbit_class

2. Classify orbits¶

Query the MPCORB table for the orbital parameters joined with the SSObject table to query for the discoverySubmissionDate for all DP1 Solar System objects.

In [5]:

query = "SELECT mpc.ssObjectId, mpc.q, mpc.e, mpc.incl, "\
        "sso.discoverySubmissionDate "\
        "FROM dp1.MPCORB as mpc "\
        "INNER JOIN dp1.SSObject as sso "\
        "ON mpc.ssObjectId = sso.ssObjectId "\
        "ORDER BY mpc.ssObjectId"

In [6]:

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 job results and assign to an astropy result table. There are 431 Solar System objects in DP1.

In [7]:

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

431

Convert the result astropy table to a pandas dataframe.

In [8]:

result_df = pd.DataFrame(result)

Add a column to result_df with each object's orbit_class.

In [9]:

result_df['orbit_class'] = calc_orbit_class(result_df['q'], result_df['e'])

Option to print the results

In [10]:

# result_df

Find the objects with non-zero discoverySubmissionDates. These are the Solar System objects discovered by Rubin. 93 of the 431 DP1 Solar System objects are discoveries.

In [11]:

discoveries = result_df[result_df['discoverySubmissionDate'] > 0]
len(discoveries)

Out[11]:

93

Plot the semimajor axes, eccentricities, and inclinations of the 431 DP1 Solar System objects as well as the 93 DP1 Solar System object discoveries, distinguishing each orbital class. The DP1 Solar System object discoveries include

In [12]:

fig, axes = plt.subplots(1, 2, figsize=(8, 4))

for ax, oo in zip(axes, (result_df, discoveries)):
    markers = "s^*+xDh"*2
    for label, orb in oo.groupby('orbit_class'):
        if len(orb) > 500:
            s, marker = 0.1, 'o'
        else:
            s = 10
            if label == 'APO':
                marker = "s"
                color = colors[0]
            elif label == 'ATE':
                marker = "^"
                color = colors[1]
            elif label == 'IMB':
                marker = "*"
                color = colors[2]
            elif label == 'MBA':
                marker = "+"
                color = colors[3]
            elif label == 'MCA':
                marker = "x"
                color = colors[4]
            elif label == 'OMB':
                marker = "D"
                color = colors[5]
            elif label == 'TJN':
                marker = "h"
                color = colors[0]
        ax.scatter((orb["q"]/(1-orb["e"])), orb["incl"],
                   s=s, marker=marker, color=color,
                   label=f"{label} ({len(orb)})")

    for handle in ax.legend().legend_handles:
        handle.set_sizes([10])
    ax.set_xlabel("semimajor axis (au)")
    ax.set_ylabel("inclination (deg)")
    ax.set_xlim([0.8, 5.4])
    ax.set_ylim([0, 40.5])

fig.suptitle("DP1 Solar System Object Detections (Left) and Discoveries (Right)")

fig, axes = plt.subplots(1, 2, figsize=(8, 4))

for ax, oo in zip(axes, (result_df, discoveries)):
    markers = "s^*+xDh"*2
    for label, orb in oo.groupby('orbit_class'):
        if len(orb) > 500:
            s, marker = 0.1, 'o'
        else:
            s = 10
            if label == 'APO':
                marker = "s"
                color = colors[0]
            elif label == 'ATE':
                marker = "^"
                color = colors[1]
            elif label == 'IMB':
                marker = "*"
                color = colors[2]
            elif label == 'MBA':
                marker = "+"
                color = colors[3]
            elif label == 'MCA':
                marker = "x"
                color = colors[4]
            elif label == 'OMB':
                marker = "D"
                color = colors[5]
            elif label == 'TJN':
                marker = "h"
                color = colors[0]
        ax.scatter((orb["q"]/(1-orb["e"])), orb["e"],
                   s=s, marker=marker, color=color,
                   label=f"{label} ({len(orb)})")

    for handle in ax.legend().legend_handles:
        handle.set_sizes([10])
    ax.set_xlabel("semimajor axis (au)")
    ax.set_ylabel("eccentricity")
    ax.set_xlim([0.8, 5.4])
    ax.set_ylim([0, 0.6])

Figure 1: Semimajor axes, eccentricities, and inclinations of the 431 DP1 Solar System objects (left) and the 93 DP1 Solar System object discoveries (right), with the colors distinguishing each orbital class. The number of objects in each orbital class is indicated in the figure legends.

Clean up.

In [13]:

del result, result_df, discoveries

In [ ]: