Omar Bohsali

Writing
Quotes
Links

Solving GPS Trilateration in SQL

GPS math is really cool. This morning I remembered a problem set I had in undergrad where we had to understand GPS math from first principles. We used either Octave or MATLAB at the time; this morning I was driving and had the idea: can I do this in postgres? If you've read this blog before, you know that getting postgres to do things it wasn't designed to do like solving crossword puzzles is fun.

Today we'll start with a raw L1 NAV GPS message, convert that to a lat/lng, and implement core GPS math... in postgres.

GPS: Basics & Trilateration

I could go into the GPS protocol and message format, but we'll leave that for a future post. For now, let's assume that anywhere in the world, you can reach at least 3 satellites via radio.

The message you get back from each satellite, the L1 NAV message, contains two critical pieces of data:

  1. Ephemeris Data: Orbital parameters to calculate the satellite's exact position.
  2. Time of Transmission (T_{tx}): The exact timestamp the signal was fired, measured in nanoseconds.

By measuring how long it takes for a signal to travel from a satellite to your receiver (at the speed of light c), we can derive distance.

1ms = 300km

For these calculations, we measure time in nanoseconds, because that accuracy really matters.

Imagine a circle with a 300km radius drawn over a map of your state. If your math is off by just 1 millisecond, your search party will have to search that entire area.

Solving IS-GPS-200 Standard in SQL

GPS satellites are not geosynchronous, meaning you don't inherently know their latitude and longitude. The Ephemeris data actually gives us Keplerian orbital elements. Converting Keplerian elements into precise 3D coordinates requires iteratively solving Kepler's equation first:

M_k = E_k - e \sin(E_k)

This is a transcendental equation. It cannot be solved algebraically for E_k. You have to guess an answer, check the error, and refine the guess.

Writing a numerical solver in C++ is easy. Writing a numerical solver in pure SQL requires weaponizing the WITH RECURSIVE CTE, forcing Postgres to loop over itself to execute fixed-point iteration.

CREATE OR REPLACE FUNCTION st_gps_ephemeris_to_wgs84(
    sqrt_a float8, e float8, m_0 float8, omega float8, 
    i_0 float8, omega_0 float8, delta_n float8, 
    omega_dot float8, i_dot float8, 
    cuc float8, cus float8, crc float8, crs float8, cic float8, cis float8, 
    toe float8, current_t float8
) 
RETURNS geometry AS $$
  WITH constants AS (
    SELECT 
      3.986005e14 AS mu,          -- Earth's gravitational constant (m^3/s^2)
      7.2921151467e-5 AS omega_e  -- Earth's rotation rate (rad/s)
  ),
  step1_prelim AS (
    SELECT 
      (sqrt_a * sqrt_a) AS a,
      (current_t - toe) AS tk,
      SQRT(c.mu / POWER(sqrt_a * sqrt_a, 3)) AS n_0,
      c.omega_e
    FROM constants c
  ),
  step2_mean_anomaly AS (
    SELECT *, (m_0 + (n_0 + delta_n) * tk) AS mk
    FROM step1_prelim
  ),
  -- The Recursive CTE to solve Kepler's Equation
  kepler_iter(iteration, ek, mk, e_val) AS (
    SELECT 0, mk, mk, e FROM step2_mean_anomaly
    UNION ALL
    SELECT iteration + 1, mk + e_val * SIN(ek), mk, e_val
    FROM kepler_iter
    WHERE iteration < 7
  ),
  kepler_solved AS (
    SELECT ek FROM kepler_iter ORDER BY iteration DESC LIMIT 1
  ),
  step3_true_anomaly AS (
    SELECT 
      s2.*, ks.ek,
      ATAN2(SQRT(1 - e * e) * SIN(ks.ek), COS(ks.ek) - e) AS nu_k
    FROM step2_mean_anomaly s2 CROSS JOIN kepler_solved ks
  ),
  step4_perturbations AS (
    SELECT *,
      (nu_k + omega) AS phi_k,
      (cus * SIN(2 * (nu_k + omega)) + cuc * COS(2 * (nu_k + omega))) AS delta_uk,
      (crs * SIN(2 * (nu_k + omega)) + crc * COS(2 * (nu_k + omega))) AS delta_rk,
      (cis * SIN(2 * (nu_k + omega)) + cic * COS(2 * (nu_k + omega))) AS delta_ik
    FROM step3_true_anomaly
  ),
  step5_corrected AS (
    SELECT *,
      (phi_k + delta_uk) AS uk,
      (a * (1 - e * COS(ek)) + delta_rk) AS rk,
      (i_0 + delta_ik + i_dot * tk) AS ik,
      (omega_0 + (omega_dot - omega_e) * tk - omega_e * toe) AS omega_k
    FROM step4_perturbations
  ),
  step6_ecef AS (
    SELECT 
      ((rk * COS(uk)) * COS(omega_k) - (rk * SIN(uk)) * COS(ik) * SIN(omega_k)) AS x_ecef,
      ((rk * COS(uk)) * SIN(omega_k) + (rk * SIN(uk)) * COS(ik) * COS(omega_k)) AS y_ecef,
      ((rk * SIN(uk)) * SIN(ik)) AS z_ecef
    FROM step5_corrected
  )
  -- Return a 3D WGS84 PointZ (Longitude, Latitude, Altitude)
  SELECT 
    ST_Transform(ST_SetSRID(ST_MakePoint(x_ecef, y_ecef, z_ecef), 4978), 4326)
  FROM step6_ecef
