Trajectories: loading data and running batch IOD

This tutorial shows how to work with TrajectorySet, the container for many objects with time‑ordered astrometric observations. You will learn how to:

import trajectories from files (MPC 80‑column and ADES),
build trajectories from in‑memory NumPy arrays,
estimate preliminary orbits for all trajectories or just one.

The heavy lifting is performed by the Rust engine; the Python API keeps things concise and composable.

Prerequisites

You will need a global environment and at least one observing site:

Register an observing site
# Observer registration example for documentation inclusion
from py_outfit import PyOutfit, Observer

env = PyOutfit("horizon:DE440", "FCCT14")

# Define a custom observing site; elevation is expressed in kilometers
obs = Observer(
    longitude=12.345,  # degrees east
    latitude=-5.0,     # degrees
    elevation=1.0,     # kilometers above MSL
    name="DemoSite",
    ra_accuracy=0.0,   # radians (optional, example value)
    dec_accuracy=0.0,  # radians (optional, example value)
)

# Register the observer in the environment
env.add_observer(obs)

# On the first use, you will see only the added custom site from below
print(env.show_observatories())

# Obtain the ZTF (Palomar) observatory by its MPC code
ztf = env.get_observer_from_mpc_code("I41")  # returns an Observer instance
print(ztf)  # human-readable summary (code, name, geodetic position)

# Listing available observatories (optional helper)
# Now that the MPC code lookup has been used, the internal registry
# has been populated with all MPC observatories, so the output is much longer.
print(env.show_observatories())  # table of currently known / registered sites

Units used in this API: angles are radians unless stated otherwise; epochs are MJD (TT, days); uncertainties may be provided in arcseconds for convenience where noted.

Import from files

MPC 80‑column

Create a set from a single MPC 80‑column file, or append into an existing set.

From MPC 80-column
# Import trajectories from an MPC 80-column file and append another
from pathlib import Path
from py_outfit import PyOutfit, TrajectorySet

# Create environment (ephemerides + error model)
env = PyOutfit("horizon:DE440", "FCCT14")

# Build from a single MPC 80-column file
mpc_path = Path("tests/data/2015AB.obs")
ts = TrajectorySet.new_from_mpc_80col(env, mpc_path)
print("n_traj=", ts.number_of_trajectories(), "total_obs=", ts.total_observations())

# Append from a second file (no de-duplication)
mpc_path2 = Path("tests/data/33803.obs")
ts.add_from_mpc_80col(env, mpc_path2)

Notes

Input parsing mirrors the Rust engine. Avoid ingesting the same file twice: no de‑duplication is performed.

ADES (JSON or XML)

When creating from ADES, you can provide global uncertainties (arcsec) if they are not specified per row.

From ADES (JSON/XML)
# Import trajectories from an ADES file (JSON or XML) and optionally append others
from pathlib import Path
from py_outfit import PyOutfit, TrajectorySet

env = PyOutfit("horizon:DE440", "FCCT14")

ades_path = Path("tests/data/example_ades.xml")
ts = TrajectorySet.new_from_ades(env, ades_path, error_ra_arcsec=0.3, error_dec_arcsec=0.3)

# Append another ADES file into the same set (avoid re-ingesting the same file)
ts.add_from_ades(env, Path("tests/data/flat_ades.xml"), 0.3, 0.3)

Build from in‑memory arrays

Two ingestion helpers are available. Use degrees/arcseconds for convenience, or supply radians for a zero‑copy path.

Degrees + arcseconds (converted once to radians)

From NumPy (degrees + arcsec)
# Build a TrajectorySet from degrees/arcseconds and MJD(TT)
import numpy as np
from py_outfit import PyOutfit, TrajectorySet, Observer

env = PyOutfit("horizon:DE440", "FCCT14")
obs = Observer(0.0, 0.0, 1.0, "DemoSite", np.deg2rad(0.3/3600.0), np.deg2rad(0.3/3600.0))
env.add_observer(obs)

trajectory_id = np.array([10, 10, 10, 11, 11, 11], dtype=np.uint32)
ra_deg        = np.array([10.0, 10.01, 10.02, 180.0, 180.02, 180.05])
dec_deg       = np.array([ 5.0,  5.01,  5.015, -10.0, -10.02, -10.03])
mjd_tt        = np.array([60000.0, 60000.01, 60000.03, 60000.0, 60000.02, 60000.05])

# Performs one conversion to radians under the hood
ts = TrajectorySet.from_numpy_degrees(
    env,
    trajectory_id,
    ra_deg,
    dec_deg,
    error_ra_arcsec=0.3,
    error_dec_arcsec=0.3,
    mjd_tt=mjd_tt,
    observer=obs,
)

