Skip to content

The SQL Advantage

The previous pages compared pg_orbit with specific tools — Skyfield, JPL Horizons, GMAT, Radio Jupiter Pro. This page steps back and looks at the thing they all have in common: none of them are SQL.

That sounds obvious, but the implications are deeper than “you can write queries instead of scripts.” SQL brings a set of patterns — compositional, declarative, and optimized by decades of database engine development — that change what’s practical to compute.

Each pattern below includes a concrete pg_orbit example.

generate_series: time series without loops

Section titled “generate_series: time series without loops”

The most common pattern in orbital computation is “evaluate this function at regular time intervals.” In Python, that’s a for-loop. In SQL, it’s generate_series.

-- Mars elevation every 15 minutes for a week
SELECT t,
topo_elevation(planet_observe(4, '40.0N 105.3W 1655m'::observer, t)) AS el
FROM generate_series(
now(),
now() + interval '7 days',
interval '15 minutes'
) AS t;

672 data points. One statement. Change the interval to '1 minute' and you get 10,080 points — same statement, same structure. The query planner decides how to execute it; you describe what you want.

This replaces the boilerplate that dominates astronomical computation scripts: initializing time arrays, iterating, collecting results into lists, converting units. Here, the time array is the generate_series call, and the computation is inline.

generate_series handles regular intervals. For irregular grids — say, the timestamps from an observation log — use a subquery or VALUES list:

-- Compute Jupiter position at each recorded observation time
SELECT o.obs_id,
o.obs_time,
topo_azimuth(planet_observe(5, o.location, o.obs_time)) AS az,
topo_elevation(planet_observe(5, o.location, o.obs_time)) AS el
FROM observations o
WHERE o.target = 'Jupiter'
AND o.obs_time > now() - interval '30 days';

The function evaluates at whatever timestamps exist in your data. No pre-generating a time grid and interpolating.

When you need to evaluate a function across every combination of multiple parameters, SQL’s CROSS JOIN (or implicit comma-join) generates the Cartesian product.

The canonical example: departure date x arrival date for Lambert transfers.

SELECT dep::date AS departure,
arr::date AS arrival,
round(c3_departure::numeric, 2) AS c3_km2s2
FROM generate_series(
'2028-08-01'::timestamptz, '2029-01-01'::timestamptz,
interval '2 days') AS dep,
generate_series(
'2029-03-01'::timestamptz, '2029-10-01'::timestamptz,
interval '2 days') AS arr,
LATERAL lambert_transfer(3, 4, dep, arr) AS xfer
WHERE tof_days > 90;

Two generate_series calls, one LATERAL function. The database generates every combination, evaluates the Lambert solver for each, and returns the results. No nested loops, no progress bars, no managing iteration state.

Which observer has the best view of each planet right now?

SELECT body_id,
obs_name,
round(topo_elevation(
planet_observe(body_id, location, now())
)::numeric, 1) AS el
FROM generate_series(1, 8) AS body_id,
(VALUES
('Boulder', '40.0N 105.3W 1655m'::observer),
('Mauna Kea', '19.8N 155.5W 4205m'::observer),
('Paranal', '24.6S 70.4W 2635m'::observer),
('Tenerife', '28.3N 16.5W 2390m'::observer)
) AS obs(obs_name, location)
WHERE topo_elevation(planet_observe(body_id, location, now())) > 0
ORDER BY body_id, el DESC;

8 planets times 4 observers = 32 evaluations. Filtered to only above-horizon results. Sorted so the best observer for each planet appears first.

This is the pattern that no standalone computation tool can replicate. When orbital data lives in the same database as your other data, correlation is a JOIN — not an export-import-match pipeline.

-- Which satellites are visible AND transmitting on frequencies
-- our receiver can handle?
SELECT s.norad_id,
s.name,
f.downlink_mhz,
f.mode,
round(topo_elevation(
observe(s.tle, '40.0N 105.3W 1655m'::observer, now())
)::numeric, 1) AS el,
round(topo_range(
observe(s.tle, '40.0N 105.3W 1655m'::observer, now())
)::numeric, 0) AS range_km
FROM satellites s
JOIN freq_assignments f ON f.norad_id = s.norad_id
WHERE f.downlink_mhz BETWEEN 144.0 AND 146.0 -- 2m band
AND topo_elevation(observe(s.tle, '40.0N 105.3W 1655m'::observer, now())) > 10
ORDER BY el DESC;

