From Skyfield to SQL

Skyfield is an excellent Python library for positional astronomy — well-documented, well-tested, and built on the same JPL ephemeris data used by spacecraft navigation teams. If you already use Skyfield, you’ll recognize the computations pg_orrery performs. The difference is where they happen.

This page shows the same tasks done both ways. Not to argue one is better than the other — they make different trade-offs — but to help you decide which fits your workflow.

Observing a planet

The most common starting point: where is Jupiter from my location, right now?

Skyfield (Python)

from skyfield.api import load, Topos

ts = load.timescale()          # downloads finals2000A.all
eph = load('de421.bsp')        # downloads 17MB BSP file

observer = Topos('40.0 N', '105.3 W', elevation_m=1655)
t = ts.now()

astrometric = (eph['earth'] + observer).at(t).observe(eph['jupiter barycenter'])
alt, az, distance = astrometric.apparent().altaz()

print(f"Az: {az.degrees:.2f}  El: {alt.degrees:.2f}  Dist: {distance.au:.4f} AU")

Before this runs, Skyfield downloads two files:

• de421.bsp (17 MB) — JPL planetary ephemeris
• finals2000A.all (3.5 MB) — Earth orientation parameters

These files expire. finals2000A.all needs refreshing every few months. The BSP file itself is stable, but managing file paths across environments (local dev, CI, production) adds friction.

pg_orrery (SQL)

SELECT topo_azimuth(t) AS az,
       topo_elevation(t) AS el,
       topo_range(t) / 149597870.7 AS distance_au
FROM planet_observe(5, '40.0N 105.3W 1655m'::observer, now()) t;

No files to download. No timescale object. The VSOP87 coefficients are compiled into the extension — they ship with CREATE EXTENSION pg_orrery and never expire.

Batch satellite observation

Observe many satellites at the same timestamp. This is where the architectural difference starts to matter.

Skyfield (Python)

from skyfield.api import load, EarthSatellite, Topos

ts = load.timescale()
observer = Topos('40.0 N', '105.3 W', elevation_m=1655)
t = ts.now()

# Load TLE file
with open('catalog.tle') as f:
    lines = f.readlines()

results = []
for i in range(0, len(lines), 3):
    name = lines[i].strip()
    sat = EarthSatellite(lines[i+1], lines[i+2], name, ts)
    topo = (sat - observer).at(t)
    alt, az, dist = topo.altaz()
    if alt.degrees > 0:
        results.append({
            'name': name,
            'az': az.degrees,
            'el': alt.degrees,
            'range_km': dist.km
        })

# Now you have results in a Python list.
# To correlate with other data, you need to:
#   1. Load that data from your database
#   2. Match by satellite name or NORAD ID
#   3. Merge in Python (pandas, dict lookup, etc.)

The results live in Python memory. If you need to correlate with operator contact schedules, frequency assignments, or historical observation logs that live in PostgreSQL, you have to bridge two systems.

pg_orrery (SQL)

-- TLEs already in your database? Just JOIN.
SELECT s.norad_id,
       s.name,
       topo_azimuth(observe(s.tle, obs.location, now())) AS az,
       topo_elevation(observe(s.tle, obs.location, now())) AS el,
       topo_range(observe(s.tle, obs.location, now())) AS range_km
FROM satellites s,
     (VALUES ('40.0N 105.3W 1655m'::observer)) AS obs(location)
WHERE topo_elevation(observe(s.tle, obs.location, now())) > 0;

-- Correlate with frequency assignments in the same query
SELECT s.norad_id,
       s.name,
       f.downlink_mhz,
       topo_elevation(observe(s.tle, obs.location, now())) AS el
FROM satellites s
JOIN freq_assignments f ON f.norad_id = s.norad_id,
     (VALUES ('40.0N 105.3W 1655m'::observer)) AS obs(location)
WHERE topo_elevation(observe(s.tle, obs.location, now())) > 10
ORDER BY el DESC;

The computation and the correlation happen in the same process. No data transfer between Python and PostgreSQL. The query planner can parallelize across cores when scanning large catalogs.

Time series generation

Generate positions over a time range — for plotting an elevation profile, building a ground track, or analyzing visibility windows.

Skyfield (Python)

from skyfield.api import load, Topos
import numpy as np

ts = load.timescale()
eph = load('de421.bsp')
observer = Topos('40.0 N', '105.3 W', elevation_m=1655)

# Generate 144 timestamps over 24 hours
t0 = ts.now()
t1 = ts.tt_jd(t0.tt + 1.0)
times = ts.linspace(t0, t1, 144)

earth_obs = eph['earth'] + observer
jupiter = eph['jupiter barycenter']