Radians (zero‑copy)

From NumPy (radians, zero-copy)
# Build a TrajectorySet from radians (zero-copy) and MJD(TT)
import numpy as np
from py_outfit import PyOutfit, TrajectorySet, Observer

env = PyOutfit("horizon:DE440", "FCCT14")
obs = Observer(0.0, 0.0, 1.0, "DemoSite", np.deg2rad(0.3/3600.0), np.deg2rad(0.3/3600.0))
env.add_observer(obs)

trajectory_id = np.array([10, 10, 10, 11, 11, 11], dtype=np.uint32)
ra_deg        = np.array([10.0, 10.01, 10.02, 180.0, 180.02, 180.05])
dec_deg       = np.array([ 5.0,  5.01,  5.015, -10.0, -10.02, -10.03])
ra_rad        = np.deg2rad(ra_deg)
dec_rad       = np.deg2rad(dec_deg)
mjd_tt        = np.array([60000.0, 60000.01, 60000.03, 60000.0, 60000.02, 60000.05])

# Zero-copy path when inputs are already radians
ts = TrajectorySet.from_numpy_radians(
    env,
    trajectory_id,
    ra_rad,
    dec_rad,
    error_ra_rad=np.deg2rad(0.3 / 3600.0),
    error_dec_rad=np.deg2rad(0.3 / 3600.0),
    mjd_tt=mjd_tt,
    observer=obs,
)

Estimate orbits

You can estimate preliminary orbits for all trajectories in a set, or for a single trajectory.

Batch over all trajectories

Batch IOD across the set
# Batch orbit estimation across all trajectories
from py_outfit import PyOutfit, TrajectorySet, IODParams, Observer
import numpy as np
from astropy.time import Time

env = PyOutfit("horizon:DE440", "FCCT14")
obs = Observer(
    0.0, 0.0, 1.0, "DemoSite", np.deg2rad(0.3 / 3600.0), np.deg2rad(0.3 / 3600.0)
)
env.add_observer(obs)

# Minimal synthetic data (single trajectory)
trajectory_id = np.array(
    [0, 1, 2, 1, 2, 1, 0, 0, 0, 1, 2, 1, 1, 0, 2, 2, 0, 2, 2],
    dtype=np.uint32,
)

ra_deg = np.array(
    [
        20.9191548,
        33.4247141,
        32.1435128,
        33.4159091,
        32.1347282,
        33.3829299,
        20.6388309,
        20.6187259,
        20.6137886,
        32.7525147,
        31.4874917,
        32.4518231,
        32.4495403,
        19.8927380,
        30.6416348,
        30.0938936,
        18.2218784,
        28.3859403,
        28.3818327,
    ],
    dtype=np.float64,
)

dec_deg = np.array(
    [
        20.0550441,
        23.5516817,
        26.5139615,
        23.5525348,
        26.5160622,
        23.5555991,
        20.1218532,
        20.1264229,
        20.1275173,
        23.6064063,
        26.6622284,
        23.6270392,
        23.6272157,
        20.2977473,
        26.8303010,
        26.9256271,
        20.7096409,
        27.1602652,
        27.1606420,
    ],
    dtype=np.float64,
)

jd_utc = np.array(
    [
        2458789.6362963,
        2458789.6381250,
        2458789.6381250,
        2458789.6663773,
        2458789.6663773,
        2458789.7706481,
        2458790.6995023,
        2458790.7733333,
        2458790.7914120,
        2458791.8445602,
        2458791.8445602,
        2458792.8514699,
        2458792.8590741,
        2458793.6896759,
        2458794.7996759,
        2458796.7965162,
        2458801.7863426,
        2458803.7699537,
        2458803.7875231,
    ],
    dtype=np.float64,
)
# Convert times to MJD (TT) using astropy
t_utc = Time(jd_utc, format="jd", scale="utc")
mjd_tt = t_utc.tt.mjd.astype(np.float64)


ts = TrajectorySet.from_numpy_degrees(
    env,
    trajectory_id,
    ra_deg,
    dec_deg,
    error_ra_arcsec=0.3,
    error_dec_arcsec=0.3,
    mjd_tt=mjd_tt,
    observer=obs,
)

params = IODParams.builder().max_triplets(100).do_sequential().build()
ok, errors = ts.estimate_all_orbits(env, params, seed=42)

print("ok keys:", list(ok.keys()))
print("errors:", errors)