The satellite catalog, frequency database, and orbit propagation all participate in the same query. No intermediate files, no API calls, no data format conversion.

-- Predict passes for our constellation and check against existing schedule
SELECT s.name,
pass_aos(p) AS rise,
pass_max_el(p) AS max_el,
pass_los(p) AS set,
CASE WHEN cs.id IS NOT NULL THEN 'SCHEDULED'
ELSE 'AVAILABLE'
END AS status
FROM satellites s,
LATERAL predict_passes(
s.tle, '40.0N 105.3W 1655m'::observer,
now(), now() + interval '24 hours', 15.0
) p
LEFT JOIN contact_schedule cs
ON cs.norad_id = s.norad_id
AND cs.start_time < pass_los(p)
AND cs.end_time > pass_aos(p)
WHERE s.constellation = 'ORBCOMM'
ORDER BY pass_aos(p);

Every predicted pass, annotated with whether it overlaps an existing scheduled contact. The LEFT JOIN means unscheduled windows show up as ‘AVAILABLE’ — these are the gaps you can fill.

-- Jupiter burst windows, filtered by weather forecast
SELECT t AT TIME ZONE 'America/Denver' AS local_time,
round(jupiter_burst_probability(
io_phase_angle(t),
jupiter_cml('40.0N 105.3W 1655m'::observer, t)
)::numeric, 3) AS burst_prob,
w.cloud_cover_pct,
w.precipitation_mm
FROM generate_series(
now(), now() + interval '12 hours', interval '15 minutes'
) AS t
LEFT JOIN weather_forecast w
ON w.station = 'KBDU'
AND w.forecast_time = date_trunc('hour', t)
WHERE topo_elevation(planet_observe(5, '40.0N 105.3W 1655m'::observer, t)) > 10
AND jupiter_burst_probability(
io_phase_angle(t),
jupiter_cml('40.0N 105.3W 1655m'::observer, t)
) > 0.2
ORDER BY t;

Radio observation is affected by weather differently than optical — rain matters less than ionospheric conditions — but the pattern is the same. Whatever data you have that affects observing decisions, JOIN it in.

All pg_orbit computation functions are declared PARALLEL SAFE. This means PostgreSQL’s query planner can distribute work across multiple CPU cores without any explicit threading or multiprocessing code.

-- Observe all 12,000+ satellites in a catalog
-- PostgreSQL will parallelize this across available cores
EXPLAIN ANALYZE
SELECT s.norad_id,
topo_elevation(observe(s.tle, '40.0N 105.3W 1655m'::observer, now())) AS el
FROM satellites s;

The execution plan will show Parallel Seq Scan or Gather nodes when the planner decides parallelism is worthwhile. You don’t request it, configure worker pools, or manage thread safety. The database handles it.

For explicit control when testing:

-- Force parallel execution with 4 workers
SET max_parallel_workers_per_gather = 4;
SET parallel_tuple_cost = 0;
SET parallel_setup_cost = 0;
SELECT s.norad_id,
topo_elevation(observe(s.tle, '40.0N 105.3W 1655m'::observer, now())) AS el
FROM satellites s
WHERE topo_elevation(observe(s.tle, '40.0N 105.3W 1655m'::observer, now())) > 0;

CREATE MATERIALIZED VIEW: cache expensive computations

Section titled “CREATE MATERIALIZED VIEW: cache expensive computations”

Some computations are expensive to run repeatedly — a 30-day burst calendar, a full catalog observation, a pork chop plot survey. Materialized views store the result and let you query it like a table.

-- Cache tonight's satellite visibility for all observers
CREATE MATERIALIZED VIEW tonight_visibility AS
SELECT s.norad_id,
s.name,
obs.obs_name,
pass_aos(p) AS rise,
pass_max_el(p) AS max_el,
pass_los(p) AS set
FROM satellites s,
(VALUES
('Boulder', '40.0N 105.3W 1655m'::observer),
('London', '51.5N 0.1W 11m'::observer),
('Tokyo', '35.7N 139.7E 40m'::observer)
) AS obs(obs_name, location),
LATERAL predict_passes(
s.tle, obs.location,
now(), now() + interval '24 hours', 10.0
) p;
-- Query it instantly
SELECT * FROM tonight_visibility
WHERE obs_name = 'Boulder'
AND max_el > 45
ORDER BY rise;
-- Refresh when TLEs update
REFRESH MATERIALIZED VIEW tonight_visibility;