elevations = []
for t in times:
    alt, az, dist = earth_obs.at(t).observe(jupiter).apparent().altaz()
    elevations.append(alt.degrees)

# Plot with matplotlib
import matplotlib.pyplot as plt
plt.plot(range(len(elevations)), elevations)
plt.ylabel('Elevation (degrees)')
plt.show()

The loop is explicit. For 144 points this is fast, but the pattern doesn’t parallelize automatically. For larger sweeps (thousands of satellites, days of 1-minute resolution), you manage the iteration yourself.

pg_orrery (SQL)

SELECT t AS time,
       topo_azimuth(obs) AS az,
       topo_elevation(obs) AS el
FROM generate_series(
       now(),
       now() + interval '24 hours',
       interval '10 minutes'
     ) AS t,
LATERAL planet_observe(5, '40.0N 105.3W 1655m'::observer, t) AS obs;

generate_series replaces the Python loop. To change the resolution from 10 minutes to 1 minute, change one parameter. The same pattern works for any observable — planets, satellites, moons, stars.

Export for plotting:

COPY (
  SELECT t,
         topo_elevation(planet_observe(5, '40.0N 105.3W 1655m'::observer, t)) AS el
  FROM generate_series(now(), now() + interval '24 hours', interval '10 minutes') t
) TO '/tmp/jupiter_elevation.csv' WITH CSV HEADER;

Then plot with whatever tool you prefer — gnuplot, matplotlib, Observable, a spreadsheet. pg_orrery produces data; visualization is a separate concern.

Pass prediction

Predict when a satellite will be visible from a location. This is where Skyfield and pg_orrery take genuinely different approaches.

Skyfield (Python)

from skyfield.api import load, EarthSatellite, Topos

ts = load.timescale()
observer = Topos('40.0 N', '105.3 W', elevation_m=1655)

line1 = '1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9025'
line2 = '2 25544  51.6400 208.9163 0006703  30.1694  61.7520 15.50100486 00001'
sat = EarthSatellite(line1, line2, 'ISS', ts)

t0 = ts.now()
t1 = ts.tt_jd(t0.tt + 1.0)

# find_events returns (times, event_types)
# event_type: 0=rise, 1=culminate, 2=set
times, events = sat.find_events(observer, t0, t1, altitude_degrees=10.0)

for ti, event in zip(times, events):
    name = ('rise', 'culminate', 'set')[event]
    print(f"{ti.utc_iso()} — {name}")

Skyfield’s find_events uses root-finding to locate the exact moments when elevation crosses the threshold. This gives sub-second precision for AOS and LOS times.

pg_orrery (SQL)

WITH iss AS (
  SELECT '1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9025
2 25544  51.6400 208.9163 0006703  30.1694  61.7520 15.50100486 00001'::tle AS tle
)
SELECT pass_aos(p) AS rise_time,
       pass_tca(p) AS max_el_time,
       pass_max_el(p) AS max_elevation,
       pass_los(p) AS set_time,
       pass_aos_az(p) AS rise_azimuth,
       pass_los_az(p) AS set_azimuth
FROM iss,
     predict_passes(tle, '40.0N 105.3W 1655m'::observer,
                    now(), now() + interval '24 hours', 10.0) p;

predict_passes returns structured pass_event records. Batch prediction across many satellites is a JOIN:

SELECT s.name,
       pass_aos(p) AS rise,
       pass_max_el(p) AS max_el,
       pass_los(p) AS set
FROM satellites s,
LATERAL predict_passes(s.tle, '40.0N 105.3W 1655m'::observer,
                       now(), now() + interval '24 hours', 10.0) p
WHERE s.constellation = 'Iridium'
ORDER BY pass_aos(p);

Every Iridium pass in the next 24 hours, filtered by constellation, sorted chronologically. Adding JOIN schedules or JOIN ground_contacts keeps the correlation inside the database.

Where Skyfield wins

Skyfield has real advantages

pg_orrery does not replace Skyfield for all use cases. Be clear about where the trade-offs fall.

Apparent-position corrections. Since v0.9.0, pg_orrery applies light-time correction, and since v0.10.0, annual stellar aberration (~20 arcsec) in all _apparent() functions — closing the two largest gaps with Skyfield. However, Skyfield still uses the full IAU 2000A nutation model (~9 arcsec more accurate than pg_orrery’s IAU 1976 precession), polar motion corrections, and delta-T from IERS data. These matter for precision apparent-coordinate work like occultation timing or sub-arcsecond astrometry.

Rise/set finding. find_events() uses numerical root-finding to pinpoint the exact moment a body crosses an elevation threshold. pg_orrery’s predict_passes uses a step-and-refine approach that’s fast for batches but less precise for individual events.

