Overview¶
Global Navigation Satellite Systems (GNSS) are space-based systems that provide precise positioning, navigation, and timing (PNT) information to users anywhere on or near the Earth. GPS (Global Positioning System), developed and operated by the United States, is the most widely used GNSS, but it is part of a broader family that also includes Russia’s GLONASS, Europe’s Galileo, and China’s BeiDou.
GNSS works by measuring the travel time of radio signals transmitted from multiple satellites to a receiver on the ground. By combining these measurements with accurate satellite orbit and clock information, a receiver can determine its position, velocity, and time with high accuracy. Modern GNSS applications range from everyday uses such as smartphone navigation and time synchronization to advanced scientific, aviation, maritime, and geodetic applications requiring centimeter-level precision.
GNSS positioning problems can be encoded as factor graphs in GTSAM; this notebook shows an example of identifying a GNSS receiver’s position using pseudoranges from a RINEX file.
Prerequisites¶
We’ll need to install a couple python packages to help with today’s example. We’ll be using python bindings for the RTKLIB library to parse RINEX files from an example receiver source. Numpy is also installed here for easy array manipulation.
Premise¶
Single-point positioning (SPP) is the most basic GNSS positioning technique, in which a receiver determines its position using pseudorange measurements from multiple satellites and broadcast satellite orbit and clock information. It requires only a single receiver and does not rely on external corrections or reference stations, making it simple and widely accessible. While SPP typically provides meter-level accuracy, its performance is limited by errors such as satellite clock and orbit uncertainties, ionospheric and tropospheric delays, and measurement noise.
We’ll be looking at a receiver mounted on top of an FAA air traffic control center; ZOA1 in Fremont, California. It records pseudorange measurements 24/7 as part of a nationwide NOAA Continuously Operating Reference Stations (CORS) network that provides high-accuracy positioning data for geodesy, surveying, and navigation applications. Data from the CORS network can be accessed through the NOAA CORS Network (NCN) data portal, where we’ll download ZOA1’s observation and navigation files for January 18th, 2026.
import tempfile
import urllib.request
import gzip
from pathlib import Path
import shutil
from IPython.display import IFrame
import numpy as np
from pyproj import Transformer
import pyrtklib as rtklib
import gtsam
from gtsam.symbol_shorthand import B, X
# Create a temporary directory (auto-deleted when closed):
with tempfile.TemporaryDirectory() as tmpdir:
tmpdir = Path(tmpdir)
print(f"Using temporary directory: {str(tmpdir)}")
navigation_compressed = tmpdir / "brdc0180.26n.gz"
navigation_file = tmpdir / "brdc0180.26n"
observable_compressed = tmpdir / "zoa10180.26o.gz"
observable_file = tmpdir / "zoa10180.26o"
# Download ZOA1's files using built-in HTTP library:
print("Downloading receiver files...")
with urllib.request.urlopen("https://noaa-cors-pds.s3.amazonaws.com/rinex/2026/018/brdc0180.26n.gz") as response, open(navigation_compressed, "wb") as f:
f.write(response.read())
with urllib.request.urlopen("https://noaa-cors-pds.s3.amazonaws.com/rinex/2026/018/zoa1/zoa10180.26o.gz") as response, open(observable_compressed, "wb") as f:
f.write(response.read())
# Extract files:
print("Extracting files...")
with gzip.open(navigation_compressed, "rb") as gz, open(navigation_file, "wb") as out:
shutil.copyfileobj(gz, out)
with gzip.open(observable_compressed, "rb") as gz, open(observable_file, "wb") as out:
shutil.copyfileobj(gz, out)
print("Parsing RINEX files...")
obs = rtklib.obs_t()
nav = rtklib.nav_t()
sta = rtklib.sta_t()
load_status = rtklib.readrnx(str(navigation_file), 1, "", obs, nav, sta)
assert load_status == 1, "Loading broadcast ephemeris RINEX"
load_status = rtklib.readrnx(str(observable_file), 1, "", obs, nav, sta)
assert load_status == 1, "Loading receiver observation RINEX"
# Temporary directory and contents are cleaned up automatically
print("Loading done!")
print(f"Found {obs.n} observations")
print(f"Found {nav.n} navigation messages")Compute Satellite Positions¶
Satellite positions must be computed first before receiver position can be solved. GNSS broadcast ephemerides are orbit and clock parameters transmitted directly by navigation satellites as part of their navigation message. They allow receivers to compute satellite positions and clock corrections in real time using a simplified orbital model. While readily available and sufficient for standard navigation, broadcast ephemerides are less accurate than precise ephemerides due to prediction errors and limited modeling of perturbations.
The brdc0180.26n RINEX file contains orbital parameters we need to compute satellite positions. We won’t cover how orbits are stored and modeled in the RINEX files, but you can check out the RTCM website for more information. For now, we’ll simply have RTKLIB’s satposs() function do it for us.
n = 5 # number of observations to process.
rs = rtklib.Arr1Ddouble(6*n) # Satellite positions and velocities:
# rs [(0:2)+i*6] = obs[i] sat position {x,y,z} (m)
# rs [(3:5)+i*6] = obs[i] sat velocity {vx,vy,vz} (m/s)
dts = rtklib.Arr1Ddouble(2*n) # Satellite clock drift:
# dts[(0:1)+i*2] = obs[i] sat clock {bias,drift} (s|s/s)
var = rtklib.Arr1Ddouble(1*n) # Satellite position and clock variances:
# var[i] = obs[i] sat position and clock error variance (m^2)
svh = rtklib.Arr1Dint(1*n) # Satellite health flag (-1:correction not available):
# svh[i] = obs[i] sat health flag
teph = obs.data[0].time
# Compute satellite positions:
rtklib.satposs(
teph, # gtime_t teph I time to select ephemeris (gpst)
obs.data.ptr, # obsd_t *obs I observation data
n, # int n I number of observation data
nav, # nav_t *nav I navigation data
0, # int ephopt I ephemeris option (EPHOPT_???) (EPHOPT_BRDC == 0)
rs, # double *rs O satellite positions and velocities (ecef)
dts, # double *dts O satellite clocks
var, # double *var O sat position and clock error variances (m^2)
svh # int *svh O sat health flag (-1:correction not available)
)Solve Receiver Position and Time¶
Finally! We’re ready for factor graphs. Using the PseudorangeFactor, we can build a simple system that estimates receiver clock bias and 3D ECEF position. The graph topology is fairly straight-forward: One node for clock bias, another node for 3D position, all connected to a series of psuedorange factors that encapsulates satellite positions.
# Prepare factor graph:
cb_key = B(0) # Receiver clock bias variable key.
pos_key = X(0) # Receiver position variable key.
nm = gtsam.noiseModel.Diagonal.Sigmas(np.array([1.0]))
graph = gtsam.NonlinearFactorGraph()
# Add a PseudorangeFactor for each observation:
for i in range(n):
obsd = obs.data[i]
sat_bias = dts[2*i+0]
sat_drift = dts[2*i+1]
sat_pos = np.array([rs[6*i+0], rs[6*i+1], rs[6*i+2]])
# Identify pseudorange code CODE_L1C:
for j, code in enumerate(obsd.code):
if code == rtklib.CODE_L1C:
pseudorange = obsd.P[j]
break
else:
# If no CODE_L1C pseudorange found, skip this observation:
continue
graph.add(
gtsam.PseudorangeFactor(pos_key, cb_key, pseudorange, sat_pos, sat_bias, nm)
)
print(f"=== Satellite {obsd.sat} ===")
print(f"Receiver pseudorange: {pseudorange} meters")
print(f"Observation time (unixtime): {rtklib.time_str(obsd.time, 3)} ({obsd.time.time})")
print(f"Satellite position (ECEF): {1e-3*sat_pos} km")
print(f"Satellite clock drift: bias {sat_bias} seconds, drift {sat_drift} seconds-per-second")
print()We initialize receiver position at origin, and then use a LM solver to optimize the system.
# Solve the system:
initial_values = gtsam.Values()
initial_values.insert(cb_key, 0.0)
initial_values.insert(pos_key, np.zeros(3))
result = gtsam.LevenbergMarquardtOptimizer(graph, initial_values).optimize()
# Process the results:
clock_bias = result.atDouble(cb_key)
receiver_position = result.atVector(pos_key)
print(f"Final clock bias {clock_bias} seconds and receiver position {receiver_position} (ECEF meters)")Nice! We see that the solver converged at a reasonable solution. ECEF is difficult to visualize, however. So we’ll convert the result to geodetic coordinates and plot it on a map:
# Create a transformer from ECEF (EPSG:4978) to WGS84 LLA (EPSG:4326)
ecef2lla = Transformer.from_crs("epsg:4978", "epsg:4326", always_xy=True)
# Convert coordinates coordinates (latitude, longitude)
lon, lat, alt = ecef2lla.transform(receiver_position[0], receiver_position[1], receiver_position[2])
print(f"Latitude: {lat} degrees")
print(f"Longitude: {lon} degrees")
print(f"Altitude: {alt} meters")
# OpenStreetMap embed URL
osm_url = f"https://www.openstreetmap.org/export/embed.html?bbox={lon-0.01}%2C{lat-0.01}%2C{lon+0.01}%2C{lat+0.01}&layer=map&marker={lat}%2C{lon}"
# Embed in notebook
IFrame(osm_url, width=700, height=500)Validate Results¶
Our results appear to be roughly in the right place, but CORS stations actually provide ground-truth coordinates we can use as a measure of accuracy in our solutions. ZOA1’s full precision-surveyed coordinates are published in their data portal, but we’ll examine a just a snippet copied here:
| ITRF2020 POSITION (EPOCH 2020.0) |
| Computed in Apr 2025 using data through gpswk 2237. |
| X = -2684436.824 m latitude = 37 32 34.99617 N |
| Y = -4293336.957 m longitude = 122 00 57.42005 W |
| Z = 3865351.638 m ellipsoid height = -3.960 m |We can compute a simple error metric by subtracting the two coordinates as follows:
ground_truth = np.array([-2684436.824, -4293336.957, 3865351.638])
error = ground_truth - receiver_position
print(f"Total positioning error: {np.linalg.norm(error)} meters")Conclusions¶
33 meters of error puts us in the same building as the receiver, which is a good starting point for simple positioning. As mentioned above, however, decades of advancements have vastly improved GNSS positioning accuracy down to sub-meter or even centimeter-scale precision and error. The simple PseudorangeFactor does not account for atmospheric effects, carrier-phase measurements, nor differential GNSS corrections. So there’s a lot more potential for improvement, but this introductory example provides a foundation for expanding GNSS measurements into the broader GTSAM ecosystem.