The initial computation might take seconds for a large catalog. Subsequent queries against the materialized view are instant — it’s just reading a table. Refresh it when the underlying data changes (new TLEs, new day).

For production systems where you don’t want to block readers during refresh:

CREATE UNIQUE INDEX ON tonight_visibility (norad_id, obs_name, rise);
REFRESH MATERIALIZED VIEW CONCURRENTLY tonight_visibility;

The CONCURRENTLY option requires a unique index but allows queries to continue reading the old data while the refresh runs.

pg_orbit produces structured data. Getting it out of PostgreSQL and into other tools is a COPY statement.

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;

CSV feeds into gnuplot, matplotlib, Excel, R, Observable, or any visualization tool. JSON feeds into web applications, APIs, or document databases. The computation stays in PostgreSQL; rendering happens wherever you prefer.

GiST INDEX: spatial queries on orbital elements

Section titled “GiST INDEX: spatial queries on orbital elements”

pg_orbit’s TLE type supports GiST indexing. This enables spatial-style queries over orbital elements — finding satellites that share similar orbits or screening for conjunction risks.

-- Create a GiST index on the satellite catalog
CREATE INDEX idx_satellites_tle ON satellites USING gist (tle);
-- Find satellites with orbital overlap (similar altitude, inclination, RAAN)
SELECT a.name AS sat_a,
b.name AS sat_b,
a.tle <-> b.tle AS altitude_separation_km
FROM satellites a, satellites b
WHERE a.norad_id < b.norad_id -- Avoid duplicate pairs
AND a.tle && b.tle -- Orbital overlap (GiST accelerated)
AND a.tle <-> b.tle < 50 -- Within 50 km altitude
ORDER BY a.tle <-> b.tle
LIMIT 100;

The && operator tests for orbital overlap — whether two TLEs describe orbits in the same region of space. The <-> operator returns altitude separation in kilometers. Both are accelerated by the GiST index, meaning the database can prune the search space before evaluating expensive propagation.

For a catalog of 12,000 satellites, a full cross-product would be 72 million pairs. The GiST index reduces this to the pairs that are actually in the same orbital regime.

-- Catalog-wide conjunction screening for the next 24 hours
-- GiST index pre-filters to nearby orbital regimes
SELECT a.norad_id AS sat_a,
b.norad_id AS sat_b,
a.tle <-> b.tle AS alt_sep_km
FROM satellites a, satellites b
WHERE a.norad_id < b.norad_id
AND a.tle && b.tle
AND a.tle <-> b.tle < 10
ORDER BY a.tle <-> b.tle;

This is a screening filter, not a precision conjunction analysis. It identifies pairs worth investigating further — the ones where orbital elements suggest close approaches. Detailed conjunction assessment would then propagate those specific pairs at high time resolution.

Window functions: tracking changes over time

Section titled “Window functions: tracking changes over time”

SQL window functions let you compute values relative to neighboring rows — previous values, running averages, ranks within groups — without self-joins or subqueries.