Nutation and polar motion. pg_orrery uses IAU 1976 precession, which introduces a ~9 arcsecond gap compared to Skyfield’s IAU 2000A nutation model. Skyfield also applies polar motion corrections from IERS data. For work requiring sub-arcsecond accuracy (occultation timing, precise astrometric reductions), this gap still matters.

Visualization integration. Skyfield works directly with matplotlib, numpy, and the rest of the Python scientific stack. pg_orrery produces rows and columns — you export to CSV or JSON and then plot separately.

Extensibility. Skyfield handles arbitrary BSP kernels — Pluto, spacecraft, asteroids with precise ephemerides. pg_orrery’s body catalog is fixed at compile time.

Where pg_orrery wins

No file management. No BSP kernels, no timescale data files, no expiring Earth orientation parameters. The computation ships with the extension.

Batch operations at database speed. Observing 12,000 satellites in 17ms. 22,500 Lambert solves in 8.3 seconds. These aren’t optimized benchmarks — they’re SELECT statements running on commodity hardware.

Data correlation. The computation happens where your data lives. JOIN orbital results with contact schedules, frequency assignments, observation logs, or any other table. No ETL pipeline between Python and PostgreSQL.

Automatic parallelism. PostgreSQL’s query planner distributes PARALLEL SAFE functions across available cores. You don’t manage threads or multiprocessing pools.

Reproducibility. A SQL query is a complete, self-contained specification of a computation. No virtual environment, no package versions, no file paths. The same query produces the same result on any PostgreSQL instance with pg_orrery installed.

Orbit determination. Skyfield is a propagation and observation library — it does not fit orbits from observations. Since v0.4.0, pg_orrery can fit TLEs from ECI, topocentric, or angles-only observations via differential correction, entirely in SQL. This closes the loop: observe a satellite, fit its orbit, and propagate the result — all within the database.

-- The closed-loop workflow Skyfield can't do:
-- observe → fit (with range rate + weighted obs) → catalog → predict
WITH obs_data AS (
    -- Collect observations from tonight's tracking pass
    SELECT array_agg(observe(seed_tle, station, t) ORDER BY t) AS topo_obs,
           array_agg(t ORDER BY t) AS times
    FROM generate_series(
             '2024-01-01 12:00+00'::timestamptz,
             '2024-01-01 12:10+00'::timestamptz,
             interval '30 seconds') t
),
fit AS (
    SELECT fitted_tle, rms_final, condition_number, status
    FROM obs_data,
         tle_from_topocentric(
             obs_data.topo_obs, obs_data.times,
             '40.0N 105.3W 1655m'::observer,
             seed := seed_tle,
             fit_range_rate := true,   -- v0.6.0: use Doppler data
             weights := obs_weights    -- v0.6.0: per-obs weighting
         )
    WHERE status = 'converged'
)
-- Immediately predict passes with the freshly-fitted TLE
SELECT pass_aos(p) AS rise,
       pass_max_el(p) AS max_el,
       pass_los(p) AS set
FROM fit,
     predict_passes(fit.fitted_tle,
                    '40.0N 105.3W 1655m'::observer,
                    now(), now() + interval '48 hours', 10.0) p;

Migrating gradually

You don’t have to choose one or the other. A practical migration path:

1. Keep Skyfield for high-precision apparent positions. Anything requiring polar motion, nutation at IAU 2000A level (~9 arcsec beyond pg_orrery’s IAU 1976 precession), or custom BSP kernels stays in Python. For geometric accuracy and apparent-position work including light-time and aberration, pg_orrery v0.10.0 with DE enabled is now competitive.

2. Move batch observation to SQL. If you’re computing positions for hundreds of objects to filter or correlate with database records, pg_orrery eliminates the Python-to-PostgreSQL round trip.

3. Move scheduling to SQL. Pass prediction and visibility windows over time ranges are natural generate_series + predict_passes queries.

4. Move reporting to SQL. “What was above 20 degrees from each of our 5 observers last night?” is a single query with a CROSS JOIN, not a Python loop over observers and timestamps.

5. Move orbit determination to SQL. If you’re fitting orbits with separate Python tools (find_orb, OREKit bindings, custom scripts), pg_orrery’s OD solver can do it alongside your existing satellite catalog. Observations in, TLE out — no data pipeline between Python and PostgreSQL.

6. Enable DE when accuracy matters. If you find VSOP87’s ~1 arcsecond isn’t enough for a specific use case, configure a DE file and switch to _de() function variants. Same SQL patterns, same parameters — just add _de to the function name.