$$ LANGUAGE sql IMMUTABLE;

Once we get the Earth-Centered, Earth-Fixed (X_ecef, Y_ecef, Z_ecef) coordinates, converting those back to a usable Latitude/Longitude usually requires another nightmare mathematical process called Bowring's algorithm. But because we are in PostGIS, we just define the raw 3D point as EPSG:4978 (the spatial standard for ECEF space) and ask PostGIS to ST_Transform it to EPSG:4326 (WGS84). PostGIS silently handles the geodetic calculus.

When we pass a satellite's raw ephemeris data into our new function, it hands us back a perfect 3D spatial geometry:

SELECT ST_AsText(
  st_gps_ephemeris_to_wgs84(
    5153.654321, 0.0054321, 1.54321, 1.01234, 0.95876, 0.51234, 1.123e-9, -1.234e-8, 1.123e-10, 
    1.2e-6, 1.3e-6, 1.4e-6, 1.5e-6, 1.6e-6, 1.7e-6, 259200.0, 262800.0
  )
) AS satellite_wgs84;

-- Output: POINT Z (-78.123 38.123 20200000)

2D Trilateration

Assuming you hold a perfectly synchronized atomic clock, your receiver notes the exact Time of Reception (T_{rx}). We calculate the distance d to satellite i:

d_i = c(T_{rx} - T_{tx,i})

In a 3D coordinate system, the intersection of these distances is represented by a system of spheres:

(x - x_i)^2 + (y - y_i)^2 + (z - z_i)^2 = d_i^2

To solve this purely on a 2D surface using PostGIS, we cast satellite coordinates to the geography type and draw distance buffers (ST_Buffer). We then find the geometric center (ST_Centroid) of their intersection (ST_Intersection). For simplicity in these early steps, we'll feed raw Latitude/Longitude directly into the CTE:

WITH constants AS (
  SELECT 
    299792458.0 AS c,                     -- Speed of light (m/s)
    1715284872.000000000 AS clock_time    -- Perfect receiver time
),
l1_nav_messages AS (
  SELECT 0 AS id, 38.123 AS lat, -78.123 AS lng, 1715284871.932000000 AS emission_time, 5.0 AS error_margin UNION ALL
  SELECT 1 AS id, 39.456 AS lat, -77.456 AS lng, 1715284871.921000000 AS emission_time, 5.0 AS error_margin UNION ALL
  SELECT 2 AS id, 37.789 AS lat, -79.789 AS lng, 1715284871.945000000 AS emission_time, 5.0 AS error_margin
)
SELECT 
  ST_AsText(
    ST_Centroid(
      ST_Intersection(
        ST_Intersection(
          ST_Buffer(
            ST_MakePoint((SELECT lng FROM l1_nav_messages WHERE id = 0), (SELECT lat FROM l1_nav_messages WHERE id = 0))::geography,
            ((clock_time - (SELECT emission_time FROM l1_nav_messages WHERE id = 0)) * c) + (SELECT error_margin FROM l1_nav_messages WHERE id = 0)
           )::geometry,

          ST_Buffer(
            ST_MakePoint((SELECT lng FROM l1_nav_messages WHERE id = 1), (SELECT lat FROM l1_nav_messages WHERE id = 1))::geography,
            ((clock_time - (SELECT emission_time FROM l1_nav_messages WHERE id = 1)) * c) + (SELECT error_margin FROM l1_nav_messages WHERE id = 1)
           )::geometry
        ),
        ST_Buffer(
          ST_MakePoint((SELECT lng FROM l1_nav_messages WHERE id = 2), (SELECT lat FROM l1_nav_messages WHERE id = 2))::geography,
          ((clock_time - (SELECT emission_time FROM l1_nav_messages WHERE id = 2)) * c) + (SELECT error_margin FROM l1_nav_messages WHERE id = 2)
         )::geometry
      )
    )
  ) AS estimated_position