-- Track ISS range and range rate, flagging closest approach
WITH iss_track AS (
SELECT t,
topo_range(observe(iss.tle, '40.0N 105.3W 1655m'::observer, t)) AS range_km,
topo_range_rate(observe(iss.tle, '40.0N 105.3W 1655m'::observer, t)) AS range_rate,
topo_elevation(observe(iss.tle, '40.0N 105.3W 1655m'::observer, t)) AS el
FROM (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 iss(tle),
generate_series(now(), now() + interval '2 hours', interval '30 seconds') AS t
WHERE topo_elevation(observe(iss.tle, '40.0N 105.3W 1655m'::observer, t)) > 0
)
SELECT t,
round(range_km::numeric, 1) AS range_km,
round(range_rate::numeric, 3) AS range_rate_kms,
round(el::numeric, 1) AS elevation,
CASE
WHEN range_rate < 0 AND lead(range_rate) OVER (ORDER BY t) >= 0
THEN 'CLOSEST APPROACH'
ELSE ''
END AS event
FROM iss_track
ORDER BY t;

The lead() window function looks at the next row’s range rate. When range rate crosses from negative to positive, the satellite has passed closest approach. No separate analysis step — it’s computed inline.

-- Which planet reaches the highest elevation each night this month?
WITH hourly AS (
SELECT body_id,
t::date AS night,
topo_elevation(planet_observe(body_id, '40.0N 105.3W 1655m'::observer, t)) AS el
FROM generate_series(1, 8) AS body_id,
generate_series(now(), now() + interval '30 days', interval '1 hour') AS t
)
SELECT night,
body_id,
CASE body_id
WHEN 1 THEN 'Mercury' WHEN 2 THEN 'Venus' WHEN 3 THEN 'Earth'
WHEN 4 THEN 'Mars' WHEN 5 THEN 'Jupiter' WHEN 6 THEN 'Saturn'
WHEN 7 THEN 'Uranus' WHEN 8 THEN 'Neptune'
END AS planet,
round(max_el::numeric, 1) AS peak_el
FROM (
SELECT *,
max(el) OVER (PARTITION BY body_id, night) AS max_el,
ROW_NUMBER() OVER (PARTITION BY night ORDER BY max(el) OVER (PARTITION BY body_id, night) DESC) AS rn
FROM hourly
WHERE el > 0
) sub
WHERE rn = 1
GROUP BY night, body_id, max_el
ORDER BY night;

For each night, this finds which planet reaches the highest elevation — useful for deciding what to observe. The window function handles the ranking within each night without a correlated subquery.

Composition: building complex queries from simple parts

Section titled “Composition: building complex queries from simple parts”

The real power of SQL is that these patterns compose. A single query can use generate_series for time steps, CROSS JOIN for parameter sweeps, JOIN for data correlation, window functions for change detection, and COPY TO for export — all in one statement.

-- Complete observation planning query:
-- For each of 3 observers, for each visible planet tonight,
-- find the 2-hour window with the highest average elevation,
-- export to CSV
COPY (
WITH planet_track AS (
SELECT obs_name, body_id, t,
topo_elevation(planet_observe(body_id, location, t)) AS el
FROM (VALUES
('Boulder', '40.0N 105.3W 1655m'::observer),
('Mauna Kea','19.8N 155.5W 4205m'::observer),
('Paranal', '24.6S 70.4W 2635m'::observer)
) AS obs(obs_name, location),
generate_series(1, 8) AS body_id,
generate_series(
date_trunc('day', now()) + interval '1 hour',
date_trunc('day', now()) + interval '13 hours',
interval '15 minutes'
) AS t
WHERE topo_elevation(planet_observe(body_id, location, t)) > 10
),
windowed AS (
SELECT obs_name, body_id, t,
el,
avg(el) OVER (
PARTITION BY obs_name, body_id
ORDER BY t
ROWS BETWEEN 4 PRECEDING AND 4 FOLLOWING
) AS rolling_avg_el
FROM planet_track
)
SELECT DISTINCT ON (obs_name, body_id)
obs_name,
body_id,
CASE body_id
WHEN 1 THEN 'Mercury' WHEN 2 THEN 'Venus'
WHEN 3 THEN 'Earth' WHEN 4 THEN 'Mars'
WHEN 5 THEN 'Jupiter' WHEN 6 THEN 'Saturn'
WHEN 7 THEN 'Uranus' WHEN 8 THEN 'Neptune'
END AS planet,
t AS best_window_center,
round(rolling_avg_el::numeric, 1) AS avg_el_in_window
FROM windowed
ORDER BY obs_name, body_id, rolling_avg_el DESC
) TO '/tmp/observation_plan.csv' WITH CSV HEADER;

This query:

  1. Generates time steps across the night (generate_series)
  2. Evaluates all 8 planets from 3 observers (CROSS JOIN)
  3. Filters to above-horizon results (WHERE)
  4. Computes a rolling 2-hour average elevation (window function)
  5. Selects the best window for each observer/planet pair (DISTINCT ON)
  6. Exports to CSV (COPY TO)

In a traditional workflow, each of these steps would be a separate script, a separate data file, and a separate tool. In SQL, they compose into a single declarative statement that the database engine optimizes and parallelizes.

That’s the advantage. Not that SQL is a better programming language — it isn’t. But for the specific pattern of “evaluate a function over structured parameter spaces and correlate the results with existing data,” SQL is exactly the right tool. And pg_orbit puts the functions inside the tool.