FROM constants;

Taking It to the Third Dimension

The 2D solution assumes we are at sea level. In reality, 3D distances form intersecting spheres, not flat circles. Because PostGIS lacks a native ST_3DSphereIntersection function, we must use a 3D Voxel Grid Search.

We generate a vertical line of thousands of points (voxels) at our rough 2D location, stretching into the atmosphere. We then use a spatial cost function (ST_3DDistance) to find the point that minimizes the error between our guessed distances and the actual time-of-flight distances.

WITH constants AS (
  SELECT 299792458.0 AS c, 1715284872.000000000 AS clock_time
),
l1_nav_messages AS (
  SELECT 0 AS id, 38.123 AS lat, -78.123 AS lng, 20200000 AS alt, 1715284871.932000000 AS emission_time UNION ALL
  SELECT 1 AS id, 39.456 AS lat, -77.456 AS lng, 20200000 AS alt, 1715284871.921000000 AS emission_time UNION ALL
  SELECT 2 AS id, 37.789 AS lat, -79.789 AS lng, 20200000 AS alt, 1715284871.945000000 AS emission_time
),
target_distances AS (
  SELECT
    id,
    ST_MakePoint(lng, lat, alt) AS geom_3d,
    ((constants.clock_time - emission_time) * constants.c) AS true_distance
  FROM l1_nav_messages, constants
),
rough_2d_fix AS (
  SELECT ST_SetSRID(ST_MakePoint(-78.503, 38.033), 4326) AS geom 
),
altitude_guesses AS (
  -- Voxel grid from -100m to 10,000m at 1m intervals
  SELECT ST_MakePoint(ST_X(geom), ST_Y(geom), generate_series(-100, 10000, 1)) AS geom_z
  FROM rough_2d_fix
)
SELECT 
  ST_AsText(geom_z) AS estimated_3d_position,
  ST_Z(geom_z) AS calculated_altitude_meters,
  (
    ABS(ST_3DDistance(geom_z, (SELECT geom_3d FROM target_distances WHERE id = 0)) - (SELECT true_distance FROM target_distances WHERE id = 0)) +
    ABS(ST_3DDistance(geom_z, (SELECT geom_3d FROM target_distances WHERE id = 1)) - (SELECT true_distance FROM target_distances WHERE id = 1)) +
    ABS(ST_3DDistance(geom_z, (SELECT geom_3d FROM target_distances WHERE id = 2)) - (SELECT true_distance FROM target_distances WHERE id = 2))
  ) AS total_error
FROM altitude_guesses ag
ORDER BY total_error ASC LIMIT 1;

What if you don't have an atomic clock?

A standard GPS receiver uses a cheap quartz oscillator. At the speed of light, a 1-millisecond clock drift introduces 300 kilometers of error. We must solve for a fourth variable: Receiver Clock Bias (t_{bias}). This requires n=4 satellites.

\sqrt{(X - x_i)^2 + (Y - y_i)^2 + (Z - z_i)^2} = c(T_{rx} - T_{tx,i} - t_{bias})

We upgrade to a 4D Grid Search, executing a Cartesian product (CROSS JOIN) between altitude guesses and nanosecond clock bias guesses.

WITH constants AS (SELECT 299792458.0 AS c),
l1_nav_messages AS (
  SELECT 0 AS id, 38.123 AS lat, -78.123 AS lng, 20200000 AS alt, 0.068 AS time_diff UNION ALL
  SELECT 1 AS id, 39.456 AS lat, -77.456 AS lng, 20200000 AS alt, 0.079 AS time_diff UNION ALL
  SELECT 2 AS id, 37.789 AS lat, -79.789 AS lng, 20200000 AS alt, 0.055 AS time_diff UNION ALL
  SELECT 3 AS id, 36.123 AS lat, -80.123 AS lng, 20200000 AS alt, 0.062 AS time_diff
),
target_sats AS (
  SELECT id, ST_MakePoint(lng, lat, alt) AS geom_3d, time_diff FROM l1_nav_messages
),
altitude_guesses AS (
  SELECT ST_MakePoint(-78.503, 38.033, generate_series(-100, 10000, 10)) AS geom_z
),
clock_bias_guesses AS (
  SELECT generate_series(-5000, 5000, 10) AS t_bias_ns
)
SELECT 
  ST_AsText(ag.geom_z) AS estimated_3d_position,
  cbg.t_bias_ns AS calculated_clock_bias_ns,
  (
    ABS(ST_3DDistance(ag.geom_z, (SELECT geom_3d FROM target_sats WHERE id=0)) - (constants.c * ((SELECT time_diff FROM target_sats WHERE id=0) - (cbg.t_bias_ns / 1e9)))) +
    ABS(ST_3DDistance(ag.geom_z, (SELECT geom_3d FROM target_sats WHERE id=1)) - (constants.c * ((SELECT time_diff FROM target_sats WHERE id=1) - (cbg.t_bias_ns / 1e9)))) +
    ABS(ST_3DDistance(ag.geom_z, (SELECT geom_3d FROM target_sats WHERE id=2)) - (constants.c * ((SELECT time_diff FROM target_sats WHERE id=2) - (cbg.t_bias_ns / 1e9)))) +
    ABS(ST_3DDistance(ag.geom_z, (SELECT geom_3d FROM target_sats WHERE id=3)) - (constants.c * ((SELECT time_diff FROM target_sats WHERE id=3) - (cbg.t_bias_ns / 1e9))))
  ) AS total_error
FROM altitude_guesses ag
CROSS JOIN clock_bias_guesses cbg, constants
ORDER BY total_error ASC LIMIT 1;

Einstein and Spacetime Curvature

Our math is now perfectly aligned, but our physical location would drift by 11.4 kilometers every day.

Combined, satellite clocks run fast by a net 38 \mu s/day. The DoD solves the bulk of this in hardware, deliberately tuning satellite clocks down to 10.22999999543 MHz prior to launch so they speed up to exactly 10.23 MHz in orbit.

However, because orbits are slightly elliptical, we must dynamically adjust the time of transmission (T_{tx}) in software using the relativistic drift (\Delta t_r) broadcasted by the satellite:

T_{tx,corrected} = T_{tx} + \Delta t_r

For the grand finale, we are bringing it all together. We will use our custom st_gps_ephemeris_to_wgs84 UDF to parse raw orbital vectors into 3D positions dynamically, correct the time of flight using relativity, and brute-force the remaining 4 dimensions.

WITH constants AS (SELECT 299792458.0 AS c),
-- Raw L1 Data intercepted by our antenna
l1_nav_messages AS (
  SELECT 
    0 AS id, 0.068 AS time_diff, 0.000000012 AS t_relativity,
    5153.654 AS sqrt_a, 0.005 AS e, 1.543 AS m_0, 1.012 AS omega, 
    0.958 AS i_0, 0.512 AS omega_0, 1.12e-9 AS delta_n, -1.23e-8 AS omega_dot, 
    1.12e-10 AS i_dot, 1.2e-6 AS cuc, 1.3e-6 AS cus, 1.4e-6 AS crc, 1.5e-6 AS crs, 
    1.6e-6 AS cic, 1.7e-6 AS cis, 259200.0 AS toe, 262800.0 AS current_t
  UNION ALL
  SELECT 1, 0.079, 0.000000015, 5153.111, 0.004, 1.1, 1.2, 0.9, 0.5, 1e-9, -1e-8, 1e-10, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 259200.0, 262800.0
  UNION ALL
  SELECT 2, 0.055, 0.000000009, 5153.222, 0.003, 1.2, 1.3, 0.8, 0.6, 1e-9, -1e-8, 1e-10, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 259200.0, 262800.0
  UNION ALL
  SELECT 3, 0.062, 0.000000011, 5153.333, 0.002, 1.3, 1.4, 0.7, 0.7, 1e-9, -1e-8, 1e-10, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 259200.0, 262800.0
),
-- Map the raw parameters through our custom UDF
target_sats AS (
  SELECT 
    id, time_diff, t_relativity,
    st_gps_ephemeris_to_wgs84(
      sqrt_a, e, m_0, omega, i_0, omega_0, delta_n, omega_dot, i_dot, 
      cuc, cus, crc, crs, cic, cis, toe, current_t
    ) AS geom_3d
  FROM l1_nav_messages
),
altitude_guesses AS (
  SELECT ST_MakePoint(-78.503, 38.033, generate_series(-100, 10000, 10)) AS geom_z
),
clock_bias_guesses AS (
  SELECT generate_series(-5000, 5000, 10) AS t_bias_ns
)
SELECT 
  ST_AsText(ag.geom_z) AS estimated_3d_position,
  cbg.t_bias_ns AS calculated_clock_bias_ns,
  (
    ABS(ST_3DDistance(ag.geom_z, (SELECT geom_3d FROM target_sats WHERE id=0)) - (constants.c * (((SELECT time_diff FROM target_sats WHERE id=0) - (SELECT t_relativity FROM target_sats WHERE id=0)) - (cbg.t_bias_ns / 1e9)))) +
    ABS(ST_3DDistance(ag.geom_z, (SELECT geom_3d FROM target_sats WHERE id=1)) - (constants.c * (((SELECT time_diff FROM target_sats WHERE id=1) - (SELECT t_relativity FROM target_sats WHERE id=1)) - (cbg.t_bias_ns / 1e9)))) +
    ABS(ST_3DDistance(ag.geom_z, (SELECT geom_3d FROM target_sats WHERE id=2)) - (constants.c * (((SELECT time_diff FROM target_sats WHERE id=2) - (SELECT t_relativity FROM target_sats WHERE id=2)) - (cbg.t_bias_ns / 1e9)))) +
    ABS(ST_3DDistance(ag.geom_z, (SELECT geom_3d FROM target_sats WHERE id=3)) - (constants.c * (((SELECT time_diff FROM target_sats WHERE id=3) - (SELECT t_relativity FROM target_sats WHERE id=3)) - (cbg.t_bias_ns / 1e9))))
  ) AS total_error
FROM altitude_guesses ag
CROSS JOIN clock_bias_guesses cbg, constants
ORDER BY total_error ASC LIMIT 1;

How Real Receivers Actually Do It

Real navigation chips don't brute-force grids. They use continuous, iterative numerical methods:

1. Gauss-Newton (Least Squares) Iteration: The receiver constructs a Jacobian matrix (J) of partial derivatives, calculating the error (\Delta \rho) between its current guess and actual measurements. It updates its position iteratively:

\Delta x = (J^T J)^{-1} J^T \Delta \rho

2. Extended Kalman Filter (EKF): Instead of calculating a static point, an EKF maintains a state vector of position, time bias, and momentum:

\mathbf{x} = [X, Y, Z, \dot{X}, \dot{Y}, \dot{Z}, t_{bias}, \dot{t}_{bias}]^T

It continuously predicts its future state and updates it based on incoming satellite noise.

Because we ORDER BY total_error ASC LIMIT 1, Postgres has to compute four PostGIS floating-point calculations for all 1.01 million rows, run out of work_mem, spill to the physical hard drive, sort it all on disk, and hand back one row.

It is the antithesis of how a relational database is meant to be used, and the absolute worst way to build a navigation engine. But as an exercise in software engineering? It's pure magic.

Happy Sunday and think twice before you go camping.

Get new posts by email

Get an email whenever I publish a new post.

© 2026 | RSS | @omarish